PWG1/AliGenInfoMaker.cxx

   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
  20 Origin: marian.ivanov@cern.ch
  21
  22 Generate complex MC information - used for Comparison later on
  23 How to use it?
  24
  25
  26 ---Usage outside of the analysis framework
  27
  28 gSystem->Load("libANALYSIS.so")
  29 gSystem->Load("libPWG1.so")
  30 AliGenInfoMaker *t = new AliGenInfoMaker("galice.root","genTracks.root",0,0)
  31 t->Exec();
  32
  33 */
  34
  35 #if !defined(__CINT__) || defined(__MAKECINT__)
  36 #include <stdio.h>
  37 #include <string.h>
  38 //ROOT includes
  39 #include "TROOT.h"
  40 #include "Rtypes.h"
  41 #include "TFile.h"
  42 #include "TTree.h"
  43 #include "TStopwatch.h"
  44 #include "TParticle.h"
  45
  46 //ALIROOT includes
  47 #include "AliMCEvent.h"
  48 #include "AliMCEventHandler.h"
  49
  50 #include "AliRun.h"
  51 #include "AliStack.h"
  52 #include "AliSimDigits.h"
  53 #include "AliTPCParam.h"
  54 #include "AliTPC.h"
  55 #include "AliTPCLoader.h"
  56 #include "AliTrackReference.h"
  57 #include "AliTPCParamSR.h"
  58 #include "AliTrackPointArray.h"
  59
  60 #endif
  61 #include "AliMCInfo.h"
  62 #include "AliGenV0Info.h"
  63 #include "AliGenKinkInfo.h"
  64 #include "AliGenInfoMaker.h"
  65 //
  66 //
  67
  68 ClassImp(AliGenInfoMaker)
  69
  70
  71
  72
  73
  74
  75
  76 ////////////////////////////////////////////////////////////////////////
  77 AliGenInfoMaker::AliGenInfoMaker():
  78   fGenTracksArray(0),          //clones array with filtered particles
  79   fGenKinkArray(0),            //clones array with filtered Kinks
  80   fGenV0Array(0),              //clones array with filtered V0s
  81   fDebug(0),                   //! debug flag
  82   fEventNr(0),                 //! current event number
  83   fLabel(0),                   //! track label
  84   fNEvents(0),                 //! number of events to process
  85   fFirstEventNr(0),            //! first event to process
  86   fNParticles(0),              //! number of particles in TreeK
  87   fTreeGenTracks(0),           //! output tree with generated tracks
  88   fTreeKinks(0),               //!  output tree with Kinks
  89   fTreeV0(0),                  //!  output tree with V0
  90   fFileGenTracks(0),           //! output file with stored fTreeGenTracks
  91   fLoader(0),                  //! pointer to the run loader
  92   fTreeD(0),                   //! current tree with digits
  93   fTreeTR(0),                  //! current tree with TR
  94   fStack(0),                   //! current stack
  95   fGenInfo(0),                 //! array with pointers to gen info
  96   fNInfos(0),                  //! number of tracks with infos
  97   fParamTPC(0),                //! AliTPCParam
  98   fTPCPtCut(0.1),              //  TPC pt cut
  99   fITSPtCut(0.1),              //  ITS pt cut
 100   fTRDPtCut(0.1),              //  TRD pt cut
 101   fTOFPtCut(0.1)               //  TOF pt cut
 102 {
 103   sprintf(fFnRes,"%s","genTracks.root");
 104 }
 105
 106
 107
 108
 109 ////////////////////////////////////////////////////////////////////////
 110 AliGenInfoMaker::AliGenInfoMaker(const char * fnGalice, const char* fnRes,
 111                                  Int_t nEvents, Int_t firstEvent):
 112   fGenTracksArray(0),          //clones array with filtered particles
 113   fGenKinkArray(0),            //clones array with filtered Kinks
 114   fGenV0Array(0),              //clones array with filtered V0s
 115   fDebug(0),                   //! debug flag
 116   fEventNr(0),                 //! current event number
 117   fLabel(0),                   //! track label
 118   fNEvents(0),                 //! number of events to process
 119   fFirstEventNr(0),            //! first event to process
 120   fNParticles(0),              //! number of particles in TreeK
 121   fTreeGenTracks(0),           //! output tree with generated tracks
 122   fTreeKinks(0),               //!  output tree with Kinks
 123   fTreeV0(0),                  //!  output tree with V0
 124   fFileGenTracks(0),           //! output file with stored fTreeGenTracks
 125   fLoader(0),                  //! pointer to the run loader
 126   fTreeD(0),                   //! current tree with digits
 127   fTreeTR(0),                  //! current tree with TR
 128   fStack(0),                   //! current stack
 129   fGenInfo(0),                 //! array with pointers to gen info
 130   fNInfos(0),                  //! number of tracks with infos
 131   fParamTPC(0),                //! AliTPCParam
 132   fTPCPtCut(0.1),
 133   fITSPtCut(0.1),
 134   fTRDPtCut(0.1),
 135   fTOFPtCut(0.1)
 136 {
 137   //
 138   //
 139   //
 140   fFirstEventNr = firstEvent;
 141   fEventNr = firstEvent;
 142   fNEvents = nEvents;
 143   sprintf(fFnRes,"%s",fnRes);
 144   //
 145   //
 146   //
 147   fLoader = AliRunLoader::Open(fnGalice);
 148   if (gAlice){
 149     delete gAlice->GetRunLoader();
 150     delete gAlice;
 151     gAlice = 0x0;
 152   }
 153   if (fLoader->LoadgAlice()){
 154     cerr<<"Error occured while l"<<endl;
 155   }
 156   Int_t nall = fLoader->GetNumberOfEvents();
 157   if (nEvents==0) {
 158     nEvents =nall;
 159     fNEvents=nall;
 160     fFirstEventNr=0;
 161   }
 162
 163   if (nall<=0){
 164     cerr<<"no events available"<<endl;
 165     fEventNr = 0;
 166     return;
 167   }
 168   if (firstEvent+nEvents>nall) {
 169     fEventNr = nall-firstEvent;
 170     cerr<<"restricted number of events availaible"<<endl;
 171   }
 172 }
 173
 174
 175
 176
 177 ////////////////////////////////////////////////////////////////////////
 178 AliGenInfoMaker::~AliGenInfoMaker()
 179 {
 180   //
 181   // Destructor
 182   //
 183
 184   if (fLoader){
 185     fLoader->UnloadgAlice();
 186     gAlice = 0;
 187     delete fLoader;
 188   }
 189 }
 190
 191
 192 AliMCInfo * AliGenInfoMaker::MakeInfo(UInt_t i)
 193 {
 194   //
 195   // Make info structure for given particle index
 196   //
 197   if (i<fNParticles) {
 198     if (fGenInfo[i]) return  fGenInfo[i];
 199     fGenInfo[i] = new AliMCInfo;
 200     fNInfos++;
 201     return fGenInfo[i];
 202   }
 203   else
 204     return 0;
 205 }
 206
 207
 208 Int_t  AliGenInfoMaker::ProcessEvent(AliMCEventHandler* mcinfo){
 209   //
 210   // Process MC info from the task
 211   //
 212   if (!fParamTPC) {
 213     SetIO();
 214     fParamTPC = GetTPCParam();
 215   }
 216   fStack  = mcinfo->MCEvent()->Stack();
 217   fTreeTR = mcinfo->TreeTR();
 218   fTreeD  = 0;
 219   // array with preprocessedprocessed information
 220   fGenTracksArray = new TObjArray;
 221   fGenKinkArray   = new TObjArray;
 222   fGenV0Array     = new TObjArray;
 223   //
 224   ProcessEvent();
 225   fEventNr++;
 226   return 0;
 227 }
 228
 229
 230
 231 Int_t AliGenInfoMaker::ProcessEvent(){
 232   //
 233   // Process Event
 234   //
 235   fNParticles = fStack->GetNtrack();
 236   //
 237   fGenInfo = new AliMCInfo*[fNParticles];
 238   for (UInt_t i = 0; i<fNParticles; i++) {
 239     fGenInfo[i]=0;
 240   }
 241   //
 242   cout<<"Start to process event "<<fEventNr<<endl;
 243   cout<<"\tfNParticles = "<<fNParticles<<endl;
 244   if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeTR"<<endl;
 245   if (TreeTRLoop()>0) return 1;
 246   //
 247   if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeD"<<endl;
 248   if (TreeDLoop()>0) return 1;
 249   //
 250   if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeK"<<endl;
 251   if (TreeKLoop()>0) return 1;
 252   if (fDebug>2) cout<<"\tEnd loop over TreeK"<<endl;
 253   //
 254   //  if (BuildKinkInfo()>0) return 1;
 255   if (BuildV0Info()>0) return 1;
 256   if (fDebug>2) cout<<"\tEnd loop over TreeK"<<endl;
 257   //
 258   for (UInt_t i = 0; i<fNParticles; i++) {
 259     if (fGenInfo[i]) delete fGenInfo[i];
 260   }
 261   delete []fGenInfo;
 262   return 0;
 263 }
 264
 265
 266
 267
 268
 269
 270
 271 Int_t  AliGenInfoMaker::SetIO()
 272 {
 273   //
 274   // Set IO for given event
 275   //
 276   CreateTreeGenTracks();
 277   if (!fTreeGenTracks) return 1;
 278
 279   fParamTPC = GetTPCParam();
 280   //
 281   return 0;
 282 }
 283
 284 ////////////////////////////////////////////////////////////////////////
 285 Int_t AliGenInfoMaker::SetIO(Int_t eventNr)
 286 {
 287   //
 288   //
 289   // SET INPUT
 290   fLoader->SetEventNumber(eventNr);
 291   //
 292   fLoader->LoadHeader();
 293   fLoader->LoadKinematics();
 294   fStack = fLoader->Stack();
 295   //
 296   fLoader->LoadTrackRefs();
 297   fLoader->LoadHits();
 298   fTreeTR = fLoader->TreeTR();
 299   //
 300   AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
 301   tpcl->LoadDigits();
 302   fTreeD = tpcl->TreeD();
 303   return 0;
 304 }
 305
 306 Int_t AliGenInfoMaker::CloseIOEvent()
 307 {
 308   //
 309   // Close IO for current event
 310   //
 311   fLoader->UnloadHeader();
 312   fLoader->UnloadKinematics();
 313   fLoader->UnloadTrackRefs();
 314   AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
 315   tpcl->UnloadDigits();
 316   return 0;
 317 }
 318
 319 Int_t AliGenInfoMaker::CloseIO()
 320 {
 321   fLoader->UnloadgAlice();
 322   return 0;
 323 }
 324
 325
 326
 327 ////////////////////////////////////////////////////////////////////////
 328 Int_t AliGenInfoMaker::Exec()
 329 {
 330   //
 331   // Make a comparision MC tree
 332   // Connect MC information -TPArticle - AliTrackRefernces ...
 333   //
 334   TStopwatch timer;
 335   timer.Start();
 336   Int_t status =SetIO();
 337   if (status>0) return status;
 338   //
 339   for (fEventNr = fFirstEventNr; fEventNr < fFirstEventNr+fNEvents;
 340        fEventNr++) {
 341     SetIO(fEventNr);
 342     ProcessEvent();
 343     CloseIOEvent();
 344   }
 345   //
 346   CloseIO();
 347   CloseOutputFile();
 348   cerr<<"Exec finished"<<endl;
 349   timer.Stop();
 350   timer.Print();
 351   return 0;
 352 }
 353
 354
 355
 356
 357
 358
 359 ////////////////////////////////////////////////////////////////////////
 360 void AliGenInfoMaker::CreateTreeGenTracks()
 361 {
 362   //
 363   //
 364   //
 365   fFileGenTracks = TFile::Open(fFnRes,"RECREATE");
 366   if (!fFileGenTracks) {
 367     cerr<<"Error in CreateTreeGenTracks: cannot open file "<<fFnRes<<endl;
 368     return;
 369   }
 370   fTreeGenTracks = new TTree("genTracksTree","genTracksTree");
 371   AliMCInfo * info = new AliMCInfo;
 372   fTreeGenTracks->Branch("MC","AliMCInfo",&info);
 373   delete info;
 374   //
 375   AliGenKinkInfo *kinkinfo = new AliGenKinkInfo;
 376   fTreeKinks = new TTree("genKinksTree","genKinksTree");
 377   fTreeKinks->Branch("MC","AliGenKinkInfo",&kinkinfo);
 378   delete kinkinfo;
 379   //
 380   AliGenV0Info *v0info = new AliGenV0Info;
 381   fTreeV0 = new TTree("genV0Tree","genV0Tree");
 382   fTreeV0->Branch("MC","AliGenV0Info",&v0info);
 383   delete v0info;
 384   //
 385   //
 386   //
 387   fTreeGenTracks->AutoSave();
 388   fTreeKinks->AutoSave();
 389   fTreeV0->AutoSave();
 390 }
 391
 392
 393
 394 ////////////////////////////////////////////////////////////////////////
 395 void AliGenInfoMaker::CloseOutputFile()
 396 {
 397   //
 398   // Close Output files
 399   //
 400   if (!fFileGenTracks) {
 401     cerr<<"File "<<fFnRes<<" not found as an open file."<<endl;
 402     return;
 403   }
 404   fFileGenTracks->cd();
 405   fTreeGenTracks->Write();
 406   delete fTreeGenTracks;
 407   fTreeKinks->Write();
 408   delete fTreeKinks;
 409   fTreeV0->Write();
 410   delete fTreeV0;
 411
 412   fFileGenTracks->Close();
 413   delete fFileGenTracks;
 414   return;
 415 }
 416
 417 ////////////////////////////////////////////////////////////////////////
 418 Int_t AliGenInfoMaker::TreeKLoop()
 419 {
 420 //
 421 // open the file with treeK
 422 // loop over all entries there and save information about some tracks
 423 //
 424
 425   AliStack * stack = fStack;
 426   if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
 427
 428   if (fDebug > 0) {
 429     cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
 430         <<fEventNr<<endl;
 431   }
 432   Int_t  ipdg = 0;
 433   TParticlePDG *ppdg = 0;
 434   // not all generators give primary vertex position. Take the vertex
 435   // of the particle 0 as primary vertex.
 436   TDatabasePDG  pdg; //get pdg table
 437
 438   TParticle *particle = stack->ParticleFromTreeK(0);
 439   fVPrim[0] = particle->Vx();
 440   fVPrim[1] = particle->Vy();
 441   fVPrim[2] = particle->Vz();
 442   Int_t accepted =0;
 443   for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
 444     // load only particles with TR
 445     AliMCInfo * info = GetInfo(iParticle);
 446     if (!info) continue;
 447     //////////////////////////////////////////////////////////////////////
 448     info->fLabel = iParticle;
 449     //
 450     info->fParticle = *(stack->Particle(iParticle));
 451     info->fVDist[0] = info->fParticle.Vx()-fVPrim[0];
 452     info->fVDist[1] = info->fParticle.Vy()-fVPrim[1];
 453     info->fVDist[2] = info->fParticle.Vz()-fVPrim[2];
 454     info->fVDist[3] = TMath::Sqrt(info->fVDist[0]*info->fVDist[0]+
 455                                   info->fVDist[1]*info->fVDist[1]+info->fVDist[2]*info->fVDist[2]);
 456     //
 457     //
 458     ipdg = info->fParticle.GetPdgCode();
 459     info->fPdg = ipdg;
 460     ppdg = pdg.GetParticle(ipdg);
 461     info->fEventNr = fEventNr;
 462     info->Update();
 463     if (fGenTracksArray){
 464       fGenTracksArray->AddLast(info->Clone());
 465     }
 466     accepted++;
 467   }
 468   //
 469   // write the results to the tree - if specified
 470   //
 471   TBranch * br = fTreeGenTracks->GetBranch("MC");
 472   for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
 473     // load only particles with TR
 474     AliMCInfo * info = GetInfo(iParticle);
 475     if (!info) continue;
 476     //////////////////////////////////////////////////////////////////////
 477     br->SetAddress(&info);
 478     fTreeGenTracks->Fill();
 479   }
 480   fTreeGenTracks->AutoSave();
 481   //
 482   //
 483   if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl;
 484   return 0;
 485 }
 486
 487
 488 ////////////////////////////////////////////////////////////////////////////////////
 489 Int_t  AliGenInfoMaker::BuildKinkInfo()
 490 {
 491   //
 492   // Build tree with Kink Information
 493   //
 494   AliStack * stack = fStack;
 495   if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
 496
 497   if (fDebug > 0) {
 498     cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
 499         <<fEventNr<<endl;
 500   }
 501   //  Int_t  ipdg = 0;
 502   //TParticlePDG *ppdg = 0;
 503   // not all generators give primary vertex position. Take the vertex
 504   // of the particle 0 as primary vertex.
 505   TDatabasePDG  pdg; //get pdg table
 506   //thank you very much root for this
 507   TBranch * br = fTreeKinks->GetBranch("MC");
 508   //
 509   AliGenKinkInfo * kinkinfo = new AliGenKinkInfo;
 510   //
 511   for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
 512     // load only particles with TR
 513     AliMCInfo * info = GetInfo(iParticle);
 514     if (!info) continue;
 515     if (info->fCharge==0) continue;
 516     if (info->fTRdecay.R()>500)  continue;  //R cut - decay outside barrel
 517     TParticle & particle = info->fParticle;
 518     Int_t first = particle.GetDaughter(0);
 519     if (first ==0) continue;
 520
 521     Int_t last = particle.GetDaughter(1);
 522     if (last ==0) last=first;
 523     AliMCInfo * info2 =  0;
 524     AliMCInfo * dinfo =  0;
 525
 526     for (Int_t id2=first;id2<=last;id2++){
 527       info2 = GetInfo(id2);
 528       if (!info2) continue;
 529       TParticle &p2 = info2->fParticle;
 530       Double_t r = TMath::Sqrt(p2.Vx()*p2.Vx()+p2.Vy()*p2.Vy());
 531       if ( TMath::Abs(info->fTRdecay.R()-r)>0.01) continue;
 532       if (!(p2.GetPDG())) continue;
 533       dinfo =info2;
 534
 535       kinkinfo->SetInfoMother(*info);
 536       kinkinfo->SetInfoDaughter(*dinfo);
 537       if (kinkinfo->GetMinus().GetParticle().GetPDG()==0 ||  kinkinfo->GetPlus().GetParticle().GetPDG()==0) continue;
 538       kinkinfo->Update();
 539       br->SetAddress(&kinkinfo);
 540       fTreeKinks->Fill();
 541     }
 542   }
 543   fTreeKinks->AutoSave();
 544   if (fDebug > 2) cerr<<"end of Kink Loop"<<endl;
 545   return 0;
 546 }
 547
 548
 549 ////////////////////////////////////////////////////////////////////////////////////
 550 Int_t  AliGenInfoMaker::BuildV0Info()
 551 {
 552   //
 553   // Build tree with V0 Information
 554   //
 555   AliStack * stack = fStack;
 556   if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
 557
 558   if (fDebug > 0) {
 559     cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
 560         <<fEventNr<<endl;
 561   }
 562   //
 563   TDatabasePDG  pdg; //get pdg table
 564   //thank you very much root for this
 565   TBranch * br = fTreeV0->GetBranch("MC");
 566   //
 567   AliGenV0Info * v0info = new AliGenV0Info;
 568   //
 569   //
 570   for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
 571     TParticle * mParticle = fStack->Particle(iParticle);
 572     if (!mParticle) continue;
 573     if (mParticle->GetPDG()==0) continue;
 574     if (mParticle->GetPDG()->Charge()!=0) continue;  //only neutral particle
 575     //
 576     Int_t first = mParticle->GetDaughter(0);
 577     Int_t last  = mParticle->GetDaughter(1);
 578     if (first-last==0) continue;
 579     //
 580     AliMCInfo * info0 = GetInfo(first);
 581     AliMCInfo * info1 = GetInfo(first+1);
 582     if (info0 && info1){
 583       if (info0->GetPdg()==0) continue;
 584       if (info1->GetPdg()==0) continue;
 585       v0info->SetInfoP(*info0);
 586       v0info->SetInfoM(*info1);
 587       v0info->SetMother(*mParticle);
 588       v0info->Update(fVPrim);
 589       if (fGenV0Array) fGenV0Array->AddLast(v0info);
 590       br->SetAddress(&v0info);
 591       fTreeV0->Fill();
 592     }
 593   }
 594   fTreeV0->AutoSave();
 595   if (fDebug > 2) cerr<<"end of V0 Loop"<<endl;
 596   return 0;
 597 }
 598
 599
 600
 601
 602
 603 ////////////////////////////////////////////////////////////////////////
 604 Int_t AliGenInfoMaker::TreeDLoop()
 605 {
 606   //
 607   // open the file with treeD
 608   // loop over all entries there and save information about some tracks
 609   //
 610   if (!fTreeD) return 0;
 611   Int_t nInnerSector = fParamTPC->GetNInnerSector();
 612   Int_t rowShift = 0;
 613   Int_t zero=fParamTPC->GetZeroSup()+6;
 614   //
 615   //
 616   AliSimDigits digitsAddress, *digits=&digitsAddress;
 617   fTreeD->GetBranch("Segment")->SetAddress(&digits);
 618
 619   Int_t sectorsByRows=(Int_t)fTreeD->GetEntries();
 620   if (fDebug > 1) cout<<"\tsectorsByRows = "<<sectorsByRows<<endl;
 621   for (Int_t i=0; i<sectorsByRows; i++) {
 622     if (!fTreeD->GetEvent(i)) continue;
 623     Int_t sec,row;
 624     fParamTPC->AdjustSectorRow(digits->GetID(),sec,row);
 625     if (fDebug > 1) cout<<sec<<' '<<row<<"                          \r";
 626     // here I expect that upper sectors follow lower sectors
 627     if (sec > nInnerSector) rowShift = fParamTPC->GetNRowLow();
 628     //
 629     digits->ExpandTrackBuffer();
 630     digits->First();
 631     do {
 632       Int_t iRow=digits->CurrentRow();
 633       Int_t iColumn=digits->CurrentColumn();
 634       Short_t digitValue = digits->CurrentDigit();
 635       if (digitValue >= zero) {
 636         Int_t label;
 637         for (Int_t j = 0; j<3; j++) {
 638           //      label = digits->GetTrackID(iRow,iColumn,j);
 639           label = digits->GetTrackIDFast(iRow,iColumn,j)-2;
 640           if (label >= (Int_t)fNParticles) { //don't label from bakground event
 641             continue;
 642           }
 643           if (label >= 0 && label <= (Int_t)fNParticles) {
 644             if (fDebug > 6 ) {
 645               cout<<"Inner loop: sector, iRow, iColumn, label, value, row "
 646                   <<sec<<" "
 647                   <<iRow<<" "<<iColumn<<" "<<label<<" "<<digitValue
 648                   <<" "<<row<<endl;
 649             }
 650             AliMCInfo * info = GetInfo(label);
 651             if (info){
 652               info->fTPCRow.SetRow(row+rowShift);
 653             }
 654           }
 655         }
 656       }
 657     } while (digits->Next());
 658   }
 659
 660   if (fDebug > 2) cerr<<"end of TreeDLoop"<<endl;
 661   return 0;
 662 }
 663
 664
 665
 666
 667 ////////////////////////////////////////////////////////////////////////
 668 Int_t AliGenInfoMaker::TreeTRLoop()
 669 {
 670   //
 671   // loop over TrackReferences and store the first one for each track
 672   //
 673   TTree * treeTR = fTreeTR;
 674   Int_t nPrimaries = (Int_t) treeTR->GetEntries();
 675   if (fDebug > 1) cout<<"There are "<<nPrimaries<<" entries in TreeTR"<<endl;
 676   //
 677   //
 678   //
 679   TClonesArray* runArrayTR = new TClonesArray("AliTrackReference");
 680   //
 681   if (treeTR->GetBranch("TrackReferences")) treeTR->GetBranch("TrackReferences")->SetAddress(&runArrayTR);
 682   //
 683   //
 684   //
 685   for (Int_t ipart = 0; ipart<fStack->GetNtrack(); ipart++) {
 686     TParticle * part = fStack->Particle(ipart);
 687     if (!part) continue;
 688     if (part->GetPDG()==0) continue;
 689     if (part->GetPDG()->Charge()==0) continue;
 690
 691     if (part->Pt()<fITSPtCut) continue;
 692     if (part->R()>250.) continue;
 693     if (TMath::Abs(part->Vz())>250.) continue;
 694     if (TMath::Abs(part->Pz()/part->Pt())>2.5) continue;
 695     //AliMCInfo * info = ;
 696     MakeInfo(ipart);
 697   }
 698   //
 699   //
 700
 701   for (Int_t iPrimPart = 0; iPrimPart<nPrimaries; iPrimPart++) {
 702     treeTR->GetEntry(iPrimPart);
 703     //
 704     // gettrack references
 705     //
 706     for (Int_t iTrackRef = 0; iTrackRef < runArrayTR->GetEntriesFast(); iTrackRef++) {
 707       AliTrackReference *trackRef = (AliTrackReference*)runArrayTR->At(iTrackRef);
 708       //
 709       Int_t label = trackRef->GetTrack();
 710       AliMCInfo * cinfo = GetInfo(label);
 711       if (cinfo) cinfo->CalcTPCrows(runArrayTR);
 712       //
 713       // TPC
 714       //
 715       if (trackRef->DetectorId()== AliTrackReference::kTPC){
 716         //
 717         Int_t label = trackRef->GetTrack();
 718         AliMCInfo * info = GetInfo(label);
 719         if (!info && trackRef->Pt()<fTPCPtCut) continue;
 720         if (!info) {
 721           info = MakeInfo(label);
 722         }
 723         if (!info) continue;
 724         info->fPrimPart =  iPrimPart;
 725         TClonesArray & arr = *(info->fTPCReferences);
 726         new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);
 727       }
 728       //
 729       // ITS
 730       //
 731       if (trackRef->DetectorId()== AliTrackReference::kITS){
 732         //
 733         Int_t label = trackRef->GetTrack();
 734         AliMCInfo * info = GetInfo(label);
 735         if (!info && trackRef->Pt()<fITSPtCut) continue;
 736         if (!info) info = MakeInfo(label);
 737         if (!info) continue;
 738         info->fPrimPart =  iPrimPart;
 739         TClonesArray & arr = *(info->fITSReferences);
 740         new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);
 741       }
 742       //
 743       // TRD
 744       //
 745       if (trackRef->DetectorId()== AliTrackReference::kTRD){
 746         //
 747         Int_t label = trackRef->GetTrack();
 748         AliMCInfo * info = GetInfo(label);
 749         if (!info&&trackRef->P()<fTRDPtCut) continue;
 750         if (!info) info = MakeInfo(label);
 751         if (!info) continue;
 752         info->fPrimPart =  iPrimPart;
 753         TClonesArray & arr = *(info->fTRDReferences);
 754         new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);
 755       }
 756       //
 757       // TOF
 758       //
 759       if (trackRef->DetectorId()== AliTrackReference::kTOF){
 760         //
 761         Int_t label = trackRef->GetTrack();
 762         AliMCInfo * info = GetInfo(label);
 763         if (!info&&trackRef->P()<fTOFPtCut) continue;
 764         if (!info) info = MakeInfo(label);
 765         if (!info) continue;
 766         info->fPrimPart =  iPrimPart;
 767         TClonesArray & arr = *(info->fTOFReferences);
 768         new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);
 769       }
 770       //
 771       // decay case
 772       //
 773       if (trackRef->DetectorId()== AliTrackReference::kDisappeared){
 774         //
 775         AliMCInfo * info = GetInfo(label);
 776         if (!info) continue;
 777         info->fTRdecay = *trackRef;
 778       }
 779     }
 780   }
 781   //
 782   //
 783   return 0;
 784 }
 785
 786 ////////////////////////////////////////////////////////////////////////
 787 Float_t AliGenInfoMaker::TR2LocalX(AliTrackReference *trackRef,
 788                                     AliTPCParam *paramTPC) const {
 789
 790   //
 791   // Transformation from Gloabal to local tacking system
 792   //
 793   Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
 794   Int_t index[4];
 795   paramTPC->Transform0to1(x,index);
 796   paramTPC->Transform1to2(x,index);
 797   return x[0];
 798 }
 799 ////////////////////////////////////////////////////////////////////////
 800
 801
 802
 803 AliTPCParam * AliGenInfoMaker::GetTPCParam(){
 804   //
 805   // create default TPC parameters
 806   //
 807   AliTPCParamSR * par = new AliTPCParamSR;
 808   par->Update();
 809   return par;
 810 }
 811
 812