first checkin of Claus' code with a few modifications
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / makeCorrection.C
CommitLineData
75ec0f41 1//
2// Script to make correction maps for dndeta measurements using the
3// dNdEtaCorrection class.
4//
5
6makeCorrection(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}