ROOT logo
/* $Id$ */

//
// Script to calculate PtvsEta correction map at different z of vtx.
// 

makeCorrectionPtEtaVSz(Char_t* dataDir, Int_t nRuns=20) {

  Char_t str[256];
static const Int_t  NZbin=5;
Float_t Zbin[]={-10., -3.,0., 3., 10.};

  gSystem->Load("../libPWG0base.so");
  // ########################################################
  // selection of esd tracks
  AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts();    
  esdTrackCuts->DefineHistograms(1);
  esdTrackCuts->SetMinNClustersTPC(50);
  esdTrackCuts->SetMaxChi2PerClusterTPC(3.5);
  esdTrackCuts->SetMaxCovDiagonalElements(2,2,0.5,0.5,2);
  esdTrackCuts->SetRequireTPCRefit(kTRUE);
  esdTrackCuts->SetMinNsigmaToVertex(3);
  esdTrackCuts->SetAcceptKingDaughters(kFALSE);
  AliLog::SetClassDebugLevel("AliESDtrackCuts",1);
  AliLog::SetClassDebugLevel("ESDtrackQualityCuts",1);

  // ########################################################
  // definition of PtEta correction objects


CorrectionMatrix2D *PtEtaMap[NZbin];
for(Int_t i=0; i<NZbin-1; i++)
  PtEtaMap[i] = new CorrectionMatrix2D(Form("PtvsEta%d",i),Form("%f<z_{vtx}<%f",Zbin[i],Zbin[i+1]),80,0.,10.,120,-2.,2.);

  PtEtaMap[NZbin-1] = new CorrectionMatrix2D("xxx","xxx",80,0.,10.,120,-2.,2.);


  // ########################################################
  // get the data dir  
  Char_t execDir[256];
  sprintf(execDir,gSystem->pwd());
  TSystemDirectory* baseDir = new TSystemDirectory(".",dataDir);
  TList* dirList            = baseDir->GetListOfFiles();
  Int_t nDirs               = dirList->GetEntries();
  // go back to the dir where this script is executed
  gSystem->cd(execDir);

  // ########################################################
  // definition of used pointers
  TFile* esdFile;
  TTree* esdTree;
  TBranch* esdBranch;

  AliESD* esd =0;
  // ########################################################
  // loop over runs
  Int_t nRunCounter = 0;
  for (Int_t r=0; r<nDirs; r++) {

    TSystemFile* presentDir = (TSystemFile*)dirList->At(r);
    if (!presentDir || !presentDir->IsDirectory())
      continue;
    // check that the files are there
    TString currentDataDir;
    currentDataDir.Form("%s/%s",dataDir, presentDir->GetName());
    cout << "Processing directory " << currentDataDir.Data() << endl;
    if ((!gSystem->Which(currentDataDir,"galice.root")) ||
          (!gSystem->Which(currentDataDir,"AliESDs.root")))
      continue;

    if (nRunCounter++ >= nRuns)
      break;

    cout << "run #" << nRunCounter << endl;

    // #########################################################
    // setup galice and runloader
    if (gAlice) {
      delete gAlice->GetRunLoader();
      delete gAlice;
      gAlice=0;
    }

    sprintf(str,"%s/galice.root",currentDataDir.Data());
    AliRunLoader* runLoader = AliRunLoader::Open(str);

    runLoader->LoadgAlice();
    gAlice = runLoader->GetAliRun();
    runLoader->LoadKinematics();
    runLoader->LoadHeader();

   
    // #########################################################
    // open esd file and get the tree

    sprintf(str,"%s/AliESDs.root",currentDataDir.Data());
    // close it first to avoid memory leak
    if (esdFile)
      if (esdFile->IsOpen())
        esdFile->Close();

    esdFile = TFile::Open(str);
    esdTree = (TTree*)esdFile->Get("esdTree");
    if (!esdTree)
      continue;
    esdBranch = esdTree->GetBranch("ESD");
    esdBranch->SetAddress(&esd);
    if (!esdBranch)
      continue;

    // ########################################################
    // Magnetic field
    AliTracker::SetFieldMap(gAlice->Field(),kTRUE); // kTRUE means uniform magnetic field
    //AliKalmanTrack::SetConvConst(1000/0.299792458/5.);

    // ########################################################
    // getting number of events
    Int_t nEvents = (Int_t)runLoader->GetNumberOfEvents();
    Int_t nESDEvents = esdBranch->GetEntries();
    
    if (nEvents!=nESDEvents) {
      cout << " Different number of events from runloader and esdtree!!!" << nEvents << " / " << nESDEvents << endl;
      return;
    }
    // ########################################################
    // loop over number of events
    cout << " looping over events..." << endl;
    for(Int_t i=0; i<nEvents; i++) {

      esdBranch->GetEntry(i);
      runLoader->GetEvent(i);              

      // ########################################################
      // get the EDS vertex
      AliESDVertex* vtxESD = esd->GetVertex();

      Double_t vtx[3];
      Double_t vtx_res[3];
      vtxESD->GetXYZ(vtx);          

    
      vtx_res[0] = vtxESD->GetXRes();
      vtx_res[1] = vtxESD->GetYRes();
      vtx_res[2] = vtxESD->GetZRes();

      // the vertex should be reconstructed
      if (strcmp(vtxESD->GetName(),"default")==0) 
	continue;

      // the resolution should be reasonable???
      if (vtx_res[2]==0 || vtx_res[2]>0.1)
	continue;

Int_t IZbin=NZbin-1;
	  for(Int_t l=0; l<NZbin-1; l++)
            if((vtx[3]<Zbin[l+1])&&((vtx[3]>Zbin[l])))IZbin=l;



      // ########################################################
      // loop over mc particles
      AliStack* particleStack = runLoader->Stack();
      Int_t nPrim    = particleStack->GetNprimary();

      for (Int_t i_mc=0; i_mc<nPrim; i_mc++) {
      
	TParticle* part = particleStack->Particle(i_mc);      
	if (!part || strcmp(part->GetName(),"XXX")==0) 
	  continue;
      
	TParticlePDG* pdgPart = part->GetPDG();

	Bool_t prim = kFALSE;
	// check if it's a primary particle - is there a better way ???
	if ((part->GetFirstDaughter() >= nPrim) || (part->GetFirstDaughter()==-1)) {
	  if (TMath::Abs(pdgPart->PdgCode())>10 && pdgPart->PdgCode()!=21 && strcmp(pdgPart->ParticleClass(),"Unknown")!=0)
	    prim = kTRUE;
	}
	if (!prim)
	  continue;

	if (pdgPart->Charge()==0)
	  continue;

	if (prim)
            PtEtaMap[IZbin]->FillGene(part->Pt(), part->Eta());	
	
      }// end of mc particle

      // ########################################################
      // loop over esd tracks      
      Int_t nTracks = esd->GetNumberOfTracks();            
      for (Int_t t=0; t<nTracks; t++) {
	AliESDtrack* esdTrack = esd->GetTrack(t);      
	
	// cut the esd track?
	if (!esdTrackCuts->AcceptTrack(esdTrack))
	  continue;

        AliTPCtrack* tpcTrack = new AliTPCtrack(*esdTrack);
        if (tpcTrack->GetAlpha()==0) {
          cout << " Warning esd track alpha = 0" << endl;
          continue;
        }

        Float_t theta = tpcTrack->Theta();
        Float_t eta   = -TMath::Log(TMath::Tan(theta/2.));
	Float_t pt=  = tpcTrack->Pt();


	PtEtaMap[IZbin]->FillMeas(eta, pt);	

      } // end of track loop
    } // end  of event loop
  } // end of run loop

  
  for(Int_t l=0; l<NZbin-1; l++)
  {
  PtEtaMap[l]->Divide();  
  PtEtaMap[l]->RemoveEdges();  
  }
  TFile* fout = new TFile("PtEtaCorrectionMap.root","RECREATE");
  
  esdTrackCuts->SaveHistograms("EsdTrackCuts");
  for(Int_t l=0; l<NZbin-1; l++) PtEtaMap[l]->SaveHistograms();
  
  fout->Write();
  fout->Close();

}
 makeCorrectionPtEtaVSz.C:1
 makeCorrectionPtEtaVSz.C:2
 makeCorrectionPtEtaVSz.C:3
 makeCorrectionPtEtaVSz.C:4
 makeCorrectionPtEtaVSz.C:5
 makeCorrectionPtEtaVSz.C:6
 makeCorrectionPtEtaVSz.C:7
 makeCorrectionPtEtaVSz.C:8
 makeCorrectionPtEtaVSz.C:9
 makeCorrectionPtEtaVSz.C:10
 makeCorrectionPtEtaVSz.C:11
 makeCorrectionPtEtaVSz.C:12
 makeCorrectionPtEtaVSz.C:13
 makeCorrectionPtEtaVSz.C:14
 makeCorrectionPtEtaVSz.C:15
 makeCorrectionPtEtaVSz.C:16
 makeCorrectionPtEtaVSz.C:17
 makeCorrectionPtEtaVSz.C:18
 makeCorrectionPtEtaVSz.C:19
 makeCorrectionPtEtaVSz.C:20
 makeCorrectionPtEtaVSz.C:21
 makeCorrectionPtEtaVSz.C:22
 makeCorrectionPtEtaVSz.C:23
 makeCorrectionPtEtaVSz.C:24
 makeCorrectionPtEtaVSz.C:25
 makeCorrectionPtEtaVSz.C:26
 makeCorrectionPtEtaVSz.C:27
 makeCorrectionPtEtaVSz.C:28
 makeCorrectionPtEtaVSz.C:29
 makeCorrectionPtEtaVSz.C:30
 makeCorrectionPtEtaVSz.C:31
 makeCorrectionPtEtaVSz.C:32
 makeCorrectionPtEtaVSz.C:33
 makeCorrectionPtEtaVSz.C:34
 makeCorrectionPtEtaVSz.C:35
 makeCorrectionPtEtaVSz.C:36
 makeCorrectionPtEtaVSz.C:37
 makeCorrectionPtEtaVSz.C:38
 makeCorrectionPtEtaVSz.C:39
 makeCorrectionPtEtaVSz.C:40
 makeCorrectionPtEtaVSz.C:41
 makeCorrectionPtEtaVSz.C:42
 makeCorrectionPtEtaVSz.C:43
 makeCorrectionPtEtaVSz.C:44
 makeCorrectionPtEtaVSz.C:45
 makeCorrectionPtEtaVSz.C:46
 makeCorrectionPtEtaVSz.C:47
 makeCorrectionPtEtaVSz.C:48
 makeCorrectionPtEtaVSz.C:49
 makeCorrectionPtEtaVSz.C:50
 makeCorrectionPtEtaVSz.C:51
 makeCorrectionPtEtaVSz.C:52
 makeCorrectionPtEtaVSz.C:53
 makeCorrectionPtEtaVSz.C:54
 makeCorrectionPtEtaVSz.C:55
 makeCorrectionPtEtaVSz.C:56
 makeCorrectionPtEtaVSz.C:57
 makeCorrectionPtEtaVSz.C:58
 makeCorrectionPtEtaVSz.C:59
 makeCorrectionPtEtaVSz.C:60
 makeCorrectionPtEtaVSz.C:61
 makeCorrectionPtEtaVSz.C:62
 makeCorrectionPtEtaVSz.C:63
 makeCorrectionPtEtaVSz.C:64
 makeCorrectionPtEtaVSz.C:65
 makeCorrectionPtEtaVSz.C:66
 makeCorrectionPtEtaVSz.C:67
 makeCorrectionPtEtaVSz.C:68
 makeCorrectionPtEtaVSz.C:69
 makeCorrectionPtEtaVSz.C:70
 makeCorrectionPtEtaVSz.C:71
 makeCorrectionPtEtaVSz.C:72
 makeCorrectionPtEtaVSz.C:73
 makeCorrectionPtEtaVSz.C:74
 makeCorrectionPtEtaVSz.C:75
 makeCorrectionPtEtaVSz.C:76
 makeCorrectionPtEtaVSz.C:77
 makeCorrectionPtEtaVSz.C:78
 makeCorrectionPtEtaVSz.C:79
 makeCorrectionPtEtaVSz.C:80
 makeCorrectionPtEtaVSz.C:81
 makeCorrectionPtEtaVSz.C:82
 makeCorrectionPtEtaVSz.C:83
 makeCorrectionPtEtaVSz.C:84
 makeCorrectionPtEtaVSz.C:85
 makeCorrectionPtEtaVSz.C:86
 makeCorrectionPtEtaVSz.C:87
 makeCorrectionPtEtaVSz.C:88
 makeCorrectionPtEtaVSz.C:89
 makeCorrectionPtEtaVSz.C:90
 makeCorrectionPtEtaVSz.C:91
 makeCorrectionPtEtaVSz.C:92
 makeCorrectionPtEtaVSz.C:93
 makeCorrectionPtEtaVSz.C:94
 makeCorrectionPtEtaVSz.C:95
 makeCorrectionPtEtaVSz.C:96
 makeCorrectionPtEtaVSz.C:97
 makeCorrectionPtEtaVSz.C:98
 makeCorrectionPtEtaVSz.C:99
 makeCorrectionPtEtaVSz.C:100
 makeCorrectionPtEtaVSz.C:101
 makeCorrectionPtEtaVSz.C:102
 makeCorrectionPtEtaVSz.C:103
 makeCorrectionPtEtaVSz.C:104
 makeCorrectionPtEtaVSz.C:105
 makeCorrectionPtEtaVSz.C:106
 makeCorrectionPtEtaVSz.C:107
 makeCorrectionPtEtaVSz.C:108
 makeCorrectionPtEtaVSz.C:109
 makeCorrectionPtEtaVSz.C:110
 makeCorrectionPtEtaVSz.C:111
 makeCorrectionPtEtaVSz.C:112
 makeCorrectionPtEtaVSz.C:113
 makeCorrectionPtEtaVSz.C:114
 makeCorrectionPtEtaVSz.C:115
 makeCorrectionPtEtaVSz.C:116
 makeCorrectionPtEtaVSz.C:117
 makeCorrectionPtEtaVSz.C:118
 makeCorrectionPtEtaVSz.C:119
 makeCorrectionPtEtaVSz.C:120
 makeCorrectionPtEtaVSz.C:121
 makeCorrectionPtEtaVSz.C:122
 makeCorrectionPtEtaVSz.C:123
 makeCorrectionPtEtaVSz.C:124
 makeCorrectionPtEtaVSz.C:125
 makeCorrectionPtEtaVSz.C:126
 makeCorrectionPtEtaVSz.C:127
 makeCorrectionPtEtaVSz.C:128
 makeCorrectionPtEtaVSz.C:129
 makeCorrectionPtEtaVSz.C:130
 makeCorrectionPtEtaVSz.C:131
 makeCorrectionPtEtaVSz.C:132
 makeCorrectionPtEtaVSz.C:133
 makeCorrectionPtEtaVSz.C:134
 makeCorrectionPtEtaVSz.C:135
 makeCorrectionPtEtaVSz.C:136
 makeCorrectionPtEtaVSz.C:137
 makeCorrectionPtEtaVSz.C:138
 makeCorrectionPtEtaVSz.C:139
 makeCorrectionPtEtaVSz.C:140
 makeCorrectionPtEtaVSz.C:141
 makeCorrectionPtEtaVSz.C:142
 makeCorrectionPtEtaVSz.C:143
 makeCorrectionPtEtaVSz.C:144
 makeCorrectionPtEtaVSz.C:145
 makeCorrectionPtEtaVSz.C:146
 makeCorrectionPtEtaVSz.C:147
 makeCorrectionPtEtaVSz.C:148
 makeCorrectionPtEtaVSz.C:149
 makeCorrectionPtEtaVSz.C:150
 makeCorrectionPtEtaVSz.C:151
 makeCorrectionPtEtaVSz.C:152
 makeCorrectionPtEtaVSz.C:153
 makeCorrectionPtEtaVSz.C:154
 makeCorrectionPtEtaVSz.C:155
 makeCorrectionPtEtaVSz.C:156
 makeCorrectionPtEtaVSz.C:157
 makeCorrectionPtEtaVSz.C:158
 makeCorrectionPtEtaVSz.C:159
 makeCorrectionPtEtaVSz.C:160
 makeCorrectionPtEtaVSz.C:161
 makeCorrectionPtEtaVSz.C:162
 makeCorrectionPtEtaVSz.C:163
 makeCorrectionPtEtaVSz.C:164
 makeCorrectionPtEtaVSz.C:165
 makeCorrectionPtEtaVSz.C:166
 makeCorrectionPtEtaVSz.C:167
 makeCorrectionPtEtaVSz.C:168
 makeCorrectionPtEtaVSz.C:169
 makeCorrectionPtEtaVSz.C:170
 makeCorrectionPtEtaVSz.C:171
 makeCorrectionPtEtaVSz.C:172
 makeCorrectionPtEtaVSz.C:173
 makeCorrectionPtEtaVSz.C:174
 makeCorrectionPtEtaVSz.C:175
 makeCorrectionPtEtaVSz.C:176
 makeCorrectionPtEtaVSz.C:177
 makeCorrectionPtEtaVSz.C:178
 makeCorrectionPtEtaVSz.C:179
 makeCorrectionPtEtaVSz.C:180
 makeCorrectionPtEtaVSz.C:181
 makeCorrectionPtEtaVSz.C:182
 makeCorrectionPtEtaVSz.C:183
 makeCorrectionPtEtaVSz.C:184
 makeCorrectionPtEtaVSz.C:185
 makeCorrectionPtEtaVSz.C:186
 makeCorrectionPtEtaVSz.C:187
 makeCorrectionPtEtaVSz.C:188
 makeCorrectionPtEtaVSz.C:189
 makeCorrectionPtEtaVSz.C:190
 makeCorrectionPtEtaVSz.C:191
 makeCorrectionPtEtaVSz.C:192
 makeCorrectionPtEtaVSz.C:193
 makeCorrectionPtEtaVSz.C:194
 makeCorrectionPtEtaVSz.C:195
 makeCorrectionPtEtaVSz.C:196
 makeCorrectionPtEtaVSz.C:197
 makeCorrectionPtEtaVSz.C:198
 makeCorrectionPtEtaVSz.C:199
 makeCorrectionPtEtaVSz.C:200
 makeCorrectionPtEtaVSz.C:201
 makeCorrectionPtEtaVSz.C:202
 makeCorrectionPtEtaVSz.C:203
 makeCorrectionPtEtaVSz.C:204
 makeCorrectionPtEtaVSz.C:205
 makeCorrectionPtEtaVSz.C:206
 makeCorrectionPtEtaVSz.C:207
 makeCorrectionPtEtaVSz.C:208
 makeCorrectionPtEtaVSz.C:209
 makeCorrectionPtEtaVSz.C:210
 makeCorrectionPtEtaVSz.C:211
 makeCorrectionPtEtaVSz.C:212
 makeCorrectionPtEtaVSz.C:213
 makeCorrectionPtEtaVSz.C:214
 makeCorrectionPtEtaVSz.C:215
 makeCorrectionPtEtaVSz.C:216
 makeCorrectionPtEtaVSz.C:217
 makeCorrectionPtEtaVSz.C:218
 makeCorrectionPtEtaVSz.C:219
 makeCorrectionPtEtaVSz.C:220
 makeCorrectionPtEtaVSz.C:221
 makeCorrectionPtEtaVSz.C:222
 makeCorrectionPtEtaVSz.C:223
 makeCorrectionPtEtaVSz.C:224
 makeCorrectionPtEtaVSz.C:225
 makeCorrectionPtEtaVSz.C:226
 makeCorrectionPtEtaVSz.C:227
 makeCorrectionPtEtaVSz.C:228
 makeCorrectionPtEtaVSz.C:229
 makeCorrectionPtEtaVSz.C:230
 makeCorrectionPtEtaVSz.C:231
 makeCorrectionPtEtaVSz.C:232