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