TPC/AliTPCComparisonMI.C

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16
  17 ///////////////////////////////////////////////////////////////////////////////
  18 //                                                                           //
  19 //  Time Projection Chamber                                                  //
  20 //  Comparison macro for TPC                                                 //
  21 //  responsible:
  22 //  marian.ivanov@cern.ch                                                    //
  23 //
  24 //
  25 // usage:
  26
  27 //
  28 //  TO CREATE OF TREE WITH MC INFO
  29
  30 //  .L AliTPCComparisonMI.C+
  31 //  TPCFindGenTracks *t = new TPCFindGenTracks("galice.root","genTracks.root",1,0)
  32 //  t->Exec()
  33 //
  34 //  TO CREATE COMPARISON TREE
  35 //  TPCCmpTr *t2 = new TPCCmpTr("genTracks.root","cmpTracks.root","galice.root",1,0);
  36 //  t2->Exec()
  37 //
  38 //  EXAMPLE OF COMPARISON VISUALIZATION SESSION
  39 //
  40 // .L AliTPCComparisonMI.C+
  41 // TCut cprim("cprim","MC.fVDist[3]<1");
  42 // TCut cprim("cprim","MC.fVDist[3]<1");
  43 // TCut crec("crec","fReconstructed==1");
  44 // TCut cteta1("cteta1","abs(MC.fTrackRef.Theta()/3.1415-0.5)<0.25");
  45 // TCut cpos1("cpos1","abs(MC.fParticle.fVz/sqrt(MC.fParticle.fVx*MC.fParticle.fVx+MC.fParticle.fVy*MC.fParticle.fVy))<1");
  46 // TCut csens("csens","abs(sqrt(fVDist[0]**2+fVDist[1]**2)-170)<50");
  47 // TCut cmuon("cmuon","abs(MC.fParticle.fPdgCode==-13)");
  48 //
  49 //
  50 // AliTPCComparisonDraw comp;
  51 // comp.SetIO();
  52 // (comp.EffVsPt("MC.fRowsWithDigits>120"+cteta1+cprim,"1"))->Draw()
  53 // (comp.EffVsRM("MC.fRowsWithDigits>20"+cteta1+cprim,"1"))->Draw()
  54 // (comp.EffVsRS("MC.fRowsWithDigits>20"+cteta1+cpos1,"1"))->Draw()
  55 //
  56 // (comp.ResPtvsPt("MC.fRowsWithDigits>20"+cteta1+cpos1,"1",0.15,2.))->Draw()
  57 // (comp.MeanPtvsPt("MC.fRowsWithDigits>20"+cteta1+cpos1,"1",0.15,2.))->Draw()
  58 // (comp.ResdEdxvsN("RC.fReconstructed==1&&MC.fPrim<1.5&&abs(MC.fParticle.fPdgCode)==211&&MC.fParticle.P()>0.35"+cteta1+cpos1+cprim,"1",100,160,4))->Draw()
  59
  60
  61 ///////////////////////////////////////////////////////////////////////////////
  62
  63
  64
  65
  66 #if !defined(__CINT__) || defined(__MAKECINT__)
  67 #include <stdio.h>
  68 #include <string.h>
  69 //ROOT includes
  70 #include "Rtypes.h"
  71 #include "TFile.h"
  72 #include "TTree.h"
  73 #include "TChain.h"
  74 #include "TCut.h"
  75 #include "TString.h"
  76 #include "TBenchmark.h"
  77 #include "TStopwatch.h"
  78 #include "TParticle.h"
  79 #include "TSystem.h"
  80 #include "TTimer.h"
  81 #include "TVector3.h"
  82 #include "TPad.h"
  83 #include "TCanvas.h"
  84 #include "TH1F.h"
  85 #include "TH2F.h"
  86 #include "TF1.h"
  87 #include "TText.h"
  88 #include "Getline.h"
  89 #include "TStyle.h"
  90
  91 //ALIROOT includes
  92 #include "AliRun.h"
  93 #include "AliStack.h"
  94 #include "AliTPCtrack.h"
  95 #include "AliSimDigits.h"
  96 #include "AliTPCParam.h"
  97 #include "AliTPC.h"
  98 #include "AliTPCLoader.h"
  99 #include "AliDetector.h"
 100 #include "AliTrackReference.h"
 101 #include "AliRun.h"
 102 #include "AliTPCParamSR.h"
 103 #include "AliTracker.h"
 104 #include "AliComplexCluster.h"
 105 #include "AliMagF.h"
 106 #endif
 107 #include "AliTPCComparisonMI.h"
 108
 109
 110
 111
 112
 113
 114
 115
 116 //
 117 //
 118
 119 AliTPCParam * GetTPCParam(){
 120   AliTPCParamSR * par = new AliTPCParamSR;
 121   par->Update();
 122   return par;
 123 }
 124
 125
 126 //_____________________________________________________________________________
 127 Float_t TPCBetheBloch(Float_t bg)
 128 {
 129   //
 130   // Bethe-Bloch energy loss formula
 131   //
 132   const Double_t kp1=0.76176e-1;
 133   const Double_t kp2=10.632;
 134   const Double_t kp3=0.13279e-4;
 135   const Double_t kp4=1.8631;
 136   const Double_t kp5=1.9479;
 137
 138   Double_t dbg = (Double_t) bg;
 139
 140   Double_t beta = dbg/TMath::Sqrt(1.+dbg*dbg);
 141
 142   Double_t aa = TMath::Power(beta,kp4);
 143   Double_t bb = TMath::Power(1./dbg,kp5);
 144
 145   bb=TMath::Log(kp3+bb);
 146
 147   return ((Float_t)((kp2-aa-bb)*kp1/aa));
 148 }
 149
 150
 151
 152
 153
 154 ////////////////////////////////////////////////////////////////////////
 155 AliTPCGenInfo::AliTPCGenInfo()
 156 {
 157   fReferences  = 0;
 158   fTRdecay.SetTrack(-1);
 159 }
 160
 161 AliTPCGenInfo::~AliTPCGenInfo()
 162 {
 163   if (fReferences) delete fReferences;
 164   fReferences =0;
 165
 166 }
 167
 168 void AliTPCGenInfo::Update()
 169 {
 170   //
 171   //
 172   fMCtracks =1;
 173   if (!fReferences) return;
 174   Float_t direction=1;
 175   //Float_t rlast=0;
 176   for (Int_t iref =0;iref<fReferences->GetEntriesFast();iref++){
 177     AliTrackReference * ref = (AliTrackReference *) fReferences->At(iref);
 178     //Float_t r = (ref->X()*ref->X()+ref->Y()*ref->Y());
 179     Float_t newdirection = ref->X()*ref->Px()+ref->Y()*ref->Py(); //inside or outside
 180     if (iref==0) direction = newdirection;
 181     if ( newdirection*direction<0){
 182       //changed direction
 183       direction = newdirection;
 184       fMCtracks+=1;
 185     }
 186     //rlast=r;
 187   }
 188 }
 189
 190 ////////////////////////////////////////////////////////////////////////
 191
 192 ////////////////////////////////////////////////////////////////////////
 193 //
 194 // End of implementation of the class AliTPCGenInfo
 195 //
 196 ////////////////////////////////////////////////////////////////////////
 197
 198
 199
 200 ////////////////////////////////////////////////////////////////////////
 201 digitRow::digitRow()
 202 {
 203   Reset();
 204 }
 205 ////////////////////////////////////////////////////////////////////////
 206 digitRow & digitRow::operator=(const digitRow &digOld)
 207 {
 208   for (Int_t i = 0; i<kgRowBytes; i++) fDig[i] = digOld.fDig[i];
 209   return (*this);
 210 }
 211 ////////////////////////////////////////////////////////////////////////
 212 void digitRow::SetRow(Int_t row)
 213 {
 214   if (row >= 8*kgRowBytes) {
 215     cerr<<"digitRow::SetRow: index "<<row<<" out of bounds."<<endl;
 216     return;
 217   }
 218   Int_t iC = row/8;
 219   Int_t iB = row%8;
 220   SETBIT(fDig[iC],iB);
 221 }
 222
 223 ////////////////////////////////////////////////////////////////////////
 224 Bool_t digitRow::TestRow(Int_t row)
 225 {
 226 //
 227 // return kTRUE if row is on
 228 //
 229   Int_t iC = row/8;
 230   Int_t iB = row%8;
 231   return TESTBIT(fDig[iC],iB);
 232 }
 233 ////////////////////////////////////////////////////////////////////////
 234 Int_t digitRow::RowsOn(Int_t upto)
 235 {
 236 //
 237 // returns number of rows with a digit
 238 // count only rows less equal row number upto
 239 //
 240   Int_t total = 0;
 241   for (Int_t i = 0; i<kgRowBytes; i++) {
 242     for (Int_t j = 0; j < 8; j++) {
 243       if (i*8+j > upto) return total;
 244       if (TESTBIT(fDig[i],j))  total++;
 245     }
 246   }
 247   return total;
 248 }
 249 ////////////////////////////////////////////////////////////////////////
 250 void digitRow::Reset()
 251 {
 252 //
 253 // resets all rows to zero
 254 //
 255   for (Int_t i = 0; i<kgRowBytes; i++) {
 256     fDig[i] <<= 8;
 257   }
 258 }
 259 ////////////////////////////////////////////////////////////////////////
 260 Int_t digitRow::Last()
 261 {
 262 //
 263 // returns the last row number with a digit
 264 // returns -1 if now digits
 265 //
 266   for (Int_t i = kgRowBytes-1; i>=0; i--) {
 267     for (Int_t j = 7; j >= 0; j--) {
 268       if TESTBIT(fDig[i],j) return i*8+j;
 269     }
 270   }
 271   return -1;
 272 }
 273 ////////////////////////////////////////////////////////////////////////
 274 Int_t digitRow::First()
 275 {
 276 //
 277 // returns the first row number with a digit
 278 // returns -1 if now digits
 279 //
 280   for (Int_t i = 0; i<kgRowBytes; i++) {
 281     for (Int_t j = 0; j < 8; j++) {
 282       if (TESTBIT(fDig[i],j)) return i*8+j;
 283     }
 284   }
 285   return -1;
 286 }
 287
 288 ////////////////////////////////////////////////////////////////////////
 289 //
 290 // end of implementation of a class digitRow
 291 //
 292 ////////////////////////////////////////////////////////////////////////
 293
 294 ////////////////////////////////////////////////////////////////////////
 295 TPCFindGenTracks::TPCFindGenTracks()
 296 {
 297   fMCInfo = new AliTPCGenInfo;
 298   fMCInfo->fReferences = new TClonesArray("AliTrackReference");
 299
 300   Reset();
 301 }
 302
 303 ////////////////////////////////////////////////////////////////////////
 304 TPCFindGenTracks::TPCFindGenTracks(const char * fnGalice, const char* fnRes,
 305                                    Int_t nEvents, Int_t firstEvent)
 306 {
 307   Reset();
 308   fFirstEventNr = firstEvent;
 309   fEventNr = firstEvent;
 310   fNEvents = nEvents;
 311   //   fFnRes = fnRes;
 312   sprintf(fFnRes,"%s",fnRes);
 313   //
 314   fLoader = AliRunLoader::Open(fnGalice);
 315   if (gAlice){
 316     delete gAlice->GetRunLoader();
 317     delete gAlice;
 318     gAlice = 0x0;
 319   }
 320   if (fLoader->LoadgAlice()){
 321     cerr<<"Error occured while l"<<endl;
 322   }
 323   Int_t nall = fLoader->GetNumberOfEvents();
 324   if (nEvents==0) {
 325     nEvents =nall;
 326     fNEvents=nall;
 327     fFirstEventNr=0;
 328   }
 329
 330   if (nall<=0){
 331     cerr<<"no events available"<<endl;
 332     fEventNr = 0;
 333     return;
 334   }
 335   if (firstEvent+nEvents>nall) {
 336     fEventNr = nall-firstEvent;
 337     cerr<<"restricted number of events availaible"<<endl;
 338   }
 339   AliMagF * magf = gAlice->Field();
 340   AliTracker::SetFieldMap(magf);
 341 }
 342
 343 ////////////////////////////////////////////////////////////////////////
 344 void TPCFindGenTracks::Reset()
 345 {
 346   fEventNr = 0;
 347   fNEvents = 0;
 348   fTreeGenTracks = 0;
 349   fFileGenTracks = 0;
 350   fContainerDigitRow = 0;
 351   //
 352   fReferences = 0;
 353   fReferenceIndex0 = 0;
 354   fReferenceIndex1 = 0;
 355   fDecayRef   = 0;
 356   //
 357   fDebug = 0;
 358   fVPrim[0] = -1000.;
 359   fVPrim[1] = -1000.;
 360   fVPrim[2] = -1000.;
 361   fParamTPC = 0;
 362 }
 363 ////////////////////////////////////////////////////////////////////////
 364 TPCFindGenTracks::~TPCFindGenTracks()
 365 {
 366
 367   if (fLoader){
 368     fLoader->UnloadgAlice();
 369     gAlice = 0;
 370     delete fLoader;
 371   }
 372 }
 373
 374
 375 Int_t  TPCFindGenTracks::SetIO()
 376 {
 377   //
 378   //
 379   CreateTreeGenTracks();
 380   if (!fTreeGenTracks) return 1;
 381   //  AliTracker::SetFieldFactor();
 382
 383   fParamTPC = GetTPCParam();
 384   //
 385   return 0;
 386 }
 387
 388 ////////////////////////////////////////////////////////////////////////
 389 Int_t TPCFindGenTracks::SetIO(Int_t eventNr)
 390 {
 391   //
 392   //
 393   // SET INPUT
 394   fLoader->SetEventNumber(eventNr);
 395   //
 396   fLoader->LoadHeader();
 397   fLoader->LoadKinematics();
 398   fStack = fLoader->Stack();
 399   //
 400   fLoader->LoadTrackRefs();
 401   fTreeTR = fLoader->TreeTR();
 402   //
 403   AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
 404   tpcl->LoadDigits();
 405   fTreeD = tpcl->TreeD();
 406   return 0;
 407 }
 408
 409 Int_t TPCFindGenTracks::CloseIOEvent()
 410 {
 411   fLoader->UnloadHeader();
 412   fLoader->UnloadKinematics();
 413   fLoader->UnloadTrackRefs();
 414   AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
 415   tpcl->UnloadDigits();
 416   return 0;
 417 }
 418
 419 Int_t TPCFindGenTracks::CloseIO()
 420 {
 421   fLoader->UnloadgAlice();
 422   return 0;
 423 }
 424
 425
 426
 427 ////////////////////////////////////////////////////////////////////////
 428 Int_t TPCFindGenTracks::Exec(Int_t nEvents, Int_t firstEventNr)
 429 {
 430   fNEvents = nEvents;
 431   fFirstEventNr = firstEventNr;
 432   return Exec();
 433 }
 434
 435 ////////////////////////////////////////////////////////////////////////
 436 Int_t TPCFindGenTracks::Exec()
 437 {
 438   TStopwatch timer;
 439   timer.Start();
 440   Int_t status =SetIO();
 441   if (status>0) return status;
 442   //
 443
 444   for (fEventNr = fFirstEventNr; fEventNr < fFirstEventNr+fNEvents;
 445        fEventNr++) {
 446     SetIO(fEventNr);
 447     fNParticles = fStack->GetNtrack();
 448     fContainerDigitRow = new digitRow[fNParticles];
 449     //
 450     fReferences      = new AliTrackReference[fgMaxTR];
 451     fReferenceIndex0 = new Int_t[fNParticles];
 452     fReferenceIndex1 = new Int_t[fNParticles];
 453     fDecayRef        = new AliTrackReference[fNParticles];
 454
 455     for (Int_t i = 0; i<fNParticles; i++) {
 456       fReferenceIndex0[i] = -1;
 457       fReferenceIndex1[i] = -1;
 458     }
 459
 460     cout<<"Start to process event "<<fEventNr<<endl;
 461     cout<<"\tfNParticles = "<<fNParticles<<endl;
 462     if (fDebug>2) cout<<"\tStart loop over TreeD"<<endl;
 463     if (TreeDLoop()>0) return 1;
 464     if (fDebug>2) cout<<"\tStart loop over TreeTR"<<endl;
 465     if (TreeTRLoop()>0) return 1;
 466     if (fDebug>2) cout<<"\tStart loop over TreeK"<<endl;
 467     if (TreeKLoop()>0) return 1;
 468     if (fDebug>2) cout<<"\tEnd loop over TreeK"<<endl;
 469
 470     delete [] fContainerDigitRow;
 471     delete [] fReferences;
 472     delete [] fReferenceIndex0;
 473     delete [] fReferenceIndex1;
 474     delete [] fDecayRef;
 475     CloseIOEvent();
 476   }
 477   //
 478   CloseIO();
 479   CloseOutputFile();
 480
 481   cerr<<"Exec finished"<<endl;
 482
 483   timer.Stop();
 484   timer.Print();
 485   return 0;
 486 }
 487 ////////////////////////////////////////////////////////////////////////
 488 void TPCFindGenTracks::CreateTreeGenTracks()
 489 {
 490   fFileGenTracks = TFile::Open(fFnRes,"RECREATE");
 491   if (!fFileGenTracks) {
 492     cerr<<"Error in CreateTreeGenTracks: cannot open file "<<fFnRes<<endl;
 493     return;
 494   }
 495   fTreeGenTracks = new TTree("genTracksTree","genTracksTree");
 496
 497
 498
 499   fMCInfo = new AliTPCGenInfo;
 500   fMCInfo->fReferences = new TClonesArray("AliTrackReference");
 501
 502   fTreeGenTracks->Branch("MC","AliTPCGenInfo",&fMCInfo);
 503
 504   fTreeGenTracks->Branch("fEventNr",&fEventNr,"fEventNr/I");
 505
 506   fTreeGenTracks->AutoSave();
 507 }
 508 ////////////////////////////////////////////////////////////////////////
 509 void TPCFindGenTracks::CloseOutputFile()
 510 {
 511   if (!fFileGenTracks) {
 512     cerr<<"File "<<fFnRes<<" not found as an open file."<<endl;
 513     return;
 514   }
 515   fFileGenTracks->cd();
 516   fTreeGenTracks->Write();
 517   delete fTreeGenTracks;
 518   fFileGenTracks->Close();
 519   delete fFileGenTracks;
 520   return;
 521 }
 522
 523 ////////////////////////////////////////////////////////////////////////
 524 Int_t TPCFindGenTracks::TreeKLoop()
 525 {
 526 //
 527 // open the file with treeK
 528 // loop over all entries there and save information about some tracks
 529 //
 530
 531   AliStack * stack = fStack;
 532   if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
 533
 534   if (fDebug > 0) {
 535     cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
 536         <<fEventNr<<endl;
 537   }
 538
 539   Int_t  ipdg = 0;
 540   TParticlePDG *ppdg = 0;
 541
 542
 543   // not all generators give primary vertex position. Take the vertex
 544   // of the particle 0 as primary vertex.
 545   TDatabasePDG  pdg; //get pdg table
 546
 547   //thank you very much root for this
 548   TBranch * br = fTreeGenTracks->GetBranch("MC");
 549
 550   TParticle *particle = stack->ParticleFromTreeK(0);
 551   fVPrim[0] = particle->Vx();
 552   fVPrim[1] = particle->Vy();
 553   fVPrim[2] = particle->Vz();
 554   for (Int_t iParticle = 0; iParticle < fNParticles; iParticle++) {
 555
 556 // load only particles with TR
 557     if (fReferenceIndex0[iParticle]<0) continue;
 558     //////////////////////////////////////////////////////////////////////
 559     fMCInfo->fLabel = iParticle;
 560
 561     fMCInfo->fRow = (fContainerDigitRow[iParticle]);
 562     fMCInfo->fRowsWithDigits = fMCInfo->fRow.RowsOn();
 563     if (fMCInfo->fRowsWithDigits<10) continue;
 564     fMCInfo->fRowsWithDigitsInn = fMCInfo->fRow.RowsOn(63); // 63 = number of inner rows
 565     fMCInfo->fRowsTrackLength = fMCInfo->fRow.Last() - fMCInfo->fRow.First();
 566     fMCInfo->fDigitsInSeed = 0;
 567     if (fMCInfo->fRow.TestRow(seedRow11) && fMCInfo->fRow.TestRow(seedRow12))
 568       fMCInfo->fDigitsInSeed = 1;
 569     if (fMCInfo->fRow.TestRow(seedRow21) && fMCInfo->fRow.TestRow(seedRow22))
 570       fMCInfo->fDigitsInSeed += 10;
 571     //
 572     //
 573     //
 574     fMCInfo->fParticle = *(stack->Particle(iParticle));
 575     //
 576     //
 577     fMCInfo->fLabel = iParticle;
 578     fMCInfo->fVDist[0] = fMCInfo->fParticle.Vx()-fVPrim[0];
 579     fMCInfo->fVDist[1] = fMCInfo->fParticle.Vy()-fVPrim[1];
 580     fMCInfo->fVDist[2] = fMCInfo->fParticle.Vz()-fVPrim[2];
 581     fMCInfo->fVDist[3] = TMath::Sqrt(fMCInfo->fVDist[0]*fMCInfo->fVDist[0]+
 582                                       fMCInfo->fVDist[1]*fMCInfo->fVDist[1]+fMCInfo->fVDist[2]*fMCInfo->fVDist[2]);
 583
 584     //
 585     Int_t index = fReferenceIndex0[iParticle];
 586     AliTrackReference  ref  = fReferences[index];
 587     // if (ref.GetTrack()!=iParticle)
 588     //  printf("problem2\n");
 589     fMCInfo->fTrackRef = ref;
 590
 591     Int_t rfindex =0;
 592     if (fMCInfo->fReferences !=0) delete fMCInfo->fReferences;
 593     fMCInfo->fReferences =
 594       new TClonesArray("AliTrackReference");
 595     fMCInfo->fReferences->ExpandCreateFast(fReferenceIndex1[iParticle]-fReferenceIndex0[iParticle]+1);
 596     //
 597     for (Int_t i = fReferenceIndex0[iParticle];i<=fReferenceIndex1[iParticle];i++){
 598       AliTrackReference  ref  = fReferences[i];
 599       AliTrackReference *ref2 = (AliTrackReference*) fMCInfo->fReferences->At(rfindex);
 600       if (ref.GetTrack()!=iParticle){
 601         //printf("problem5\n");
 602         continue;
 603       }
 604       *ref2 = ref;
 605       rfindex++;
 606     }
 607     //
 608     //
 609     ipdg = fMCInfo->fParticle.GetPdgCode();
 610     fMCInfo->fPdg = ipdg;
 611     ppdg = pdg.GetParticle(ipdg);
 612     // calculate primary ionization per cm
 613     if (ppdg){
 614       Float_t mass = ppdg->Mass();
 615       Float_t p = TMath::Sqrt(fMCInfo->fTrackRef.Px()*fMCInfo->fTrackRef.Px()+
 616                               fMCInfo->fTrackRef.Py()*fMCInfo->fTrackRef.Py()+
 617                               fMCInfo->fTrackRef.Pz()*fMCInfo->fTrackRef.Pz());
 618
 619       //      Float_t betagama = fMCInfo->fParticle.P()/mass;
 620       Float_t betagama = p /mass;
 621       fMCInfo->fPrim = TPCBetheBloch(betagama);
 622     }
 623     fMCInfo->fTPCdecay=kFALSE;
 624     if (fDecayRef[iParticle].GetTrack()>0){
 625       fMCInfo->fTRdecay  = fDecayRef[iParticle];
 626       fMCInfo->fDecayCoord[0] = fMCInfo->fTRdecay.X();
 627       fMCInfo->fDecayCoord[1] = fMCInfo->fTRdecay.Y();
 628       fMCInfo->fDecayCoord[2] = fMCInfo->fTRdecay.Z();
 629       if ( (fMCInfo->fTRdecay.R()<250)&&(fMCInfo->fTRdecay.R()>85) && (TMath::Abs(fMCInfo->fTRdecay.Z())<250) ){
 630         fMCInfo->fTPCdecay=kTRUE;
 631       }
 632     }
 633     else{
 634       fMCInfo->fTRdecay.SetTrack(-1);
 635       fMCInfo->fDecayCoord[0] = 0;
 636       fMCInfo->fDecayCoord[1] = 0;
 637       fMCInfo->fDecayCoord[2] = 0;
 638     }
 639     fMCInfo->Update();
 640     //////////////////////////////////////////////////////////////////////
 641
 642     br->SetAddress(&fMCInfo);
 643     fTreeGenTracks->Fill();
 644
 645   }
 646   fTreeGenTracks->AutoSave();
 647
 648   if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl;
 649
 650   return 0;
 651 }
 652
 653
 654
 655
 656 ////////////////////////////////////////////////////////////////////////
 657 Int_t TPCFindGenTracks::TreeDLoop()
 658 {
 659 //
 660 // open the file with treeD
 661 // loop over all entries there and save information about some tracks
 662 //
 663
 664   Int_t nInnerSector = fParamTPC->GetNInnerSector();
 665   Int_t rowShift = 0;
 666   //  Int_t zero=fParamTPC->GetZeroSup();
 667   Int_t zero=fParamTPC->GetZeroSup()+6;
 668
 669
 670   //char treeDName[100];
 671   //sprintf(treeDName,"TreeD_75x40_100x60_150x60_%d",fEventNr);
 672   //TTree *treeD=(TTree*)fFileTreeD->Get(treeDName);
 673   //
 674   //
 675   AliSimDigits digitsAddress, *digits=&digitsAddress;
 676   fTreeD->GetBranch("Segment")->SetAddress(&digits);
 677
 678   Int_t sectorsByRows=(Int_t)fTreeD->GetEntries();
 679   if (fDebug > 1) cout<<"\tsectorsByRows = "<<sectorsByRows<<endl;
 680   for (Int_t i=0; i<sectorsByRows; i++) {
 681 //  for (Int_t i=5720; i<sectorsByRows; i++) {
 682     if (!fTreeD->GetEvent(i)) continue;
 683     Int_t sec,row;
 684     fParamTPC->AdjustSectorRow(digits->GetID(),sec,row);
 685     if (fDebug > 1) cout<<sec<<' '<<row<<"                          \r";
 686 //    cerr<<sec<<' '<<row<<endl;
 687
 688 // here I expect that upper sectors follow lower sectors
 689     if (sec > nInnerSector) rowShift = fParamTPC->GetNRowLow();
 690     //
 691     //digits->ExpandTrackBuffer();
 692     //Int_t *tracks = digits->GetTracks();
 693     digits->First();
 694     do {
 695       Int_t iRow=digits->CurrentRow();
 696       Int_t iColumn=digits->CurrentColumn();
 697       Short_t digitValue = digits->CurrentDigit();
 698 //      cout<<"Inner loop: sector, iRow, iColumn "
 699 //        <<sec<<" "<<iRow<<" "<<iColumn<<endl;
 700       if (digitValue >= zero) {
 701         Int_t label;
 702         for (Int_t j = 0; j<3; j++) {
 703           label = digits->GetTrackID(iRow,iColumn,j);
 704           //label = digits->GetTrackIDFast(iRow,iColumn,j);
 705
 706           if (label >= fNParticles) { //don't label from bakground event
 707             continue;
 708           }
 709           if (label >= 0 && label <= fNParticles) {
 710 //        if (label >= 0 && label <= fDebug) {
 711             if (fDebug > 6 ) {
 712               cout<<"Inner loop: sector, iRow, iColumn, label, value, row "
 713                   <<sec<<" "
 714                   <<iRow<<" "<<iColumn<<" "<<label<<" "<<digitValue
 715                   <<" "<<row<<endl;
 716             }
 717             fContainerDigitRow[label].SetRow(row+rowShift);
 718           }
 719         }
 720       }
 721     } while (digits->Next());
 722   }
 723
 724   if (fDebug > 2) cerr<<"end of TreeDLoop"<<endl;
 725
 726   return 0;
 727 }
 728
 729
 730 ////////////////////////////////////////////////////////////////////////
 731 Int_t TPCFindGenTracks::TreeTRLoop()
 732 {
 733   //
 734   // loop over TrackReferences and store the first one for each track
 735   //
 736
 737   TTree * treeTR = fTreeTR;
 738   if (!treeTR) {
 739     cerr<<"TreeTR not found"<<endl;
 740     return 1;
 741   }
 742   Int_t nPrimaries = (Int_t) treeTR->GetEntries();
 743   if (fDebug > 1) cout<<"There are "<<nPrimaries<<" entries in TreeTR"<<endl;
 744   //
 745   //
 746   //track references for TPC
 747   TBranch *TPCBranchTR  = treeTR->GetBranch("TPC");
 748   if (!TPCBranchTR) {
 749     cerr<<"TPC branch in TR not found"<<endl;
 750     return 1;
 751   }
 752   TClonesArray* TPCArrayTR = new TClonesArray("AliTrackReference");
 753   TPCBranchTR->SetAddress(&TPCArrayTR);
 754   //get decay point if exist
 755   TBranch *runbranch  = treeTR->GetBranch("AliRun");
 756   if (!runbranch) {
 757     cerr<<"Run branch in TR not found"<<endl;
 758     return 1;
 759   }
 760   TClonesArray* RunArrayTR = new TClonesArray("AliTrackReference");
 761   runbranch->SetAddress(&RunArrayTR);
 762   //
 763   //
 764   //
 765   Int_t index     =  0;
 766   for (Int_t iPrimPart = 0; iPrimPart<nPrimaries; iPrimPart++) {
 767     TPCBranchTR->GetEntry(iPrimPart);
 768     Float_t ptstart = 0;
 769     for (Int_t iTrackRef = 0; iTrackRef < TPCArrayTR->GetEntriesFast(); iTrackRef++) {
 770       AliTrackReference *trackRef = (AliTrackReference*)TPCArrayTR->At(iTrackRef);
 771       //
 772       Int_t label = trackRef->GetTrack();
 773       if (label<0  || label > fNParticles) {
 774         continue;
 775       }
 776       if (fReferenceIndex0[label]<0) ptstart = trackRef->Pt();  //store pt at the TPC entrance
 777       if (ptstart<fgPtCut) continue;
 778
 779       if (index>=fgMaxTR) continue;     //restricted size of buffer for TR
 780       fReferences[index] = *trackRef;
 781       fReferenceIndex1[label] = index;  // the last ref with given label
 782       if (fReferenceIndex0[label]==-1) fReferenceIndex0[label] = index;   //the first index with label
 783       index++;
 784     }
 785     // get dacay position
 786     runbranch->GetEntry(iPrimPart);
 787     for (Int_t iTrackRef = 0; iTrackRef < RunArrayTR->GetEntriesFast(); iTrackRef++) {
 788       AliTrackReference *trackRef = (AliTrackReference*)RunArrayTR->At(iTrackRef);
 789       //
 790       if (trackRef->Pt() < fgPtCut) continue;
 791       Int_t label = trackRef->GetTrack();
 792       if (label<0  || label > fNParticles) {
 793         continue;
 794       }
 795       if (trackRef->R()>450) continue;   //not decay  in TPC
 796       if (trackRef->Z()>450) continue;   //not decay  in TPC
 797       if (!trackRef->TestBit(BIT(2))) continue;  //if not decay
 798
 799       if (label>=fgMaxTR) continue;     //restricted size of buffer for TR
 800       fDecayRef[label] = *trackRef;
 801     }
 802
 803   }
 804   TPCArrayTR->Delete();
 805   delete  TPCArrayTR;
 806   RunArrayTR->Delete();
 807   delete  RunArrayTR;
 808
 809   return 0;
 810 }
 811
 812 ////////////////////////////////////////////////////////////////////////
 813 Float_t TPCFindGenTracks::TR2LocalX(AliTrackReference *trackRef,
 814                                     AliTPCParam *paramTPC) {
 815
 816   Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
 817   Int_t index[4];
 818   paramTPC->Transform0to1(x,index);
 819   paramTPC->Transform1to2(x,index);
 820   return x[0];
 821 }
 822 ////////////////////////////////////////////////////////////////////////
 823
 824
 825
 826
 827 ////////////////////////////////////////////////////////////////////////
 828 TPCCmpTr::TPCCmpTr()
 829 {
 830   Reset();
 831 }
 832
 833 ////////////////////////////////////////////////////////////////////////
 834 TPCCmpTr::TPCCmpTr(const char* fnGenTracks,
 835                    const char* fnCmp,
 836                    const char* fnGalice,
 837                    Int_t nEvents, Int_t firstEvent)
 838 {
 839   Reset();
 840   //  fFnGenTracks = fnGenTracks;
 841   //  fFnCmp = fnCmp;
 842   sprintf(fFnGenTracks,"%s",fnGenTracks);
 843   sprintf(fFnCmp,"%s",fnCmp);
 844
 845   fFirstEventNr = firstEvent;
 846   fEventNr = firstEvent;
 847   fNEvents = nEvents;
 848
 849   //
 850   fLoader = AliRunLoader::Open(fnGalice);
 851   if (gAlice){
 852     //delete gAlice->GetRunLoader();
 853     delete gAlice;
 854     gAlice = 0x0;
 855   }
 856   if (fLoader->LoadgAlice()){
 857     cerr<<"Error occured while l"<<endl;
 858   }
 859   Int_t nall = fLoader->GetNumberOfEvents();
 860   if (nEvents==0) {
 861     nEvents =nall;
 862     fNEvents=nall;
 863     fFirstEventNr=0;
 864   }
 865
 866   if (nall<=0){
 867     cerr<<"no events available"<<endl;
 868     fEventNr = 0;
 869     return;
 870   }
 871   if (firstEvent+nEvents>nall) {
 872     fEventNr = nall-firstEvent;
 873     cerr<<"restricted number of events availaible"<<endl;
 874   }
 875   AliMagF * magf = gAlice->Field();
 876   AliTracker::SetFieldMap(magf);
 877
 878 }
 879
 880
 881 ////////////////////////////////////////////////////////////////////////
 882 TPCCmpTr::~TPCCmpTr()
 883 {
 884   if (fLoader) {
 885     delete fLoader;
 886   }
 887 }
 888
 889 //////////////////////////////////////////////////////////////
 890 Int_t TPCCmpTr::SetIO()
 891 {
 892   //
 893   //
 894   CreateTreeCmp();
 895   if (!fTreeCmp) return 1;
 896   fParamTPC = GetTPCParam();
 897   //
 898   if (!ConnectGenTree()) {
 899     cerr<<"Cannot connect tree with generated tracks"<<endl;
 900     return 1;
 901   }
 902   return 0;
 903 }
 904
 905 //////////////////////////////////////////////////////////////
 906
 907 Int_t TPCCmpTr::SetIO(Int_t eventNr)
 908 {
 909   //
 910   //
 911   // SET INPUT
 912   //gAlice->GetEvent(eventNr);
 913   fLoader->SetEventNumber(eventNr);
 914   //
 915   AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
 916   tpcl->LoadTracks();
 917   fTreeRecTracks = tpcl->TreeT();
 918   //
 919   return 0;
 920 }
 921
 922
 923
 924 ////////////////////////////////////////////////////////////////////////
 925 void TPCCmpTr::Reset()
 926 {
 927   fEventNr = 0;
 928   fNEvents = 0;
 929   fTreeCmp = 0;
 930   //  fFnCmp = "cmpTracks.root";
 931   fFileGenTracks = 0;
 932   fDebug = 0;
 933   //
 934   fParamTPC = 0;
 935   fTreeRecTracks = 0;
 936   fTreePoints =0;
 937
 938   fTPCTrack = 0;
 939   fTracks   = 0;
 940   fTrackPoints =0;
 941 }
 942
 943 ////////////////////////////////////////////////////////////////////////
 944 Int_t TPCCmpTr::Exec(Int_t nEvents, Int_t firstEventNr)
 945 {
 946   fNEvents = nEvents;
 947   fFirstEventNr = firstEventNr;
 948   return Exec();
 949 }
 950
 951 ////////////////////////////////////////////////////////////////////////
 952 Int_t TPCCmpTr::Exec()
 953 {
 954   TStopwatch timer;
 955   timer.Start();
 956
 957   if (SetIO()==1)
 958     return 1;
 959
 960   fNextTreeGenEntryToRead = 0;
 961   cerr<<"fFirstEventNr, fNEvents: "<<fFirstEventNr<<" "<<fNEvents<<endl;
 962   for (Int_t eventNr = fFirstEventNr; eventNr < fFirstEventNr+fNEvents;
 963        eventNr++) {
 964     SetIO(fEventNr);
 965     fNParticles = gAlice->GetEvent(fEventNr);
 966
 967     fIndexRecTracks = new Int_t[fNParticles*4];  //write at maximum 4 tracks corresponding to particle
 968     fFakeRecTracks = new Int_t[fNParticles];
 969     fMultiRecTracks = new Int_t[fNParticles];
 970     for (Int_t i = 0; i<fNParticles; i++) {
 971       fIndexRecTracks[4*i] = -1;
 972       fIndexRecTracks[4*i+1] = -1;
 973       fIndexRecTracks[4*i+2] = -1;
 974       fIndexRecTracks[4*i+3] = -1;
 975
 976       fFakeRecTracks[i] = 0;
 977       fMultiRecTracks[i] = 0;
 978     }
 979
 980     cout<<"Start to process event "<<fEventNr<<endl;
 981     cout<<"\tfNParticles = "<<fNParticles<<endl;
 982     if (fDebug>2) cout<<"\tStart loop over TreeT"<<endl;
 983     if (TreeTLoop(eventNr)>0) return 1;
 984
 985     if (fDebug>2) cout<<"\tStart loop over tree genTracks"<<endl;
 986     if (TreeGenLoop(eventNr)>0) return 1;
 987     if (fDebug>2) cout<<"\tEnd loop over tree genTracks"<<endl;
 988
 989     delete fTreeRecTracks;
 990     delete [] fIndexRecTracks;
 991     delete [] fFakeRecTracks;
 992     delete [] fMultiRecTracks;
 993   }
 994
 995   CloseOutputFile();
 996
 997   cerr<<"Exec finished"<<endl;
 998   timer.Stop();
 999   timer.Print();
1000   return 0;
1001
1002 }
1003 ////////////////////////////////////////////////////////////////////////
1004 Bool_t TPCCmpTr::ConnectGenTree()
1005 {
1006 //
1007 // connect all branches from the genTracksTree
1008 // use the same variables as for the new cmp tree, it may work
1009 //
1010   fFileGenTracks = TFile::Open(fFnGenTracks,"READ");
1011   if (!fFileGenTracks) {
1012     cerr<<"Error in ConnectGenTree: cannot open file "<<fFnGenTracks<<endl;
1013     return kFALSE;
1014   }
1015   fTreeGenTracks = (TTree*)fFileGenTracks->Get("genTracksTree");
1016   if (!fTreeGenTracks) {
1017     cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file "
1018         <<fFnGenTracks<<endl;
1019     return kFALSE;
1020   }
1021   //
1022   fMCInfo = new AliTPCGenInfo;
1023   fMCInfo->fReferences = new TClonesArray("AliTrackReference");
1024   fTreeGenTracks->SetBranchAddress("MC",&fMCInfo);
1025
1026   //
1027   //fTreeGenTracks->SetBranchAddress("fEventNr",&fEventNr);
1028
1029
1030   if (fDebug > 1) {
1031     cout<<"Number of gen. tracks with TR: "<<fTreeGenTracks->GetEntries()<<endl;
1032   }
1033   return kTRUE;
1034 }
1035
1036
1037 ////////////////////////////////////////////////////////////////////////
1038 void TPCCmpTr::CreateTreeCmp()
1039 {
1040   fFileCmp = TFile::Open(fFnCmp,"RECREATE");
1041   if (!fFileCmp) {
1042     cerr<<"Error in CreateTreeCmp: cannot open file "<<fFnCmp<<endl;
1043     return;
1044   }
1045
1046
1047   fTreeCmp    = new TTree("TPCcmpTracks","TPCcmpTracks");
1048   //
1049   fMCInfo = new AliTPCGenInfo;
1050   fMCInfo->fReferences = new TClonesArray("AliTrackReference");
1051   fRecInfo = new AliTPCRecInfo;
1052   //
1053   fTPCTrack = new AliTPCtrack;
1054    //
1055   fTreeCmp->Branch("MC","AliTPCGenInfo",&fMCInfo);
1056   fTreeCmp->Branch("RC","AliTPCRecInfo",&fRecInfo);
1057   fTreeCmp->Branch("fEventNr",&fEventNr,"fEventNr/I");
1058   fTreeCmp->Branch("fTPCTrack","AliTPCtrack",&fTPCTrack);
1059   //
1060   fTreeCmp->AutoSave();
1061
1062 }
1063 ////////////////////////////////////////////////////////////////////////
1064 void TPCCmpTr::CloseOutputFile()
1065 {
1066   if (!fFileCmp) {
1067     cerr<<"File "<<fFnCmp<<" not found as an open file."<<endl;
1068     return;
1069   }
1070   fFileCmp->cd();
1071   fTreeCmp->Write();
1072   delete fTreeCmp;
1073
1074   fFileCmp->Close();
1075   delete fFileCmp;
1076   return;
1077 }
1078 ////////////////////////////////////////////////////////////////////////
1079
1080 TVector3 TPCCmpTr::TR2Local(AliTrackReference *trackRef,
1081                             AliTPCParam *paramTPC) {
1082
1083   Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
1084   Int_t index[4];
1085   paramTPC->Transform0to1(x,index);
1086   paramTPC->Transform1to2(x,index);
1087   return TVector3(x);
1088 }
1089 ////////////////////////////////////////////////////////////////////////
1090
1091 Int_t TPCCmpTr::TreeTLoop(Int_t eventNr)
1092 {
1093   //
1094   // loop over all TPC reconstructed tracks and store info in memory
1095   //
1096   TStopwatch  timer;
1097   timer.Start();
1098
1099   if (!fTreeRecTracks) {
1100     cerr<<"Can't get a tree with TPC rec. tracks  "<<endl;
1101     return 1;
1102   }
1103   //fTreePoints=(TTree*)fFileRecTracks->Get("trackDebug");
1104
1105   Int_t nEntries = (Int_t) fTreeRecTracks->GetEntries();
1106   if (fDebug > 2) cout<<"Event, rec. tracks: "<<eventNr<<" "
1107                       <<nEntries<<endl;
1108   TBranch * br= fTreeRecTracks->GetBranch("tracks");
1109   br->SetAddress(&fTPCTrack);
1110   TBranch *brp = 0;
1111   if (fTreePoints) brp = fTreePoints->GetBranch("debug");
1112
1113   if (fTracks){
1114     fTracks->Delete();
1115     delete fTracks;
1116   }
1117   if (fTrackPoints){
1118     fTrackPoints->Delete();
1119     delete fTrackPoints;
1120     fTrackPoints =0;
1121   }
1122   fTracks      = new TObjArray(nEntries);
1123   if (brp){
1124     fTrackPoints = new TObjArray(nEntries);
1125   }
1126   else fTrackPoints = 0;
1127
1128   //
1129   //load tracks to the memory
1130   for (Int_t i=0; i<nEntries;i++){
1131     AliTPCtrack * track = new AliTPCtrack;
1132     br->SetAddress(&track);
1133     br->GetEntry(i);
1134     fTracks->AddAt(track,i);
1135   }
1136   //
1137   //load track points to the memory
1138   if (brp) for (Int_t i=0; i<nEntries;i++){
1139     TClonesArray * arr = new TClonesArray("AliTPCTrackPoint2");
1140     brp->SetAddress(&arr);
1141     brp->GetEntry(i);
1142     if (arr!=0)
1143       for (Int_t j=0;j<arr->GetEntriesFast();j++){
1144         AliTPCTrackPoint2 * point = (AliTPCTrackPoint2*)arr->UncheckedAt(j);
1145         if (point && point->fID>=0){
1146           fTrackPoints->AddAt(arr,point->fID);
1147           break;
1148         }
1149       }
1150   }
1151   //
1152
1153   //
1154   for (Int_t iEntry=0; iEntry<nEntries;iEntry++){
1155     //br->GetEntry(iEntry);
1156     fTPCTrack = (AliTPCtrack*)fTracks->UncheckedAt(iEntry);
1157     //
1158     Int_t label = fTPCTrack->GetLabel();
1159     Int_t absLabel = abs(label);
1160     if (absLabel < fNParticles) {
1161       //      fIndexRecTracks[absLabel] =  iEntry;
1162       if (label < 0) fFakeRecTracks[absLabel]++;
1163
1164       if (fMultiRecTracks[absLabel]>0){
1165         if (fMultiRecTracks[absLabel]<4)
1166           fIndexRecTracks[absLabel*4+fMultiRecTracks[absLabel]] =  iEntry;
1167       }
1168       else
1169         fIndexRecTracks[absLabel*4] =  iEntry;
1170       fMultiRecTracks[absLabel]++;
1171     }
1172   }
1173   printf("Time spended in TreeTLoop\n");
1174   timer.Print();
1175
1176   if (fDebug > 2) cerr<<"end of TreeTLoop"<<endl;
1177
1178   return 0;
1179 }
1180 ////////////////////////////////////////////////////////////////////////
1181 Int_t TPCCmpTr::TreeGenLoop(Int_t eventNr)
1182 {
1183 //
1184 // loop over all entries for a given event, find corresponding
1185 // rec. track and store in the fTreeCmp
1186 //
1187   TStopwatch timer;
1188   timer.Start();
1189   Int_t entry = fNextTreeGenEntryToRead;
1190   Double_t nParticlesTR = fTreeGenTracks->GetEntriesFast();
1191   cerr<<"fNParticles, nParticlesTR, fNextTreeGenEntryToRead: "<<fNParticles<<" "
1192       <<nParticlesTR<<" "<<fNextTreeGenEntryToRead<<endl;
1193   TBranch * branch = fTreeCmp->GetBranch("RC");
1194   branch->SetAddress(&fRecInfo); // set all pointers
1195
1196   while (entry < nParticlesTR) {
1197     fTreeGenTracks->GetEntry(entry);
1198     entry++;
1199     if (fEventNr < eventNr) continue;
1200     if (fEventNr > eventNr) break;
1201     fNextTreeGenEntryToRead = entry-1;
1202     if (fDebug > 2 && fMCInfo->fLabel < 10) {
1203       cerr<<"Fill track with a label "<<fMCInfo->fLabel<<endl;
1204     }
1205
1206     fRecInfo->Reset();
1207
1208
1209     if (fIndexRecTracks[fMCInfo->fLabel*4] >= 0) {
1210       fTPCTrack= (AliTPCtrack*)fTracks->UncheckedAt(fIndexRecTracks[fMCInfo->fLabel*4]);
1211
1212       //      if (nBytes > 0) {
1213       if (fTPCTrack) {
1214         //
1215         TVector3 local = TR2Local(&(fMCInfo->fTrackRef),fParamTPC);
1216         local.GetXYZ(fRecInfo->fTRLocalCoord);
1217         //
1218         // find nearest track if multifound
1219         if (fIndexRecTracks[fMCInfo->fLabel*4]+1){
1220           Float_t dz = TMath::Abs(local.Z()-fTPCTrack->GetZ());
1221           for (Int_t i=1;i<4;i++){
1222             if (fIndexRecTracks[fMCInfo->fLabel*4+i]>=0){
1223               AliTPCtrack * track = (AliTPCtrack*)fTracks->UncheckedAt(fIndexRecTracks[fMCInfo->fLabel*4+i]);
1224               if  (TMath::Abs(local.Z()-track->GetZ())<dz)
1225                 fTPCTrack = track;
1226             }
1227           }
1228         }
1229         fRecInfo->fTP=0;
1230         if (fTrackPoints){
1231           Int_t id = fTPCTrack->GetUniqueID();
1232           if (fTrackPoints->UncheckedAt(id)){
1233             fRecInfo->fTP = (TClonesArray*)fTrackPoints->UncheckedAt(id);
1234             //      fTrackPoints->AddAt(0,id);   //not owner anymore
1235           }
1236         }
1237         fRecInfo->fTPCTrack =*fTPCTrack;
1238         fRecInfo->fReconstructed = 1;
1239         fRecInfo->fFake = fFakeRecTracks[fMCInfo->fLabel];
1240         fRecInfo->fMultiple = fMultiRecTracks[fMCInfo->fLabel];
1241         fRecInfo->fdEdx = fTPCTrack->GetdEdx();
1242         //
1243         fRecInfo->fTPCTrack.PropagateTo(local.X());
1244         Double_t par[5];
1245         //
1246         Double_t localX = local.X();
1247         fTPCTrack->GetExternalParameters(localX,par);
1248         fRecInfo->fRecPhi=TMath::ASin(par[2]) + fTPCTrack->GetAlpha();
1249         if (fRecInfo->fRecPhi<0) fRecInfo->fRecPhi+=2*TMath::Pi();
1250         if (fRecInfo->fRecPhi>=2*TMath::Pi()) fRecInfo->fRecPhi-=2*TMath::Pi();
1251 //        fRecInfo->fRecPhi = (fRecInfo->fRecPhi)*kRaddeg;
1252         fRecInfo->fLambda = TMath::ATan(par[3]);
1253         fRecInfo->fRecPt_1 = TMath::Abs(par[4]);
1254       }
1255
1256     }
1257
1258     fTreeCmp->Fill();
1259   }
1260   fTreeCmp->AutoSave();
1261   fTracks->Delete();
1262   fTPCTrack =0;
1263   if (fTrackPoints){
1264     fTrackPoints->Delete();
1265     delete fTrackPoints;
1266     fTrackPoints =0;
1267   }
1268   printf("Time spended in TreeKLoop\n");
1269   timer.Print();
1270   if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl;
1271
1272   return 0;
1273 }
1274 ////////////////////////////////////////////////////////////////////////
1275 ////////////////////////////////////////////////////////////////////////
1276
1277 void AliTPCComparisonDraw::SetIO(const char *fname)
1278 {
1279   //
1280    TFile* file = TFile::Open(fname);
1281    if (!file) {
1282      printf("Could not open file  - generated new one\n");
1283      TFile* filegen = TFile::Open("genTracks.root");
1284      if (!filegen){
1285        printf("FILE with MC information is generated\n");
1286        TPCFindGenTracks *t = new TPCFindGenTracks("galice.root","genTracks.root",0);
1287        t->Exec();
1288        delete t;
1289      }
1290      filegen = TFile::Open("genTracks.root");
1291      if (!filegen){
1292        printf("COMPARISON  FILE COULDNT BE GENERATED \n");
1293        return;
1294      }
1295      printf("COMPARISON  FILE IS GENERATED \n");
1296      TPCCmpTr *t2 = new TPCCmpTr("genTracks.root","cmpTracks.root","galice.root",0);
1297      t2->Exec();
1298    }
1299    file = TFile::Open(fname);
1300    if (!file){
1301      printf("Comparison file couldn't be generated\n");
1302      return;
1303    }
1304    //
1305    fTree = (TTree*) file->Get("TPCcmpTracks");
1306    if (!fTree) {
1307     printf("no track comparison tree found\n");
1308     file->Close();
1309     delete file;
1310   }
1311 }
1312
1313 void AliTPCComparisonDraw::ResPt()
1314 {
1315  //
1316   //
1317   gStyle->SetOptFit();
1318   TCanvas *c1=new TCanvas("TPC pt resolution","TPC pt resolution",0,0,700,850);
1319   c1->Draw(); c1->cd();
1320   TPad *p1=new TPad("p1","p1",0.01,0.51,.99,.99);
1321   p1->Draw();p1->cd();
1322   p1->SetGridx(); p1->SetGridy();
1323   //
1324   c1->cd();
1325   TPad *p2=new TPad("p2","p2",0.01,0.01,.99,.49);
1326   p2->Draw();p2->cd();
1327   p2->SetGridx(); p2->SetGridy();
1328   //
1329   //Default cuts
1330   TCut cprim("cprim","MC.fVDist[3]<1");
1331   TCut cnprim("cnprim","MC.fVDist[3]>1");
1332   TCut cteta1("cteta1","abs(MC.fTrackRef.Theta()/3.1415-0.5)<0.25");
1333   TCut cpos1("cpos1","abs(MC.fParticle.fVz/sqrt(MC.fParticle.fVx*MC.fParticle.fVx+MC.fParticle.fVy*MC.fParticle.fVy))<1");
1334   //
1335   c1->cd();  p1->cd();
1336   TH1F* hisp = ResPtvsPt("MC.fRowsWithDigits>100"+cteta1+cpos1,"1",0.15,2.,6);
1337   c1->cd();
1338   c1->Draw();
1339   p1->cd();
1340   p1->Draw();
1341   hisp->Draw();
1342   //
1343   c1->cd();
1344   c1->Draw();
1345   p2->cd();
1346   p2->Draw();
1347   TH1F* his2 =  new TH1F("Ptresolution","Ptresolution",40,-5.,5.);
1348   fTree->Draw("100.*(abs(1./fTPCTrack.Get1Pt())-MC.fTrackRef.Pt())/MC.fTrackRef.Pt()>>Ptresolution","MC.fRowsWithDigits>100&&RC.fTPCTrack.fN>50&&RC.fMultiple==1"+cteta1+cpos1+cprim);
1349   AliLabelAxes(his2, "#Delta p_{t} / p_{t} [%]", "entries");
1350   his2->Fit("gaus");
1351   his2->Draw();
1352
1353 }
1354
1355 void AliTPCComparisonDraw::Eff()
1356 {
1357   //
1358   //
1359   TCanvas *c1=new TCanvas("TPC efficiency","TPC efficiency",0,0,700,850);
1360   c1->Draw(); c1->cd();
1361   TPad *p1=new TPad("p1","p1",0.01,0.51,.99,.99);
1362   p1->Draw();p1->cd();
1363   p1->SetGridx(); p1->SetGridy();
1364   //
1365   c1->cd();
1366   TPad *p2=new TPad("p2","p2",0.01,0.01,.99,.49);
1367   p2->Draw();p2->cd();
1368   p2->SetGridx(); p2->SetGridy();
1369   //
1370   //Default cuts
1371   TCut cprim("cprim","MC.fVDist[3]<1");
1372   TCut cnprim("cnprim","MC.fVDist[3]>1");
1373   TCut cteta1("cteta1","abs(MC.fTrackRef.Theta()/3.1415-0.5)<0.25");
1374   TCut cpos1("cpos1","abs(MC.fParticle.fVz/sqrt(MC.fParticle.fVx*MC.fParticle.fVx+MC.fParticle.fVy*MC.fParticle.fVy))<1");
1375   //
1376   c1->cd();  p1->cd();
1377   TH1F* hisp = EffVsPt("MC.fRowsWithDigits>100"+cteta1+cpos1+cprim,"RC.fTPCTrack.fN>50");
1378   hisp->Draw();
1379   //hisp->DrawClone();
1380   TText * text = new TText(0.25,102.,"Primary particles");
1381   text->SetTextSize(0.05);
1382   text->Draw();
1383   //
1384   c1->cd();  p2->cd();
1385   TH1F* hiss = EffVsPt("MC.fRowsWithDigits>100"+cteta1+cpos1+cnprim,"RC.fTPCTrack.fN>50");
1386   hiss->Draw();
1387   //hiss->DrawClone();
1388   text = new TText(0.25,102.,"Secondary particles");
1389   text->SetTextSize(0.05);
1390   text->Draw();
1391
1392
1393 }
1394
1395
1396 TH1F * AliTPCComparisonDraw::EffVsPt(const char* selection, const char * quality, Float_t min, Float_t max)
1397 {
1398   //
1399   //
1400   Int_t nBins = 10;
1401   Double_t* bins = CreateLogBins(nBins, min, max);
1402   TH1F* hGen = new TH1F("hGen", "gen. tracks", nBins, bins);
1403   TH1F* hRec = new TH1F("hRec", "rec. tracks", nBins, bins);
1404
1405   fTree->Draw("MC.fParticle.Pt()>>hGen", selection, "groff");
1406   char selectionRec[256];
1407   sprintf(selectionRec, "%s && RC.fReconstructed && %s", selection, quality);
1408   fTree->Draw("MC.fParticle.Pt()>>hRec", selectionRec, "groff");
1409
1410   TH1F* hEff = CreateEffHisto(hGen, hRec);
1411   AliLabelAxes(hEff, "p_{t} [GeV/c]", "#epsilon [%]");
1412
1413   delete hRec;
1414   delete hGen;
1415   delete[] bins;
1416   return hEff;
1417 }
1418
1419
1420 TH1F * AliTPCComparisonDraw::EffVsRM(const char* selection, const char * quality, Float_t min, Float_t max, Int_t nBins)
1421 {
1422   //
1423   TH1F* hGen = new TH1F("hGen", "gen. tracks", nBins, min, max);
1424   TH1F* hRec = new TH1F("hRec", "rec. tracks", nBins, min, max);
1425   //
1426   fTree->Draw("sqrt(MC.fDecayCoord[0]*MC.fDecayCoord[0] + MC.fDecayCoord[1]*MC.fDecayCoord[1])>>hGen", selection, "groff");
1427   char selectionRec[256];
1428   sprintf(selectionRec, "%s && RC.fReconstructed && %s", selection, quality);
1429   fTree->Draw("sqrt(MC.fDecayCoord[0]*MC.fDecayCoord[0] + MC.fDecayCoord[1]*MC.fDecayCoord[1])>>hRec", selectionRec, "groff");
1430   //
1431   TH1F* hEff = CreateEffHisto(hGen, hRec);
1432   AliLabelAxes(hEff, "r_{vertex} [cm]", "#epsilon [%]");
1433   //
1434   delete hRec;
1435   delete hGen;
1436   return hEff;
1437 }
1438
1439 TH1F * AliTPCComparisonDraw::EffVsRS(const char* selection, const char * quality, Float_t min, Float_t max, Int_t nBins)
1440 {
1441   //
1442   TH1F* hGen = new TH1F("hGen", "gen. tracks", nBins, min, max);
1443   TH1F* hRec = new TH1F("hRec", "rec. tracks", nBins, min, max);
1444   //
1445   fTree->Draw("sqrt(MC.fVDist[0]*MC.fVDist[0] + MC.fVDist[1]*MC.fVDist[1])>>hGen", selection, "groff");
1446   char selectionRec[256];
1447   sprintf(selectionRec, "%s && RC.fReconstructed && %s", selection, quality);
1448   fTree->Draw("sqrt(MC.fVDist[0]*MC.fVDist[0] + MC.fVDist[1]*MC.fVDist[1])>>hRec", selectionRec, "groff");
1449   //
1450   TH1F* hEff = CreateEffHisto(hGen, hRec);
1451   AliLabelAxes(hEff, "r_{vertex} [cm]", "#epsilon [%]");
1452   //
1453   delete hRec;
1454   delete hGen;
1455   return hEff;
1456
1457 }
1458
1459 TH1F * AliTPCComparisonDraw::ResPtvsPt(const char* selection, const char * quality, Float_t min, Float_t max, Int_t nBins)
1460 {
1461   //
1462   Double_t* bins = CreateLogBins(nBins, min, max);
1463   Int_t nBinsRes = 30;
1464   Double_t maxRes = 10.;
1465   TH2F* hRes2 = new TH2F("hRes2", "residuals", nBins, bins, nBinsRes, -maxRes, maxRes);
1466
1467   fTree->Draw("100.*(abs(1./fTPCTrack.Get1Pt())-MC.fTrackRef.Pt())/MC.fTrackRef.Pt():MC.fTrackRef.Pt()>>hRes2", selection, "groff");
1468
1469   TH1F* hMean=0;
1470   TH1F* hRes = CreateResHisto(hRes2, &hMean);
1471   AliLabelAxes(hRes, "p_{t} [GeV/c]", "#Delta p_{t} / p_{t} [%]");
1472   //
1473   delete hRes2;
1474   delete[] bins;
1475   return hRes;
1476
1477 }
1478
1479 TH1F * AliTPCComparisonDraw::MeanPtvsPt(const char* selection, const char * quality, Float_t min, Float_t max, Int_t nBins)
1480 {
1481   //
1482   Double_t* bins = CreateLogBins(nBins, min, max);
1483   Int_t nBinsRes = 30;
1484   Double_t maxRes = 10.;
1485   TH2F* hRes2 = new TH2F("hRes2", "residuals", nBins, bins, nBinsRes, -maxRes, maxRes);
1486
1487   fTree->Draw("100.*(1./fTPCTrack.Get1Pt()-MC.fTrackRef.Pt())/MC.fTrackRef.Pt():MC.fTrackRef.Pt()>>hRes2", selection, "groff");
1488
1489   TH1F* hMean=0;
1490   TH1F* hRes = CreateResHisto(hRes2, &hMean);
1491   AliLabelAxes(hRes, "p_{t} [GeV/c]", "mean p_{t} / p_{t} [%]");
1492   //
1493   delete hRes2;
1494   delete[] bins;
1495   if (!hMean) return 0;
1496   return hMean;
1497
1498 }
1499
1500
1501 TH1F * AliTPCComparisonDraw::ResdEdxvsN(const char* selection, const char * quality, Float_t min, Float_t max, Int_t nBins)
1502 {
1503   //
1504   Int_t nBinsRes = 15;
1505
1506   TH2F* hRes2 = new TH2F("hRes2", "residuals", nBins, min, max, nBinsRes, 34, 60);
1507
1508   fTree->Draw("RC.fTPCTrack.fdEdx/MC.fPrim:RC.fTPCTrack.fN>>hRes2", selection, "groff");
1509
1510   TH1F* hMean=0;
1511   TH1F* hRes = CreateResHisto(hRes2, &hMean);
1512   AliLabelAxes(hRes, "N points", "sigma dEdx/Nprim [%]");
1513   //
1514   delete hRes2;
1515   return hRes;
1516
1517 }
1518
1519
1520
1521 Double_t* AliTPCComparisonDraw::CreateLogBins(Int_t nBins, Double_t xMin, Double_t xMax)
1522 {
1523   Double_t* bins = new Double_t[nBins+1];
1524   bins[0] = xMin;
1525   Double_t factor = pow(xMax/xMin, 1./nBins);
1526   for (Int_t i = 1; i <= nBins; i++)
1527     bins[i] = factor * bins[i-1];
1528   return bins;
1529 }
1530
1531
1532
1533
1534 void AliTPCComparisonDraw::AliLabelAxes(TH1* histo, const char* xAxisTitle, const char* yAxisTitle)
1535 {
1536   //
1537   histo->GetXaxis()->SetTitle(xAxisTitle);
1538   histo->GetYaxis()->SetTitle(yAxisTitle);
1539 }
1540
1541
1542 TH1F* AliTPCComparisonDraw::CreateEffHisto(TH1F* hGen, TH1F* hRec)
1543 {
1544   //
1545   Int_t nBins = hGen->GetNbinsX();
1546   TH1F* hEff = (TH1F*) hGen->Clone("hEff");
1547   hEff->SetTitle("");
1548   hEff->SetStats(kFALSE);
1549   hEff->SetMinimum(0.);
1550   hEff->SetMaximum(110.);
1551   //
1552   for (Int_t iBin = 0; iBin <= nBins; iBin++) {
1553     Double_t nGen = hGen->GetBinContent(iBin);
1554     Double_t nRec = hRec->GetBinContent(iBin);
1555     if (nGen > 0) {
1556       Double_t eff = nRec/nGen;
1557       hEff->SetBinContent(iBin, 100. * eff);
1558       Double_t error = sqrt((eff*(1.-eff)+0.01) / nGen);
1559       //      if (error == 0) error = sqrt(0.1/nGen);
1560       //
1561       if (error == 0) error = 0.0001;
1562       hEff->SetBinError(iBin, 100. * error);
1563     } else {
1564       hEff->SetBinContent(iBin, 100. * 0.5);
1565       hEff->SetBinError(iBin, 100. * 0.5);
1566     }
1567   }
1568   return hEff;
1569 }
1570
1571
1572
1573 TH1F* AliTPCComparisonDraw::CreateResHisto(TH2F* hRes2, TH1F **phMean,  Bool_t drawBinFits,
1574                      Bool_t overflowBinFits)
1575 {
1576   TVirtualPad* currentPad = gPad;
1577   TAxis* axis = hRes2->GetXaxis();
1578   Int_t nBins = axis->GetNbins();
1579   TH1F* hRes, *hMean;
1580   if (axis->GetXbins()->GetSize()){
1581     hRes = new TH1F("hRes", "", nBins, axis->GetXbins()->GetArray());
1582     hMean = new TH1F("hMean", "", nBins, axis->GetXbins()->GetArray());
1583   }
1584   else{
1585     hRes = new TH1F("hRes", "", nBins, axis->GetXmin(), axis->GetXmax());
1586     hMean = new TH1F("hMean", "", nBins, axis->GetXmin(), axis->GetXmax());
1587
1588   }
1589   hRes->SetStats(false);
1590   hRes->SetOption("E");
1591   hRes->SetMinimum(0.);
1592   //
1593   hMean->SetStats(false);
1594   hMean->SetOption("E");
1595
1596   // create the fit function
1597   //TKFitGaus* fitFunc = new TKFitGaus("resFunc");
1598   //   TF1 * fitFunc = new TF1("G","[3]+[0]*exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))",-3,3);
1599   TF1 * fitFunc = new TF1("G","[0]*exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))",-3,3);
1600
1601   fitFunc->SetLineWidth(2);
1602   fitFunc->SetFillStyle(0);
1603   // create canvas for fits
1604   TCanvas* canBinFits = NULL;
1605   Int_t nPads = (overflowBinFits) ? nBins+2 : nBins;
1606   Int_t nx = Int_t(sqrt(nPads-1.));// + 1;
1607   Int_t ny = (nPads-1) / nx + 1;
1608   if (drawBinFits) {
1609     canBinFits = (TCanvas*)gROOT->FindObject("canBinFits");
1610     if (canBinFits) delete canBinFits;
1611     canBinFits = new TCanvas("canBinFits", "fits of bins", 200, 100, 500, 700);
1612     canBinFits->Divide(nx, ny);
1613   }
1614
1615   // loop over x bins and fit projection
1616   Int_t dBin = ((overflowBinFits) ? 1 : 0);
1617   for (Int_t bin = 1-dBin; bin <= nBins+dBin; bin++) {
1618     if (drawBinFits) canBinFits->cd(bin + dBin);
1619     TH1D* hBin = hRes2->ProjectionY("hBin", bin, bin);
1620     //
1621     if (hBin->GetEntries() > 5) {
1622       fitFunc->SetParameters(hBin->GetMaximum(),hBin->GetMean(),hBin->GetRMS());
1623       hBin->Fit(fitFunc,"s");
1624       Double_t sigma = TMath::Abs(fitFunc->GetParameter(2));
1625
1626       if (sigma > 0.){
1627         hRes->SetBinContent(bin, TMath::Abs(fitFunc->GetParameter(2)));
1628         hMean->SetBinContent(bin, fitFunc->GetParameter(1));
1629       }
1630       else{
1631         hRes->SetBinContent(bin, 0.);
1632         hMean->SetBinContent(bin,0);
1633       }
1634       hRes->SetBinError(bin, fitFunc->GetParError(2));
1635       hMean->SetBinError(bin, fitFunc->GetParError(1));
1636
1637       //
1638       //
1639
1640     } else {
1641       hRes->SetBinContent(bin, 0.);
1642       hRes->SetBinError(bin, 0.);
1643       hMean->SetBinContent(bin, 0.);
1644       hMean->SetBinError(bin, 0.);
1645     }
1646
1647
1648     if (drawBinFits) {
1649       char name[256];
1650       if (bin == 0) {
1651         sprintf(name, "%s < %.4g", axis->GetTitle(), axis->GetBinUpEdge(bin));
1652       } else if (bin == nBins+1) {
1653         sprintf(name, "%.4g < %s", axis->GetBinLowEdge(bin), axis->GetTitle());
1654       } else {
1655         sprintf(name, "%.4g < %s < %.4g", axis->GetBinLowEdge(bin),
1656                 axis->GetTitle(), axis->GetBinUpEdge(bin));
1657       }
1658       canBinFits->cd(bin + dBin);
1659       hBin->SetTitle(name);
1660       hBin->SetStats(kTRUE);
1661       hBin->DrawCopy("E");
1662       canBinFits->Update();
1663       canBinFits->Modified();
1664       canBinFits->Update();
1665     }
1666
1667     delete hBin;
1668   }
1669
1670   delete fitFunc;
1671   currentPad->cd();
1672   *phMean = hMean;
1673   return hRes;
1674 }
1675
1676
1677