X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PHOS%2FAliPHOSPIDv0.cxx;h=a393a4f9988cedc9b031733dd014ad773c217abb;hb=dddc33c571b008eccdc8f659de11ae42fc5a5509;hp=ac6adde79b6f16e6b159c94100282f1d31800d44;hpb=21cd0c07a367c6e6836907ecf5222aef9cb05c05;p=u%2Fmrichter%2FAliRoot.git diff --git a/PHOS/AliPHOSPIDv0.cxx b/PHOS/AliPHOSPIDv0.cxx index ac6adde79b6..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,237 +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 --- // --- 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, const Bool_t toSplit) : AliPHOSPID(headerFile, name,toSplit) +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 ; } - AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ; - if(gime->BranchExists("RecParticles") ) - 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 ) { -// Error("Exec", "RecParticles and/or PIDtMaker branch with name %s already exists", taskName.Data() ) ; -// return ; -// } - - Int_t nevents = gime->MaxEvent() ; //(Int_t) gAlice->TreeE()->GetEntries() ; - Int_t ievent ; - - for(ievent = 0; ievent < nevents; ievent++){ - gime->Event(ievent,"R") ; - Info("Exec", "event %d %d %d", ievent, gime->EmcRecPoints(), gime->TrackSegments()) ; - 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"); - Info("Exec", "took %f seconds for PID %f seconds per event", - gBenchmark->GetCpuTime("PHOSPID"), gBenchmark->GetCpuTime("PHOSPID")/nevents) ; + 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(),fToSplit) ; - 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) ; - // 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 ; @@ -298,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() ; @@ -341,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 ; @@ -358,35 +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 TString message ; - message = "=============== AliPHOSPID1 ================\n" ; + message = "=============== AliPHOSPIDv0 ================\n" ; message += "Making PID\n" ; - message += " Headers file: %s\n" ; - message += " RecPoints branch title: %s\n" ; - message += " TrackSegments Branch title: %s\n" ; - message += " RecParticles Branch title %s\n" ; message += "with parameters:\n" ; message += " Maximal EMC - CPV distance (cm) %f\n" ; - Info("Print", message.Data(), - fHeaderFileName.Data(), - fRecPointsTitle.Data(), - fTrackSegmentsTitle.Data(), - fRecParticlesTitle.Data(), - fCpvEmcDistance ); + AliInfo(Form( message.Data(), + fCpvEmcDistance )); if(fIDOptions.Contains("dis",TString::kIgnoreCase )) - Info("Print", " dispersion cut %f", fDispersion ) ; + AliInfo(Form(" dispersion cut %f", fDispersion )) ; if(fIDOptions.Contains("ell",TString::kIgnoreCase )) - Info("Print", " Eliptic cuts function: %s", fFormula->GetTitle() ) ; + AliInfo(Form(" Eliptic cuts function: %s", + fFormula->GetTitle() )) ; if(fIDOptions.Contains("tim",TString::kIgnoreCase )) - Info("Print", " Time Gate used: %f", fTimeGate) ; + 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. @@ -394,79 +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() ) ; - - 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() ; - } - -// //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 = 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 = 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() ; - - treeR->AutoSave() ; //Write(0,kOverwrite) ; - -} //____________________________________________________________________________ void AliPHOSPIDv0::PlotDispersionCuts()const @@ -502,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 @@ -538,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 ; } @@ -551,49 +409,43 @@ void AliPHOSPIDv0::PrintRecParticles(Option_t * option) { // Print table of reconstructed particles - AliPHOSGetter *gime = AliPHOSGetter::GetInstance() ; - - TString taskName(GetName()) ; - taskName.Remove(taskName.Index(Version())-1) ; - TClonesArray * recParticles = gime->RecParticles(taskName) ; - TString message ; - message = "event %d\n" ; - message += " found %d RecParticles\n" ; - Info("PrintRecParticles", message.Data(), gAlice->GetEvNumber(), recParticles->GetEntriesFast() ) ; + message = "Found %d RecParticles\n" ; + AliInfo(Form(message.Data(), + fRecParticles->GetEntriesFast() )) ; if(strstr(option,"all")) { // printing found TS - Info("PrintRecParticles"," PARTICLE Index \n" ) ; - + 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 ; } @@ -601,10 +453,8 @@ void AliPHOSPIDv0::PrintRecParticles(Option_t * option) // Int_t nprimaries; // primaries = rp->GetPrimaries(nprimaries); - Info("PrintRecParticles", " %s %d\n", particle, rp->GetIndexInList()) ; + AliInfo(Form(" %s %d", + particle, rp->GetIndexInList())) ; } } } - - -