]>
Commit | Line | Data |
---|---|---|
4ebdd20e | 1 | #include <TFile.h> |
2 | #include <TList.h> | |
3 | #include <TTree.h> | |
4 | #include <TChain.h> | |
5 | #include <TClonesArray.h> | |
6 | #include <TObjString.h> | |
7 | #include <TString.h> | |
8 | #include <TROOT.h> | |
9 | ||
10 | #include "AliHighPtDeDxData.h" | |
11 | #include "AliHighPtDeDxMc.h" | |
12 | #include "AliHighPtDeDxCalib.h" | |
13 | ||
14 | #include <AliXRDPROOFtoolkit.h> | |
15 | ||
16 | #include "DebugClasses.C" | |
17 | #include "my_tools.C" | |
18 | #include "my_functions.C" | |
19 | ||
20 | #include <iostream> | |
21 | #include <fstream> | |
22 | #include <string> | |
23 | ||
24 | using namespace std; | |
25 | ||
26 | /* | |
27 | ||
28 | In this version the start step for calibrations is 2. We do not correct for | |
29 | the Ncl dependence as this is still uncertain of this should be done or not. | |
30 | ||
31 | It is good to make a small test first because the MIP peak might have to be | |
32 | adjusted. With pass0 the default 40-60 window should be good, but | |
33 | ||
34 | To run calibrations: | |
35 | ==================== | |
36 | ||
37 | Use AliRoot because of AliXRDPROOFtoolkit: | |
38 | gSystem->AddIncludePath("-I$ALICE_ROOT/TPC") | |
39 | gSystem->AddIncludePath("-I../lib") | |
40 | gSystem->AddIncludePath("-I../grid") | |
41 | gSystem->AddIncludePath("-I../macros") | |
42 | gROOT->SetMacroPath(".:../macros:../grid:../lib/") | |
43 | .L libMyDeDxAnalysis.so | |
44 | .L my_functions.C+ | |
45 | .L my_tools.C+ | |
46 | .L DebugClasses.C+ | |
47 | .L run_code.C+ | |
48 | ||
49 | // Step 1: create calibrations | |
50 | ||
51 | CreateCalib("7tev_b.dat", kFALSE, "7tev_b.root", 0, 2) | |
52 | ||
53 | // This function is not used anymore | |
54 | // CreateCalibV0("7tev_b_v0.dat", kFALSE, "7tev_b_v0_loose.root", 0) | |
55 | ||
56 | ||
57 | // Step 2 and 3 are found in calibrate_de_dx.C | |
58 | ||
59 | ||
60 | // Step 5: create data | |
61 | CreateOutput("7tev_b.dat", kFALSE, "7tev_b.root", 0, "fitparameters/7tev_b.root") | |
62 | CreateOutputV0("7tev_b_v0.dat", kFALSE, "7tev_b_v0_loose.root", 0, "fitparameters/7tev_b.root") | |
63 | ||
64 | ||
65 | ||
66 | ||
67 | // Test examples | |
68 | // Step 1 | |
69 | CreateCalib("7tev_b_test.dat", kFALSE, "7tev_b_test.root", 0, 2) | |
70 | // Step 5 | |
71 | CreateOutput("7tev_b_test.dat", kFALSE, "7tev_b_test.root", 0, "fitparameters/7tev_b_test.root") | |
72 | CreateOutputV0("7tev_b_test_v0.dat", kFALSE, "7tev_b_test_v0_loose.root", 0, "fitparameters/7tev_b_test.root") | |
73 | ||
74 | */ | |
75 | ||
76 | ||
77 | const Int_t nPtBins = 68; | |
78 | //Double_t xBins[nPtBins+1] = {0., 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0, 1.1 , 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.2, 3.4, 3.6, 3.8, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 8.0, 9.0, 10.0}; | |
79 | Double_t xBins[nPtBins+1] = {0. , 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, | |
80 | 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, | |
81 | 1.0, 1.1 , 1.2, 1.3 , 1.4, 1.5 , 1.6, 1.7 , 1.8, 1.9 , | |
82 | 2.0, 2.2 , 2.4, 2.6 , 2.8, 3.0 , 3.2, 3.4 , 3.6, 3.8 , | |
83 | 4.0, 4.5 , 5.0, 5.5 , 6.0, 6.5 , 7.0, 8.0 , 9.0, 10.0, | |
84 | 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 18.0, 20.0, 22.0, 24.0, | |
85 | 26.0, 28.0, 30.0, 32.0, 34.0, 36.0, 40.0, 45.0, 50.0 }; | |
86 | ||
87 | TFile* dataOutFile = 0; | |
88 | TFile* mcOutFile = 0; | |
89 | ||
90 | TF1* piFunc = 0; | |
91 | TF1* kFunc = 0; | |
92 | TF1* pFunc = 0; | |
93 | TF1* eFunc = 0; | |
94 | TF1* sigmaFunc = 0; | |
95 | ||
96 | void CreateOutput(const Char_t* dataFileName, Bool_t isMc, | |
97 | const Char_t* outfilename, Int_t maxEvents, | |
98 | const Char_t* fitFileName=0); | |
99 | void AddObject(TList* list, Int_t filter, Bool_t phiCut, Int_t run, | |
100 | Bool_t analyzeMc, Bool_t etaCut = kFALSE, Bool_t etaAbs = kFALSE, | |
101 | Int_t etaLow=-8, Int_t etaHigh=8); | |
102 | void CreateOutputV0(const Char_t* dataFileName, Bool_t isMc, | |
103 | const Char_t* outFileName, Int_t maxEvents, | |
104 | const Char_t* fitFileName=0); | |
105 | void AddObjectV0(TList* list, const Char_t* baseName, Bool_t phiCut, Int_t run, | |
106 | Bool_t analyzeMc, Bool_t etaCut = kFALSE, Bool_t etaAbs = kFALSE, | |
107 | Int_t etaLow=-8, Int_t etaHigh=8); | |
108 | ||
109 | void CreateCalib(const Char_t* dataFileName, Bool_t isMc, | |
110 | const Char_t* outfilename, Int_t maxEvents, Int_t startStep); | |
111 | void AddCalibObject(TList* list, Int_t filter, Bool_t phiCut, Int_t run, | |
112 | Bool_t analyzeMc, Bool_t etaCut = kFALSE, Bool_t etaAbs = kFALSE, | |
113 | Int_t etaLow=-8, Int_t etaHigh=8); | |
114 | void CreateCalibV0(const Char_t* dataFileName, Bool_t isMc, | |
115 | const Char_t* outFileName, Int_t maxEvents); | |
116 | void AddCalibObjectV0(TList* list, const Char_t* baseName, Bool_t phiCut, | |
117 | Int_t run, Bool_t analyzeMc); | |
118 | Int_t CalculateFilter(DeDxTrack* track); | |
119 | ||
120 | //____________________________________________________________________________ | |
121 | void CreateOutput(const Char_t* dataFileName, Bool_t isMc, | |
122 | const Char_t* outFileName, Int_t maxEvents, | |
123 | const Char_t* fitFileName) | |
124 | { | |
125 | ||
126 | TF1* fDeDxVsEtaNeg = 0; | |
127 | TF1* fDeDxVsEtaPos = 0; | |
128 | TF1* fDeDxVsNcl = 0; | |
129 | ||
130 | if(fitFileName) { | |
131 | ||
132 | TFile* fitFile = FindFileFresh(fitFileName); | |
133 | if(!fitFile) | |
134 | return; | |
135 | DeDxFitInfo* fitPar = (DeDxFitInfo*)fitFile->Get("fitInfo"); | |
136 | fitPar->Print(); | |
137 | ||
138 | ||
139 | if(!fitPar->calibFileName.IsNull()) { | |
140 | ||
141 | cout << "Setting calibFile: " << fitPar->calibFileName << endl; | |
142 | TFile* calibFile = FindFileFresh(fitPar->calibFileName); | |
143 | if(!calibFile) | |
144 | return; | |
145 | AliHighPtDeDxCalib* calib = (AliHighPtDeDxCalib*)GetObject(calibFile, 1, kTRUE, 0); | |
146 | calib->Print(); | |
147 | fDeDxVsEtaNeg = calib->GetDeDxVsEtaNeg(); | |
148 | fDeDxVsEtaPos = calib->GetDeDxVsEtaPos(); | |
149 | fDeDxVsNcl = calib->GetDeDxVsNcl(); | |
150 | } | |
151 | ||
152 | fixMIP = fitPar->MIP; | |
153 | fixPlateau = fitPar->plateau; | |
154 | ||
155 | Double_t dedxPar[6] = {0, 0, 0, 0, 0, 0}; | |
156 | Double_t sigmaPar[6] = {0, 0, 0, 0, 0, 0}; | |
157 | ||
158 | dedxPar[0] = fitPar->optionDeDx; | |
159 | for(Int_t i = 0; i < fitPar->nDeDxPar; i++) { | |
160 | dedxPar[i+1] = fitPar->parDeDx[i]; | |
161 | } | |
162 | ||
163 | sigmaPar[0] = fitPar->optionSigma; | |
164 | for(Int_t i = 0; i < fitPar->nSigmaPar; i++) { | |
165 | sigmaPar[i+1] = fitPar->parSigma[i]; | |
166 | } | |
167 | ||
168 | piFunc = new TF1("piFunc", FitFunc, 0, 100, fitPar->nDeDxPar+1); | |
169 | piFunc->SetParameters(dedxPar); | |
170 | ||
171 | kFunc = new TF1("kFunc", FitFunc, 0, 100, fitPar->nDeDxPar+1); | |
172 | kFunc->SetParameters(dedxPar); | |
173 | kFunc->SetParameter(0, kFunc->GetParameter(0)+10); | |
174 | ||
175 | pFunc = new TF1("pFunc", FitFunc, 0, 100, fitPar->nDeDxPar+1); | |
176 | pFunc->SetParameters(dedxPar); | |
177 | pFunc->SetParameter(0, pFunc->GetParameter(0)+20); | |
178 | ||
179 | eFunc = new TF1("eFunc", FitFunc, 0, 100, fitPar->nDeDxPar+1); | |
180 | eFunc->SetParameters(dedxPar); | |
181 | eFunc->SetParameter(0, eFunc->GetParameter(0)+30); | |
182 | ||
183 | sigmaFunc = new TF1("sigmaFunc", SigmaFunc, 0, 100, fitPar->nSigmaPar+1); | |
184 | sigmaFunc->SetParameters(sigmaPar); | |
185 | } | |
186 | ||
187 | CreateDir("data"); | |
188 | dataOutFile = new TFile(Form("data/%s", outFileName), "RECREATE"); | |
189 | if(isMc) { | |
190 | ||
191 | CreateDir("mc"); | |
192 | mcOutFile = new TFile(Form("mc/%s", outFileName), "RECREATE"); | |
193 | } | |
194 | ||
195 | // Create output objects | |
196 | dataOutFile->cd(); | |
197 | TList* runList = new TList(); | |
198 | runList->SetOwner(kTRUE); | |
199 | runList->SetBit(TObject::kSingleKey); | |
200 | ||
201 | TList* analysisList = new TList(); | |
202 | analysisList->SetOwner(kFALSE); | |
203 | ||
204 | // TList filter, phi cut, run, isMc | |
205 | AddObject(analysisList, 1, kTRUE, 0, isMc); | |
206 | if(kTRUE) { | |
207 | AddObject(analysisList, 1, kTRUE, 0, isMc, kTRUE, kFALSE, -8, -6); | |
208 | AddObject(analysisList, 1, kTRUE, 0, isMc, kTRUE, kFALSE, -6, -4); | |
209 | AddObject(analysisList, 1, kTRUE, 0, isMc, kTRUE, kFALSE, -4, -2); | |
210 | AddObject(analysisList, 1, kTRUE, 0, isMc, kTRUE, kFALSE, -2, 0); | |
211 | AddObject(analysisList, 1, kTRUE, 0, isMc, kTRUE, kFALSE, 0, 2); | |
212 | AddObject(analysisList, 1, kTRUE, 0, isMc, kTRUE, kFALSE, 2, 4); | |
213 | AddObject(analysisList, 1, kTRUE, 0, isMc, kTRUE, kFALSE, 4, 6); | |
214 | AddObject(analysisList, 1, kTRUE, 0, isMc, kTRUE, kFALSE, 6, 8); | |
215 | } | |
216 | AddObject(analysisList, 1, kTRUE, 0, isMc, kTRUE, kFALSE, -8, 0); | |
217 | AddObject(analysisList, 1, kTRUE, 0, isMc, kTRUE, kFALSE, 0, 8); | |
218 | AddObject(analysisList, 1, kTRUE, 0, isMc, kFALSE, kTRUE, 0, 4); | |
219 | AddObject(analysisList, 1, kTRUE, 0, isMc, kFALSE, kTRUE, 4, 8); | |
220 | AddObject(analysisList, 1, kFALSE, 0, isMc); | |
221 | AddObject(analysisList, 2, kTRUE, 0, isMc); | |
222 | AddObject(analysisList, 2, kFALSE, 0, isMc); | |
223 | ||
224 | TTree* Tree = 0; | |
225 | ||
226 | if(strstr(dataFileName, ".dat")) { | |
227 | ||
228 | AliXRDPROOFtoolkit tool; | |
229 | TChain* chain = tool.MakeChain(dataFileName,"tree", 0, 1000); | |
230 | // chain->Lookup(); | |
231 | Tree = chain; | |
232 | } else { | |
233 | TFile* dataFile = FindFileFresh(dataFileName); | |
234 | if(!dataFile) | |
235 | return; | |
236 | ||
237 | Tree = (TTree*)dataFile->Get("tree"); | |
238 | } | |
239 | ||
240 | TClonesArray* trackArray = 0; | |
241 | TClonesArray* mcTrackArray = 0; | |
242 | DeDxEvent* event = 0; | |
243 | Tree->SetBranchAddress("event", &event); | |
244 | Tree->SetBranchAddress("trackGlobalPar" , &trackArray); | |
245 | if(isMc) | |
246 | Tree->SetBranchAddress("trackMC" , &mcTrackArray); | |
247 | ||
248 | Int_t nEvents = Tree->GetEntries(); | |
249 | cout << "Number of events: " << nEvents << endl; | |
250 | ||
251 | if(maxEvents>0 && maxEvents < nEvents) { | |
252 | ||
253 | nEvents = maxEvents; | |
254 | cout << "N events was reduced to: " << maxEvents << endl; | |
255 | } | |
256 | ||
257 | Int_t currentRun = 0; | |
258 | Int_t nBad = 0; | |
259 | TIter* iter = new TIter(analysisList); | |
260 | ||
261 | for(Int_t n = 0; n < nEvents; n++) { | |
262 | ||
263 | Tree->GetEntry(n); | |
264 | ||
265 | if((n+1)%1000000==0) | |
266 | cout << "Event: " << n+1 << "/" << nEvents << endl; | |
267 | ||
268 | // A few bad entries | |
269 | if(event->run == -1) { | |
270 | nBad++; | |
271 | continue; | |
272 | } | |
273 | ||
274 | if(event->run != currentRun) { | |
275 | ||
276 | cout << "New run: " << event->run << endl; | |
277 | currentRun = event->run; | |
278 | ||
279 | // Check if run objects exist | |
280 | TObjString* runString = new TObjString(Form("%d", currentRun)); | |
281 | if(!runList->FindObject(runString->GetString().Data())) { | |
282 | ||
283 | runList->Add(runString); | |
284 | ||
285 | // TList filter, phi cut, run, isMc | |
286 | AddObject(analysisList, 1, kTRUE, currentRun, isMc); | |
287 | AddObject(analysisList, 1, kFALSE, currentRun, isMc); | |
288 | // AddObject(analysisList, 2, kTRUE, currentRun, isMc); | |
289 | // AddObject(analysisList, 2, kFALSE, currentRun, isMc); | |
290 | ||
291 | // Is this really necessary? | |
292 | delete iter; | |
293 | iter = new TIter(analysisList); | |
294 | ||
295 | } else { | |
296 | ||
297 | delete runString; | |
298 | } | |
299 | } | |
300 | ||
301 | // iterate over analysis list | |
302 | iter->Reset(); | |
303 | AliHighPtDeDxBase* analysis = 0; | |
304 | while ((analysis = dynamic_cast<AliHighPtDeDxBase*> (iter->Next()))) { | |
305 | ||
306 | // First we set the global properties | |
307 | // If I want to use a narrower vtx I have here to redefine my | |
308 | // vtxstatus according to the new vtx range | |
309 | analysis->SetEventRun(event->run); | |
310 | analysis->SetEventMag(event->mag); | |
311 | analysis->SetEventVtxStatus(event->vtxstatus); | |
312 | analysis->SetEventTrigger(event->trig); | |
313 | ||
314 | Int_t vtxstatusmc = 0; | |
315 | if(TMath::Abs(event->zvtxMC) < 10.0) | |
316 | vtxstatusmc = 1; | |
317 | analysis->SetEventVtxStatusMc(vtxstatusmc); | |
318 | ||
319 | if(!analysis->EventAccepted()) // only checks for runs currently | |
320 | continue; | |
321 | ||
322 | // There is a small prolem in making this really nice, because we need | |
323 | // also event info from events that are rejected by the vtx cut to | |
324 | // correct our data. That is why we for bnow only have the run cut there | |
325 | ||
326 | analysis->FillEventInfo(); | |
327 | ||
328 | // First we do the special MC part | |
329 | if(isMc) { | |
330 | ||
331 | AliHighPtDeDxMc* mcAnalysis = dynamic_cast<AliHighPtDeDxMc*> (analysis); | |
332 | if(mcAnalysis) { | |
333 | if(vtxstatusmc) { | |
334 | ||
335 | // loop over mc tracks | |
336 | const Int_t nMcTracks = mcTrackArray->GetEntries(); | |
337 | ||
338 | for(Int_t i = 0; i < nMcTracks; i++) { | |
339 | ||
340 | DeDxTrackMC* trackMC = (DeDxTrackMC*)mcTrackArray->At(i); | |
341 | ||
342 | mcAnalysis->SetTrackPtMc(trackMC->ptMC); | |
343 | mcAnalysis->SetTrackEtaMc(trackMC->etaMC); | |
344 | mcAnalysis->SetTrackChargeMc(trackMC->qMC); | |
345 | mcAnalysis->SetTrackPidMc(trackMC->pidMC); | |
346 | if(mcAnalysis->TrackAcceptedMc()) { | |
347 | // if(trackMC->ptMC<2.0) | |
348 | // mcAnalysis->FillTrackInfoMc(100.0); | |
349 | // else | |
350 | mcAnalysis->FillTrackInfoMc(1.0); | |
351 | } | |
352 | } | |
353 | } | |
354 | } | |
355 | } | |
356 | ||
357 | // The trig==1 is always true for real data, but not for MC data | |
358 | if(!event->trig) | |
359 | continue; | |
360 | ||
361 | if(event->vtxstatus<1) // only fill tracks for events with vtx inside cuts | |
362 | continue; | |
363 | ||
364 | const Int_t nTracks = trackArray->GetEntries(); | |
365 | ||
366 | for(Int_t i = 0; i < nTracks; i++) { | |
367 | ||
368 | DeDxTrack* track = (DeDxTrack*)trackArray->At(i); | |
369 | ||
370 | Double_t eta = track->eta; | |
371 | Double_t dedx = track->dedx; | |
372 | Int_t ncl = track->ncl; | |
373 | if(fDeDxVsEtaNeg) { | |
374 | ||
375 | dedx *= 50.0 / fDeDxVsNcl->Eval(ncl); | |
376 | ||
377 | if(eta < 0) | |
378 | dedx *= 50.0 / fDeDxVsEtaNeg->Eval(eta); | |
379 | else | |
380 | dedx *= 50.0 / fDeDxVsEtaPos->Eval(eta); | |
381 | } | |
382 | ||
383 | analysis->SetTrackCharge(track->q); | |
384 | analysis->SetTrackP(track->p); | |
385 | analysis->SetTrackPt(track->pt); | |
386 | // | |
387 | Int_t filter = CalculateFilter(track); | |
388 | analysis->SetTrackFilter(filter); | |
389 | analysis->SetTrackPhi(track->phi); | |
390 | analysis->SetTrackDeDx(dedx); | |
391 | analysis->SetTrackNcl(ncl); | |
392 | analysis->SetTrackEta(eta); | |
393 | analysis->SetTrackPidMc(track->pid); | |
394 | ||
395 | if(analysis->TrackAccepted()) { | |
396 | // if(track->pt<2.0) | |
397 | // analysis->FillTrackInfo(100.0); | |
398 | // else | |
399 | analysis->FillTrackInfo(1.0); | |
400 | } | |
401 | } | |
402 | } | |
403 | } | |
404 | ||
405 | if(isMc) { | |
406 | ||
407 | TList* runListMc = (TList*)runList->Clone(); | |
408 | mcOutFile->cd(); | |
409 | runListMc->Write(); | |
410 | ||
411 | // I need to somehow specify autowrite, but I am not sure how, so for now | |
412 | // this a bit stupid workaround... | |
413 | iter->Reset(); | |
414 | AliHighPtDeDxBase* analysis = 0; | |
415 | while ((analysis = dynamic_cast<AliHighPtDeDxBase*> (iter->Next()))) { | |
416 | ||
417 | AliHighPtDeDxMc* mcAnalysis = dynamic_cast<AliHighPtDeDxMc*> (analysis); | |
418 | if(mcAnalysis) { | |
419 | mcAnalysis->Write(); | |
420 | } | |
421 | } | |
422 | mcOutFile->Close(); | |
423 | delete mcOutFile; | |
424 | mcOutFile = 0; | |
425 | } | |
426 | ||
427 | dataOutFile->cd(); | |
428 | runList->Write("runList"); | |
429 | iter->Reset(); | |
430 | AliHighPtDeDxBase* analysis = 0; | |
431 | while ((analysis = dynamic_cast<AliHighPtDeDxBase*> (iter->Next()))) { | |
432 | ||
433 | AliHighPtDeDxData* dataAnalysis = dynamic_cast<AliHighPtDeDxData*> (analysis); | |
434 | if(dataAnalysis) { | |
435 | dataAnalysis->Write(); | |
436 | } | |
437 | } | |
438 | dataOutFile->Close(); | |
439 | delete dataOutFile; | |
440 | dataOutFile = 0; | |
441 | ||
442 | cout << "Nbad (runno == -1) : " << nBad << endl; | |
443 | } | |
444 | ||
445 | //___________________________________________________________________________ | |
446 | void AddObject(TList* list, Int_t filter, Bool_t phiCut, Int_t run, | |
447 | Bool_t analyzeMc, Bool_t etaCut, Bool_t etaAbs, | |
448 | Int_t etaLow, Int_t etaHigh) | |
449 | { | |
450 | if(etaCut && etaAbs) { | |
451 | ||
452 | Fatal("AddObject", "You cannot have a cut on both abs and signed eta"); | |
453 | return; | |
454 | } | |
455 | TString objectName("filter"); | |
456 | objectName += filter; | |
457 | if(phiCut) | |
458 | objectName += "phicut"; | |
459 | if(run) { | |
460 | objectName += "_"; | |
461 | objectName += run; | |
462 | } | |
463 | if(etaCut) { | |
464 | objectName += "eta"; | |
465 | objectName += etaLow; | |
466 | objectName += etaHigh; | |
467 | } | |
468 | if(etaAbs) { | |
469 | objectName += "etaabs"; | |
470 | objectName += etaLow; | |
471 | objectName += etaHigh; | |
472 | } | |
473 | ||
474 | dataOutFile->cd(); | |
475 | AliHighPtDeDxData* data = new AliHighPtDeDxData(objectName.Data(), "Analysis data"); | |
476 | if(run) { | |
477 | data->SetUseRunCut(kTRUE); | |
478 | data->SetRun(run); | |
479 | } | |
480 | ||
481 | list->Add(data); | |
482 | data->SetIsMc(analyzeMc); | |
483 | data->SetUseFilterCut(kTRUE); | |
484 | data->SetFilter(filter); | |
485 | data->SetUsePhiCut(phiCut); | |
486 | if(phiCut) { | |
487 | data->SetPhiCutLow(data->GetStandardPhiCutLow()); | |
488 | data->SetPhiCutHigh(data->GetStandardPhiCutHigh()); | |
489 | } | |
490 | if(etaCut) | |
491 | data->SetUseEtaCut(kTRUE); | |
492 | if(etaAbs) | |
493 | data->SetUseEtaCutAbs(kTRUE); | |
494 | if(etaCut || etaAbs) { | |
495 | data->SetEtaLow(Double_t(etaLow)/10.0); | |
496 | data->SetEtaHigh(Double_t(etaHigh)/10.0); | |
497 | } | |
498 | ||
499 | data->Init(nPtBins, xBins); | |
500 | data->Print(); | |
501 | ||
502 | if(piFunc && kFunc && pFunc && eFunc && sigmaFunc) { | |
503 | ||
504 | data->SetPionDeDxFunction(piFunc); | |
505 | data->SetKaonDeDxFunction(kFunc); | |
506 | data->SetProtonDeDxFunction(pFunc); | |
507 | data->SetElectronDeDxFunction(eFunc); | |
508 | data->SetSigmaDeDxFunction(sigmaFunc); | |
509 | } | |
510 | ||
511 | if(analyzeMc) { | |
512 | // create the correction class also | |
513 | ||
514 | mcOutFile->cd(); | |
515 | AliHighPtDeDxMc* mc = new AliHighPtDeDxMc(objectName.Data(), "Analysis mc"); | |
516 | if(run) { | |
517 | mc->SetUseRunCut(kTRUE); | |
518 | mc->SetRun(run); | |
519 | } | |
520 | ||
521 | list->Add(mc); | |
522 | mc->SetUseFilterCut(kTRUE); | |
523 | mc->SetFilter(filter); | |
524 | mc->SetUsePhiCut(phiCut); | |
525 | if(phiCut) { | |
526 | mc->SetPhiCutLow(mc->GetStandardPhiCutLow()); | |
527 | mc->SetPhiCutHigh(mc->GetStandardPhiCutHigh()); | |
528 | } | |
529 | mc->Init(nPtBins, xBins); | |
530 | mc->Print(); | |
531 | } | |
532 | } | |
533 | ||
534 | //____________________________________________________________________________ | |
535 | void CreateOutputV0(const Char_t* dataFileName, Bool_t isMc, | |
536 | const Char_t* outFileName, Int_t maxEvents, | |
537 | const Char_t* fitFileName) | |
538 | { | |
539 | ||
540 | TF1* fDeDxVsEtaNeg = 0; | |
541 | TF1* fDeDxVsEtaPos = 0; | |
542 | TF1* fDeDxVsNcl = 0; | |
543 | ||
544 | if(fitFileName) { | |
545 | ||
546 | TFile* fitFile = FindFileFresh(fitFileName); | |
547 | if(!fitFile) | |
548 | return; | |
549 | DeDxFitInfo* fitPar = (DeDxFitInfo*)fitFile->Get("fitInfo"); | |
550 | fitPar->Print(); | |
551 | ||
552 | ||
553 | if(!fitPar->calibFileName.IsNull()) { | |
554 | ||
555 | cout << "Setting calibFile: " << fitPar->calibFileName << endl; | |
556 | TFile* calibFile = FindFileFresh(fitPar->calibFileName); | |
557 | if(!calibFile) | |
558 | return; | |
559 | AliHighPtDeDxCalib* calib = (AliHighPtDeDxCalib*)GetObject(calibFile, 1, kTRUE, 0); | |
560 | calib->Print(); | |
561 | fDeDxVsEtaNeg = calib->GetDeDxVsEtaNeg(); | |
562 | fDeDxVsEtaPos = calib->GetDeDxVsEtaPos(); | |
563 | fDeDxVsNcl = calib->GetDeDxVsNcl(); | |
564 | } | |
565 | ||
566 | fixMIP = fitPar->MIP; | |
567 | fixPlateau = fitPar->plateau; | |
568 | ||
569 | Double_t dedxPar[6] = {0, 0, 0, 0, 0, 0}; | |
570 | Double_t sigmaPar[6] = {0, 0, 0, 0, 0, 0}; | |
571 | ||
572 | dedxPar[0] = fitPar->optionDeDx; | |
573 | for(Int_t i = 0; i < fitPar->nDeDxPar; i++) { | |
574 | dedxPar[i+1] = fitPar->parDeDx[i]; | |
575 | } | |
576 | ||
577 | sigmaPar[0] = fitPar->optionSigma; | |
578 | for(Int_t i = 0; i < fitPar->nSigmaPar; i++) { | |
579 | sigmaPar[i+1] = fitPar->parSigma[i]; | |
580 | } | |
581 | ||
582 | piFunc = new TF1("piFunc", FitFunc, 0, 100, fitPar->nDeDxPar+1); | |
583 | piFunc->SetParameters(dedxPar); | |
584 | ||
585 | kFunc = new TF1("kFunc", FitFunc, 0, 100, fitPar->nDeDxPar+1); | |
586 | kFunc->SetParameters(dedxPar); | |
587 | kFunc->SetParameter(0, kFunc->GetParameter(0)+10); | |
588 | ||
589 | pFunc = new TF1("pFunc", FitFunc, 0, 100, fitPar->nDeDxPar+1); | |
590 | pFunc->SetParameters(dedxPar); | |
591 | pFunc->SetParameter(0, pFunc->GetParameter(0)+20); | |
592 | ||
593 | eFunc = new TF1("eFunc", FitFunc, 0, 100, fitPar->nDeDxPar+1); | |
594 | eFunc->SetParameters(dedxPar); | |
595 | eFunc->SetParameter(0, eFunc->GetParameter(0)+30); | |
596 | ||
597 | sigmaFunc = new TF1("sigmaFunc", SigmaFunc, 0, 100, fitPar->nSigmaPar+1); | |
598 | sigmaFunc->SetParameters(sigmaPar); | |
599 | } | |
600 | ||
601 | CreateDir("data"); | |
602 | dataOutFile = new TFile(Form("data/%s", outFileName), "RECREATE"); | |
603 | ||
604 | // Create output objects | |
605 | dataOutFile->cd(); | |
606 | TList* runList = new TList(); | |
607 | runList->SetOwner(kTRUE); | |
608 | runList->SetBit(TObject::kSingleKey); | |
609 | ||
610 | TList* analysisList = new TList(); | |
611 | analysisList->SetOwner(kFALSE); | |
612 | ||
613 | // TList v0type, phi cut, run, isMc | |
614 | AddObjectV0(analysisList, "lambda", kTRUE, 0, isMc); | |
615 | AddObjectV0(analysisList, "lambda", kTRUE, 0, isMc, kTRUE, kFALSE, -8, 0); | |
616 | AddObjectV0(analysisList, "lambda", kTRUE, 0, isMc, kTRUE, kFALSE, 0, 8); | |
617 | AddObjectV0(analysisList, "lambda", kTRUE, 0, isMc, kFALSE, kTRUE, 0, 4); | |
618 | AddObjectV0(analysisList, "lambda", kTRUE, 0, isMc, kFALSE, kTRUE, 4, 8); | |
619 | AddObjectV0(analysisList, "kaon", kTRUE, 0, isMc); | |
620 | AddObjectV0(analysisList, "kaon", kTRUE, 0, isMc, kTRUE, kFALSE, -8, 0); | |
621 | AddObjectV0(analysisList, "kaon", kTRUE, 0, isMc, kTRUE, kFALSE, 0, 8); | |
622 | AddObjectV0(analysisList, "kaon", kTRUE, 0, isMc, kFALSE, kTRUE, 0, 4); | |
623 | AddObjectV0(analysisList, "kaon", kTRUE, 0, isMc, kFALSE, kTRUE, 4, 8); | |
624 | ||
625 | TTree* Tree = 0; | |
626 | ||
627 | if(strstr(dataFileName, ".dat")) { | |
628 | ||
629 | AliXRDPROOFtoolkit tool; | |
630 | TChain* chain = tool.MakeChain(dataFileName,"tree", 0, 1000); | |
631 | // chain->Lookup(); | |
632 | Tree = chain; | |
633 | } else { | |
634 | TFile* dataFile = FindFileFresh(dataFileName); | |
635 | if(!dataFile) | |
636 | return; | |
637 | ||
638 | Tree = (TTree*)dataFile->Get("tree"); | |
639 | } | |
640 | ||
641 | TClonesArray* v0Array = 0; | |
642 | DeDxEvent* event = 0; | |
643 | Tree->SetBranchAddress("event", &event); | |
644 | Tree->SetBranchAddress("v0GlobalPar" , &v0Array); | |
645 | ||
646 | Int_t nEvents = Tree->GetEntries(); | |
647 | cout << "Number of events: " << nEvents << endl; | |
648 | ||
649 | if(maxEvents>0 && maxEvents < nEvents) { | |
650 | ||
651 | nEvents = maxEvents; | |
652 | cout << "N events was reduced to: " << maxEvents << endl; | |
653 | } | |
654 | ||
655 | Int_t currentRun = 0; | |
656 | Int_t nBad = 0; | |
657 | TIter* iter = new TIter(analysisList); | |
658 | ||
659 | for(Int_t n = 0; n < nEvents; n++) { | |
660 | ||
661 | Tree->GetEntry(n); | |
662 | ||
663 | if((n+1)%1000000==0) | |
664 | cout << "Event: " << n+1 << "/" << nEvents << endl; | |
665 | ||
666 | if(event->run == -1) { | |
667 | nBad++; | |
668 | continue; | |
669 | } | |
670 | // if(event->run == 126437) | |
671 | // continue; | |
672 | ||
673 | if(event->run != currentRun) { | |
674 | ||
675 | cout << "New run: " << event->run << endl; | |
676 | currentRun = event->run; | |
677 | ||
678 | // Check if run objects exist | |
679 | TObjString* runString = new TObjString(Form("%d", currentRun)); | |
680 | if(!runList->FindObject(runString->GetString().Data())) { | |
681 | ||
682 | runList->Add(runString); | |
683 | ||
684 | // TList v0type, phi cut, run, isMc | |
685 | AddObjectV0(analysisList, "lambda", kTRUE, currentRun, isMc); | |
686 | AddObjectV0(analysisList, "kaon", kTRUE, currentRun, isMc); | |
687 | ||
688 | // Is this really necessary? | |
689 | delete iter; | |
690 | iter = new TIter(analysisList); | |
691 | ||
692 | } else { | |
693 | ||
694 | delete runString; | |
695 | } | |
696 | } | |
697 | ||
698 | // iterate over analysis list | |
699 | iter->Reset(); | |
700 | AliHighPtDeDxBase* analysis = 0; | |
701 | while ((analysis = dynamic_cast<AliHighPtDeDxBase*> (iter->Next()))) { | |
702 | ||
703 | // First we set the global properties | |
704 | // If I want to use a narrower vtx I have here to redefine my | |
705 | // vtxstatus according to the new vtx range | |
706 | analysis->SetEventRun(event->run); | |
707 | analysis->SetEventMag(event->mag); | |
708 | analysis->SetEventVtxStatus(event->vtxstatus); | |
709 | analysis->SetEventTrigger(event->trig); | |
710 | ||
711 | Int_t vtxstatusmc = 0; | |
712 | if(TMath::Abs(event->zvtxMC) < 10.0) | |
713 | vtxstatusmc = 1; | |
714 | analysis->SetEventVtxStatusMc(vtxstatusmc); | |
715 | ||
716 | if(!analysis->EventAccepted()) // only checks for runs currently | |
717 | continue; | |
718 | // There is a small prolem in making this really nice, because we need | |
719 | // also event info from events that are rejected by the vtx cut to | |
720 | // correct our data. That is why we for bnow only have the run cut there | |
721 | ||
722 | analysis->FillEventInfo(); | |
723 | ||
724 | ||
725 | // The trig==1 is always true for real data, but not for MC data | |
726 | if(!event->trig) | |
727 | continue; | |
728 | ||
729 | if(event->vtxstatus<1) // only fill tracks for events with vtx inside cuts | |
730 | continue; | |
731 | ||
732 | const Int_t nV0s = v0Array->GetEntries(); | |
733 | ||
734 | for(Int_t i = 0; i < nV0s; i++) { | |
735 | ||
736 | DeDxV0* v0 = (DeDxV0*)v0Array->At(i); | |
737 | ||
738 | // if(v0->ptrack.filter != 2) | |
739 | // continue; | |
740 | ||
741 | if(v0->status != 0) | |
742 | continue; | |
743 | ||
744 | if(v0->dmassG < 0.1) | |
745 | continue; | |
746 | ||
747 | // if(v0->decayr < 5.0) | |
748 | // continue; | |
749 | ||
750 | // if(v0->decayr > 40.0) | |
751 | // continue; | |
752 | ||
753 | const Double_t dmassK = TMath::Abs(v0->dmassK0); | |
754 | const Double_t dmassL = TMath::Abs(v0->dmassL); | |
755 | const Double_t dmassAL = TMath::Abs(v0->dmassAL); | |
756 | ||
757 | Bool_t fillPos = kFALSE; | |
758 | Bool_t fillNeg = kFALSE; | |
759 | ||
760 | ||
761 | if(strstr(analysis->GetName(), "lambda")) { | |
762 | ||
763 | if(dmassK<0.01) | |
764 | continue; | |
765 | ||
766 | if(dmassL<0.01&&dmassL>0.01) | |
767 | fillPos = kTRUE; | |
768 | ||
769 | if(dmassAL<0.01&&dmassL) | |
770 | fillNeg = kTRUE; | |
771 | } else { // kaons | |
772 | ||
773 | if(dmassL<0.01 || dmassAL<0.01) | |
774 | continue; | |
775 | ||
776 | if(dmassK<0.01) { | |
777 | fillPos = kTRUE; | |
778 | fillNeg = kTRUE; | |
779 | } | |
780 | } | |
781 | ||
782 | for(Int_t j = 0; j < 2; j++) { | |
783 | ||
784 | DeDxTrack* track = 0; | |
785 | ||
786 | if(j==0) { | |
787 | ||
788 | if(fillNeg) | |
789 | track = &(v0->ntrack); | |
790 | else | |
791 | continue; | |
792 | } else { | |
793 | ||
794 | if(fillPos) | |
795 | track = &(v0->ptrack); | |
796 | else | |
797 | continue; | |
798 | } | |
799 | ||
800 | Double_t eta = track->eta; | |
801 | Double_t dedx = track->dedx; | |
802 | Int_t ncl = track->ncl; | |
803 | if(fDeDxVsEtaNeg) { | |
804 | ||
805 | dedx *= 50.0 / fDeDxVsNcl->Eval(ncl); | |
806 | ||
807 | if(eta < 0) | |
808 | dedx *= 50.0 / fDeDxVsEtaNeg->Eval(eta); | |
809 | else | |
810 | dedx *= 50.0 / fDeDxVsEtaPos->Eval(eta); | |
811 | } | |
812 | ||
813 | analysis->SetTrackCharge(track->q); | |
814 | analysis->SetTrackP(track->p); | |
815 | analysis->SetTrackPt(track->pt); | |
816 | // NB! Not used | |
817 | analysis->SetTrackFilter(track->filter); | |
818 | analysis->SetTrackPhi(track->phi); | |
819 | analysis->SetTrackDeDx(dedx); | |
820 | analysis->SetTrackNcl(ncl); | |
821 | analysis->SetTrackBeta(track->beta); | |
822 | analysis->SetTrackPidMc(track->pid); | |
823 | analysis->SetTrackEta(eta); | |
824 | ||
825 | if(analysis->TrackAccepted()) { | |
826 | analysis->FillTrackInfo(1.0); | |
827 | } | |
828 | } | |
829 | } | |
830 | } | |
831 | } | |
832 | ||
833 | dataOutFile->cd(); | |
834 | runList->Write("runList"); | |
835 | iter->Reset(); | |
836 | AliHighPtDeDxBase* analysis = 0; | |
837 | while ((analysis = dynamic_cast<AliHighPtDeDxBase*> (iter->Next()))) { | |
838 | ||
839 | AliHighPtDeDxData* dataAnalysis = dynamic_cast<AliHighPtDeDxData*> (analysis); | |
840 | if(dataAnalysis) { | |
841 | dataAnalysis->Write(); | |
842 | } | |
843 | } | |
844 | dataOutFile->Close(); | |
845 | delete dataOutFile; | |
846 | dataOutFile = 0; | |
847 | ||
848 | cout << "Nbad (runno == -1) : " << nBad << endl; | |
849 | } | |
850 | ||
851 | ||
852 | //___________________________________________________________________________ | |
853 | void AddObjectV0(TList* list, const Char_t* baseName, Bool_t phiCut, Int_t run, | |
854 | Bool_t analyzeMc, Bool_t etaCut, Bool_t etaAbs, | |
855 | Int_t etaLow, Int_t etaHigh) | |
856 | { | |
857 | if(etaCut && etaAbs) { | |
858 | ||
859 | Fatal("AddObject", "You cannot have a cut on both abs and signed eta"); | |
860 | return; | |
861 | } | |
862 | ||
863 | TString objectName(baseName); | |
864 | if(phiCut) | |
865 | objectName += "phicut"; | |
866 | if(run) { | |
867 | objectName += "_"; | |
868 | objectName += run; | |
869 | } | |
870 | if(etaCut) { | |
871 | objectName += "eta"; | |
872 | objectName += etaLow; | |
873 | objectName += etaHigh; | |
874 | } | |
875 | if(etaAbs) { | |
876 | objectName += "etaabs"; | |
877 | objectName += etaLow; | |
878 | objectName += etaHigh; | |
879 | } | |
880 | dataOutFile->cd(); | |
881 | AliHighPtDeDxData* data = new AliHighPtDeDxData(objectName.Data(), "Analysis data"); | |
882 | if(run) { | |
883 | data->SetUseRunCut(kTRUE); | |
884 | data->SetRun(run); | |
885 | } | |
886 | ||
887 | list->Add(data); | |
888 | data->SetIsMc(analyzeMc); | |
889 | data->SetUseFilterCut(kFALSE); | |
890 | ||
891 | data->SetUsePhiCut(phiCut); | |
892 | if(phiCut) { | |
893 | data->SetPhiCutLow(data->GetStandardPhiCutLow()); | |
894 | data->SetPhiCutHigh(data->GetStandardPhiCutHigh()); | |
895 | } | |
896 | if(etaCut) | |
897 | data->SetUseEtaCut(kTRUE); | |
898 | if(etaAbs) | |
899 | data->SetUseEtaCutAbs(kTRUE); | |
900 | if(etaCut || etaAbs) { | |
901 | data->SetEtaLow(Double_t(etaLow)/10.0); | |
902 | data->SetEtaHigh(Double_t(etaHigh)/10.0); | |
903 | } | |
904 | ||
905 | data->Init(nPtBins, xBins); | |
906 | data->Print(); | |
907 | ||
908 | if(piFunc && kFunc && pFunc && eFunc && sigmaFunc) { | |
909 | ||
910 | data->SetPionDeDxFunction(piFunc); | |
911 | data->SetKaonDeDxFunction(kFunc); | |
912 | data->SetProtonDeDxFunction(pFunc); | |
913 | data->SetElectronDeDxFunction(eFunc); | |
914 | data->SetSigmaDeDxFunction(sigmaFunc); | |
915 | } | |
916 | } | |
917 | ||
918 | //____________________________________________________________________________ | |
919 | void CreateCalib(const Char_t* dataFileName, Bool_t isMc, | |
920 | const Char_t* outFileName, Int_t maxEvents, Int_t startStep) | |
921 | { | |
922 | if(startStep < 1 || startStep > 3) { | |
923 | ||
924 | cout << "Start step should be 1,2, or 3 - 2 is recommended" << endl; | |
925 | } | |
926 | ||
927 | if(startStep == 1) { | |
928 | ||
929 | for(Int_t i = 0; i < 10; i++) | |
930 | cout << "Step 1 calibration is depreceated! - 2 is recommended" << endl; | |
931 | } | |
932 | ||
933 | if(startStep<3) { | |
934 | ||
935 | CreateDir("calib_eta"); | |
936 | dataOutFile = new TFile(Form("calib_eta/%s", outFileName), "RECREATE"); | |
937 | } else { | |
938 | ||
939 | CreateDir("no_calib_eta"); | |
940 | dataOutFile = new TFile(Form("no_calib_eta/%s", outFileName), "RECREATE"); | |
941 | } | |
942 | ||
943 | // Create output objects | |
944 | dataOutFile->cd(); | |
945 | TList* runList = new TList(); | |
946 | runList->SetOwner(kTRUE); | |
947 | runList->SetBit(TObject::kSingleKey); | |
948 | ||
949 | TList* calibList = new TList(); | |
950 | calibList->SetOwner(kFALSE); | |
951 | ||
952 | // TList filter, phi cut, run, isMc | |
953 | AddCalibObject(calibList, 1, kTRUE, 0, isMc); | |
954 | AddCalibObject(calibList, 2, kTRUE, 0, isMc); | |
955 | AddCalibObject(calibList, 1, kTRUE, 0, isMc, kTRUE, kFALSE, -8, 0); | |
956 | AddCalibObject(calibList, 1, kTRUE, 0, isMc, kTRUE, kFALSE, 0, 8); | |
957 | AddCalibObject(calibList, 1, kTRUE, 0, isMc, kFALSE, kTRUE, 0, 4); | |
958 | AddCalibObject(calibList, 1, kTRUE, 0, isMc, kFALSE, kTRUE, 4, 8); | |
959 | // AddCalibObject(calibList, 1, kFALSE, 0, isMc); | |
960 | // AddCalibObject(calibList, 2, kTRUE, 0, isMc); | |
961 | // AddCalibObject(calibList, 2, kFALSE, 0, isMc); | |
962 | ||
963 | TTree* Tree = 0; | |
964 | ||
965 | if(strstr(dataFileName, ".dat")) { | |
966 | ||
967 | AliXRDPROOFtoolkit tool; | |
968 | TChain* chain = tool.MakeChain(dataFileName,"tree", 0, 1000); | |
969 | // chain->Lookup(); | |
970 | Tree = chain; | |
971 | } else { | |
972 | TFile* dataFile = FindFileFresh(dataFileName); | |
973 | if(!dataFile) | |
974 | return; | |
975 | ||
976 | Tree = (TTree*)dataFile->Get("tree"); | |
977 | } | |
978 | ||
979 | TClonesArray* trackArray = 0; | |
980 | TClonesArray* mcTrackArray = 0; | |
981 | DeDxEvent* event = 0; | |
982 | Tree->SetBranchAddress("event", &event); | |
983 | Tree->SetBranchAddress("trackGlobalPar" , &trackArray); | |
984 | if(isMc) | |
985 | Tree->SetBranchAddress("trackMC" , &mcTrackArray); | |
986 | ||
987 | Int_t nEvents = Tree->GetEntries(); | |
988 | cout << "Number of events: " << nEvents << endl; | |
989 | ||
990 | if(maxEvents>0 && maxEvents < nEvents) { | |
991 | ||
992 | nEvents = maxEvents; | |
993 | cout << "N events was reduced to: " << maxEvents << endl; | |
994 | } | |
995 | ||
996 | Int_t currentRun = 0; | |
997 | Int_t nBad = 0; | |
998 | TIter* iter = new TIter(calibList); | |
999 | AliHighPtDeDxCalib* calib = 0; | |
1000 | ||
1001 | for(Int_t step = startStep; step <= 3; step++) { | |
1002 | ||
1003 | iter->Reset(); | |
1004 | while ((calib = dynamic_cast<AliHighPtDeDxCalib*> (iter->Next()))) { | |
1005 | ||
1006 | calib->SetStep(step); | |
1007 | // Here you can set the MIP interval | |
1008 | // default is 40-60 | |
1009 | if(step==startStep) { | |
1010 | calib->SetDeDxMIPMin(25); | |
1011 | calib->SetDeDxMIPMax(55); | |
1012 | } else { | |
1013 | calib->SetDeDxMIPMin(40); | |
1014 | calib->SetDeDxMIPMax(60); | |
1015 | } | |
1016 | calib->Init(step, nPtBins, xBins); | |
1017 | } | |
1018 | ||
1019 | for(Int_t n = 0; n < nEvents; n++) { | |
1020 | ||
1021 | Tree->GetEntry(n); | |
1022 | ||
1023 | if((n+1)%1000000==0) | |
1024 | cout << "Event: " << n+1 << "/" << nEvents << endl; | |
1025 | ||
1026 | if(event->run == -1) { | |
1027 | nBad++; | |
1028 | continue; | |
1029 | } | |
1030 | // if(event->run == 126437) | |
1031 | // continue; | |
1032 | ||
1033 | if(event->run != currentRun) { | |
1034 | ||
1035 | cout << "New run: " << event->run << endl; | |
1036 | currentRun = event->run; | |
1037 | ||
1038 | // Check if run objects exist | |
1039 | TObjString* runString = new TObjString(Form("%d", currentRun)); | |
1040 | if(!runList->FindObject(runString->GetString().Data())) { | |
1041 | ||
1042 | runList->Add(runString); | |
1043 | ||
1044 | // TList filter, phi cut, run, isMc | |
1045 | AddCalibObject(calibList, 1, kTRUE, currentRun, isMc); | |
1046 | // AddCalibObject(calibList, 1, kFALSE, currentRun, isMc); | |
1047 | // AddCalibObject(calibList, 2, kTRUE, currentRun, isMc); | |
1048 | // AddCalibObject(calibList, 2, kFALSE, currentRun, isMc); | |
1049 | ||
1050 | // Is this really necessary? | |
1051 | delete iter; | |
1052 | iter = new TIter(calibList); | |
1053 | ||
1054 | iter->Reset(); | |
1055 | while ((calib = dynamic_cast<AliHighPtDeDxCalib*> (iter->Next()))) { | |
1056 | ||
1057 | calib->SetStep(step); | |
1058 | if(step==startStep) { | |
1059 | calib->SetDeDxMIPMin(25); | |
1060 | calib->SetDeDxMIPMax(55); | |
1061 | } else { | |
1062 | calib->SetDeDxMIPMin(40); | |
1063 | calib->SetDeDxMIPMax(60); | |
1064 | } | |
1065 | calib->Init(step, nPtBins, xBins); | |
1066 | } | |
1067 | ||
1068 | } else { | |
1069 | ||
1070 | delete runString; | |
1071 | } | |
1072 | } | |
1073 | ||
1074 | // iterate over calib list | |
1075 | iter->Reset(); | |
1076 | while ((calib = dynamic_cast<AliHighPtDeDxCalib*> (iter->Next()))) { | |
1077 | ||
1078 | // First we set the global properties | |
1079 | // If I want to use a narrower vtx I have here to redefine my | |
1080 | // vtxstatus according to the new vtx range | |
1081 | calib->SetEventRun(event->run); | |
1082 | calib->SetEventMag(event->mag); | |
1083 | calib->SetEventVtxStatus(event->vtxstatus); | |
1084 | calib->SetEventTrigger(event->trig); | |
1085 | ||
1086 | if(!calib->EventAccepted()) // only checks for runs currently | |
1087 | continue; | |
1088 | ||
1089 | calib->FillEventInfo(); | |
1090 | ||
1091 | // The trig==1 is always true for real data, but not for MC data | |
1092 | if(!event->trig) | |
1093 | continue; | |
1094 | ||
1095 | if(event->vtxstatus<1) // only fill tracks for events with vtx inside cuts | |
1096 | continue; | |
1097 | ||
1098 | const Int_t nTracks = trackArray->GetEntries(); | |
1099 | ||
1100 | for(Int_t i = 0; i < nTracks; i++) { | |
1101 | ||
1102 | DeDxTrack* track = (DeDxTrack*)trackArray->At(i); | |
1103 | ||
1104 | calib->SetTrackCharge(track->q); | |
1105 | calib->SetTrackP(track->p); | |
1106 | calib->SetTrackPt(track->pt); | |
1107 | Int_t filter = CalculateFilter(track); | |
1108 | calib->SetTrackFilter(filter); | |
1109 | calib->SetTrackPhi(track->phi); | |
1110 | calib->SetTrackDeDx(track->dedx); | |
1111 | calib->SetTrackNcl(track->ncl); | |
1112 | calib->SetTrackBeta(track->beta); | |
1113 | calib->SetTrackPidMc(track->pid); | |
1114 | calib->SetTrackEta(track->eta); | |
1115 | ||
1116 | if(calib->TrackAccepted()) { | |
1117 | // if(track->pt<2.0) | |
1118 | // calib->FillTrackInfo(100.0); | |
1119 | // else | |
1120 | calib->FillTrackInfo(1.0); | |
1121 | } | |
1122 | } | |
1123 | } | |
1124 | } | |
1125 | ||
1126 | if(step == 2) { // perform eta calibration | |
1127 | ||
1128 | cout << "Starting Eta cal" << endl; | |
1129 | iter->Reset(); | |
1130 | while ((calib = dynamic_cast<AliHighPtDeDxCalib*> (iter->Next()))) { | |
1131 | ||
1132 | calib->PerformEtaCal(); | |
1133 | } | |
1134 | cout << "Finishing Eta cal" << endl; | |
1135 | } else if(step == 1) { // perform ncl calibration | |
1136 | ||
1137 | ||
1138 | cout << "Starting Ncl cal" << endl; | |
1139 | iter->Reset(); | |
1140 | while ((calib = dynamic_cast<AliHighPtDeDxCalib*> (iter->Next()))) { | |
1141 | ||
1142 | calib->PerformNclCal(); | |
1143 | } | |
1144 | cout << "Finished Ncl cal" << endl; | |
1145 | } | |
1146 | } | |
1147 | ||
1148 | cout << "Nbad (runno == -1) : " << nBad << endl; | |
1149 | ||
1150 | dataOutFile->cd(); | |
1151 | runList->Write("runList"); | |
1152 | iter->Reset(); | |
1153 | cout << "Writing calibration data" << endl; | |
1154 | while ((calib = dynamic_cast<AliHighPtDeDxCalib*> (iter->Next()))) { | |
1155 | ||
1156 | cout << "Writing calibration object " << calib->GetName() << endl; | |
1157 | calib->Write(); | |
1158 | } | |
1159 | cout << "Finsihed writing calibration data" << endl; | |
1160 | dataOutFile->Close(); | |
1161 | delete dataOutFile; | |
1162 | dataOutFile = 0; | |
1163 | } | |
1164 | ||
1165 | //___________________________________________________________________________ | |
1166 | void AddCalibObject(TList* list, Int_t filter, Bool_t phiCut, Int_t run, | |
1167 | Bool_t analyzeMc, Bool_t etaCut, Bool_t etaAbs, | |
1168 | Int_t etaLow, Int_t etaHigh) | |
1169 | { | |
1170 | if(etaCut && etaAbs) { | |
1171 | ||
1172 | Fatal("AddCalibObject", "You cannot have a cut on both abs and signed eta"); | |
1173 | return; | |
1174 | } | |
1175 | TString objectName("filter"); | |
1176 | objectName += filter; | |
1177 | if(phiCut) | |
1178 | objectName += "phicut"; | |
1179 | if(run) { | |
1180 | objectName += "_"; | |
1181 | objectName += run; | |
1182 | } | |
1183 | if(etaCut) { | |
1184 | objectName += "eta"; | |
1185 | objectName += etaLow; | |
1186 | objectName += etaHigh; | |
1187 | } | |
1188 | if(etaAbs) { | |
1189 | objectName += "etaabs"; | |
1190 | objectName += etaLow; | |
1191 | objectName += etaHigh; | |
1192 | } | |
1193 | ||
1194 | // dataOutFile->cd(); | |
1195 | AliHighPtDeDxCalib* data = new AliHighPtDeDxCalib(objectName.Data(), "Calib data"); | |
1196 | if(run) { | |
1197 | data->SetUseRunCut(kTRUE); | |
1198 | data->SetRun(run); | |
1199 | } | |
1200 | ||
1201 | list->Add(data); | |
1202 | data->SetIsMc(analyzeMc); | |
1203 | data->SetUseFilterCut(kTRUE); | |
1204 | data->SetFilter(filter); | |
1205 | data->SetUsePhiCut(phiCut); | |
1206 | if(phiCut) { | |
1207 | data->SetPhiCutLow(data->GetStandardPhiCutLow()); | |
1208 | data->SetPhiCutHigh(data->GetStandardPhiCutHigh()); | |
1209 | } | |
1210 | if(etaCut) | |
1211 | data->SetUseEtaCut(kTRUE); | |
1212 | if(etaAbs) | |
1213 | data->SetUseEtaCutAbs(kTRUE); | |
1214 | if(etaCut || etaAbs) { | |
1215 | data->SetEtaLow(Double_t(etaLow)/10.0); | |
1216 | data->SetEtaHigh(Double_t(etaHigh)/10.0); | |
1217 | } | |
1218 | if(run) { | |
1219 | data->SetUseRunCut(kTRUE); | |
1220 | data->SetRun(run); | |
1221 | } | |
1222 | ||
1223 | data->Init(nPtBins, xBins); | |
1224 | data->Print(); | |
1225 | } | |
1226 | ||
1227 | //____________________________________________________________________________ | |
1228 | void CreateCalibV0(const Char_t* dataFileName, Bool_t isMc, | |
1229 | const Char_t* outFileName, Int_t maxEvents) | |
1230 | { | |
1231 | CreateDir("no_calib_eta"); | |
1232 | dataOutFile = new TFile(Form("no_calib_eta/%s", outFileName), "RECREATE"); | |
1233 | ||
1234 | // Create output objects | |
1235 | dataOutFile->cd(); | |
1236 | TList* runList = new TList(); | |
1237 | runList->SetOwner(kTRUE); | |
1238 | runList->SetBit(TObject::kSingleKey); | |
1239 | ||
1240 | TList* calibList = new TList(); | |
1241 | calibList->SetOwner(kFALSE); | |
1242 | ||
1243 | // TList filter, phi cut, run, isMc | |
1244 | AddCalibObjectV0(calibList, "lambda", kTRUE, 0, isMc); | |
1245 | // AddCalibObject(calibList, 1, kFALSE, 0, isMc); | |
1246 | AddCalibObjectV0(calibList, "kaon", kTRUE, 0, isMc); | |
1247 | // AddCalibObject(calibList, 2, kFALSE, 0, isMc); | |
1248 | ||
1249 | TTree* Tree = 0; | |
1250 | ||
1251 | if(strstr(dataFileName, ".dat")) { | |
1252 | ||
1253 | AliXRDPROOFtoolkit tool; | |
1254 | TChain* chain = tool.MakeChain(dataFileName,"tree", 0, 1000); | |
1255 | // chain->Lookup(); | |
1256 | Tree = chain; | |
1257 | } else { | |
1258 | TFile* dataFile = FindFileFresh(dataFileName); | |
1259 | if(!dataFile) | |
1260 | return; | |
1261 | ||
1262 | Tree = (TTree*)dataFile->Get("tree"); | |
1263 | } | |
1264 | ||
1265 | TClonesArray* v0Array = 0; | |
1266 | DeDxEvent* event = 0; | |
1267 | Tree->SetBranchAddress("event", &event); | |
1268 | Tree->SetBranchAddress("v0" , &v0Array); | |
1269 | ||
1270 | Int_t nEvents = Tree->GetEntries(); | |
1271 | cout << "Number of events: " << nEvents << endl; | |
1272 | ||
1273 | if(maxEvents>0 && maxEvents < nEvents) { | |
1274 | ||
1275 | nEvents = maxEvents; | |
1276 | cout << "N events was reduced to: " << maxEvents << endl; | |
1277 | } | |
1278 | ||
1279 | Int_t currentRun = 0; | |
1280 | Int_t nBad = 0; | |
1281 | TIter* iter = new TIter(calibList); | |
1282 | AliHighPtDeDxCalib* calib = 0; | |
1283 | ||
1284 | iter->Reset(); | |
1285 | while ((calib = dynamic_cast<AliHighPtDeDxCalib*> (iter->Next()))) { | |
1286 | ||
1287 | calib->SetStep(3); | |
1288 | } | |
1289 | ||
1290 | for(Int_t n = 0; n < nEvents; n++) { | |
1291 | ||
1292 | Tree->GetEntry(n); | |
1293 | ||
1294 | if((n+1)%1000000==0) | |
1295 | cout << "Event: " << n+1 << "/" << nEvents << endl; | |
1296 | ||
1297 | if(event->run == -1) { | |
1298 | nBad++; | |
1299 | continue; | |
1300 | } | |
1301 | // if(event->run == 126437) | |
1302 | // continue; | |
1303 | ||
1304 | if(event->run != currentRun) { | |
1305 | ||
1306 | cout << "New run: " << event->run << endl; | |
1307 | currentRun = event->run; | |
1308 | ||
1309 | // Check if run objects exist | |
1310 | TObjString* runString = new TObjString(Form("%d", currentRun)); | |
1311 | if(!runList->FindObject(runString->GetString().Data())) { | |
1312 | ||
1313 | runList->Add(runString); | |
1314 | ||
1315 | // TList filter, phi cut, run, isMc | |
1316 | AddCalibObjectV0(calibList, "lambda", kTRUE, currentRun, isMc); | |
1317 | AddCalibObjectV0(calibList, "kaon", kTRUE, currentRun, isMc); | |
1318 | ||
1319 | // Is this really necessary? | |
1320 | delete iter; | |
1321 | iter = new TIter(calibList); | |
1322 | ||
1323 | iter->Reset(); | |
1324 | while ((calib = dynamic_cast<AliHighPtDeDxCalib*> (iter->Next()))) { | |
1325 | ||
1326 | calib->SetStep(3); | |
1327 | } | |
1328 | ||
1329 | } else { | |
1330 | ||
1331 | delete runString; | |
1332 | } | |
1333 | } | |
1334 | ||
1335 | // iterate over calib list | |
1336 | iter->Reset(); | |
1337 | while ((calib = dynamic_cast<AliHighPtDeDxCalib*> (iter->Next()))) { | |
1338 | ||
1339 | // First we set the global properties | |
1340 | // If I want to use a narrower vtx I have here to redefine my | |
1341 | // vtxstatus according to the new vtx range | |
1342 | calib->SetEventRun(event->run); | |
1343 | calib->SetEventMag(event->mag); | |
1344 | calib->SetEventVtxStatus(event->vtxstatus); | |
1345 | calib->SetEventTrigger(event->trig); | |
1346 | ||
1347 | if(!calib->EventAccepted()) // only checks for runs currently | |
1348 | continue; | |
1349 | ||
1350 | calib->FillEventInfo(); | |
1351 | ||
1352 | // The trig==1 is always true for real data, but not for MC data | |
1353 | if(!event->trig) | |
1354 | continue; | |
1355 | ||
1356 | if(event->vtxstatus<1) // only fill tracks for events with vtx inside cuts | |
1357 | continue; | |
1358 | ||
1359 | const Int_t nV0s = v0Array->GetEntries(); | |
1360 | ||
1361 | for(Int_t i = 0; i < nV0s; i++) { | |
1362 | ||
1363 | DeDxV0* v0 = (DeDxV0*)v0Array->At(i); | |
1364 | ||
1365 | // if(v0->ptrack.filter != 2) | |
1366 | // continue; | |
1367 | ||
1368 | if(v0->status != 0) | |
1369 | continue; | |
1370 | ||
1371 | if(v0->dmassG < 0.1) | |
1372 | continue; | |
1373 | ||
1374 | // if(v0->decayr < 5.0) | |
1375 | // continue; | |
1376 | ||
1377 | // if(v0->decayr > 40.0) | |
1378 | // continue; | |
1379 | ||
1380 | const Double_t dmassK = TMath::Abs(v0->dmassK0); | |
1381 | const Double_t dmassL = TMath::Abs(v0->dmassL); | |
1382 | const Double_t dmassAL = TMath::Abs(v0->dmassAL); | |
1383 | ||
1384 | Bool_t fillPos = kFALSE; | |
1385 | Bool_t fillNeg = kFALSE; | |
1386 | ||
1387 | ||
1388 | if(strstr(calib->GetName(), "lambda")) { | |
1389 | ||
1390 | if(dmassK<0.01) | |
1391 | continue; | |
1392 | ||
1393 | if(dmassL<0.01) | |
1394 | fillPos = kTRUE; | |
1395 | ||
1396 | if(dmassAL<0.01) | |
1397 | fillNeg = kTRUE; | |
1398 | } else { // kaons | |
1399 | ||
1400 | if(dmassL<0.01 || dmassAL<0.01) | |
1401 | continue; | |
1402 | ||
1403 | if(dmassK<0.01) { | |
1404 | fillPos = kTRUE; | |
1405 | fillNeg = kTRUE; | |
1406 | } | |
1407 | } | |
1408 | ||
1409 | for(Int_t j = 0; j < 2; j++) { | |
1410 | ||
1411 | DeDxTrack* track = 0; | |
1412 | ||
1413 | if(j==0) { | |
1414 | ||
1415 | if(fillNeg) | |
1416 | track = &(v0->ntrack); | |
1417 | else | |
1418 | continue; | |
1419 | } else { | |
1420 | ||
1421 | if(fillPos) | |
1422 | track = &(v0->ptrack); | |
1423 | else | |
1424 | continue; | |
1425 | } | |
1426 | ||
1427 | calib->SetTrackCharge(track->q); | |
1428 | calib->SetTrackP(track->p); | |
1429 | calib->SetTrackPt(track->pt); | |
1430 | // NB! Filter is not used for these tracks! | |
1431 | calib->SetTrackFilter(track->filter); | |
1432 | calib->SetTrackPhi(track->phi); | |
1433 | calib->SetTrackDeDx(track->dedx); | |
1434 | calib->SetTrackNcl(track->ncl); | |
1435 | calib->SetTrackBeta(track->beta); | |
1436 | calib->SetTrackPidMc(track->pid); | |
1437 | calib->SetTrackEta(track->eta); | |
1438 | ||
1439 | if(calib->TrackAccepted()) { | |
1440 | calib->FillTrackInfo(1.0); | |
1441 | } | |
1442 | } | |
1443 | } | |
1444 | } | |
1445 | } | |
1446 | dataOutFile->cd(); | |
1447 | runList->Write("runList"); | |
1448 | iter->Reset(); | |
1449 | while ((calib = dynamic_cast<AliHighPtDeDxCalib*> (iter->Next()))) { | |
1450 | ||
1451 | calib->Write(); | |
1452 | } | |
1453 | dataOutFile->Close(); | |
1454 | delete dataOutFile; | |
1455 | dataOutFile = 0; | |
1456 | ||
1457 | cout << "Nbad (runno == -1) : " << nBad << endl; | |
1458 | } | |
1459 | ||
1460 | //___________________________________________________________________________ | |
1461 | void AddCalibObjectV0(TList* list, const Char_t* baseName, Bool_t phiCut, | |
1462 | Int_t run, Bool_t analyzeMc) | |
1463 | { | |
1464 | TString objectName(baseName); | |
1465 | if(phiCut) | |
1466 | objectName += "phicut"; | |
1467 | if(run) { | |
1468 | objectName += "_"; | |
1469 | objectName += run; | |
1470 | } | |
1471 | dataOutFile->cd(); | |
1472 | AliHighPtDeDxCalib* data = new AliHighPtDeDxCalib(objectName.Data(), "Calib data"); | |
1473 | if(run) { | |
1474 | data->SetUseRunCut(kTRUE); | |
1475 | data->SetRun(run); | |
1476 | } | |
1477 | ||
1478 | list->Add(data); | |
1479 | data->SetIsMc(analyzeMc); | |
1480 | data->SetUseFilterCut(kFALSE); | |
1481 | data->SetUsePhiCut(phiCut); | |
1482 | if(phiCut) { | |
1483 | data->SetPhiCutLow(data->GetStandardPhiCutLow()); | |
1484 | data->SetPhiCutHigh(data->GetStandardPhiCutHigh()); | |
1485 | } | |
1486 | data->Init(nPtBins, xBins); | |
1487 | data->Print(); | |
1488 | } | |
1489 | ||
1490 | //___________________________________________________________________________ | |
1491 | Int_t CalculateFilter(DeDxTrack* track) | |
1492 | { | |
1493 | Int_t filter = 0; | |
1494 | if(track->filterset3) | |
1495 | filter += 1; | |
1496 | if(track->filterset2 && !track->filterset3) | |
1497 | filter += 2; | |
1498 | if(track->filterset1) | |
1499 | filter += 4; | |
1500 | return filter; | |
1501 | } |