]>
Commit | Line | Data |
---|---|---|
75ec0f41 | 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 | } |