From: schutz Date: Wed, 8 Jan 2003 17:20:44 +0000 (+0000) Subject: classes for the PID construction X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=commitdiff_plain;h=c9b2f1047ec9ad1ced6bfff1f944be6ddbd26f44;ds=sidebyside classes for the PID construction --- diff --git a/EMCAL/AliEMCALPID.cxx b/EMCAL/AliEMCALPID.cxx new file mode 100644 index 00000000000..64e4c0f50bd --- /dev/null +++ b/EMCAL/AliEMCALPID.cxx @@ -0,0 +1,70 @@ +/************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +/* $Id$ */ + +//_________________________________________________________________________ +// Algorithm class for the identification of particles detected in EMCAL +// base class of identified particle +// Why should I put meaningless comments +// just to satisfy +// the code checker + +// +//*-- Author: Yves Schutz (SUBATECH) & Dmitri Peressounko + + +// --- ROOT system --- +#include "TGeometry.h" +#include "TDirectory.h" +#include "TFile.h" +#include "TTree.h" + +// --- Standard library --- +#include + + +// --- AliRoot header files --- +#include "AliRun.h" +#include "AliEMCALPID.h" +#include "AliHeader.h" + +ClassImp(AliEMCALPID) + +//____________________________________________________________________________ + AliEMCALPID::AliEMCALPID():TTask("","") +{ + // ctor + fSplitFile= 0 ; + +} + + +//____________________________________________________________________________ +AliEMCALPID::AliEMCALPID(const char* headerFile, const char * name, const Bool_t toSplit):TTask(name, headerFile) +{ + // ctor + + fToSplit = toSplit ; + fSplitFile= 0 ; +} + +//____________________________________________________________________________ +AliEMCALPID::~AliEMCALPID() +{ + // dtor + + fSplitFile = 0 ; +} diff --git a/EMCAL/AliEMCALPID.h b/EMCAL/AliEMCALPID.h new file mode 100644 index 00000000000..d35f0ad917b --- /dev/null +++ b/EMCAL/AliEMCALPID.h @@ -0,0 +1,55 @@ +#ifndef ALIEMCALPID_H +#define ALIEMCALPID_H +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +/* $Id$ */ + +//_________________________________________________________________________ +// Algorithm class for the identification of particles detected in EMCAL +// base class +// of identified particles +//*-- Author: Yves Schutz (SUBATECH) + +// --- ROOT system --- + +#include "TTask.h" +class TFormula ; +class TClonesArray ; + +// --- Standard library --- + +// --- AliRoot header files --- + +class AliEMCALGeometry ; +class AliEMCALClusterizer ; +class AliEMCALTrackSegmentMaker ; + +class AliEMCALPID : public TTask { + + public: + + AliEMCALPID() ; // ctor + AliEMCALPID(const char* headerFile,const char * name, const Bool_t toSplit) ; + virtual ~AliEMCALPID() ; // dtor + + virtual void Exec(Option_t * option) { Warning("Exec", "not defined" ) ; } + virtual const Int_t GetRecParticlesInRun() const { Warning("GetRecParticlesInRun", "not defined" ) ; return 0 ;} + virtual void Print(Option_t * option) const { Warning("Print", "not defined" ) ;} + virtual void SetCpvtoEmcDistanceCut(Float_t Cluster_En, TString Eff_Pur,Float_t cut ) { Warning("SetCpvtoEmcDistanceCut", "not defined" ) ;} + virtual void SetTimeGate(Float_t Cluster_En, TString Eff_Pur, Float_t gate) { Warning("SetTimeGate", "not defined" ) ; } + virtual const char * Version() const { Warning("Version", "not defined" ) ; return 0 ; } + virtual void WriteRecParticles(Int_t event) { Warning("WriteRecParticles", "not defined" ) ; } + +private: + virtual void Init() { Warning("Init", "not defined" ) ; } + +protected: + + TFile * fSplitFile ; //! file in which RecParticles will eventually be stored + Bool_t fToSplit ; //! do we in the split mode + ClassDef(AliEMCALPID,1) // Particle Identifier algorithm (base class) + +} ; + +#endif // ALIEMCALPID_H diff --git a/EMCAL/AliEMCALPIDv1.cxx b/EMCAL/AliEMCALPIDv1.cxx new file mode 100644 index 00000000000..7184d7b1c78 --- /dev/null +++ b/EMCAL/AliEMCALPIDv1.cxx @@ -0,0 +1,386 @@ +/************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +/* $Id$ */ + +//_________________________________________________________________________ +// Implementation version v1 of the EMCAL particle identifier + +//*-- Author: Yves Schutz (SUBATECH) +// +// --- ROOT system --- +#include "TROOT.h" +#include "TTree.h" +#include "TFile.h" +#include "TF2.h" +#include "TCanvas.h" +#include "TFolder.h" +#include "TSystem.h" +#include "TBenchmark.h" + +#include "TSystem.h" + +// --- Standard library --- + +#include + +// --- AliRoot header files --- + +#include "AliRun.h" +#include "AliGenerator.h" +#include "AliEMCAL.h" +#include "AliEMCALPIDv1.h" +#include "AliEMCALClusterizerv1.h" +#include "AliEMCALTrackSegment.h" +#include "AliEMCALTrackSegmentMakerv1.h" +#include "AliEMCALRecParticle.h" +#include "AliEMCALGeometry.h" +#include "AliEMCALGetter.h" + +ClassImp( AliEMCALPIDv1) + +//____________________________________________________________________________ +AliEMCALPIDv1::AliEMCALPIDv1():AliEMCALPID() +{ + // default ctor + + InitParameters() ; + fDefaultInit = kTRUE ; + +} + +//____________________________________________________________________________ +AliEMCALPIDv1::AliEMCALPIDv1(const char * headerFile,const char * name, const Bool_t toSplit) +:AliEMCALPID(headerFile, name,toSplit) +{ + //ctor with the indication on where to look for the track segments + + InitParameters() ; + + Init() ; + fDefaultInit = kFALSE ; + +} + +//____________________________________________________________________________ +AliEMCALPIDv1::~AliEMCALPIDv1() +{ + // dtor + + if (!fDefaultInit) { + fSplitFile = 0 ; + } +} + +//____________________________________________________________________________ +const TString AliEMCALPIDv1::BranchName() const +{ + TString branchName(GetName() ) ; + branchName.Remove(branchName.Index(Version())-1) ; + return branchName ; +} + +//____________________________________________________________________________ +void AliEMCALPIDv1::Init() +{ + // Make all memory allocations that are not possible in default constructor + // Add the PID task to the list of EMCAL tasks + + if ( strcmp(GetTitle(), "") == 0 ) + SetTitle("galice.root") ; + + TString branchname(GetName()) ; + branchname.Remove(branchname.Index(Version())-1) ; + AliEMCALGetter * gime = AliEMCALGetter::GetInstance(GetTitle(),branchname.Data(),fToSplit ) ; + + // gime->SetRecParticlesTitle(BranchName()) ; + if ( gime == 0 ) { + Error("Init", "Could not obtain the Getter object !" ) ; + return ; + } + + fSplitFile = 0 ; + if(fToSplit){ + //First - extract full path if necessary + TString fileName(GetTitle()) ; + Ssiz_t islash = fileName.Last('/') ; + if(islash(gROOT->GetFile(fileName.Data())); + if(!fSplitFile) + fSplitFile = TFile::Open(fileName.Data(),"update") ; + } + + gime->PostPID(this) ; + gime->PostRecParticles(branchname) ; + +} + +//____________________________________________________________________________ +void AliEMCALPIDv1::InitParameters() +{ + + fRecParticlesInRun = 0 ; + fNEvent = 0 ; + fRecParticlesInRun = 0 ; + TString pidName( GetName()) ; + if (pidName.IsNull() ) + pidName = "Default" ; + pidName.Append(":") ; + pidName.Append(Version()) ; + SetName(pidName) ; + fPi0Analysis = kFALSE ; +} + +//____________________________________________________________________________ + +void AliEMCALPIDv1::Exec(Option_t * option) +{ + //Steering method + + if( strcmp(GetName(), "")== 0 ) + Init() ; + + if(strstr(option,"tim")) + gBenchmark->Start("EMCALPID"); + + if(strstr(option,"print")) { + Print("") ; + return ; + } + AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ; + if(gime->BranchExists("RecParticles") ) + return ; + Int_t nevents = gime->MaxEvent() ; //(Int_t) gAlice->TreeE()->GetEntries() ; + Int_t ievent ; + + + for(ievent = 0; ievent < nevents; ievent++){ + gime->Event(ievent,"R") ; + + MakeRecParticles() ; + + WriteRecParticles(ievent); + + if(strstr(option,"deb")) + PrintRecParticles(option) ; + + //increment the total number of rec particles per run + fRecParticlesInRun += gime->RecParticles(BranchName())->GetEntriesFast() ; + + } + + if(strstr(option,"tim")){ + gBenchmark->Stop("EMCALPID"); + Info("Exec", "took %f seconds for PID %f seconds per event", + gBenchmark->GetCpuTime("EMCALPID"), + gBenchmark->GetCpuTime("EMCALPID")/nevents) ; + } +} + +//____________________________________________________________________________ +void AliEMCALPIDv1::MakeRecParticles(){ + + // Makes a RecParticle out of a TrackSegment + + AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ; + TObjArray * aECRecPoints = gime->ECALRecPoints() ; + TObjArray * aPRRecPoints = gime->PRERecPoints() ; + TObjArray * aHCRecPoints = gime->HCALRecPoints() ; + TClonesArray * trackSegments = gime->TrackSegments() ; + if ( !aECRecPoints || !aPRRecPoints || !aHCRecPoints || !trackSegments ) { + Fatal("MakeRecParticles", "RecPoints or TrackSegments not found !") ; + } + TClonesArray * recParticles = gime->RecParticles() ; + recParticles->Clear(); + + TIter next(trackSegments) ; + AliEMCALTrackSegment * ts ; + Int_t index = 0 ; + AliEMCALRecParticle * rp ; + while ( (ts = (AliEMCALTrackSegment *)next()) ) { + + new( (*recParticles)[index] ) AliEMCALRecParticle() ; + rp = (AliEMCALRecParticle *)recParticles->At(index) ; + rp->SetTrackSegment(index) ; + rp->SetIndexInList(index) ; + + AliEMCALTowerRecPoint * ecal = 0 ; + if(ts->GetECIndex()>=0) + ecal = dynamic_cast(aECRecPoints->At(ts->GetECIndex())) ; + + AliEMCALTowerRecPoint * pre = 0 ; + if(ts->GetPREIndex()>=0) + pre = dynamic_cast(aPRRecPoints->At(ts->GetPREIndex())) ; + + AliEMCALTowerRecPoint * hcal = 0 ; + if(ts->GetHCIndex()>=0) + hcal = dynamic_cast(aHCRecPoints->At(ts->GetHCIndex())) ; + + // Now set type (reconstructed) of the particle + + // Choose the cluster energy range + + if (!ecal) { + Fatal("MakeRecParticles", "-> emcal(%d) = %d", ts->GetECIndex(), ecal ) ; + } + + Float_t e = ecal->GetEnergy() ; + + Float_t lambda[2] ; + ecal->GetElipsAxis(lambda) ; + + if((lambda[0]>0.01) && (lambda[1]>0.01)){ + // Looking PCA. Define and calculate the data (X), + // introduce in the function X2P that gives the components (P). + + Float_t Spher = 0. ; + Float_t Emaxdtotal = 0. ; + + if((lambda[0]+lambda[1])!=0) + Spher=fabs(lambda[0]-lambda[1])/(lambda[0]+lambda[1]); + + Emaxdtotal=ecal->GetMaximalEnergy()/ecal->GetEnergy(); + } + + // Float_t time = ecal->GetTime() ; + + //Set momentum, energy and other parameters + Float_t encal = e; + TVector3 dir(0., 0., 0.) ; + dir.SetMag(encal) ; + rp->SetMomentum(dir.X(),dir.Y(),dir.Z(),encal) ; + rp->SetCalcMass(0); + rp->Name(); //If photon sets the particle pdg name to gamma + rp->SetProductionVertex(0,0,0,0); + rp->SetFirstMother(-1); + rp->SetLastMother(-1); + rp->SetFirstDaughter(-1); + rp->SetLastDaughter(-1); + rp->SetPolarisation(0,0,0); + index++ ; + } + +} + +//____________________________________________________________________________ +void AliEMCALPIDv1:: Print() +{ + // Print the parameters used for the particle type identification + + TString message ; + message = "\n=============== AliEMCALPID1 ================\n" ; + message += "Making PID\n"; + message += " Pricipal analysis file from 0.5 to 100 %s\n" ; + message += " Name of parameters file %s\n" ; + message += " Matrix of Parameters: 9x4\n" ; + message += " RCPV 2x3 rows x and z, columns function cut parameters\n" ; + message += " TOF 1x3 [High Eff-Low Pur,Medium Eff-Pur, Low Eff-High Pur]\n" ; + message += " PCA 5x4 [5 ellipse parametres and 4 parametres to calculate them: A/Sqrt(E) + B* E + C * E^2 + D]\n" ; + message += " Energy Calibration 1x3 [3 parametres to calibrate energy: A + B* E + C * E^2]\n" ; + Info("Print", message.Data(), fFileName.Data(), fFileNamePar.Data() ) ; +} + +//____________________________________________________________________________ +void AliEMCALPIDv1::WriteRecParticles(Int_t event) +{ + + AliEMCALGetter *gime = AliEMCALGetter::GetInstance() ; + + TClonesArray * recParticles = gime->RecParticles() ; + recParticles->Expand(recParticles->GetEntriesFast() ) ; + TTree * treeR ; + + if(fToSplit){ + if(!fSplitFile) + return ; + fSplitFile->cd() ; + char name[10] ; + sprintf(name,"%s%d", "TreeR",event) ; + treeR = dynamic_cast(fSplitFile->Get(name)); + } + else{ + treeR = gAlice->TreeR(); + } + + if(!treeR){ + gAlice->MakeTree("R", fSplitFile); + treeR = gAlice->TreeR() ; + } + + //First rp + Int_t bufferSize = 32000 ; + TBranch * rpBranch = treeR->Branch("EMCALRP",&recParticles,bufferSize); + rpBranch->SetTitle(BranchName()); + + + //second, pid + Int_t splitlevel = 0 ; + AliEMCALPIDv1 * pid = this ; + TBranch * pidBranch = treeR->Branch("AliEMCALPID","AliEMCALPIDv1",&pid,bufferSize,splitlevel); + pidBranch->SetTitle(BranchName()); + + rpBranch->Fill() ; + pidBranch->Fill() ; + + treeR->AutoSave() ; //Write(0,kOverwrite) ; + if(gAlice->TreeR()!=treeR){ + treeR->Delete(); + } +} + +//____________________________________________________________________________ +void AliEMCALPIDv1::PrintRecParticles(Option_t * option) +{ + // Print table of reconstructed particles + + AliEMCALGetter *gime = AliEMCALGetter::GetInstance() ; + + TClonesArray * recParticles = gime->RecParticles(BranchName()) ; + + TString message ; + message = "\nevent " ; + message += gAlice->GetEvNumber() ; + message += " found " ; + message += recParticles->GetEntriesFast(); + message += " RecParticles\n" ; + + if(strstr(option,"all")) { // printing found TS + message += "\n PARTICLE Index \n" ; + + Int_t index ; + for (index = 0 ; index < recParticles->GetEntries() ; index++) { + AliEMCALRecParticle * rp = (AliEMCALRecParticle * ) recParticles->At(index) ; + message += "\n" ; + message += rp->Name().Data() ; + message += " " ; + message += rp->GetIndexInList() ; + message += " " ; + message += rp->GetType() ; + } + } + Info("Print", message.Data() ) ; +} + + + diff --git a/EMCAL/AliEMCALPIDv1.h b/EMCAL/AliEMCALPIDv1.h new file mode 100644 index 00000000000..fc89e4451c4 --- /dev/null +++ b/EMCAL/AliEMCALPIDv1.h @@ -0,0 +1,77 @@ +#ifndef ALIEMCALPIDV1_H +#define ALIEMCALPIDV1_H +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +/* $Id$ */ + + +//_________________________________________________________________________ +// Implementation version v1 of the EMCAL particle identifier +// Identification is based on information from CPV and EMC +// Oh yeah +//*-- Author: Yves Schutz (SUBATECH), Gustavo Conesa. + +// --- ROOT system --- +//class TFormula ; +class TVector3 ; +class TMatrixD ; +class TPrincipal ; + +// --- Standard library --- + +// --- AliRoot header files --- +class AliEMCALEmcRecPoint ; +class AliEMCALRecPoint ; + +#include "AliEMCALPID.h" + +class AliEMCALPIDv1 : public AliEMCALPID { + + public: + + AliEMCALPIDv1() ; // ctor + AliEMCALPIDv1(const char* headerFile, const char * tsBranch = "Default", const Bool_t toSplit=kFALSE) ; + + virtual ~AliEMCALPIDv1() ; // dtor + + virtual void Exec(Option_t * option) ; + // virtual char * GetRecParticlesBranch()const {return (char*) fRecParticlesTitle.Data() ;} + // virtual char * GetTrackSegmentsBranch()const{return (char*) fTrackSegmentsTitle.Data(); } + virtual const Int_t GetRecParticlesInRun() const {return fRecParticlesInRun ;} + + virtual void Print(Option_t * option) const {} + void Print() ; + + + //To turn on or off the Pi0 analysis + const Bool_t GetPi0Analysis(){return fPi0Analysis;} + void SetPi0Analysis(Bool_t turnonoff){ fPi0Analysis = turnonoff; } + + virtual const char * Version() const { return "pid-v1" ; } + + private: + + const TString BranchName() const ; + virtual void Init() ; + virtual void InitParameters() ; + void MakeRecParticles(void ) ; + void PrintRecParticles(Option_t * option) ; + virtual void WriteRecParticles(Int_t event) ; + + + + private: + + Bool_t fDefaultInit; //! Says if the task was created by defaut ctor (only parameters are initialized) + TString fFileName ; // File that contains the Principal file for analysis + TString fFileNamePar ;// File that contains the parameters for analysis + Int_t fNEvent ; //! current event number + Bool_t fPi0Analysis; //! Pi0 analysis on or off + Int_t fRecParticlesInRun ; //! Total number of recparticles in one run + + ClassDef( AliEMCALPIDv1,6) // Particle identifier implementation version 1 + +}; + +#endif // AliEMCALPIDV1_H