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