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