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
  33 #include "TSystem.h"
  34
  35 // --- Standard library ---
  36
  37 #include <Riostream.h>
  38
  39 // --- AliRoot header files ---
  40
  41 #include "AliRun.h"
  42 #include "AliGenerator.h"
  43 #include "AliEMCAL.h"
  44 #include "AliEMCALPIDv1.h"
  45 #include "AliEMCALClusterizerv1.h"
  46 #include "AliEMCALTrackSegment.h"
  47 #include "AliEMCALTrackSegmentMakerv1.h"
  48 #include "AliEMCALRecParticle.h"
  49 #include "AliEMCALGeometry.h"
  50 #include "AliEMCALGetter.h"
  51
  52 ClassImp( AliEMCALPIDv1)
  53
  54 //____________________________________________________________________________
  55 AliEMCALPIDv1::AliEMCALPIDv1():AliEMCALPID()
  56 {
  57   // default ctor
  58
  59   InitParameters() ;
  60   fDefaultInit = kTRUE ;
  61
  62 }
  63
  64 //____________________________________________________________________________
  65 AliEMCALPIDv1::AliEMCALPIDv1(const char * headerFile,const char * name, const Bool_t toSplit)
  66 :AliEMCALPID(headerFile, name,toSplit)
  67 {
  68   //ctor with the indication on where to look for the track segments
  69
  70   InitParameters() ;
  71
  72   Init() ;
  73   fDefaultInit = kFALSE ;
  74
  75 }
  76
  77 //____________________________________________________________________________
  78 AliEMCALPIDv1::~AliEMCALPIDv1()
  79 {
  80   // dtor
  81
  82   if (!fDefaultInit) {
  83     fSplitFile = 0 ;
  84   }
  85 }
  86
  87 //____________________________________________________________________________
  88 const TString AliEMCALPIDv1::BranchName() const
  89 {
  90   TString branchName(GetName() ) ;
  91   branchName.Remove(branchName.Index(Version())-1) ;
  92   return branchName ;
  93 }
  94
  95 //____________________________________________________________________________
  96 void AliEMCALPIDv1::Init()
  97 {
  98   // Make all memory allocations that are not possible in default constructor
  99   // Add the PID task to the list of EMCAL tasks
 100
 101   if ( strcmp(GetTitle(), "") == 0 )
 102     SetTitle("galice.root") ;
 103
 104   TString branchname(GetName()) ;
 105   branchname.Remove(branchname.Index(Version())-1) ;
 106   AliEMCALGetter * gime = AliEMCALGetter::GetInstance(GetTitle(),branchname.Data(),fToSplit ) ;
 107
 108   //  gime->SetRecParticlesTitle(BranchName()) ;
 109   if ( gime == 0 ) {
 110     Error("Init", "Could not obtain the Getter object !" ) ;
 111     return ;
 112   }
 113
 114   fSplitFile = 0 ;
 115   if(fToSplit){
 116     //First - extract full path if necessary
 117     TString fileName(GetTitle()) ;
 118     Ssiz_t islash = fileName.Last('/') ;
 119     if(islash<fileName.Length())
 120       fileName.Remove(islash+1,fileName.Length()) ;
 121     else
 122       fileName="" ;
 123     fileName+="EMCAL.RecData." ;
 124     if((strcmp(branchname.Data(),"Default")!=0)&&(strcmp(branchname.Data(),"")!=0)){
 125       fileName+=branchname.Data() ;
 126       fileName+="." ;
 127     }
 128     fileName+="root" ;
 129     fSplitFile = static_cast<TFile*>(gROOT->GetFile(fileName.Data()));
 130     if(!fSplitFile)
 131       fSplitFile =  TFile::Open(fileName.Data(),"update") ;
 132   }
 133
 134   gime->PostPID(this) ;
 135   gime->PostRecParticles(branchname) ;
 136
 137 }
 138
 139 //____________________________________________________________________________
 140 void AliEMCALPIDv1::InitParameters()
 141 {
 142
 143   fRecParticlesInRun = 0 ;
 144   fNEvent            = 0 ;
 145   fRecParticlesInRun = 0 ;
 146   TString pidName( GetName()) ;
 147   if (pidName.IsNull() )
 148     pidName = "Default" ;
 149   pidName.Append(":") ;
 150   pidName.Append(Version()) ;
 151   SetName(pidName) ;
 152   fPi0Analysis = kFALSE ;
 153 }
 154
 155 //____________________________________________________________________________
 156
 157 void  AliEMCALPIDv1::Exec(Option_t * option)
 158 {
 159   //Steering method
 160
 161   if( strcmp(GetName(), "")== 0 )
 162     Init() ;
 163
 164   if(strstr(option,"tim"))
 165     gBenchmark->Start("EMCALPID");
 166
 167   if(strstr(option,"print")) {
 168     Print("") ;
 169     return ;
 170   }
 171   AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
 172   if(gime->BranchExists("RecParticles") )
 173     return ;
 174   Int_t nevents = gime->MaxEvent() ;       //(Int_t) gAlice->TreeE()->GetEntries() ;
 175   Int_t ievent ;
 176
 177
 178   for(ievent = 0; ievent < nevents; ievent++){
 179     gime->Event(ievent,"R") ;
 180
 181     MakeRecParticles() ;
 182
 183     WriteRecParticles(ievent);
 184
 185     if(strstr(option,"deb"))
 186       PrintRecParticles(option) ;
 187
 188     //increment the total number of rec particles per run
 189     fRecParticlesInRun += gime->RecParticles(BranchName())->GetEntriesFast() ;
 190
 191   }
 192
 193   if(strstr(option,"tim")){
 194     gBenchmark->Stop("EMCALPID");
 195     Info("Exec", "took %f seconds for PID %f seconds per event",
 196          gBenchmark->GetCpuTime("EMCALPID"),
 197          gBenchmark->GetCpuTime("EMCALPID")/nevents) ;
 198   }
 199 }
 200
 201 //____________________________________________________________________________
 202 void  AliEMCALPIDv1::MakeRecParticles(){
 203
 204   // Makes a RecParticle out of a TrackSegment
 205
 206   AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
 207   TObjArray * aECRecPoints = gime->ECALRecPoints() ;
 208   TObjArray * aPRRecPoints = gime->PRERecPoints() ;
 209   TObjArray * aHCRecPoints = gime->HCALRecPoints() ;
 210   TClonesArray * trackSegments = gime->TrackSegments() ;
 211   if ( !aECRecPoints || !aPRRecPoints || !aHCRecPoints || !trackSegments ) {
 212     Fatal("MakeRecParticles", "RecPoints or TrackSegments not found !") ;
 213   }
 214   TClonesArray * recParticles  = gime->RecParticles() ;
 215   recParticles->Clear();
 216
 217   TIter next(trackSegments) ;
 218   AliEMCALTrackSegment * ts ;
 219   Int_t index = 0 ;
 220   AliEMCALRecParticle * rp ;
 221   while ( (ts = (AliEMCALTrackSegment *)next()) ) {
 222
 223     new( (*recParticles)[index] ) AliEMCALRecParticle() ;
 224     rp = (AliEMCALRecParticle *)recParticles->At(index) ;
 225     rp->SetTrackSegment(index) ;
 226     rp->SetIndexInList(index) ;
 227
 228     AliEMCALTowerRecPoint * ecal = 0 ;
 229     if(ts->GetECIndex()>=0)
 230       ecal = dynamic_cast<AliEMCALTowerRecPoint *>(aECRecPoints->At(ts->GetECIndex())) ;
 231
 232     AliEMCALTowerRecPoint    * pre = 0 ;
 233     if(ts->GetPREIndex()>=0)
 234       pre = dynamic_cast<AliEMCALTowerRecPoint *>(aPRRecPoints->At(ts->GetPREIndex())) ;
 235
 236     AliEMCALTowerRecPoint    * hcal = 0 ;
 237     if(ts->GetHCIndex()>=0)
 238       hcal = dynamic_cast<AliEMCALTowerRecPoint *>(aHCRecPoints->At(ts->GetHCIndex())) ;
 239
 240     // Now set type (reconstructed) of the particle
 241
 242     // Choose the cluster energy range
 243
 244     if (!ecal) {
 245       Fatal("MakeRecParticles", "-> emcal(%d) = %d", ts->GetECIndex(), ecal ) ;
 246     }
 247
 248     Float_t    e = ecal->GetEnergy() ;
 249
 250     Float_t  lambda[2] ;
 251     ecal->GetElipsAxis(lambda) ;
 252
 253     if((lambda[0]>0.01) && (lambda[1]>0.01)){
 254       // Looking PCA. Define and calculate the data (X),
 255       // introduce in the function X2P that gives the components (P).
 256
 257       Float_t  Spher = 0. ;
 258       Float_t  Emaxdtotal = 0. ;
 259
 260       if((lambda[0]+lambda[1])!=0)
 261         Spher=fabs(lambda[0]-lambda[1])/(lambda[0]+lambda[1]);
 262
 263       Emaxdtotal=ecal->GetMaximalEnergy()/ecal->GetEnergy();
 264     }
 265
 266     //    Float_t time = ecal->GetTime() ;
 267
 268     //Set momentum, energy and other parameters
 269     Float_t  encal = e;
 270     TVector3 dir(0., 0., 0.) ;
 271     dir.SetMag(encal) ;
 272     rp->SetMomentum(dir.X(),dir.Y(),dir.Z(),encal) ;
 273     rp->SetCalcMass(0);
 274     rp->Name(); //If photon sets the particle pdg name to gamma
 275     rp->SetProductionVertex(0,0,0,0);
 276     rp->SetFirstMother(-1);
 277     rp->SetLastMother(-1);
 278     rp->SetFirstDaughter(-1);
 279     rp->SetLastDaughter(-1);
 280     rp->SetPolarisation(0,0,0);
 281     index++ ;
 282   }
 283
 284 }
 285
 286 //____________________________________________________________________________
 287 void  AliEMCALPIDv1:: Print()
 288 {
 289   // Print the parameters used for the particle type identification
 290
 291   TString message ;
 292     message  = "\n=============== AliEMCALPID1 ================\n" ;
 293     message += "Making PID\n";
 294     message += "    Pricipal analysis file from 0.5 to 100 %s\n" ;
 295     message += "    Name of parameters file     %s\n" ;
 296     message += "    Matrix of Parameters: 9x4\n" ;
 297     message += "        RCPV 2x3 rows x and z, columns function cut parameters\n" ;
 298     message += "        TOF  1x3 [High Eff-Low Pur,Medium Eff-Pur, Low Eff-High Pur]\n" ;
 299     message += "        PCA  5x4 [5 ellipse parametres and 4 parametres to calculate them: A/Sqrt(E) + B* E + C * E^2 + D]\n" ;
 300     message += "        Energy Calibration  1x3 [3 parametres to calibrate energy: A + B* E + C * E^2]\n" ;
 301     Info("Print", message.Data(), fFileName.Data(), fFileNamePar.Data() ) ;
 302 }
 303
 304 //____________________________________________________________________________
 305 void  AliEMCALPIDv1::WriteRecParticles(Int_t event)
 306 {
 307
 308   AliEMCALGetter *gime = AliEMCALGetter::GetInstance() ;
 309
 310   TClonesArray * recParticles = gime->RecParticles() ;
 311   recParticles->Expand(recParticles->GetEntriesFast() ) ;
 312   TTree * treeR ;
 313
 314   if(fToSplit){
 315     if(!fSplitFile)
 316       return ;
 317     fSplitFile->cd() ;
 318     char name[10] ;
 319     sprintf(name,"%s%d", "TreeR",event) ;
 320     treeR = dynamic_cast<TTree*>(fSplitFile->Get(name));
 321   }
 322   else{
 323     treeR = gAlice->TreeR();
 324   }
 325
 326   if(!treeR){
 327     gAlice->MakeTree("R", fSplitFile);
 328     treeR = gAlice->TreeR() ;
 329   }
 330
 331   //First rp
 332   Int_t bufferSize = 32000 ;
 333   TBranch * rpBranch = treeR->Branch("EMCALRP",&recParticles,bufferSize);
 334   rpBranch->SetTitle(BranchName());
 335
 336
 337   //second, pid
 338   Int_t splitlevel = 0 ;
 339   AliEMCALPIDv1 * pid = this ;
 340   TBranch * pidBranch = treeR->Branch("AliEMCALPID","AliEMCALPIDv1",&pid,bufferSize,splitlevel);
 341   pidBranch->SetTitle(BranchName());
 342
 343   rpBranch->Fill() ;
 344   pidBranch->Fill() ;
 345
 346   treeR->AutoSave() ; //Write(0,kOverwrite) ;
 347   if(gAlice->TreeR()!=treeR){
 348     treeR->Delete();
 349   }
 350 }
 351
 352 //____________________________________________________________________________
 353 void AliEMCALPIDv1::PrintRecParticles(Option_t * option)
 354 {
 355   // Print table of reconstructed particles
 356
 357   AliEMCALGetter *gime = AliEMCALGetter::GetInstance() ;
 358
 359   TClonesArray * recParticles = gime->RecParticles(BranchName()) ;
 360
 361   TString message ;
 362   message  = "\nevent " ;
 363   message += gAlice->GetEvNumber() ;
 364   message += "       found " ;
 365   message += recParticles->GetEntriesFast();
 366   message += " RecParticles\n" ;
 367
 368   if(strstr(option,"all")) {  // printing found TS
 369     message += "\n  PARTICLE         Index    \n" ;
 370
 371     Int_t index ;
 372     for (index = 0 ; index < recParticles->GetEntries() ; index++) {
 373       AliEMCALRecParticle * rp = (AliEMCALRecParticle * ) recParticles->At(index) ;
 374       message += "\n" ;
 375       message += rp->Name().Data() ;
 376       message += " " ;
 377       message += rp->GetIndexInList() ;
 378       message += " " ;
 379       message += rp->GetType()  ;
 380     }
 381   }
 382   Info("Print", message.Data() ) ;
 383 }
 384
 385
 386