Removed ; after ClassImp(dNdEtaAnalysis)
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / makeCorrection.C
1 //
2 // Script to make correction maps for dndeta measurements using the
3 // dNdEtaCorrection class.
4 // 
5
6 makeCorrection(Char_t* dataDir, Int_t nRuns=20) {
7
8   Char_t str[256];
9
10   gSystem->Load("../esdTrackCuts/libESDtrackQuality.so");
11   gSystem->Load("libdNdEta.so");
12
13   // ########################################################
14   // selection of esd tracks
15   ESDtrackQualityCuts* esdTrackCuts = new ESDtrackQualityCuts();    
16   esdTrackCuts->DefineHistograms(1);
17   
18   esdTrackCuts->SetMinNClustersTPC(50);
19   esdTrackCuts->SetMaxChi2PerClusterTPC(3.5);
20   esdTrackCuts->SetMaxCovDiagonalElements(2,2,0.5,0.5,2);
21   esdTrackCuts->SetRequireTPCRefit(kTRUE);
22   
23   esdTrackCuts->SetMinNsigmaToVertex(3);
24   esdTrackCuts->SetAcceptKingDaughters(kFALSE);
25
26   AliLog::SetClassDebugLevel("ESDtrackQualityCuts",1);
27
28   // ########################################################
29   // definition of dNdEta correction object
30
31   dNdEtaCorrection* dNdEtaMap = new dNdEtaCorrection();
32
33   // ########################################################
34   // get the data dir  
35   Char_t execDir[256];
36   sprintf(execDir,gSystem->pwd());
37   TSystemDirectory* baseDir = new TSystemDirectory(".",dataDir);
38   TList* dirList            = baseDir->GetListOfFiles();
39   Int_t nDirs               = dirList->GetEntries();
40   // go back to the dir where this script is executed
41   gSystem->cd(execDir);
42
43   // ########################################################
44   // definition of used pointers
45   TFile* esdFile;
46   TTree* esdTree;
47   TBranch* esdBranch;
48
49   AliESD* esd =0;
50
51   // ########################################################
52   // loop over runs
53   Int_t nRunCounter = 0;
54   for (Int_t r=0; r<nDirs; r++) {
55
56     TSystemFile* presentDir = (TSystemFile*)dirList->At(r);
57     if (!presentDir || !presentDir->IsDirectory())
58       continue;
59     // check that the files are there
60     TString currentDataDir;
61     currentDataDir.Form("%s/%s",dataDir, presentDir->GetName());
62     cout << "Processing directory " << currentDataDir.Data() << endl;
63     if ((!gSystem->Which(currentDataDir,"galice.root")) ||
64           (!gSystem->Which(currentDataDir,"AliESDs.root"))) 
65       continue;
66     
67     if (nRunCounter++ >= nRuns)
68       break;    
69     
70     cout << "run #" << nRunCounter << endl;
71
72     // #########################################################
73     // setup galice and runloader
74     if (gAlice) {
75       delete gAlice->GetRunLoader();
76       delete gAlice;
77       gAlice=0;
78     }
79
80     sprintf(str,"%s/galice.root",currentDataDir.Data());
81     AliRunLoader* runLoader = AliRunLoader::Open(str);
82
83     runLoader->LoadgAlice();
84     gAlice = runLoader->GetAliRun();
85     runLoader->LoadKinematics();
86     runLoader->LoadHeader();
87
88     // #########################################################
89     // open esd file and get the tree
90
91     sprintf(str,"%s/AliESDs.root",currentDataDir.Data());
92     // close it first to avoid memory leak
93     if (esdFile)
94       if (esdFile->IsOpen())
95         esdFile->Close();
96
97     esdFile = TFile::Open(str);
98     esdTree = (TTree*)esdFile->Get("esdTree");
99     if (!esdTree)
100       continue;
101     esdBranch = esdTree->GetBranch("ESD");
102     esdBranch->SetAddress(&esd);
103     if (!esdBranch)
104       continue;
105
106     // ########################################################
107     // Magnetic field
108     AliTracker::SetFieldMap(gAlice->Field(),kTRUE); // kTRUE means uniform magnetic field
109     //AliKalmanTrack::SetConvConst(1000/0.299792458/5.);
110
111     // ########################################################
112     // getting number of events
113     Int_t nEvents = (Int_t)runLoader->GetNumberOfEvents();
114     Int_t nESDEvents = esdBranch->GetEntries();
115     
116     if (nEvents!=nESDEvents) {
117       cout << " Different number of events from runloader and esdtree!!!" << nEvents << " / " << nESDEvents << endl;
118       return;
119     }
120     // ########################################################
121     // loop over number of events
122     cout << " looping over events..." << endl;
123     for(Int_t i=0; i<nEvents; i++) {
124
125       esdBranch->GetEntry(i);
126       runLoader->GetEvent(i);              
127
128       // ########################################################
129       // get the EDS vertex
130       AliESDVertex* vtxESD = esd->GetVertex();
131
132       Double_t vtx[3];
133       Double_t vtx_res[3];
134       vtxESD->GetXYZ(vtx);          
135     
136       vtx_res[0] = vtxESD->GetXRes();
137       vtx_res[1] = vtxESD->GetYRes();
138       vtx_res[2] = vtxESD->GetZRes();
139
140       // the vertex should be reconstructed
141       if (strcmp(vtxESD->GetName(),"default")==0) 
142         continue;
143
144       // the resolution should be reasonable???
145       if (vtx_res[2]==0 || vtx_res[2]>0.1)
146         continue;
147
148
149       // ########################################################
150       // get the MC vertex
151       AliGenPythiaEventHeader* genHeader =
152         (AliGenPythiaEventHeader*)runLoader->GetHeader()->GenEventHeader();
153       
154       TArrayF vtxMC(3);
155       genHeader->PrimaryVertex(vtxMC);
156       Double_t vtx[3];
157       vtx[0] = vtxMC[0];
158       vtx[1] = vtxMC[1];
159       vtx[2] = vtxMC[2];
160
161
162       // ########################################################
163       // loop over mc particles
164       AliStack* particleStack = runLoader->Stack();
165       Int_t nPrim    = particleStack->GetNprimary();
166
167       for (Int_t i_mc=0; i_mc<nPrim; i_mc++) {
168       
169         TParticle* part = particleStack->Particle(i_mc);      
170         if (!part || strcmp(part->GetName(),"XXX")==0) 
171           continue;
172       
173         TParticlePDG* pdgPart = part->GetPDG();
174
175         Bool_t prim = kFALSE;
176         // check if it's a primary particle - is there a better way ???
177         if ((part->GetFirstDaughter() >= nPrim) || (part->GetFirstDaughter()==-1)) {
178           if (TMath::Abs(pdgPart->PdgCode())>10 && pdgPart->PdgCode()!=21 && strcmp(pdgPart->ParticleClass(),"Unknown")!=0)
179             prim = kTRUE;
180         }
181         if (!prim)
182           continue;
183
184         if (pdgPart->Charge()==0)
185           continue;     
186
187         Float_t eta = part->Eta();
188
189         if (prim)
190           dNdEtaMap->FillGene(vtx[2], eta);     
191         
192       }// end of mc particle
193
194       // ########################################################
195       // loop over esd tracks      
196       Int_t nTracks = esd->GetNumberOfTracks();            
197       for (Int_t t=0; t<nTracks; t++) {
198         AliESDtrack* esdTrack = esd->GetTrack(t);      
199         
200         // cut the esd track?
201         if (!esdTrackCuts->AcceptTrack(esdTrack))
202           continue;
203         
204         AliTPCtrack* tpcTrack = new AliTPCtrack(*esdTrack);     
205         if (tpcTrack->GetAlpha()==0) {
206           cout << " Warning esd track alpha = 0" << endl;
207           continue;
208         }
209
210         Float_t theta = tpcTrack->Theta();
211         Float_t eta   = -TMath::Log(TMath::Tan(theta/2.));
212
213         // using the eta of the mc particle
214         Int_t label = TMath::Abs(esdTrack->GetLabel());
215         if (label<0) {
216           cout << " cannot find corresponding mc part !!! " << label << endl;
217           continue;
218         }
219         TParticle* mcPart = particleStack->Particle(label);     
220         eta = mcPart->Eta();
221
222         dNdEtaMap->FillMeas(vtx[2], eta);       
223
224       } // end of track loop
225     } // end  of event loop
226   } // end of run loop
227
228   dNdEtaMap->Finish();  
229
230   TFile* fout = new TFile("correction_map.root","RECREATE");
231   
232   esdTrackCuts->SaveHistograms("esd_track_cuts");
233   dNdEtaMap->SaveHistograms();
234   
235   fout->Write();
236   fout->Close();
237   
238   dNdEtaMap->DrawHistograms();
239
240 }