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