X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PHOS%2FAliPHOSPIDv0.cxx;h=a393a4f9988cedc9b031733dd014ad773c217abb;hb=dddc33c571b008eccdc8f659de11ae42fc5a5509;hp=d1e44f4fa0681838444a3bb35c13ed0b75b83052;hpb=a278df550f3b5c4fcab0b857805d970be80252d9;p=u%2Fmrichter%2FAliRoot.git diff --git a/PHOS/AliPHOSPIDv0.cxx b/PHOS/AliPHOSPIDv0.cxx index d1e44f4fa06..a393a4f9988 100644 --- a/PHOS/AliPHOSPIDv0.cxx +++ b/PHOS/AliPHOSPIDv0.cxx @@ -15,6 +15,20 @@ /* $Id$ */ +/* History of cvs commits: + * + * $Log$ + * Revision 1.15 2007/03/06 06:57:05 kharlov + * DP:calculation of distance to CPV done in TSM + * + * Revision 1.14 2006/09/07 18:31:08 kharlov + * Effective c++ corrections (T.Pocheptsov) + * + * Revision 1.13 2005/05/28 14:19:04 schutz + * Compilation warnings fixed by T.P. + * + */ + //_________________________________________________________________________ // Implementation version v0 of the PHOS particle identifier // Particle identification based on the @@ -57,220 +71,156 @@ // Completely redesined by Dmitri Peressounko, March 2001 // --- ROOT system --- -#include "TROOT.h" #include "TTree.h" -#include "TFile.h" #include "TF2.h" #include "TFormula.h" #include "TCanvas.h" -#include "TFolder.h" -#include "TSystem.h" +#include "TClonesArray.h" + #include "TBenchmark.h" // --- Standard library --- -#include -#include - // --- AliRoot header files --- - -#include "AliRun.h" -#include "AliGenerator.h" -#include "AliPHOS.h" +#include "AliLog.h" #include "AliPHOSPIDv0.h" -#include "AliPHOSClusterizerv1.h" +#include "AliPHOSEmcRecPoint.h" #include "AliPHOSTrackSegment.h" -#include "AliPHOSTrackSegmentMakerv1.h" #include "AliPHOSRecParticle.h" #include "AliPHOSGeometry.h" -#include "AliPHOSGetter.h" ClassImp( AliPHOSPIDv0) //____________________________________________________________________________ -AliPHOSPIDv0::AliPHOSPIDv0():AliPHOSPID() +AliPHOSPIDv0::AliPHOSPIDv0(): + fIDOptions("dis time"), + fClusterizer(0), + fTSMaker(0), + fFormula(0), + fDispersion(0.f), + fCpvEmcDistance(0.f), + fTimeGate(2.e-9f) { // default ctor - fFormula = 0 ; - fDispersion = 0. ; - fCpvEmcDistance = 0 ; - fTimeGate = 2.e-9 ; - fHeaderFileName = "" ; - fTrackSegmentsTitle= "" ; - fRecPointsTitle = "" ; - fRecParticlesTitle = "" ; - fIDOptions = "dis time" ; - fRecParticlesInRun = 0 ; - fClusterizer = 0; - fTSMaker = 0; } //____________________________________________________________________________ -AliPHOSPIDv0::AliPHOSPIDv0(const char * headerFile,const char * name) : AliPHOSPID(headerFile, name) +AliPHOSPIDv0::AliPHOSPIDv0(AliPHOSGeometry *geom) : + AliPHOSPID(geom), + fIDOptions("dis time"), + fClusterizer(0), + fTSMaker(0), + fFormula(new TFormula("LambdaCuts","(x>1)*(x<2.5)*(y>0)*(y1)*(x<2.5)*(y>0)*(yPHOSGeometry() ; - TVector3 vecEmc ; - TVector3 vecCpv ; - - emc->GetLocalPosition(vecEmc) ; - cpv->GetLocalPosition(vecCpv) ; - if(emc->GetPHOSMod() == cpv->GetPHOSMod()){ - - // Correct to difference in CPV and EMC position due to different distance to center. - // we assume, that particle moves from center - Float_t dCPV = geom->GetIPtoOuterCoverDistance(); - Float_t dEMC = geom->GetIPtoCrystalSurface() ; - dEMC = dEMC / dCPV ; - vecCpv = dEMC * vecCpv - vecEmc ; - if (Axis == "X") return vecCpv.X(); - if (Axis == "Y") return vecCpv.Y(); - if (Axis == "Z") return vecCpv.Z(); - if (Axis == "R") return vecCpv.Mag(); - } - - return 100000000 ; + //Empty dtor, fFormula leaks } +//DP +////____________________________________________________________________________ +//Float_t AliPHOSPIDv0::GetDistance(AliPHOSEmcRecPoint * emc,AliPHOSRecPoint * cpv, Option_t * Axis)const +//{ +// // Calculates the distance between the EMC RecPoint and the PPSD RecPoint +// +// const AliPHOSGeometry * geom = AliPHOSLoader::GetPHOSGeometry() ; +// TVector3 vecEmc ; +// TVector3 vecCpv ; +// +// emc->GetLocalPosition(vecEmc) ; +// cpv->GetLocalPosition(vecCpv) ; +// if(emc->GetPHOSMod() == cpv->GetPHOSMod()){ +// +// // Correct to difference in CPV and EMC position due to different distance to center. +// // we assume, that particle moves from center +// Float_t dCPV = geom->GetIPtoOuterCoverDistance(); +// Float_t dEMC = geom->GetIPtoCrystalSurface() ; +// dEMC = dEMC / dCPV ; +// vecCpv = dEMC * vecCpv - vecEmc ; +// if (Axis == "X") return vecCpv.X(); +// if (Axis == "Y") return vecCpv.Y(); +// if (Axis == "Z") return vecCpv.Z(); +// if (Axis == "R") return vecCpv.Mag(); +// } +// +// return 100000000 ; +//} + //____________________________________________________________________________ -void AliPHOSPIDv0::Exec(Option_t * option) +void AliPHOSPIDv0::TrackSegments2RecParticles(Option_t * option) { //Steering method - - if( strcmp(GetName(), "")== 0 ) - Init() ; - + if(strstr(option,"tim")) gBenchmark->Start("PHOSPID"); if(strstr(option,"print")) { - Print("") ; + Print() ; return ; } - gAlice->GetEvent(0) ; - //check, if the branch with name of this" already exits? - TObjArray * lob = (TObjArray*)gAlice->TreeR()->GetListOfBranches() ; - TIter next(lob) ; - TBranch * branch = 0 ; - Bool_t phospidfound = kFALSE, pidfound = kFALSE ; - - TString taskName(GetName()) ; - taskName.Remove(taskName.Index(Version())-1) ; + AliInfo(Form("%d emc clusters, %d track segments", + fEMCRecPoints->GetEntriesFast(), + fTrackSegments->GetEntriesFast())) ; - while ( (branch = (TBranch*)next()) && (!phospidfound || !pidfound) ) { - if ( (strcmp(branch->GetName(), "PHOSPID")==0) && (strcmp(branch->GetTitle(), taskName.Data())==0) ) - phospidfound = kTRUE ; - - else if ( (strcmp(branch->GetName(), "AliPHOSPID")==0) && (strcmp(branch->GetTitle(), taskName.Data())==0) ) - pidfound = kTRUE ; - } - - if ( phospidfound || pidfound ) { - cerr << "WARNING: AliPHOSPIDv0::Exec -> RecParticles and/or PIDtMaker branch with name " - << taskName.Data() << " already exits" << endl ; - return ; - } - - Int_t nevents = (Int_t) gAlice->TreeE()->GetEntries() ; - Int_t ievent ; - AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ; - - for(ievent = 0; ievent < nevents; ievent++){ - gime->Event(ievent,"R") ; - - MakeRecParticles() ; - - WriteRecParticles(ievent); + MakeRecParticles() ; - if(strstr(option,"deb")) - PrintRecParticles(option) ; + if(strstr(option,"deb")) + PrintRecParticles(option) ; - //increment the total number of rec particles per run - fRecParticlesInRun += gime->RecParticles()->GetEntriesFast() ; - - } - if(strstr(option,"tim")){ gBenchmark->Stop("PHOSPID"); - cout << "AliPHOSPID:" << endl ; - cout << " took " << gBenchmark->GetCpuTime("PHOSPID") << " seconds for PID " - << gBenchmark->GetCpuTime("PHOSPID")/nevents << " seconds per event " << endl ; - cout << endl ; + AliInfo(Form("took %f seconds for PID", + gBenchmark->GetCpuTime("PHOSPID"))); } } + //____________________________________________________________________________ -void AliPHOSPIDv0::Init() +void AliPHOSPIDv0::MakeRecParticles() { - // Make all memory allocations that are not possible in default constructor - // Add the PID task to the list of PHOS tasks - - if ( strcmp(GetTitle(), "") == 0 ) - SetTitle("galice.root") ; - - TString taskName(GetName()) ; - taskName.Remove(taskName.Index(Version())-1) ; - - AliPHOSGetter * gime = AliPHOSGetter::GetInstance(GetTitle(), taskName.Data()) ; - if ( gime == 0 ) { - cerr << "ERROR: AliPHOSPIDv0::Init -> Could not obtain the Getter object !" << endl ; - return ; - } - - gime->PostPID(this) ; - // create a folder on the white board //YSAlice/WhiteBoard/RecParticles/PHOS/recparticlesName - gime->PostRecParticles(taskName.Data() ) ; - -} + // Reconstructs the particles from the tracksegments -//____________________________________________________________________________ -void AliPHOSPIDv0::MakeRecParticles(){ - - // Makes a RecParticle out of a TrackSegment - TString taskName(GetName()) ; - taskName.Remove(taskName.Index(Version())-1) ; - - AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ; - TObjArray * emcRecPoints = gime->EmcRecPoints(taskName) ; - TObjArray * cpvRecPoints = gime->CpvRecPoints(taskName) ; - TClonesArray * trackSegments = gime->TrackSegments(taskName) ; - TClonesArray * recParticles = gime->RecParticles(taskName) ; - recParticles->Clear(); + fRecParticles->Clear(); - TIter next(trackSegments) ; + TIter next(fTrackSegments) ; AliPHOSTrackSegment * ts ; Int_t index = 0 ; AliPHOSRecParticle * rp ; @@ -281,18 +231,20 @@ void AliPHOSPIDv0::MakeRecParticles(){ while ( (ts = (AliPHOSTrackSegment *)next()) ) { - new( (*recParticles)[index] ) AliPHOSRecParticle() ; - rp = (AliPHOSRecParticle *)recParticles->At(index) ; + new( (*fRecParticles)[index] ) AliPHOSRecParticle() ; + rp = (AliPHOSRecParticle *)fRecParticles->At(index) ; rp->SetTrackSegment(index) ; rp->SetIndexInList(index) ; AliPHOSEmcRecPoint * emc = 0 ; if(ts->GetEmcIndex()>=0) - emc = (AliPHOSEmcRecPoint *) emcRecPoints->At(ts->GetEmcIndex()) ; + emc = (AliPHOSEmcRecPoint *) fEMCRecPoints->At(ts->GetEmcIndex()) ; + else + AliFatal(Form("ts->GetEmcIndex() return illegal index %d",ts->GetEmcIndex())); AliPHOSRecPoint * cpv = 0 ; if(ts->GetCpvIndex()>=0) - cpv = (AliPHOSRecPoint *) cpvRecPoints->At(ts->GetCpvIndex()) ; + cpv = (AliPHOSRecPoint *) fCPVRecPoints->At(ts->GetCpvIndex()) ; //set momentum and energy first Float_t e = emc->GetEnergy() ; @@ -324,7 +276,7 @@ void AliPHOSPIDv0::MakeRecParticles(){ // Looking at the CPV detector Int_t cpvdetector= 0 ; //1 hit and 0 no hit if(cpv) - if(GetDistance(emc, cpv, "R") < fCpvEmcDistance) + if(ts->GetCpvDistance("R") < fCpvEmcDistance) cpvdetector = 1 ; Int_t type = showerprofile + 2 * slow + 4 * cpvdetector ; @@ -341,30 +293,28 @@ void AliPHOSPIDv0::MakeRecParticles(){ } //____________________________________________________________________________ -void AliPHOSPIDv0:: Print(Option_t * option) const +void AliPHOSPIDv0:: Print(const Option_t *) const { // Print the parameters used for the particle type identification - cout << "=============== AliPHOSPID1 ================" << endl ; - cout << "Making PID "<< endl ; - cout << " Headers file: " << fHeaderFileName.Data() << endl ; - cout << " RecPoints branch title: " << fRecPointsTitle.Data() << endl ; - cout << " TrackSegments Branch title: " << fTrackSegmentsTitle.Data() << endl ; - cout << " RecParticles Branch title " << fRecParticlesTitle.Data() << endl; - cout << "with parameters: " << endl ; - cout << " Maximal EMC - CPV distance (cm) " << fCpvEmcDistance << endl ; - if(fIDOptions.Contains("dis",TString::kIgnoreCase )) - cout << " dispersion cut " << fDispersion << endl ; - if(fIDOptions.Contains("ell",TString::kIgnoreCase )){ - cout << " Eliptic cuts function: " << endl ; - cout << fFormula->GetTitle() << endl ; - } - if(fIDOptions.Contains("tim",TString::kIgnoreCase )) - cout << " Time Gate uzed: " << fTimeGate << endl ; - cout << "============================================" << endl ; + TString message ; + message = "=============== AliPHOSPIDv0 ================\n" ; + message += "Making PID\n" ; + message += "with parameters:\n" ; + message += " Maximal EMC - CPV distance (cm) %f\n" ; + AliInfo(Form( message.Data(), + fCpvEmcDistance )); + + if(fIDOptions.Contains("dis",TString::kIgnoreCase )) + AliInfo(Form(" dispersion cut %f", fDispersion )) ; + if(fIDOptions.Contains("ell",TString::kIgnoreCase )) + AliInfo(Form(" Eliptic cuts function: %s", + fFormula->GetTitle() )) ; + if(fIDOptions.Contains("tim",TString::kIgnoreCase )) + AliInfo(Form(" Time Gate used: %f", fTimeGate)) ; } //____________________________________________________________________________ -void AliPHOSPIDv0::SetShowerProfileCut(char * formula) +void AliPHOSPIDv0::SetShowerProfileCut(const char * formula) { //set shape of the cut on the axis of ellipce, drown around shouer //shower considered "narrow" if Formula(lambda[0],lambda[1]) > 0. @@ -372,61 +322,6 @@ void AliPHOSPIDv0::SetShowerProfileCut(char * formula) delete fFormula; fFormula = new TFormula("Lambda Cut",formula) ; } -//____________________________________________________________________________ -void AliPHOSPIDv0::WriteRecParticles(Int_t event) -{ - - AliPHOSGetter *gime = AliPHOSGetter::GetInstance() ; - TString taskName(GetName()) ; - taskName.Remove(taskName.Index(Version())-1) ; - TClonesArray * recParticles = gime->RecParticles(taskName) ; - recParticles->Expand(recParticles->GetEntriesFast() ) ; - - //Make branch in TreeR for RecParticles - char * filename = 0; - if(gSystem->Getenv("CONFIG_SPLIT_FILE")!=0){ //generating file name - filename = new char[strlen(gAlice->GetBaseFile())+20] ; - sprintf(filename,"%s/PHOS.Reco.root",gAlice->GetBaseFile()) ; - } - - TDirectory *cwd = gDirectory; - - //First rp - Int_t bufferSize = 32000 ; - TBranch * rpBranch = gAlice->TreeR()->Branch("PHOSRP",&recParticles,bufferSize); - rpBranch->SetTitle(fRecParticlesTitle); - if (filename) { - rpBranch->SetFile(filename); - TIter next( rpBranch->GetListOfBranches()); - TBranch * sb ; - while ((sb=(TBranch*)next())) { - sb->SetFile(filename); - } - cwd->cd(); - } - - //second, pid - Int_t splitlevel = 0 ; - AliPHOSPIDv0 * pid = this ; - TBranch * pidBranch = gAlice->TreeR()->Branch("AliPHOSPID","AliPHOSPIDv0",&pid,bufferSize,splitlevel); - pidBranch->SetTitle(fRecParticlesTitle.Data()); - if (filename) { - pidBranch->SetFile(filename); - TIter next( pidBranch->GetListOfBranches()); - TBranch * sb ; - while ((sb=(TBranch*)next())) { - sb->SetFile(filename); - } - cwd->cd(); - } - - rpBranch->Fill() ; - pidBranch->Fill() ; - - gAlice->TreeR()->Write(0,kOverwrite) ; - - delete [] filename ; -} //____________________________________________________________________________ void AliPHOSPIDv0::PlotDispersionCuts()const @@ -462,7 +357,7 @@ void AliPHOSPIDv0::PlotDispersionCuts()const } //____________________________________________________________________________ -TVector3 AliPHOSPIDv0::GetMomentumDirection(AliPHOSEmcRecPoint * emc, AliPHOSRecPoint * cpv)const +TVector3 AliPHOSPIDv0::GetMomentumDirection(AliPHOSEmcRecPoint * emc, AliPHOSRecPoint * )const { // Calculates the momentum direction: // 1. if only a EMC RecPoint, direction is given by IP and this RecPoint @@ -498,11 +393,14 @@ TVector3 AliPHOSPIDv0::GetMomentumDirection(AliPHOSEmcRecPoint * emc, AliPHOSRec dir.SetZ( -dir.Z() ) ; // why ? dir.SetMag(1.) ; + // One can not access MC information in the reconstruction!! + // PLEASE FIT IT, EITHER BY TAKING 0,0,0 OR ACCESSING THE + // VERTEX DIAMOND FROM CDB GRP FOLDER. //account correction to the position of IP - Float_t xo,yo,zo ; //Coordinates of the origin - gAlice->Generator()->GetOrigin(xo,yo,zo) ; - TVector3 origin(xo,yo,zo); - dir = dir - origin ; + // Float_t xo,yo,zo ; //Coordinates of the origin + // gAlice->Generator()->GetOrigin(xo,yo,zo) ; + // TVector3 origin(xo,yo,zo); + // dir = dir - origin ; return dir ; } @@ -511,54 +409,43 @@ void AliPHOSPIDv0::PrintRecParticles(Option_t * option) { // Print table of reconstructed particles - AliPHOSGetter *gime = AliPHOSGetter::GetInstance() ; + TString message ; + message = "Found %d RecParticles\n" ; + AliInfo(Form(message.Data(), + fRecParticles->GetEntriesFast() )) ; - TString taskName(GetName()) ; - taskName.Remove(taskName.Index(Version())-1) ; - TClonesArray * recParticles = gime->RecParticles(taskName) ; - - cout << "AliPHOSPIDv0: event "<GetEvNumber() << endl ; - cout << " found " << recParticles->GetEntriesFast() << " RecParticles " << endl ; - if(strstr(option,"all")) { // printing found TS - - cout << " PARTICLE " - << " Index " << endl ; - // << " X " - // << " Y " - // << " Z " - // << " # of primaries " - // << " Primaries list " << endl; - + AliInfo(" PARTICLE Index" ) ; + Int_t index ; - for (index = 0 ; index < recParticles->GetEntries() ; index++) { - AliPHOSRecParticle * rp = (AliPHOSRecParticle * ) recParticles->At(index) ; + for (index = 0 ; index < fRecParticles->GetEntries() ; index++) { + AliPHOSRecParticle * rp = (AliPHOSRecParticle * ) fRecParticles->At(index) ; - Text_t particle[11]; + Text_t particle[100]; switch(rp->GetType()) { case AliPHOSFastRecParticle::kNEUTRALEMFAST: - strcpy( particle, "NEUTRAL EM FAST"); + strncpy( particle, "NEUTRAL EM FAST",100); break; case AliPHOSFastRecParticle::kNEUTRALHAFAST: - strcpy(particle, "NEUTRAL HA FAST"); + strncpy(particle, "NEUTRAL HA FAST",100); break; case AliPHOSFastRecParticle::kNEUTRALEMSLOW: - strcpy(particle, "NEUTRAL EM SLOW"); + strncpy(particle, "NEUTRAL EM SLOW",100); break ; case AliPHOSFastRecParticle::kNEUTRALHASLOW: - strcpy(particle, "NEUTRAL HA SLOW"); + strncpy(particle, "NEUTRAL HA SLOW",100); break ; case AliPHOSFastRecParticle::kCHARGEDEMFAST: - strcpy(particle, "CHARGED EM FAST") ; + strncpy(particle, "CHARGED EM FAST",100); break ; case AliPHOSFastRecParticle::kCHARGEDHAFAST: - strcpy(particle, "CHARGED HA FAST") ; + strncpy(particle, "CHARGED HA FAST",100); break ; case AliPHOSFastRecParticle::kCHARGEDEMSLOW: - strcpy(particle, "CHARGEDEMSLOW") ; + strncpy(particle, "CHARGEDEMSLOW",100); break ; case AliPHOSFastRecParticle::kCHARGEDHASLOW: - strcpy(particle, "CHARGED HA SLOW") ; + strncpy(particle, "CHARGED HA SLOW",100); break ; } @@ -566,17 +453,8 @@ void AliPHOSPIDv0::PrintRecParticles(Option_t * option) // Int_t nprimaries; // primaries = rp->GetPrimaries(nprimaries); - cout << setw(10) << particle << " " - << setw(5) << rp->GetIndexInList() << " " ; - // << setw(4) << nprimaries << " "; - // for (Int_t iprimary=0; iprimaryGetIndexInList())) ; } - cout << "-------------------------------------------" << endl ; - } - + } } - - -