EMCAL/AliEMCALPIDv1.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 /* $Id$ */
  17
  18 //_________________________________________________________________________
  19 // Implementation version v1 of the EMCAL particle identifier
  20
  21 //*-- Author: Yves Schutz (SUBATECH)
  22 //
  23 // --- ROOT system ---
  24 //#include "TROOT.h"
  25 #include "TTree.h"
  26 #include "TBenchmark.h"
  27 #include "TSystem.h"
  28
  29   // --- Standard library ---
  30
  31   // --- AliRoot header files ---
  32 #include "AliEMCALPIDv1.h"
  33 #include "AliEMCALTrackSegment.h"
  34 #include "AliEMCALRecParticle.h"
  35 #include "AliEMCALGetter.h"
  36
  37   ClassImp( AliEMCALPIDv1)
  38
  39 //____________________________________________________________________________
  40 AliEMCALPIDv1::AliEMCALPIDv1():AliEMCALPID()
  41 {
  42   // default ctor
  43
  44   InitParameters() ;
  45   fDefaultInit = kTRUE ;
  46
  47 }
  48
  49 //____________________________________________________________________________
  50 AliEMCALPIDv1::AliEMCALPIDv1(const AliEMCALPIDv1 & pid ):AliEMCALPID(pid)
  51 {
  52   // ctor
  53   InitParameters() ;
  54   Init() ;
  55
  56 }
  57
  58 //____________________________________________________________________________
  59 AliEMCALPIDv1::AliEMCALPIDv1(const TString alirunFileName, const TString eventFolderName):AliEMCALPID(alirunFileName, eventFolderName)
  60 {
  61   //ctor with the indication on where to look for the track segments
  62
  63   InitParameters() ;
  64   Init() ;
  65   fDefaultInit = kFALSE ;
  66
  67 }
  68
  69 //____________________________________________________________________________
  70 AliEMCALPIDv1::~AliEMCALPIDv1()
  71 {
  72   // dtor
  73 }
  74
  75 //____________________________________________________________________________
  76 const TString AliEMCALPIDv1::BranchName() const
  77 {
  78
  79   return GetName() ;
  80 }
  81
  82 //____________________________________________________________________________
  83 void AliEMCALPIDv1::Init()
  84 {
  85   // Make all memory allocations that are not possible in default constructor
  86   // Add the PID task to the list of EMCAL tasks
  87
  88
  89   AliEMCALGetter * gime = AliEMCALGetter::Instance(GetTitle(), fEventFolderName.Data()) ;
  90
  91   if ( !gime->PID() )
  92     gime->PostPID(this) ;
  93 }
  94
  95 //____________________________________________________________________________
  96 void AliEMCALPIDv1::InitParameters()
  97 {
  98   // Initialize the parameters
  99   fRecParticlesInRun = 0 ;
 100   fNEvent            = 0 ;
 101   fRecParticlesInRun = 0 ;
 102 }
 103
 104 //____________________________________________________________________________
 105
 106 void  AliEMCALPIDv1::Exec(Option_t * option)
 107 {
 108   //Steering method
 109
 110   if(strstr(option,"tim"))
 111     gBenchmark->Start("EMCALPID");
 112
 113   if(strstr(option,"print")) {
 114     Print("") ;
 115     return ;
 116   }
 117   AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
 118
 119   Int_t nevents = gime->MaxEvent() ;
 120   Int_t ievent ;
 121
 122
 123   for(ievent = 0; ievent < nevents; ievent++){
 124     gime->Event(ievent,"TR") ;
 125     if(gime->TrackSegments() && //Skip events, where no track segments made
 126        gime->TrackSegments()->GetEntriesFast()) {
 127       MakeRecParticles() ;
 128       WriteRecParticles();
 129       if(strstr(option,"deb"))
 130         PrintRecParticles(option) ;
 131       //increment the total number of rec particles per run
 132       fRecParticlesInRun += gime->RecParticles()->GetEntriesFast() ;
 133     }
 134   }
 135   if(strstr(option,"tim")){
 136     gBenchmark->Stop("EMCALPID");
 137     Info("Exec", "took %f seconds for PID %f seconds per event",
 138          gBenchmark->GetCpuTime("EMCALPID"),
 139          gBenchmark->GetCpuTime("EMCALPID")/nevents) ;
 140   }
 141
 142   Unload();
 143 }
 144
 145 //____________________________________________________________________________
 146 void  AliEMCALPIDv1::MakeRecParticles(){
 147
 148   // Makes a RecParticle out of a TrackSegment
 149
 150   AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
 151   TObjArray * aECARecPoints = gime->ECARecPoints() ;
 152   TObjArray * aPRERecPoints = gime->PRERecPoints() ;
 153   TObjArray * aHCARecPoints = gime->HCARecPoints() ;
 154   TClonesArray * trackSegments = gime->TrackSegments() ;
 155   if ( !aECARecPoints || !aPRERecPoints || !aHCARecPoints || !trackSegments ) {
 156     Fatal("MakeRecParticles", "RecPoints or TrackSegments not found !") ;
 157   }
 158   TClonesArray * recParticles  = gime->RecParticles() ;
 159   recParticles->Clear();
 160
 161   TIter next(trackSegments) ;
 162   AliEMCALTrackSegment * ts ;
 163   Int_t index = 0 ;
 164   AliEMCALRecParticle * rp ;
 165   while ( (ts = (AliEMCALTrackSegment *)next()) ) {
 166
 167     new( (*recParticles)[index] ) AliEMCALRecParticle() ;
 168     rp = (AliEMCALRecParticle *)recParticles->At(index) ;
 169     rp->SetTrackSegment(index) ;
 170     rp->SetIndexInList(index) ;
 171
 172     AliEMCALTowerRecPoint * eca = 0 ;
 173     if(ts->GetECAIndex()>=0)
 174       eca = dynamic_cast<AliEMCALTowerRecPoint *>(aECARecPoints->At(ts->GetECAIndex())) ;
 175
 176     AliEMCALTowerRecPoint * pre = 0 ;
 177     if(ts->GetPREIndex()>=0)
 178       pre = dynamic_cast<AliEMCALTowerRecPoint *>(aPRERecPoints->At(ts->GetPREIndex())) ;
 179
 180     AliEMCALTowerRecPoint * hca = 0 ;
 181     if(ts->GetHCAIndex()>=0)
 182       hca = dynamic_cast<AliEMCALTowerRecPoint *>(aHCARecPoints->At(ts->GetHCAIndex())) ;
 183
 184     // Now set type (reconstructed) of the particle
 185
 186     // Choose the cluster energy range
 187
 188     if (!eca) {
 189       Fatal("MakeRecParticles", "-> emcal(%d) = %d", ts->GetECAIndex(), eca ) ;
 190     }
 191
 192     Float_t    e = eca->GetEnergy() ;
 193
 194     Float_t  lambda[2] ;
 195     eca->GetElipsAxis(lambda) ;
 196
 197     if((lambda[0]>0.01) && (lambda[1]>0.01)){
 198       // Looking PCA. Define and calculate the data (X),
 199       // introduce in the function X2P that gives the components (P).
 200
 201       Float_t  spher = 0. ;
 202       Float_t  emaxdtotal = 0. ;
 203
 204       if((lambda[0]+lambda[1])!=0)
 205         spher=fabs(lambda[0]-lambda[1])/(lambda[0]+lambda[1]);
 206
 207       emaxdtotal=eca->GetMaximalEnergy()/eca->GetEnergy();
 208     }
 209
 210     //    Float_t time = ecal->GetTime() ;
 211
 212     //Set momentum, energy and other parameters
 213     Float_t  enca = e;
 214     TVector3 dir(0., 0., 0.) ;
 215     dir.SetMag(enca) ;
 216     rp->SetMomentum(dir.X(),dir.Y(),dir.Z(),enca) ;
 217     rp->SetCalcMass(0);
 218     rp->Name(); //If photon sets the particle pdg name to gamma
 219     rp->SetProductionVertex(0,0,0,0);
 220     rp->SetFirstMother(-1);
 221     rp->SetLastMother(-1);
 222     rp->SetFirstDaughter(-1);
 223     rp->SetLastDaughter(-1);
 224     rp->SetPolarisation(0,0,0);
 225     index++ ;
 226   }
 227
 228 }
 229
 230 //____________________________________________________________________________
 231 void  AliEMCALPIDv1:: Print(Option_t * /*option*/) const
 232 {
 233   // Print the parameters used for the particle type identification
 234
 235     Info("Print", "=============== AliEMCALPID1 ================") ;
 236     printf("Making PID\n");
 237     printf("    Pricipal analysis file from 0.5 to 100 %s\n", fFileName.Data() ) ;
 238     printf("    Name of parameters file     %s\n", fFileNamePar.Data() )  ;
 239 }
 240
 241 //____________________________________________________________________________
 242 void AliEMCALPIDv1::PrintRecParticles(Option_t * option)
 243 {
 244   // Print table of reconstructed particles
 245
 246   AliEMCALGetter *gime = AliEMCALGetter::Instance() ;
 247
 248   TClonesArray * recParticles = gime->RecParticles() ;
 249
 250   TString message ;
 251   message  = "\nevent " ;
 252   message += gAlice->GetEvNumber() ;
 253   message += "       found " ;
 254   message += recParticles->GetEntriesFast();
 255   message += " RecParticles\n" ;
 256
 257   if(strstr(option,"all")) {  // printing found TS
 258     message += "\n  PARTICLE         Index    \n" ;
 259
 260     Int_t index ;
 261     for (index = 0 ; index < recParticles->GetEntries() ; index++) {
 262       AliEMCALRecParticle * rp = (AliEMCALRecParticle * ) recParticles->At(index) ;
 263       message += "\n" ;
 264       message += rp->Name().Data() ;
 265       message += " " ;
 266       message += rp->GetIndexInList() ;
 267       message += " " ;
 268       message += rp->GetType()  ;
 269     }
 270   }
 271   Info("Print", message.Data() ) ;
 272 }
 273
 274 //____________________________________________________________________________
 275 void AliEMCALPIDv1::Unload()
 276 {
 277   // Unloads RecPoints, TrackSegments and RecParticles from the folder
 278   AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
 279   gime->EmcalLoader()->UnloadRecPoints() ;
 280   gime->EmcalLoader()->UnloadTracks() ;
 281   gime->EmcalLoader()->UnloadRecParticles() ;
 282 }
 283
 284 //____________________________________________________________________________
 285 void  AliEMCALPIDv1::WriteRecParticles()
 286 {
 287   // Write RecParticles array to file
 288   AliEMCALGetter *gime = AliEMCALGetter::Instance() ;
 289
 290   TClonesArray * recParticles = gime->RecParticles() ;
 291   recParticles->Expand(recParticles->GetEntriesFast() ) ;
 292   TTree * treeP =  gime->TreeP() ;
 293
 294
 295
 296   //First rp
 297   Int_t bufferSize = 32000 ;
 298   TBranch * rpBranch = treeP->Branch("EMCALRP",&recParticles,bufferSize);
 299   rpBranch->SetTitle(BranchName());
 300
 301   rpBranch->Fill() ;
 302
 303   gime->WriteRecParticles("OVERWRITE");
 304   gime->WritePID("OVERWRITE");
 305
 306 }
 307