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