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