TPC/AliTPCComparisonMI.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2 #include "iostream.h"
   3 #include <stdio.h>
   4 #include <string.h>
   5 #include "Rtypes.h"
   6 #include "TFile.h"
   7 #include "TTree.h"
   8 #include "TString.h"
   9 #include "TBenchmark.h"
  10 #include "TStopwatch.h"
  11 #include "TParticle.h"
  12 #include "AliRun.h"
  13 #include "AliStack.h"
  14 #include "AliTPCtrack.h"
  15 #include "AliSimDigits.h"
  16 #include "AliTPCParam.h"
  17 #include "TParticle.h"
  18 #include "AliTPC.h"
  19 #include "AliDetector.h"
  20 #include "AliTrackReference.h"
  21 #include "TSystem.h"
  22 #include "TTimer.h"
  23 #include "TVector3.h"
  24 //
  25 #include "iostream.h"
  26 #include "Rtypes.h"
  27 #include "TSystem.h"
  28 #include "TTimer.h"
  29 #include "Getline.h"
  30 #include "TChain.h"
  31 #include "TString.h"
  32 #include "TFile.h"
  33 #include "AliRun.h"
  34 #include "AliTPCComparisonMI.h"
  35 #include "AliTPCParamSR.h"
  36 #include "AliTracker.h"
  37 #include "AliComplexCluster.h"
  38
  39
  40 #endif
  41
  42 //
  43 //
  44 Bool_t ImportgAlice(TFile *file);
  45 TFile* OpenAliceFile(char *fn, Bool_t importgAlice = kFALSE, char *mode = "read");
  46
  47 AliTPCParam * GetTPCParam(){
  48   AliTPCParamSR * par = new AliTPCParamSR;
  49   par->Update();
  50   return par;
  51 }
  52
  53
  54
  55 ////////////////////////////////////////////////////////////////////////
  56 Bool_t ImportgAlice(TFile *file) {
  57 // read in gAlice object from the file
  58   gAlice = (AliRun*)file->Get("gAlice");
  59   if (!gAlice)  return kFALSE;
  60   return kTRUE;
  61 }
  62 ////////////////////////////////////////////////////////////////////////
  63 TFile* OpenAliceFile(char *fn, Bool_t importgAlice, char *mode) {
  64   TFile *file = TFile::Open(fn,mode);
  65   if (!file->IsOpen()) {
  66     cerr<<"OpenAliceFile: cannot open file "<<fn<<" in mode "
  67         <<mode<<endl;
  68     return 0;
  69   }
  70   if (!importgAlice) return file;
  71   if (ImportgAlice(file)) return file;
  72   return 0;
  73 }
  74 ////////////////////////////////////////////////////////////////////////
  75
  76 //_____________________________________________________________________________
  77 Float_t TPCBetheBloch(Float_t bg)
  78 {
  79   //
  80   // Bethe-Bloch energy loss formula
  81   //
  82   const Double_t kp1=0.76176e-1;
  83   const Double_t kp2=10.632;
  84   const Double_t kp3=0.13279e-4;
  85   const Double_t kp4=1.8631;
  86   const Double_t kp5=1.9479;
  87
  88   Double_t dbg = (Double_t) bg;
  89
  90   Double_t beta = dbg/TMath::Sqrt(1.+dbg*dbg);
  91
  92   Double_t aa = TMath::Power(beta,kp4);
  93   Double_t bb = TMath::Power(1./dbg,kp5);
  94
  95   bb=TMath::Log(kp3+bb);
  96
  97   return ((Float_t)((kp2-aa-bb)*kp1/aa));
  98 }
  99
 100
 101
 102
 103
 104 ////////////////////////////////////////////////////////////////////////
 105 AliTPCGenInfo::AliTPCGenInfo()
 106 {
 107   fReferences  = 0;
 108 }
 109
 110 AliTPCGenInfo::~AliTPCGenInfo()
 111 {
 112   if (fReferences) delete fReferences;
 113   fReferences =0;
 114 }
 115 ////////////////////////////////////////////////////////////////////////
 116
 117 ////////////////////////////////////////////////////////////////////////
 118 //
 119 // End of implementation of the class AliTPCGenInfo
 120 //
 121 ////////////////////////////////////////////////////////////////////////
 122
 123
 124
 125 ////////////////////////////////////////////////////////////////////////
 126 digitRow::digitRow()
 127 {
 128   Reset();
 129 }
 130 ////////////////////////////////////////////////////////////////////////
 131 digitRow & digitRow::operator=(const digitRow &digOld)
 132 {
 133   for (Int_t i = 0; i<kgRowBytes; i++) fDig[i] = digOld.fDig[i];
 134   return (*this);
 135 }
 136 ////////////////////////////////////////////////////////////////////////
 137 void digitRow::SetRow(Int_t row)
 138 {
 139   if (row >= 8*kgRowBytes) {
 140     cerr<<"digitRow::SetRow: index "<<row<<" out of bounds."<<endl;
 141     return;
 142   }
 143   Int_t iC = row/8;
 144   Int_t iB = row%8;
 145   SETBIT(fDig[iC],iB);
 146 }
 147
 148 ////////////////////////////////////////////////////////////////////////
 149 Bool_t digitRow::TestRow(Int_t row)
 150 {
 151 //
 152 // return kTRUE if row is on
 153 //
 154   Int_t iC = row/8;
 155   Int_t iB = row%8;
 156   return TESTBIT(fDig[iC],iB);
 157 }
 158 ////////////////////////////////////////////////////////////////////////
 159 Int_t digitRow::RowsOn(Int_t upto)
 160 {
 161 //
 162 // returns number of rows with a digit
 163 // count only rows less equal row number upto
 164 //
 165   Int_t total = 0;
 166   for (Int_t i = 0; i<kgRowBytes; i++) {
 167     for (Int_t j = 0; j < 8; j++) {
 168       if (i*8+j > upto) return total;
 169       if (TESTBIT(fDig[i],j))  total++;
 170     }
 171   }
 172   return total;
 173 }
 174 ////////////////////////////////////////////////////////////////////////
 175 void digitRow::Reset()
 176 {
 177 //
 178 // resets all rows to zero
 179 //
 180   for (Int_t i = 0; i<kgRowBytes; i++) {
 181     fDig[i] <<= 8;
 182   }
 183 }
 184 ////////////////////////////////////////////////////////////////////////
 185 Int_t digitRow::Last()
 186 {
 187 //
 188 // returns the last row number with a digit
 189 // returns -1 if now digits
 190 //
 191   for (Int_t i = kgRowBytes-1; i>=0; i--) {
 192     for (Int_t j = 7; j >= 0; j--) {
 193       if TESTBIT(fDig[i],j) return i*8+j;
 194     }
 195   }
 196   return -1;
 197 }
 198 ////////////////////////////////////////////////////////////////////////
 199 Int_t digitRow::First()
 200 {
 201 //
 202 // returns the first row number with a digit
 203 // returns -1 if now digits
 204 //
 205   for (Int_t i = 0; i<kgRowBytes; i++) {
 206     for (Int_t j = 0; j < 8; j++) {
 207       if (TESTBIT(fDig[i],j)) return i*8+j;
 208     }
 209   }
 210   return -1;
 211 }
 212
 213 ////////////////////////////////////////////////////////////////////////
 214 //
 215 // end of implementation of a class digitRow
 216 //
 217 ////////////////////////////////////////////////////////////////////////
 218
 219 ////////////////////////////////////////////////////////////////////////
 220 TPCFindGenTracks::TPCFindGenTracks()
 221 {
 222   fMCInfo = new AliTPCGenInfo;
 223   fMCInfo->fReferences = new TClonesArray("AliTrackReference");
 224   Reset();
 225 }
 226
 227 ////////////////////////////////////////////////////////////////////////
 228 TPCFindGenTracks::TPCFindGenTracks(char * fnGalice, char* fnRes,
 229                                    Int_t nEvents, Int_t firstEvent)
 230 {
 231   Reset();
 232   fFirstEventNr = firstEvent;
 233   fEventNr = firstEvent;
 234   fNEvents = nEvents;
 235   fFnRes = fnRes;
 236   //
 237   fLoader = AliRunLoader::Open(fnGalice);
 238   if (gAlice){
 239     delete gAlice->GetRunLoader();
 240     delete gAlice;
 241     gAlice = 0x0;
 242   }
 243   if (fLoader->LoadgAlice()){
 244     cerr<<"Error occured while l"<<endl;
 245   }
 246   Int_t nall = fLoader->GetNumberOfEvents();
 247   if (nall<=0){
 248     cerr<<"no events available"<<endl;
 249     fEventNr = 0;
 250     return;
 251   }
 252   if (firstEvent+nEvents>nall) {
 253     fEventNr = nall-firstEvent;
 254     cerr<<"restricted number of events availaible"<<endl;
 255   }
 256 }
 257
 258 ////////////////////////////////////////////////////////////////////////
 259 void TPCFindGenTracks::Reset()
 260 {
 261   fEventNr = 0;
 262   fNEvents = 0;
 263   fTreeGenTracks = 0;
 264   fFnRes = "genTracks.root";
 265   fFileGenTracks = 0;
 266   fContainerDigitRow = 0;
 267   //
 268   fReferences = 0;
 269   fReferenceIndex0 = 0;
 270   fReferenceIndex1 = 0;
 271   //
 272   fDebug = 0;
 273   fVPrim[0] = -1000.;
 274   fVPrim[1] = -1000.;
 275   fVPrim[2] = -1000.;
 276   fParamTPC = 0;
 277 }
 278 ////////////////////////////////////////////////////////////////////////
 279 TPCFindGenTracks::~TPCFindGenTracks()
 280 {
 281   ;
 282 }
 283
 284
 285 Int_t  TPCFindGenTracks::SetIO()
 286 {
 287   //
 288   //
 289   CreateTreeGenTracks();
 290   if (!fTreeGenTracks) return 1;
 291   AliTracker::SetFieldFactor();
 292   fParamTPC = GetTPCParam();
 293   //
 294   return 0;
 295 }
 296
 297 ////////////////////////////////////////////////////////////////////////
 298 Int_t TPCFindGenTracks::SetIO(Int_t eventNr)
 299 {
 300   //
 301   //
 302   // SET INPUT
 303   gAlice->GetEvent(eventNr);
 304   fLoader->SetEventNumber(eventNr);
 305   //
 306   fLoader->LoadHeader();
 307   fLoader->LoadKinematics();
 308   fStack = fLoader->Stack();
 309   //
 310   fLoader->LoadTrackRefs();
 311   fTreeTR = fLoader->TreeTR();
 312   //
 313   AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
 314   tpcl->LoadDigits();
 315   fTreeD = tpcl->TreeD();
 316
 317   return 0;
 318 }
 319
 320 ////////////////////////////////////////////////////////////////////////
 321 Int_t TPCFindGenTracks::Exec(Int_t nEvents, Int_t firstEventNr)
 322 {
 323   fNEvents = nEvents;
 324   fFirstEventNr = firstEventNr;
 325   return Exec();
 326 }
 327
 328 ////////////////////////////////////////////////////////////////////////
 329 Int_t TPCFindGenTracks::Exec()
 330 {
 331   TStopwatch timer;
 332   timer.Start();
 333   Int_t status =SetIO();
 334   if (status>0) return status;
 335   //
 336
 337   for (fEventNr = fFirstEventNr; fEventNr < fFirstEventNr+fNEvents;
 338        fEventNr++) {
 339     SetIO(fEventNr);
 340     fNParticles = fStack->GetNtrack();
 341     fContainerDigitRow = new digitRow[fNParticles];
 342     //
 343     fReferences = new AliTrackReference[fgMaxTR];
 344     fReferenceIndex0 = new Int_t[fNParticles];
 345     fReferenceIndex1 = new Int_t[fNParticles];
 346
 347     for (Int_t i = 0; i<fNParticles; i++) {
 348       fReferenceIndex0[i] = -1;
 349       fReferenceIndex1[i] = -1;
 350     }
 351
 352     cout<<"Start to process event "<<fEventNr<<endl;
 353     cout<<"\tfNParticles = "<<fNParticles<<endl;
 354     if (fDebug>2) cout<<"\tStart loop over TreeD"<<endl;
 355     if (TreeDLoop()>0) return 1;
 356     if (fDebug>2) cout<<"\tStart loop over TreeTR"<<endl;
 357     if (TreeTRLoop()>0) return 1;
 358     if (fDebug>2) cout<<"\tStart loop over TreeK"<<endl;
 359     if (TreeKLoop()>0) return 1;
 360     if (fDebug>2) cout<<"\tEnd loop over TreeK"<<endl;
 361
 362     delete [] fContainerDigitRow;
 363     delete [] fReferences;
 364     delete [] fReferenceIndex0;
 365     delete [] fReferenceIndex1;
 366   }
 367
 368   CloseOutputFile();
 369
 370   cerr<<"Exec finished"<<endl;
 371   delete gAlice;
 372
 373   timer.Stop();
 374   timer.Print();
 375   return 0;
 376 }
 377 ////////////////////////////////////////////////////////////////////////
 378 void TPCFindGenTracks::CreateTreeGenTracks()
 379 {
 380   fFileGenTracks = TFile::Open(fFnRes,"RECREATE");
 381   if (!fFileGenTracks) {
 382     cerr<<"Error in CreateTreeGenTracks: cannot open file "<<fFnRes<<endl;
 383     return;
 384   }
 385   fTreeGenTracks = new TTree("genTracksTree","genTracksTree");
 386
 387
 388
 389   fMCInfo = new AliTPCGenInfo;
 390   fMCInfo->fReferences = new TClonesArray("AliTrackReference");
 391
 392   fTreeGenTracks->Branch("MC","AliTPCGenInfo",&fMCInfo);
 393
 394   fTreeGenTracks->Branch("fEventNr",&fEventNr,"fEventNr/I");
 395
 396   fTreeGenTracks->AutoSave();
 397 }
 398 ////////////////////////////////////////////////////////////////////////
 399 void TPCFindGenTracks::CloseOutputFile()
 400 {
 401   if (!fFileGenTracks) {
 402     cerr<<"File "<<fFnRes<<" not found as an open file."<<endl;
 403     return;
 404   }
 405   fFileGenTracks->cd();
 406   fTreeGenTracks->Write();
 407   delete fTreeGenTracks;
 408   fFileGenTracks->Close();
 409   delete fFileGenTracks;
 410   return;
 411 }
 412
 413 ////////////////////////////////////////////////////////////////////////
 414 Int_t TPCFindGenTracks::TreeKLoop()
 415 {
 416 //
 417 // open the file with treeK
 418 // loop over all entries there and save information about some tracks
 419 //
 420
 421   AliStack * stack = fStack;
 422   if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
 423
 424   if (fDebug > 0) {
 425     cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
 426         <<fEventNr<<endl;
 427   }
 428
 429   Int_t  ipdg = 0;
 430   TParticlePDG *ppdg = 0;
 431
 432
 433   // not all generators give primary vertex position. Take the vertex
 434   // of the particle 0 as primary vertex.
 435   TDatabasePDG  pdg; //get pdg table
 436
 437   //thank you very much root for this
 438   TBranch * br = fTreeGenTracks->GetBranch("MC");
 439
 440   TParticle *particle = stack->ParticleFromTreeK(0);
 441   fVPrim[0] = particle->Vx();
 442   fVPrim[1] = particle->Vy();
 443   fVPrim[2] = particle->Vz();
 444   for (Int_t iParticle = 0; iParticle < fNParticles; iParticle++) {
 445
 446 // load only particles with TR
 447     if (fReferenceIndex0[iParticle]<0) continue;
 448     //////////////////////////////////////////////////////////////////////
 449     fMCInfo->fLabel = iParticle;
 450
 451     fMCInfo->fRow = (fContainerDigitRow[iParticle]);
 452     fMCInfo->fRowsWithDigits = fMCInfo->fRow.RowsOn();
 453     if (fMCInfo->fRowsWithDigits<10) continue;
 454     fMCInfo->fRowsWithDigitsInn = fMCInfo->fRow.RowsOn(63); // 63 = number of inner rows
 455     fMCInfo->fRowsTrackLength = fMCInfo->fRow.Last() - fMCInfo->fRow.First();
 456     fMCInfo->fDigitsInSeed = 0;
 457     if (fMCInfo->fRow.TestRow(seedRow11) && fMCInfo->fRow.TestRow(seedRow12))
 458       fMCInfo->fDigitsInSeed = 1;
 459     if (fMCInfo->fRow.TestRow(seedRow21) && fMCInfo->fRow.TestRow(seedRow22))
 460       fMCInfo->fDigitsInSeed += 10;
 461     //
 462     //
 463     //
 464     fMCInfo->fParticle = *(stack->Particle(iParticle));
 465     if (fMCInfo->fParticle.GetLastDaughter()>0&&fMCInfo->fParticle.GetLastDaughter()<fNParticles){
 466       TParticle * dparticle0 = 0;
 467       if (fMCInfo->fParticle.GetDaughter(0)>0) dparticle0 = stack->Particle(fMCInfo->fParticle.GetDaughter(0));
 468       TParticle * dparticle1 = 0;
 469       if (fMCInfo->fParticle.GetDaughter(1)>0) dparticle1 = stack->Particle(fMCInfo->fParticle.GetDaughter(1));
 470       if ((dparticle1)&&(dparticle1->P()>0.15)) {
 471         fMCInfo->fDecayCoord[0] = dparticle1->Vx();
 472         fMCInfo->fDecayCoord[1] = dparticle1->Vy();
 473         fMCInfo->fDecayCoord[2] = dparticle1->Vz();
 474       }
 475       else
 476         fMCInfo->fDecayCoord[0]=1000000;
 477     }else
 478       fMCInfo->fDecayCoord[0]=1000000;
 479     fMCInfo->fParticle = *(stack->ParticleFromTreeK(iParticle));
 480     //
 481     //
 482     //
 483     fMCInfo->fLabel = iParticle;
 484     fMCInfo->fVDist[0] = fMCInfo->fParticle.Vx()-fVPrim[0];
 485     fMCInfo->fVDist[1] = fMCInfo->fParticle.Vy()-fVPrim[1];
 486     fMCInfo->fVDist[2] = fMCInfo->fParticle.Vz()-fVPrim[2];
 487     fMCInfo->fVDist[3] = TMath::Sqrt(fMCInfo->fVDist[0]*fMCInfo->fVDist[0]+
 488                                       fMCInfo->fVDist[1]*fMCInfo->fVDist[1]+fMCInfo->fVDist[2]*fMCInfo->fVDist[2]);
 489
 490     //
 491     Int_t index = fReferenceIndex0[iParticle];
 492     AliTrackReference  ref  = fReferences[index];
 493     // if (ref.GetTrack()!=iParticle)
 494     //  printf("problem2\n");
 495     fMCInfo->fTrackRef = ref;
 496
 497     Int_t rfindex =0;
 498     if (fMCInfo->fReferences !=0) delete fMCInfo->fReferences;
 499     fMCInfo->fReferences =
 500       new TClonesArray("AliTrackReference");
 501     fMCInfo->fReferences->ExpandCreateFast(fReferenceIndex1[iParticle]-fReferenceIndex0[iParticle]+1);
 502     //
 503     for (Int_t i = fReferenceIndex0[iParticle];i<=fReferenceIndex1[iParticle];i++){
 504       AliTrackReference  ref  = fReferences[i];
 505       AliTrackReference *ref2 = (AliTrackReference*) fMCInfo->fReferences->At(rfindex);
 506       if (ref.GetTrack()!=iParticle)
 507         printf("problem5\n");
 508       *ref2 = ref;
 509       rfindex++;
 510     }
 511     //
 512     //
 513     ipdg = fMCInfo->fParticle.GetPdgCode();
 514     ppdg = pdg.GetParticle(ipdg);
 515     // calculate primary ionization per cm
 516     if (ppdg){
 517       Float_t mass = ppdg->Mass();
 518       Float_t p = TMath::Sqrt(fMCInfo->fTrackRef.Px()*fMCInfo->fTrackRef.Px()+
 519                               fMCInfo->fTrackRef.Py()*fMCInfo->fTrackRef.Py()+
 520                               fMCInfo->fTrackRef.Pz()*fMCInfo->fTrackRef.Pz());
 521
 522       //      Float_t betagama = fMCInfo->fParticle.P()/mass;
 523       Float_t betagama = p /mass;
 524       fMCInfo->fPrim = TPCBetheBloch(betagama);
 525     }
 526     //////////////////////////////////////////////////////////////////////
 527     br->SetAddress(&fMCInfo);
 528     fTreeGenTracks->Fill();
 529
 530   }
 531   fTreeGenTracks->AutoSave();
 532
 533   if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl;
 534
 535   return 0;
 536 }
 537
 538
 539
 540
 541 ////////////////////////////////////////////////////////////////////////
 542 Int_t TPCFindGenTracks::TreeDLoop()
 543 {
 544 //
 545 // open the file with treeD
 546 // loop over all entries there and save information about some tracks
 547 //
 548
 549   Int_t nInnerSector = fParamTPC->GetNInnerSector();
 550   Int_t rowShift = 0;
 551   //  Int_t zero=fParamTPC->GetZeroSup();
 552   Int_t zero=fParamTPC->GetZeroSup()+6;
 553
 554
 555   //char treeDName[100];
 556   //sprintf(treeDName,"TreeD_75x40_100x60_150x60_%d",fEventNr);
 557   //TTree *treeD=(TTree*)fFileTreeD->Get(treeDName);
 558   //
 559   //
 560   AliSimDigits digitsAddress, *digits=&digitsAddress;
 561   fTreeD->GetBranch("Segment")->SetAddress(&digits);
 562
 563   Int_t sectorsByRows=(Int_t)fTreeD->GetEntries();
 564   if (fDebug > 1) cout<<"\tsectorsByRows = "<<sectorsByRows<<endl;
 565   for (Int_t i=0; i<sectorsByRows; i++) {
 566 //  for (Int_t i=5720; i<sectorsByRows; i++) {
 567     if (!fTreeD->GetEvent(i)) continue;
 568     Int_t sec,row;
 569     fParamTPC->AdjustSectorRow(digits->GetID(),sec,row);
 570     if (fDebug > 1) cout<<sec<<' '<<row<<"                          \r";
 571 //    cerr<<sec<<' '<<row<<endl;
 572
 573 // here I expect that upper sectors follow lower sectors
 574     if (sec > nInnerSector) rowShift = fParamTPC->GetNRowLow();
 575     //
 576     //digits->ExpandTrackBuffer();
 577     //Int_t *tracks = digits->GetTracks();
 578     digits->First();
 579     do {
 580       Int_t iRow=digits->CurrentRow();
 581       Int_t iColumn=digits->CurrentColumn();
 582       Short_t digitValue = digits->CurrentDigit();
 583 //      cout<<"Inner loop: sector, iRow, iColumn "
 584 //        <<sec<<" "<<iRow<<" "<<iColumn<<endl;
 585       if (digitValue >= zero) {
 586         Int_t label;
 587         for (Int_t j = 0; j<3; j++) {
 588           label = digits->GetTrackID(iRow,iColumn,j);
 589           //label = digits->GetTrackIDFast(iRow,iColumn,j);
 590
 591           if (label >= fNParticles) { //don't label from bakground event
 592             continue;
 593           }
 594           if (label >= 0 && label <= fNParticles) {
 595 //        if (label >= 0 && label <= fDebug) {
 596             if (fDebug > 6 ) {
 597               cout<<"Inner loop: sector, iRow, iColumn, label, value, row "
 598                   <<sec<<" "
 599                   <<iRow<<" "<<iColumn<<" "<<label<<" "<<digitValue
 600                   <<" "<<row<<endl;
 601             }
 602             fContainerDigitRow[label].SetRow(row+rowShift);
 603           }
 604         }
 605       }
 606     } while (digits->Next());
 607   }
 608
 609   if (fDebug > 2) cerr<<"end of TreeDLoop"<<endl;
 610
 611   return 0;
 612 }
 613
 614
 615 ////////////////////////////////////////////////////////////////////////
 616 Int_t TPCFindGenTracks::TreeTRLoop()
 617 {
 618   //
 619   // loop over TrackReferences and store the first one for each track
 620   //
 621
 622   TTree * treeTR = fTreeTR;
 623   if (!treeTR) {
 624     cerr<<"TreeTR not found"<<endl;
 625     return 1;
 626   }
 627   Int_t nPrimaries = (Int_t) treeTR->GetEntries();
 628   if (fDebug > 1) cout<<"There are "<<nPrimaries<<" entries in TreeTR"<<endl;
 629   //
 630   //
 631   TBranch *TPCBranchTR  = treeTR->GetBranch("TPC");
 632   if (!TPCBranchTR) {
 633     cerr<<"TPC branch in TR not found"<<endl;
 634     return 1;
 635   }
 636   TClonesArray* TPCArrayTR = new TClonesArray("AliTrackReference");
 637   TPCBranchTR->SetAddress(&TPCArrayTR);
 638   //
 639   //
 640   //
 641   Int_t index     =  0;
 642   for (Int_t iPrimPart = 0; iPrimPart<nPrimaries; iPrimPart++) {
 643     TPCBranchTR->GetEntry(iPrimPart);
 644     for (Int_t iTrackRef = 0; iTrackRef < TPCArrayTR->GetEntriesFast(); iTrackRef++) {
 645       AliTrackReference *trackRef = (AliTrackReference*)TPCArrayTR->At(iTrackRef);
 646       //
 647       if (trackRef->Pt() < fgPtCut) continue;
 648       Int_t label = trackRef->GetTrack();
 649       if (label<0  || label > fNParticles) {
 650         continue;
 651       }
 652       if (index>=fgMaxTR) continue;     //restricted size of buffer for TR
 653       fReferences[index] = *trackRef;
 654       fReferenceIndex1[label] = index;  // the last ref with given label
 655       if (fReferenceIndex0[label]==-1) fReferenceIndex0[label] = index;   //the first index with label
 656       index++;
 657
 658     }
 659   }
 660   return 0;
 661 }
 662
 663 ////////////////////////////////////////////////////////////////////////
 664 Float_t TPCFindGenTracks::TR2LocalX(AliTrackReference *trackRef,
 665                                     AliTPCParam *paramTPC) {
 666
 667   Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
 668   Int_t index[4];
 669   paramTPC->Transform0to1(x,index);
 670   paramTPC->Transform1to2(x,index);
 671   return x[0];
 672 }
 673 ////////////////////////////////////////////////////////////////////////
 674
 675
 676
 677
 678 ////////////////////////////////////////////////////////////////////////
 679 TPCCmpTr::TPCCmpTr()
 680 {
 681   Reset();
 682 }
 683
 684 ////////////////////////////////////////////////////////////////////////
 685 TPCCmpTr::TPCCmpTr(char* fnGenTracks,
 686                    char* fnCmp,
 687                    char* fnGalice,
 688                    Int_t nEvents, Int_t firstEvent)
 689 {
 690   Reset();
 691   fFnGenTracks = fnGenTracks;
 692   fFnCmp = fnCmp;
 693   fFirstEventNr = firstEvent;
 694   fEventNr = firstEvent;
 695   fNEvents = nEvents;
 696
 697   //
 698   fLoader = AliRunLoader::Open(fnGalice);
 699   if (gAlice){
 700     delete gAlice->GetRunLoader();
 701     delete gAlice;
 702     gAlice = 0x0;
 703   }
 704   if (fLoader->LoadgAlice()){
 705     cerr<<"Error occured while l"<<endl;
 706   }
 707   Int_t nall = fLoader->GetNumberOfEvents();
 708   if (nall<=0){
 709     cerr<<"no events available"<<endl;
 710     fEventNr = 0;
 711     return;
 712   }
 713   if (firstEvent+nEvents>nall) {
 714     fEventNr = nall-firstEvent;
 715     cerr<<"restricted number of events availaible"<<endl;
 716   }
 717 }
 718
 719 //////////////////////////////////////////////////////////////
 720 Int_t TPCCmpTr::SetIO()
 721 {
 722   //
 723   //
 724   CreateTreeCmp();
 725   if (!fTreeCmp) return 1;
 726   AliTracker::SetFieldFactor();
 727   fParamTPC = GetTPCParam();
 728   //
 729   if (!ConnectGenTree()) {
 730     cerr<<"Cannot connect tree with generated tracks"<<endl;
 731     return 1;
 732   }
 733   return 0;
 734 }
 735
 736 //////////////////////////////////////////////////////////////
 737
 738 Int_t TPCCmpTr::SetIO(Int_t eventNr)
 739 {
 740   //
 741   //
 742   // SET INPUT
 743   gAlice->GetEvent(eventNr);
 744   fLoader->SetEventNumber(eventNr);
 745   //
 746   AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
 747   tpcl->LoadTracks();
 748   fTreeRecTracks = tpcl->TreeT();
 749   //
 750   return 0;
 751 }
 752
 753
 754
 755 ////////////////////////////////////////////////////////////////////////
 756 void TPCCmpTr::Reset()
 757 {
 758   fEventNr = 0;
 759   fNEvents = 0;
 760   fTreeCmp = 0;
 761   fFnCmp = "cmpTracks.root";
 762   fFileGenTracks = 0;
 763   fDebug = 0;
 764   //
 765   fParamTPC = 0;
 766   fTreeRecTracks = 0;
 767   fTreePoints =0;
 768
 769   fTPCTrack = 0;
 770   fTracks   = 0;
 771   fTrackPoints =0;
 772 }
 773 ////////////////////////////////////////////////////////////////////////
 774 TPCCmpTr::~TPCCmpTr()
 775 {
 776   ;
 777 }
 778 ////////////////////////////////////////////////////////////////////////
 779 Int_t TPCCmpTr::Exec(Int_t nEvents, Int_t firstEventNr)
 780 {
 781   fNEvents = nEvents;
 782   fFirstEventNr = firstEventNr;
 783   return Exec();
 784 }
 785
 786 ////////////////////////////////////////////////////////////////////////
 787 Int_t TPCCmpTr::Exec()
 788 {
 789   TStopwatch timer;
 790   timer.Start();
 791
 792   if (SetIO()==1)
 793     return 1;
 794
 795   fNextTreeGenEntryToRead = 0;
 796   cerr<<"fFirstEventNr, fNEvents: "<<fFirstEventNr<<" "<<fNEvents<<endl;
 797   for (Int_t eventNr = fFirstEventNr; eventNr < fFirstEventNr+fNEvents;
 798        eventNr++) {
 799     SetIO(fEventNr);
 800     fNParticles = gAlice->GetEvent(fEventNr);
 801
 802     fIndexRecTracks = new Int_t[fNParticles*4];  //write at maximum 4 tracks corresponding to particle
 803     fFakeRecTracks = new Int_t[fNParticles];
 804     fMultiRecTracks = new Int_t[fNParticles];
 805     for (Int_t i = 0; i<fNParticles; i++) {
 806       fIndexRecTracks[4*i] = -1;
 807       fIndexRecTracks[4*i+1] = -1;
 808       fIndexRecTracks[4*i+2] = -1;
 809       fIndexRecTracks[4*i+3] = -1;
 810
 811       fFakeRecTracks[i] = 0;
 812       fMultiRecTracks[i] = 0;
 813     }
 814
 815     cout<<"Start to process event "<<fEventNr<<endl;
 816     cout<<"\tfNParticles = "<<fNParticles<<endl;
 817     if (fDebug>2) cout<<"\tStart loop over TreeT"<<endl;
 818     if (TreeTLoop(eventNr)>0) return 1;
 819
 820     if (fDebug>2) cout<<"\tStart loop over tree genTracks"<<endl;
 821     if (TreeGenLoop(eventNr)>0) return 1;
 822     if (fDebug>2) cout<<"\tEnd loop over tree genTracks"<<endl;
 823
 824     delete fTreeRecTracks;
 825     delete [] fIndexRecTracks;
 826     delete [] fFakeRecTracks;
 827     delete [] fMultiRecTracks;
 828   }
 829
 830   CloseOutputFile();
 831
 832   cerr<<"Exec finished"<<endl;
 833   delete gAlice;
 834
 835   timer.Stop();
 836   timer.Print();
 837   return 0;
 838 }
 839 ////////////////////////////////////////////////////////////////////////
 840 Bool_t TPCCmpTr::ConnectGenTree()
 841 {
 842 //
 843 // connect all branches from the genTracksTree
 844 // use the same variables as for the new cmp tree, it may work
 845 //
 846   fFileGenTracks = TFile::Open(fFnGenTracks,"READ");
 847   if (!fFileGenTracks) {
 848     cerr<<"Error in ConnectGenTree: cannot open file "<<fFnGenTracks<<endl;
 849     return kFALSE;
 850   }
 851   fTreeGenTracks = (TTree*)fFileGenTracks->Get("genTracksTree");
 852   if (!fTreeGenTracks) {
 853     cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file "
 854         <<fFnGenTracks<<endl;
 855     return kFALSE;
 856   }
 857   //
 858   fMCInfo = new AliTPCGenInfo;
 859   fMCInfo->fReferences = new TClonesArray("AliTrackReference");
 860   fTreeGenTracks->SetBranchAddress("MC",&fMCInfo);
 861
 862   //
 863   //fTreeGenTracks->SetBranchAddress("fEventNr",&fEventNr);
 864
 865
 866   if (fDebug > 1) {
 867     cout<<"Number of gen. tracks with TR: "<<fTreeGenTracks->GetEntries()<<endl;
 868   }
 869   return kTRUE;
 870 }
 871
 872
 873 ////////////////////////////////////////////////////////////////////////
 874 void TPCCmpTr::CreateTreeCmp()
 875 {
 876   fFileCmp = TFile::Open(fFnCmp,"RECREATE");
 877   if (!fFileCmp) {
 878     cerr<<"Error in CreateTreeCmp: cannot open file "<<fFnCmp<<endl;
 879     return;
 880   }
 881
 882
 883   fTreeCmp    = new TTree("TPCcmpTracks","TPCcmpTracks");
 884   //
 885   fMCInfo = new AliTPCGenInfo;
 886   fMCInfo->fReferences = new TClonesArray("AliTrackReference");
 887   fRecInfo = new AliTPCRecInfo;
 888   //
 889   fTPCTrack = new AliTPCtrack;
 890    //
 891   fTreeCmp->Branch("MC","AliTPCGenInfo",&fMCInfo);
 892   fTreeCmp->Branch("RC","AliTPCRecInfo",&fRecInfo);
 893   fTreeCmp->Branch("fEventNr",&fEventNr,"fEventNr/I");
 894   fTreeCmp->Branch("fTPCTrack","AliTPCtrack",&fTPCTrack);
 895   //
 896   fTreeCmp->AutoSave();
 897
 898 }
 899 ////////////////////////////////////////////////////////////////////////
 900 void TPCCmpTr::CloseOutputFile()
 901 {
 902   if (!fFileCmp) {
 903     cerr<<"File "<<fFnCmp<<" not found as an open file."<<endl;
 904     return;
 905   }
 906   fFileCmp->cd();
 907   fTreeCmp->Write();
 908   delete fTreeCmp;
 909
 910   fFileCmp->Close();
 911   delete fFileCmp;
 912   return;
 913 }
 914 ////////////////////////////////////////////////////////////////////////
 915
 916 TVector3 TPCCmpTr::TR2Local(AliTrackReference *trackRef,
 917                             AliTPCParam *paramTPC) {
 918
 919   Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
 920   Int_t index[4];
 921   paramTPC->Transform0to1(x,index);
 922   paramTPC->Transform1to2(x,index);
 923   return TVector3(x);
 924 }
 925 ////////////////////////////////////////////////////////////////////////
 926
 927 Int_t TPCCmpTr::TreeTLoop(Int_t eventNr)
 928 {
 929   //
 930   // loop over all TPC reconstructed tracks and store info in memory
 931   //
 932
 933
 934   if (!fTreeRecTracks) {
 935     cerr<<"Can't get a tree with TPC rec. tracks  "<<endl;
 936     return 1;
 937   }
 938   //fTreePoints=(TTree*)fFileRecTracks->Get("trackDebug");
 939
 940   Int_t nEntries = (Int_t) fTreeRecTracks->GetEntries();
 941   if (fDebug > 2) cout<<"Event, rec. tracks: "<<eventNr<<" "
 942                       <<nEntries<<endl;
 943   TBranch * br= fTreeRecTracks->GetBranch("tracks");
 944   br->SetAddress(&fTPCTrack);
 945   TBranch *brp = 0;
 946   if (fTreePoints) brp = fTreePoints->GetBranch("debug");
 947
 948   if (fTracks){
 949     fTracks->Delete();
 950     delete fTracks;
 951   }
 952   if (fTrackPoints){
 953     fTrackPoints->Delete();
 954     delete fTrackPoints;
 955     fTrackPoints =0;
 956   }
 957   fTracks      = new TObjArray(nEntries);
 958   if (brp){
 959     fTrackPoints = new TObjArray(nEntries);
 960   }
 961   else fTrackPoints = 0;
 962
 963   //
 964   //load tracks to the memory
 965   for (Int_t i=0; i<nEntries;i++){
 966     AliTPCtrack * track = new AliTPCtrack;
 967     br->SetAddress(&track);
 968     br->GetEntry(i);
 969     fTracks->AddAt(track,i);
 970   }
 971   //
 972   //load track points to the memory
 973   if (brp) for (Int_t i=0; i<nEntries;i++){
 974     TClonesArray * arr = new TClonesArray("AliTPCTrackPoint2");
 975     brp->SetAddress(&arr);
 976     brp->GetEntry(i);
 977     if (arr!=0)
 978       for (Int_t j=0;j<arr->GetEntriesFast();j++){
 979         AliTPCTrackPoint2 * point = (AliTPCTrackPoint2*)arr->UncheckedAt(j);
 980         if (point && point->fID>=0){
 981           fTrackPoints->AddAt(arr,point->fID);
 982           break;
 983         }
 984       }
 985   }
 986   //
 987
 988   //
 989   for (Int_t iEntry=0; iEntry<nEntries;iEntry++){
 990     //br->GetEntry(iEntry);
 991     fTPCTrack = (AliTPCtrack*)fTracks->UncheckedAt(iEntry);
 992     //
 993     Int_t label = fTPCTrack->GetLabel();
 994     Int_t absLabel = abs(label);
 995     if (absLabel < fNParticles) {
 996       //      fIndexRecTracks[absLabel] =  iEntry;
 997       if (label < 0) fFakeRecTracks[absLabel]++;
 998
 999       if (fMultiRecTracks[absLabel]>0){
1000         if (fMultiRecTracks[absLabel]<4)
1001           fIndexRecTracks[absLabel*4+fMultiRecTracks[absLabel]] =  iEntry;
1002       }
1003       else
1004         fIndexRecTracks[absLabel*4] =  iEntry;
1005       fMultiRecTracks[absLabel]++;
1006     }
1007   }
1008
1009   if (fDebug > 2) cerr<<"end of TreeTLoop"<<endl;
1010
1011   return 0;
1012 }
1013 ////////////////////////////////////////////////////////////////////////
1014 Int_t TPCCmpTr::TreeGenLoop(Int_t eventNr)
1015 {
1016 //
1017 // loop over all entries for a given event, find corresponding
1018 // rec. track and store in the fTreeCmp
1019 //
1020
1021   Int_t entry = fNextTreeGenEntryToRead;
1022   Double_t nParticlesTR = fTreeGenTracks->GetEntriesFast();
1023   cerr<<"fNParticles, nParticlesTR, fNextTreeGenEntryToRead: "<<fNParticles<<" "
1024       <<nParticlesTR<<" "<<fNextTreeGenEntryToRead<<endl;
1025   TBranch * branch = fTreeCmp->GetBranch("RC");
1026   while (entry < nParticlesTR) {
1027     fTreeGenTracks->GetEntry(entry);
1028     entry++;
1029     if (fEventNr < eventNr) continue;
1030     if (fEventNr > eventNr) break;
1031     fNextTreeGenEntryToRead = entry-1;
1032     if (fDebug > 2 && fMCInfo->fLabel < 10) {
1033       cerr<<"Fill track with a label "<<fMCInfo->fLabel<<endl;
1034     }
1035
1036     fRecInfo->Reset();
1037
1038
1039     if (fIndexRecTracks[fMCInfo->fLabel*4] >= 0) {
1040       fTPCTrack= (AliTPCtrack*)fTracks->UncheckedAt(fIndexRecTracks[fMCInfo->fLabel*4]);
1041
1042       //      if (nBytes > 0) {
1043       if (fTPCTrack) {
1044         //
1045         TVector3 local = TR2Local(&(fMCInfo->fTrackRef),fParamTPC);
1046         local.GetXYZ(fRecInfo->fTRLocalCoord);
1047         //
1048         // find nearest track if multifound
1049         if (fIndexRecTracks[fMCInfo->fLabel*4]+1){
1050           Float_t dz = TMath::Abs(local.Z()-fTPCTrack->GetZ());
1051           for (Int_t i=1;i<4;i++){
1052             if (fIndexRecTracks[fMCInfo->fLabel*4+i]>=0){
1053               AliTPCtrack * track = (AliTPCtrack*)fTracks->UncheckedAt(fIndexRecTracks[fMCInfo->fLabel*4+i]);
1054               if  (TMath::Abs(local.Z()-track->GetZ())<dz)
1055                 fTPCTrack = track;
1056             }
1057           }
1058         }
1059         fRecInfo->fTP=0;
1060         if (fTrackPoints){
1061           Int_t id = fTPCTrack->GetUniqueID();
1062           if (fTrackPoints->UncheckedAt(id)){
1063             fRecInfo->fTP = (TClonesArray*)fTrackPoints->UncheckedAt(id);
1064             //      fTrackPoints->AddAt(0,id);   //not owner anymore
1065           }
1066         }
1067         fRecInfo->fTPCTrack =*fTPCTrack;
1068         fRecInfo->fReconstructed = 1;
1069         fRecInfo->fFake = fFakeRecTracks[fMCInfo->fLabel];
1070         fRecInfo->fMultiple = fMultiRecTracks[fMCInfo->fLabel];
1071         fRecInfo->fdEdx = fTPCTrack->GetdEdx();
1072         //
1073         fRecInfo->fTPCTrack.PropagateTo(local.X());
1074         Double_t par[5];
1075         //
1076         Double_t localX = local.X();
1077         fTPCTrack->GetExternalParameters(localX,par);
1078         fRecInfo->fRecPhi=TMath::ASin(par[2]) + fTPCTrack->GetAlpha();
1079         if (fRecInfo->fRecPhi<0) fRecInfo->fRecPhi+=2*kPI;
1080         if (fRecInfo->fRecPhi>=2*kPI) fRecInfo->fRecPhi-=2*kPI;
1081 //        fRecInfo->fRecPhi = (fRecInfo->fRecPhi)*kRaddeg;
1082         fRecInfo->fLambda = TMath::ATan(par[3]);
1083         fRecInfo->fRecPt_1 = TMath::Abs(par[4]);
1084       }
1085
1086     }
1087
1088     branch->SetAddress(&fRecInfo); // set all pointers
1089     fTreeCmp->Fill();
1090
1091   }
1092   fTreeCmp->AutoSave();
1093   fTracks->Delete();
1094   fTPCTrack =0;
1095   if (fTrackPoints){
1096     fTrackPoints->Delete();
1097     delete fTrackPoints;
1098     fTrackPoints =0;
1099   }
1100   if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl;
1101
1102   return 0;
1103 }
1104 ////////////////////////////////////////////////////////////////////////