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 }