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