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