Updates to Trains. create a job-script to help
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / IdentifiedHighPt / macros / run_code.C
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 }