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 gSystem->Load("libPWG1.so")
  26 AliGenInfoMaker *t = new AliGenInfoMaker("galice.root","genTracks.root",0,0)
  27 t->Exec();
  28
  29 */
  30
  31 #if !defined(__CINT__) || defined(__MAKECINT__)
  32 #include <stdio.h>
  33 #include <string.h>
  34 //ROOT includes
  35 #include "TROOT.h"
  36 #include "Rtypes.h"
  37 #include "TFile.h"
  38 #include "TTree.h"
  39 #include "TChain.h"
  40 #include "TCut.h"
  41 #include "TString.h"
  42 #include "TStopwatch.h"
  43 #include "TParticle.h"
  44 #include "TSystem.h"
  45 #include "TCanvas.h"
  46 #include "TGeometry.h"
  47 #include "TPolyLine3D.h"
  48
  49 //ALIROOT includes
  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 "AliDetector.h"
  57 #include "AliTrackReference.h"
  58 #include "AliTPCParamSR.h"
  59 #include "AliTracker.h"
  60 #include "AliMagF.h"
  61 #include "AliHelix.h"
  62 #include "AliTrackPointArray.h"
  63
  64 #endif
  65 #include "AliMCInfo.h"
  66 #include "AliGenV0Info.h"
  67 #include "AliGenKinkInfo.h"
  68 #include "AliGenInfoMaker.h"
  69 //
  70 //
  71
  72 ClassImp(AliGenInfoMaker)
  73
  74
  75
  76
  77
  78
  79
  80 ////////////////////////////////////////////////////////////////////////
  81 AliGenInfoMaker::AliGenInfoMaker():
  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   fTreeHitLines(0),            //! tree with hit lines
  92   fFileGenTracks(0),           //! output file with stored fTreeGenTracks
  93   fLoader(0),                  //! pointer to the run loader
  94   fTreeD(0),                   //! current tree with digits
  95   fTreeTR(0),                  //! current tree with TR
  96   fStack(0),                   //! current stack
  97   fGenInfo(0),                 //! array with pointers to gen info
  98   fNInfos(0),                  //! number of tracks with infos
  99   fParamTPC(0),                //! AliTPCParam
 100   fTPCPtCut(0.03),
 101   fITSPtCut(0.1),
 102   fTRDPtCut(0.1),
 103   fTOFPtCut(0.1)
 104 {
 105 }
 106
 107 ////////////////////////////////////////////////////////////////////////
 108 AliGenInfoMaker::AliGenInfoMaker(const char * fnGalice, const char* fnRes,
 109                                  Int_t nEvents, Int_t firstEvent):
 110   fDebug(0),                   //! debug flag
 111   fEventNr(0),                 //! current event number
 112   fLabel(0),                   //! track label
 113   fNEvents(0),                 //! number of events to process
 114   fFirstEventNr(0),            //! first event to process
 115   fNParticles(0),              //! number of particles in TreeK
 116   fTreeGenTracks(0),           //! output tree with generated tracks
 117   fTreeKinks(0),               //!  output tree with Kinks
 118   fTreeV0(0),                  //!  output tree with V0
 119   fTreeHitLines(0),            //! tree with hit lines
 120   fFileGenTracks(0),           //! output file with stored fTreeGenTracks
 121   fLoader(0),                  //! pointer to the run loader
 122   fTreeD(0),                   //! current tree with digits
 123   fTreeTR(0),                  //! current tree with TR
 124   fStack(0),                   //! current stack
 125   fGenInfo(0),                 //! array with pointers to gen info
 126   fNInfos(0),                  //! number of tracks with infos
 127   fParamTPC(0),                //! AliTPCParam
 128   fTPCPtCut(0.03),
 129   fITSPtCut(0.1),
 130   fTRDPtCut(0.1),
 131   fTOFPtCut(0.1)
 132
 133 {
 134   //
 135   //
 136   //
 137   fFirstEventNr = firstEvent;
 138   fEventNr = firstEvent;
 139   fNEvents = nEvents;
 140   sprintf(fFnRes,"%s",fnRes);
 141   //
 142   fLoader = AliRunLoader::Open(fnGalice);
 143   if (gAlice){
 144     delete gAlice->GetRunLoader();
 145     delete gAlice;
 146     gAlice = 0x0;
 147   }
 148   if (fLoader->LoadgAlice()){
 149     cerr<<"Error occured while l"<<endl;
 150   }
 151   Int_t nall = fLoader->GetNumberOfEvents();
 152   if (nEvents==0) {
 153     nEvents =nall;
 154     fNEvents=nall;
 155     fFirstEventNr=0;
 156   }
 157
 158   if (nall<=0){
 159     cerr<<"no events available"<<endl;
 160     fEventNr = 0;
 161     return;
 162   }
 163   if (firstEvent+nEvents>nall) {
 164     fEventNr = nall-firstEvent;
 165     cerr<<"restricted number of events availaible"<<endl;
 166   }
 167   AliMagF * magf = gAlice->Field();
 168   AliTracker::SetFieldMap(magf,0);
 169 }
 170
 171
 172 AliMCInfo * AliGenInfoMaker::MakeInfo(UInt_t i)
 173 {
 174   //
 175   if (i<fNParticles) {
 176     if (fGenInfo[i]) return  fGenInfo[i];
 177     fGenInfo[i] = new AliMCInfo;
 178     fNInfos++;
 179     return fGenInfo[i];
 180   }
 181   else
 182     return 0;
 183 }
 184
 185 ////////////////////////////////////////////////////////////////////////
 186 AliGenInfoMaker::~AliGenInfoMaker()
 187 {
 188
 189   if (fLoader){
 190     fLoader->UnloadgAlice();
 191     gAlice = 0;
 192     delete fLoader;
 193   }
 194 }
 195
 196 Int_t  AliGenInfoMaker::SetIO()
 197 {
 198   //
 199   //
 200   CreateTreeGenTracks();
 201   if (!fTreeGenTracks) return 1;
 202   //  AliTracker::SetFieldFactor();
 203
 204   fParamTPC = GetTPCParam();
 205   //
 206   return 0;
 207 }
 208
 209 ////////////////////////////////////////////////////////////////////////
 210 Int_t AliGenInfoMaker::SetIO(Int_t eventNr)
 211 {
 212   //
 213   //
 214   // SET INPUT
 215   fLoader->SetEventNumber(eventNr);
 216   //
 217   fLoader->LoadHeader();
 218   fLoader->LoadKinematics();
 219   fStack = fLoader->Stack();
 220   //
 221   fLoader->LoadTrackRefs();
 222   fLoader->LoadHits();
 223   fTreeTR = fLoader->TreeTR();
 224   //
 225   AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
 226   tpcl->LoadDigits();
 227   fTreeD = tpcl->TreeD();
 228   return 0;
 229 }
 230
 231 Int_t AliGenInfoMaker::CloseIOEvent()
 232 {
 233   fLoader->UnloadHeader();
 234   fLoader->UnloadKinematics();
 235   fLoader->UnloadTrackRefs();
 236   AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
 237   tpcl->UnloadDigits();
 238   return 0;
 239 }
 240
 241 Int_t AliGenInfoMaker::CloseIO()
 242 {
 243   fLoader->UnloadgAlice();
 244   return 0;
 245 }
 246
 247
 248
 249 ////////////////////////////////////////////////////////////////////////
 250 Int_t AliGenInfoMaker::Exec(Int_t nEvents, Int_t firstEventNr)
 251 {
 252   fNEvents = nEvents;
 253   fFirstEventNr = firstEventNr;
 254   return Exec();
 255 }
 256
 257 ////////////////////////////////////////////////////////////////////////
 258 Int_t AliGenInfoMaker::Exec()
 259 {
 260   TStopwatch timer;
 261   timer.Start();
 262   Int_t status =SetIO();
 263   if (status>0) return status;
 264   //
 265
 266   for (fEventNr = fFirstEventNr; fEventNr < fFirstEventNr+fNEvents;
 267        fEventNr++) {
 268     SetIO(fEventNr);
 269     fNParticles = fStack->GetNtrack();
 270     //
 271     fGenInfo = new AliMCInfo*[fNParticles];
 272     for (UInt_t i = 0; i<fNParticles; i++) {
 273       fGenInfo[i]=0;
 274     }
 275     //
 276     cout<<"Start to process event "<<fEventNr<<endl;
 277     cout<<"\tfNParticles = "<<fNParticles<<endl;
 278     if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeTR"<<endl;
 279     if (TreeTRLoop()>0) return 1;
 280     //
 281     if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeD"<<endl;
 282     if (TreeDLoop()>0) return 1;
 283     //
 284     if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeK"<<endl;
 285     if (TreeKLoop()>0) return 1;
 286     if (fDebug>2) cout<<"\tEnd loop over TreeK"<<endl;
 287     //
 288     if (BuildKinkInfo()>0) return 1;
 289     if (BuildV0Info()>0) return 1;
 290     //if (BuildHitLines()>0) return 1;
 291     if (fDebug>2) cout<<"\tEnd loop over TreeK"<<endl;
 292     //
 293     for (UInt_t i = 0; i<fNParticles; i++) {
 294       if (fGenInfo[i]) delete fGenInfo[i];
 295     }
 296     delete []fGenInfo;
 297     CloseIOEvent();
 298   }
 299   //
 300   CloseIO();
 301   CloseOutputFile();
 302
 303   cerr<<"Exec finished"<<endl;
 304
 305   timer.Stop();
 306   timer.Print();
 307   return 0;
 308 }
 309
 310 ////////////////////////////////////////////////////////////////////////
 311 void AliGenInfoMaker::CreateTreeGenTracks()
 312 {
 313   fFileGenTracks = TFile::Open(fFnRes,"RECREATE");
 314   if (!fFileGenTracks) {
 315     cerr<<"Error in CreateTreeGenTracks: cannot open file "<<fFnRes<<endl;
 316     return;
 317   }
 318   fTreeGenTracks = new TTree("genTracksTree","genTracksTree");
 319   AliMCInfo * info = new AliMCInfo;
 320   fTreeGenTracks->Branch("MC","AliMCInfo",&info);
 321   delete info;
 322   //
 323   AliGenKinkInfo *kinkinfo = new AliGenKinkInfo;
 324   fTreeKinks = new TTree("genKinksTree","genKinksTree");
 325   fTreeKinks->Branch("MC","AliGenKinkInfo",&kinkinfo);
 326   delete kinkinfo;
 327   //
 328   AliGenV0Info *v0info = new AliGenV0Info;
 329   fTreeV0 = new TTree("genV0Tree","genV0Tree");
 330   fTreeV0->Branch("MC","AliGenV0Info",&v0info);
 331   delete v0info;
 332   //
 333   //
 334   AliTrackPointArray * points0 = new AliTrackPointArray;
 335   AliTrackPointArray * points1 = new AliTrackPointArray;
 336   AliTrackPointArray * points2 = new AliTrackPointArray;
 337   fTreeHitLines = new TTree("HitLines","HitLines");
 338   fTreeHitLines->Branch("TPC.","AliTrackPointArray",&points0,32000,10);
 339   fTreeHitLines->Branch("TRD.","AliTrackPointArray",&points1,32000,10);
 340   fTreeHitLines->Branch("ITS.","AliTrackPointArray",&points2,32000,10);
 341   Int_t label=0;
 342   fTreeHitLines->Branch("Label",&label,"label/I");
 343   //
 344   fTreeGenTracks->AutoSave();
 345   fTreeKinks->AutoSave();
 346   fTreeV0->AutoSave();
 347   fTreeHitLines->AutoSave();
 348 }
 349 ////////////////////////////////////////////////////////////////////////
 350 void AliGenInfoMaker::CloseOutputFile()
 351 {
 352   if (!fFileGenTracks) {
 353     cerr<<"File "<<fFnRes<<" not found as an open file."<<endl;
 354     return;
 355   }
 356   fFileGenTracks->cd();
 357   fTreeGenTracks->Write();
 358   delete fTreeGenTracks;
 359   fTreeKinks->Write();
 360   delete fTreeKinks;
 361   fTreeV0->Write();
 362   delete fTreeV0;
 363
 364   fFileGenTracks->Close();
 365   delete fFileGenTracks;
 366   return;
 367 }
 368
 369 ////////////////////////////////////////////////////////////////////////
 370 Int_t AliGenInfoMaker::TreeKLoop()
 371 {
 372 //
 373 // open the file with treeK
 374 // loop over all entries there and save information about some tracks
 375 //
 376
 377   AliStack * stack = fStack;
 378   if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
 379
 380   if (fDebug > 0) {
 381     cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
 382         <<fEventNr<<endl;
 383   }
 384   Int_t  ipdg = 0;
 385   TParticlePDG *ppdg = 0;
 386   // not all generators give primary vertex position. Take the vertex
 387   // of the particle 0 as primary vertex.
 388   TDatabasePDG  pdg; //get pdg table
 389   //thank you very much root for this
 390   TBranch * br = fTreeGenTracks->GetBranch("MC");
 391   TParticle *particle = stack->ParticleFromTreeK(0);
 392   fVPrim[0] = particle->Vx();
 393   fVPrim[1] = particle->Vy();
 394   fVPrim[2] = particle->Vz();
 395   for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
 396     // load only particles with TR
 397     AliMCInfo * info = GetInfo(iParticle);
 398     if (!info) continue;
 399     //////////////////////////////////////////////////////////////////////
 400     info->fLabel = iParticle;
 401     //
 402     info->fParticle = *(stack->Particle(iParticle));
 403     info->fVDist[0] = info->fParticle.Vx()-fVPrim[0];
 404     info->fVDist[1] = info->fParticle.Vy()-fVPrim[1];
 405     info->fVDist[2] = info->fParticle.Vz()-fVPrim[2];
 406     info->fVDist[3] = TMath::Sqrt(info->fVDist[0]*info->fVDist[0]+
 407                                   info->fVDist[1]*info->fVDist[1]+info->fVDist[2]*info->fVDist[2]);
 408     //
 409     //
 410     ipdg = info->fParticle.GetPdgCode();
 411     info->fPdg = ipdg;
 412     ppdg = pdg.GetParticle(ipdg);
 413     info->fEventNr = fEventNr;
 414     info->Update();
 415     //////////////////////////////////////////////////////////////////////
 416     br->SetAddress(&info);
 417     fTreeGenTracks->Fill();
 418   }
 419   fTreeGenTracks->AutoSave();
 420   if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl;
 421   return 0;
 422 }
 423
 424
 425 ////////////////////////////////////////////////////////////////////////////////////
 426 Int_t  AliGenInfoMaker::BuildKinkInfo()
 427 {
 428   //
 429   // Build tree with Kink Information
 430   //
 431   AliStack * stack = fStack;
 432   if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
 433
 434   if (fDebug > 0) {
 435     cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
 436         <<fEventNr<<endl;
 437   }
 438   //  Int_t  ipdg = 0;
 439   //TParticlePDG *ppdg = 0;
 440   // not all generators give primary vertex position. Take the vertex
 441   // of the particle 0 as primary vertex.
 442   TDatabasePDG  pdg; //get pdg table
 443   //thank you very much root for this
 444   TBranch * br = fTreeKinks->GetBranch("MC");
 445   //
 446   AliGenKinkInfo * kinkinfo = new AliGenKinkInfo;
 447   //
 448   for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
 449     // load only particles with TR
 450     AliMCInfo * info = GetInfo(iParticle);
 451     if (!info) continue;
 452     if (info->fCharge==0) continue;
 453     if (info->fTRdecay.P()<0.13) continue;  //momenta cut
 454     if (info->fTRdecay.R()>500)  continue;  //R cut - decay outside barrel
 455     TParticle & particle = info->fParticle;
 456     Int_t first = particle.GetDaughter(0);
 457     if (first ==0) continue;
 458
 459     Int_t last = particle.GetDaughter(1);
 460     if (last ==0) last=first;
 461     AliMCInfo * info2 =  0;
 462     AliMCInfo * dinfo =  0;
 463
 464     for (Int_t id2=first;id2<=last;id2++){
 465       info2 = GetInfo(id2);
 466       if (!info2) continue;
 467       TParticle &p2 = info2->fParticle;
 468       Double_t r = TMath::Sqrt(p2.Vx()*p2.Vx()+p2.Vy()*p2.Vy());
 469       if ( TMath::Abs(info->fTRdecay.R()-r)>0.01) continue;
 470       if (!(p2.GetPDG())) continue;
 471       dinfo =info2;
 472
 473       kinkinfo->SetInfoMother(*info);
 474       kinkinfo->SetInfoDaughter(*dinfo);
 475       if (kinkinfo->GetMinus().GetParticle().GetPDG()==0 ||  kinkinfo->GetPlus().GetParticle().GetPDG()==0) continue;
 476       kinkinfo->Update();
 477       br->SetAddress(&kinkinfo);
 478       fTreeKinks->Fill();
 479     }
 480     /*
 481     if (dinfo){
 482       kinkinfo->fMCm = (*info);
 483       kinkinfo->GetPlus() = (*dinfo);
 484       kinkinfo->Update();
 485       br->SetAddress(&kinkinfo);
 486       fTreeKinks->Fill();
 487     }
 488     */
 489   }
 490   fTreeGenTracks->AutoSave();
 491   if (fDebug > 2) cerr<<"end of Kink Loop"<<endl;
 492   return 0;
 493 }
 494
 495
 496 ////////////////////////////////////////////////////////////////////////////////////
 497 Int_t  AliGenInfoMaker::BuildV0Info()
 498 {
 499   //
 500   // Build tree with V0 Information
 501   //
 502   AliStack * stack = fStack;
 503   if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
 504
 505   if (fDebug > 0) {
 506     cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
 507         <<fEventNr<<endl;
 508   }
 509   //  Int_t  ipdg = 0;
 510   //TParticlePDG *ppdg = 0;
 511   // not all generators give primary vertex position. Take the vertex
 512   // of the particle 0 as primary vertex.
 513   TDatabasePDG  pdg; //get pdg table
 514   //thank you very much root for this
 515   TBranch * br = fTreeV0->GetBranch("MC");
 516   //
 517   AliGenV0Info * v0info = new AliGenV0Info;
 518   //
 519   //
 520   for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
 521     // load only particles with TR
 522     AliMCInfo * info = GetInfo(iParticle);
 523     if (!info) continue;
 524     if (info->fCharge==0) continue;
 525     //
 526     //
 527     TParticle & particle = info->fParticle;
 528     Int_t mother = particle.GetMother(0);
 529     if (mother <=0) continue;
 530     //
 531     TParticle * motherparticle = stack->Particle(mother);
 532     if (!motherparticle) continue;
 533     //
 534     Int_t last = motherparticle->GetDaughter(1);
 535     if (last==(int)iParticle) continue;
 536     AliMCInfo * info2 =  info;
 537     AliMCInfo * dinfo =  GetInfo(last);
 538     if (!dinfo) continue;
 539     if (!dinfo->fParticle.GetPDG()) continue;
 540     if (!info2->fParticle.GetPDG()) continue;
 541     if (dinfo){
 542       v0info->SetInfoP(*info);
 543       v0info->SetInfoM(*dinfo);
 544       v0info->SetMother(*motherparticle);
 545       v0info->Update(fVPrim);
 546       br->SetAddress(&v0info);
 547       fTreeV0->Fill();
 548     }
 549   }
 550   fTreeV0->AutoSave();
 551   if (fDebug > 2) cerr<<"end of V0 Loop"<<endl;
 552   return 0;
 553 }
 554
 555
 556
 557
 558 ////////////////////////////////////////////////////////////////////////
 559 Int_t AliGenInfoMaker::BuildHitLines()
 560 {
 561
 562 //
 563 // open the file with treeK
 564 // loop over all entries there and save information about some tracks
 565 //
 566
 567   AliStack * stack = fStack;
 568   if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
 569
 570   if (fDebug > 0) {
 571     cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
 572         <<fEventNr<<endl;
 573   }
 574 //   Int_t  ipdg = 0;
 575 //   // TParticlePDG *ppdg = 0;
 576 //   // not all generators give primary vertex position. Take the vertex
 577 //   // of the particle 0 as primary vertex.
 578 //   TDatabasePDG  pdg; //get pdg table
 579 //   //thank you very much root for this
 580 //   AliTrackPointArray *tpcp = new AliTrackPointArray;
 581 //   AliTrackPointArray *trdp = new AliTrackPointArray;
 582 //   AliTrackPointArray *itsp = new AliTrackPointArray;
 583 //   Int_t label =0;
 584 //   //
 585 //   TBranch * brtpc = fTreeHitLines->GetBranch("TPC.");
 586 //   TBranch * brtrd = fTreeHitLines->GetBranch("TRD.");
 587 //   TBranch * brits = fTreeHitLines->GetBranch("ITS.");
 588 //   TBranch * brlabel = fTreeHitLines->GetBranch("Label");
 589 //   brlabel->SetAddress(&label);
 590 //   brtpc->SetAddress(&tpcp);
 591 //   brtrd->SetAddress(&trdp);
 592 //   brits->SetAddress(&itsp);
 593 //   //
 594 //   AliDetector *dtpc = gAlice->GetDetector("TPC");
 595 //   AliDetector *dtrd = gAlice->GetDetector("TRD");
 596 //   AliDetector *dits = gAlice->GetDetector("ITS");
 597
 598 //   for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
 599 //     // load only particles with TR
 600 //     AliMCInfo * info = GetInfo(iParticle);
 601 //     if (!info) continue;
 602 //     Int_t primpart = info->fPrimPart;
 603 //     ipdg = info->fParticle.GetPdgCode();
 604 //     label = iParticle;
 605 //     //
 606 //     gAlice->ResetHits();
 607 //     tpcp->Reset();
 608 //     itsp->Reset();
 609 //     trdp->Reset();
 610 //     tpcp->fLabel1 = ipdg;
 611 //     trdp->fLabel1 = ipdg;
 612 //     itsp->fLabel1 = ipdg;
 613 //     if (dtpc->TreeH()->GetEvent(primpart)){
 614 //       dtpc->LoadPoints(primpart);
 615 //       tpcp->Reset(dtpc,iParticle);
 616 //     }
 617 //     if (dtrd->TreeH()->GetEvent(primpart)){
 618 //       dtrd->LoadPoints(primpart);
 619 //       trdp->Reset(dtrd,iParticle);
 620 //     }
 621 //     if (dits->TreeH()->GetEvent(primpart)){
 622 //       dits->LoadPoints(primpart);
 623 //       itsp->Reset(dits,iParticle);
 624 //     }
 625 //     //
 626 //     fTreeHitLines->Fill();
 627 //     dtpc->ResetPoints();
 628 //     dtrd->ResetPoints();
 629 //     dits->ResetPoints();
 630 //   }
 631 //   delete tpcp;
 632 //   delete trdp;
 633 //   delete itsp;
 634 //   fTreeHitLines->AutoSave();
 635 //   if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl;
 636   return 0;
 637 }
 638
 639
 640 ////////////////////////////////////////////////////////////////////////
 641 Int_t AliGenInfoMaker::TreeDLoop()
 642 {
 643   //
 644   // open the file with treeD
 645   // loop over all entries there and save information about some tracks
 646   //
 647
 648   Int_t nInnerSector = fParamTPC->GetNInnerSector();
 649   Int_t rowShift = 0;
 650   Int_t zero=fParamTPC->GetZeroSup()+6;
 651   //
 652   //
 653   AliSimDigits digitsAddress, *digits=&digitsAddress;
 654   fTreeD->GetBranch("Segment")->SetAddress(&digits);
 655
 656   Int_t sectorsByRows=(Int_t)fTreeD->GetEntries();
 657   if (fDebug > 1) cout<<"\tsectorsByRows = "<<sectorsByRows<<endl;
 658   for (Int_t i=0; i<sectorsByRows; i++) {
 659     if (!fTreeD->GetEvent(i)) continue;
 660     Int_t sec,row;
 661     fParamTPC->AdjustSectorRow(digits->GetID(),sec,row);
 662     if (fDebug > 1) cout<<sec<<' '<<row<<"                          \r";
 663     // here I expect that upper sectors follow lower sectors
 664     if (sec > nInnerSector) rowShift = fParamTPC->GetNRowLow();
 665     //
 666     digits->ExpandTrackBuffer();
 667     digits->First();
 668     do {
 669       Int_t iRow=digits->CurrentRow();
 670       Int_t iColumn=digits->CurrentColumn();
 671       Short_t digitValue = digits->CurrentDigit();
 672       if (digitValue >= zero) {
 673         Int_t label;
 674         for (Int_t j = 0; j<3; j++) {
 675           //      label = digits->GetTrackID(iRow,iColumn,j);
 676           label = digits->GetTrackIDFast(iRow,iColumn,j)-2;
 677           if (label >= (Int_t)fNParticles) { //don't label from bakground event
 678             continue;
 679           }
 680           if (label >= 0 && label <= (Int_t)fNParticles) {
 681             if (fDebug > 6 ) {
 682               cout<<"Inner loop: sector, iRow, iColumn, label, value, row "
 683                   <<sec<<" "
 684                   <<iRow<<" "<<iColumn<<" "<<label<<" "<<digitValue
 685                   <<" "<<row<<endl;
 686             }
 687             AliMCInfo * info = GetInfo(label);
 688             if (info){
 689               info->fTPCRow.SetRow(row+rowShift);
 690             }
 691           }
 692         }
 693       }
 694     } while (digits->Next());
 695   }
 696
 697   if (fDebug > 2) cerr<<"end of TreeDLoop"<<endl;
 698   return 0;
 699 }
 700
 701
 702 ////////////////////////////////////////////////////////////////////////
 703 Int_t AliGenInfoMaker::TreeTRLoop()
 704 {
 705   //
 706   // loop over TrackReferences and store the first one for each track
 707   //
 708   TTree * treeTR = fTreeTR;
 709   Int_t nPrimaries = (Int_t) treeTR->GetEntries();
 710   if (fDebug > 1) cout<<"There are "<<nPrimaries<<" entries in TreeTR"<<endl;
 711   //
 712   //
 713   //track references for TPC
 714   TClonesArray* tpcArrayTR = new TClonesArray("AliTrackReference");
 715   TClonesArray* itsArrayTR = new TClonesArray("AliTrackReference");
 716   TClonesArray* trdArrayTR = new TClonesArray("AliTrackReference");
 717   TClonesArray* tofArrayTR = new TClonesArray("AliTrackReference");
 718   TClonesArray* runArrayTR = new TClonesArray("AliTrackReference");
 719   //
 720   if (treeTR->GetBranch("TPC"))    treeTR->GetBranch("TPC")->SetAddress(&tpcArrayTR);
 721   if (treeTR->GetBranch("ITS"))    treeTR->GetBranch("ITS")->SetAddress(&itsArrayTR);
 722   if (treeTR->GetBranch("TRD"))    treeTR->GetBranch("TRD")->SetAddress(&trdArrayTR);
 723   if (treeTR->GetBranch("TOF"))    treeTR->GetBranch("TOF")->SetAddress(&tofArrayTR);
 724   if (treeTR->GetBranch("AliRun")) treeTR->GetBranch("AliRun")->SetAddress(&runArrayTR);
 725   //
 726   //
 727   //
 728   for (Int_t iPrimPart = 0; iPrimPart<nPrimaries; iPrimPart++) {
 729     treeTR->GetEntry(iPrimPart);
 730     //
 731     // Loop over TPC references
 732     //
 733     for (Int_t iTrackRef = 0; iTrackRef < tpcArrayTR->GetEntriesFast(); iTrackRef++) {
 734       AliTrackReference *trackRef = (AliTrackReference*)tpcArrayTR->At(iTrackRef);
 735       //
 736       if (trackRef->TestBit(BIT(2))){
 737         //if decay
 738         if (trackRef->P()<fTPCPtCut) continue;
 739         Int_t label = trackRef->GetTrack();
 740         AliMCInfo * info = GetInfo(label);
 741         if (!info) info = MakeInfo(label);
 742         info->fTRdecay = *trackRef;
 743       }
 744       //
 745       if (trackRef->P()<fTPCPtCut) continue;
 746       Int_t label = trackRef->GetTrack();
 747       AliMCInfo * info = GetInfo(label);
 748       if (!info) info = MakeInfo(label);
 749       if (!info) continue;
 750       info->fPrimPart =  iPrimPart;
 751       TClonesArray & arr = *(info->fTPCReferences);
 752       new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);
 753     }
 754     //
 755     // Loop over ITS references
 756     //
 757     for (Int_t iTrackRef = 0; iTrackRef < itsArrayTR->GetEntriesFast(); iTrackRef++) {
 758       AliTrackReference *trackRef = (AliTrackReference*)itsArrayTR->At(iTrackRef);
 759       //
 760       //
 761       if (trackRef->P()<fTPCPtCut) continue;
 762       Int_t label = trackRef->GetTrack();
 763       AliMCInfo * info = GetInfo(label);
 764       if ( (!info) && trackRef->Pt()<fITSPtCut) continue;
 765       if (!info) info = MakeInfo(label);
 766       if (!info) continue;
 767       info->fPrimPart =  iPrimPart;
 768       TClonesArray & arr = *(info->fITSReferences);
 769       new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);
 770     }
 771     //
 772     // Loop over TRD references
 773     //
 774     for (Int_t iTrackRef = 0; iTrackRef < trdArrayTR->GetEntriesFast(); iTrackRef++) {
 775       AliTrackReference *trackRef = (AliTrackReference*)trdArrayTR->At(iTrackRef);
 776       //
 777       if (trackRef->P()<fTRDPtCut) continue;
 778       Int_t label = trackRef->GetTrack();
 779       AliMCInfo * info = GetInfo(label);
 780       if ( (!info) && trackRef->Pt()<fTRDPtCut) continue;
 781       if (!info) info = MakeInfo(label);
 782       if (!info) continue;
 783       info->fPrimPart =  iPrimPart;
 784       TClonesArray & arr = *(info->fTRDReferences);
 785       new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);
 786     }
 787     //
 788     // Loop over TOF references
 789     //
 790     for (Int_t iTrackRef = 0; iTrackRef < tofArrayTR->GetEntriesFast(); iTrackRef++) {
 791       AliTrackReference *trackRef = (AliTrackReference*)tofArrayTR->At(iTrackRef);
 792       Int_t label = trackRef->GetTrack();
 793       AliMCInfo * info = GetInfo(label);
 794       if (!info){
 795         if (trackRef->Pt()<fTPCPtCut) continue;
 796         if ( (!info) && trackRef->Pt()<fTOFPtCut) continue;
 797       }
 798       if (!info) info = MakeInfo(label);
 799       if (!info) continue;
 800       info->fPrimPart =  iPrimPart;
 801       TClonesArray & arr = *(info->fTOFReferences);
 802       new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);
 803     }
 804     //
 805     //
 806     // get dacay position
 807     //
 808     for (Int_t iTrackRef = 0; iTrackRef < runArrayTR->GetEntriesFast(); iTrackRef++) {
 809       AliTrackReference *trackRef = (AliTrackReference*)runArrayTR->At(iTrackRef);
 810       //
 811       Int_t label = trackRef->GetTrack();
 812       AliMCInfo * info = GetInfo(label);
 813       if (!info) continue;
 814       if (!trackRef->TestBit(BIT(2))) continue;  //if not decay
 815       //      if (TMath::Abs(trackRef.X());
 816       info->fTRdecay = *trackRef;
 817     }
 818   }
 819   //
 820   tpcArrayTR->Delete();
 821   delete  tpcArrayTR;
 822   trdArrayTR->Delete();
 823   delete  trdArrayTR;
 824   tofArrayTR->Delete();
 825   delete  tofArrayTR;
 826   itsArrayTR->Delete();
 827   delete  itsArrayTR;
 828   runArrayTR->Delete();
 829   delete  runArrayTR;
 830   //
 831   return TreeTRLoopNew();
 832 }
 833
 834
 835
 836
 837
 838 ////////////////////////////////////////////////////////////////////////
 839 Int_t AliGenInfoMaker::TreeTRLoopNew()
 840 {
 841   //
 842   // loop over TrackReferences and store the first one for each track
 843   //
 844   TTree * treeTR = fTreeTR;
 845   Int_t nPrimaries = (Int_t) treeTR->GetEntries();
 846   if (fDebug > 1) cout<<"There are "<<nPrimaries<<" entries in TreeTR"<<endl;
 847   //
 848   //
 849   //
 850   TClonesArray* runArrayTR = new TClonesArray("AliTrackReference");
 851   //
 852   if (treeTR->GetBranch("TrackReferences")) treeTR->GetBranch("TrackReferences")->SetAddress(&runArrayTR);
 853   //
 854   //
 855   //
 856   for (Int_t ipart = 0; ipart<fStack->GetNtrack(); ipart++) {
 857     TParticle * part = fStack->Particle(ipart);
 858     if (!part) continue;
 859     if (part->Pt()<fITSPtCut) continue;
 860     if (part->R()>250.) continue;
 861     if (TMath::Abs(part->Vz())>250.) continue;
 862     if (TMath::Abs(part->Pz()/part->Pt())>2.5) continue;
 863     MakeInfo(ipart);
 864   }
 865   //
 866   //
 867
 868   for (Int_t iPrimPart = 0; iPrimPart<nPrimaries; iPrimPart++) {
 869     treeTR->GetEntry(iPrimPart);
 870     //
 871     // gettrack references
 872     //
 873     for (Int_t iTrackRef = 0; iTrackRef < runArrayTR->GetEntriesFast(); iTrackRef++) {
 874       AliTrackReference *trackRef = (AliTrackReference*)runArrayTR->At(iTrackRef);
 875       //
 876       Int_t label = trackRef->GetTrack();
 877       //
 878       // TPC
 879       //
 880       if (trackRef->DetectorId()== AliTrackReference::kTPC){
 881         //
 882         Int_t label = trackRef->GetTrack();
 883         AliMCInfo * info = GetInfo(label);
 884         if (!info && trackRef->Pt()<fTPCPtCut) continue;
 885         if (!info) info = MakeInfo(label);
 886         if (!info) continue;
 887         info->fPrimPart =  iPrimPart;
 888         TClonesArray & arr = *(info->fTPCReferences);
 889         new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);
 890       }
 891       //
 892       // ITS
 893       //
 894       if (trackRef->DetectorId()== AliTrackReference::kITS){
 895         //
 896         Int_t label = trackRef->GetTrack();
 897         AliMCInfo * info = GetInfo(label);
 898         if (!info && trackRef->Pt()<fITSPtCut) continue;
 899         if (!info) info = MakeInfo(label);
 900         if (!info) continue;
 901         info->fPrimPart =  iPrimPart;
 902         TClonesArray & arr = *(info->fITSReferences);
 903         new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);
 904       }
 905       //
 906       // TRD
 907       //
 908       if (trackRef->DetectorId()== AliTrackReference::kTRD){
 909         //
 910         Int_t label = trackRef->GetTrack();
 911         AliMCInfo * info = GetInfo(label);
 912         if (!info&&trackRef->P()<fTRDPtCut) continue;
 913         if (!info) info = MakeInfo(label);
 914         if (!info) continue;
 915         info->fPrimPart =  iPrimPart;
 916         TClonesArray & arr = *(info->fTRDReferences);
 917         new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);
 918       }
 919       //
 920       // TOF
 921       //
 922       if (trackRef->DetectorId()== AliTrackReference::kTOF){
 923         //
 924         Int_t label = trackRef->GetTrack();
 925         AliMCInfo * info = GetInfo(label);
 926         if (!info&&trackRef->P()<fTOFPtCut) continue;
 927         if (!info) info = MakeInfo(label);
 928         if (!info) continue;
 929         info->fPrimPart =  iPrimPart;
 930         TClonesArray & arr = *(info->fTOFReferences);
 931         new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);
 932       }
 933       //
 934       // decay case
 935       //
 936       if (trackRef->DetectorId()== AliTrackReference::kDisappeared){
 937         //
 938         AliMCInfo * info = GetInfo(label);
 939         if (!info) continue;
 940         //if (!trackRef->TestBit(BIT(2))) continue;  //if not decay
 941         info->fTRdecay = *trackRef;
 942       }
 943
 944     }
 945   }
 946   //
 947   //
 948   return 0;
 949 }
 950
 951 ////////////////////////////////////////////////////////////////////////
 952 Float_t AliGenInfoMaker::TR2LocalX(AliTrackReference *trackRef,
 953                                     AliTPCParam *paramTPC) const {
 954
 955   Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
 956   Int_t index[4];
 957   paramTPC->Transform0to1(x,index);
 958   paramTPC->Transform1to2(x,index);
 959   return x[0];
 960 }
 961 ////////////////////////////////////////////////////////////////////////
 962
 963
 964
 965 AliTPCParam * AliGenInfoMaker::GetTPCParam(){
 966   AliTPCParamSR * par = new AliTPCParamSR;
 967   par->Update();
 968   return par;
 969 }
 970
 971