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   if (fLastEvent == -1)
 120     fLastEvent = gime->MaxEvent() - 1 ;
 121   else
 122     fLastEvent = TMath::Min(fLastEvent,gime->MaxEvent());
 123   Int_t nEvents   = fLastEvent - fFirstEvent + 1;
 124
 125   Int_t ievent ;
 126
 127
 128   for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
 129     gime->Event(ievent,"TR") ;
 130     if(gime->TrackSegments() && //Skip events, where no track segments made
 131        gime->TrackSegments()->GetEntriesFast()) {
 132       MakeRecParticles() ;
 133       WriteRecParticles();
 134       if(strstr(option,"deb"))
 135         PrintRecParticles(option) ;
 136       //increment the total number of rec particles per run
 137       fRecParticlesInRun += gime->RecParticles()->GetEntriesFast() ;
 138     }
 139   }
 140   if(strstr(option,"tim")){
 141     gBenchmark->Stop("EMCALPID");
 142     printf("Exec: took %f seconds for PID %f seconds per event",
 143          gBenchmark->GetCpuTime("EMCALPID"),
 144          gBenchmark->GetCpuTime("EMCALPID")/nEvents) ;
 145   }
 146
 147   Unload();
 148 }
 149
 150 //____________________________________________________________________________
 151 void  AliEMCALPIDv1::MakeRecParticles(){
 152
 153   // Makes a RecParticle out of a TrackSegment
 154
 155   AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
 156   TObjArray * aECARecPoints = gime->ECARecPoints() ;
 157   TClonesArray * trackSegments = gime->TrackSegments() ;
 158   if ( !aECARecPoints || !trackSegments ) {
 159     Fatal("MakeRecParticles", "RecPoints or TrackSegments not found !") ;
 160   }
 161   TClonesArray * recParticles  = gime->RecParticles() ;
 162   recParticles->Clear();
 163
 164   TIter next(trackSegments) ;
 165   AliEMCALTrackSegment * ts ;
 166   Int_t index = 0 ;
 167   AliEMCALRecParticle * rp ;
 168   while ( (ts = (AliEMCALTrackSegment *)next()) ) {
 169
 170     new( (*recParticles)[index] ) AliEMCALRecParticle() ;
 171     rp = (AliEMCALRecParticle *)recParticles->At(index) ;
 172     rp->SetTrackSegment(index) ;
 173     rp->SetIndexInList(index) ;
 174
 175     AliEMCALTowerRecPoint * eca = 0 ;
 176     if(ts->GetECAIndex()>=0)
 177       eca = dynamic_cast<AliEMCALTowerRecPoint *>(aECARecPoints->At(ts->GetECAIndex())) ;
 178
 179     // Now set type (reconstructed) of the particle
 180
 181     // Choose the cluster energy range
 182
 183     if (!eca) {
 184       Fatal("MakeRecParticles", "-> emcal(%d) = %d", ts->GetECAIndex(), eca ) ;
 185     }
 186
 187     Float_t    e = eca->GetEnergy() ;
 188
 189     Float_t  lambda[2] ;
 190     eca->GetElipsAxis(lambda) ;
 191
 192     if((lambda[0]>0.01) && (lambda[1]>0.01)){
 193       // Looking PCA. Define and calculate the data (X),
 194       // introduce in the function X2P that gives the components (P).
 195
 196       Float_t  spher = 0. ;
 197       Float_t  emaxdtotal = 0. ;
 198
 199       if((lambda[0]+lambda[1])!=0)
 200         spher=fabs(lambda[0]-lambda[1])/(lambda[0]+lambda[1]);
 201
 202       emaxdtotal=eca->GetMaximalEnergy()/eca->GetEnergy();
 203     }
 204
 205     //    Float_t time = ecal->GetTime() ;
 206
 207     //Set momentum, energy and other parameters
 208     Float_t  enca = e;
 209     TVector3 dir(0., 0., 0.) ;
 210     dir.SetMag(enca) ;
 211     rp->SetMomentum(dir.X(),dir.Y(),dir.Z(),enca) ;
 212     rp->SetCalcMass(0);
 213     rp->Name(); //If photon sets the particle pdg name to gamma
 214     rp->SetProductionVertex(0,0,0,0);
 215     rp->SetFirstMother(-1);
 216     rp->SetLastMother(-1);
 217     rp->SetFirstDaughter(-1);
 218     rp->SetLastDaughter(-1);
 219     rp->SetPolarisation(0,0,0);
 220     index++ ;
 221   }
 222
 223 }
 224
 225 //____________________________________________________________________________
 226 void  AliEMCALPIDv1:: Print(Option_t * /*option*/) const
 227 {
 228   // Print the parameters used for the particle type identification
 229
 230     printf("Print: =============== AliEMCALPID1 ================") ;
 231     printf("Making PID\n");
 232     printf("    Pricipal analysis file from 0.5 to 100 %s\n", fFileName.Data() ) ;
 233     printf("    Name of parameters file     %s\n", fFileNamePar.Data() )  ;
 234 }
 235
 236 //____________________________________________________________________________
 237 void  AliEMCALPIDv1::Print() const
 238 {
 239   // Print the parameters used for the particle type identification
 240
 241     Info("Print", "=============== AliEMCALPIDv1 ================") ;
 242 }
 243
 244 //____________________________________________________________________________
 245 void AliEMCALPIDv1::PrintRecParticles(Option_t * option)
 246 {
 247   // Print table of reconstructed particles
 248
 249   AliEMCALGetter *gime = AliEMCALGetter::Instance() ;
 250
 251   TClonesArray * recParticles = gime->RecParticles() ;
 252
 253   printf("\nevent %i", gAlice->GetEvNumber());
 254   printf("       found %i", recParticles->GetEntriesFast());
 255   printf(" RecParticles\n");
 256
 257   if(strstr(option,"all")) {  // printing found TS
 258     printf("\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       printf("\n");
 264       printf(rp->Name().Data());
 265       printf(" %i", rp->GetIndexInList());
 266       printf(" %i", rp->GetType());
 267     }
 268   }
 269 }
 270
 271 //____________________________________________________________________________
 272 void AliEMCALPIDv1::Unload()
 273 {
 274   // Unloads RecPoints, TrackSegments and RecParticles from the folder
 275   AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
 276   gime->EmcalLoader()->UnloadRecPoints() ;
 277   gime->EmcalLoader()->UnloadTracks() ;
 278   gime->EmcalLoader()->UnloadRecParticles() ;
 279 }
 280
 281 //____________________________________________________________________________
 282 void  AliEMCALPIDv1::WriteRecParticles()
 283 {
 284   // Write RecParticles array to file
 285   AliEMCALGetter *gime = AliEMCALGetter::Instance() ;
 286
 287   TClonesArray * recParticles = gime->RecParticles() ;
 288   recParticles->Expand(recParticles->GetEntriesFast() ) ;
 289   TTree * treeP =  gime->TreeP() ;
 290
 291
 292
 293   //First rp
 294   Int_t bufferSize = 32000 ;
 295   TBranch * rpBranch = treeP->Branch("EMCALRP",&recParticles,bufferSize);
 296   rpBranch->SetTitle(BranchName());
 297
 298   rpBranch->Fill() ;
 299
 300   gime->WriteRecParticles("OVERWRITE");
 301   gime->WritePID("OVERWRITE");
 302
 303 }
 304