X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=blobdiff_plain;f=PHOS%2FAliPHOSAnalyze.cxx;h=100d5119cb3977256e1f7f2aceeb570c85aabf56;hp=5bc8d0188a444aea140885de7818b6dfcfc60360;hb=a675b8d6e4d875046946750ab129103310114951;hpb=3fc07a8086ded4667f992b0731c9a6ac1800e16b diff --git a/PHOS/AliPHOSAnalyze.cxx b/PHOS/AliPHOSAnalyze.cxx index 5bc8d0188a4..100d5119cb3 100644 --- a/PHOS/AliPHOSAnalyze.cxx +++ b/PHOS/AliPHOSAnalyze.cxx @@ -13,66 +13,120 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ +/* $Id$ */ + //_________________________________________________________________________ -// Algorythm class to analyze PHOS events -//*-- Y. Schutz : SUBATECH +// Algorythm class to analyze PHOS events. In this class we demostrate, +// how to handle reconstructed objects with AliPHSOIndexToObject. +// As an example we propose sulotions for four most frequently used tasks: +// DrawRecon(...) - to draw reconstructed objects in the PHOS plane, +// very usefull in the debuging +// InvarianMass(...) - to calculate "REAL" and "MIXED" photon pairs +// invariant mass distributions +// EnergyResoluition(...) -\ Energy and position resolutions of the +// PositionResolution(...)-/ reconstructed photons +// Contamination(...) - calculates contamination of the photon spectrum and +// pobability of reconstruction of several primaries as +// kGAMMA,kELECTRON etc. +//// User Case: +// root [0] AliPHOSAnalyze * a = new AliPHOSAnalyze("galice.root") +// // set the file you want to analyse +// root [1] a->DrawRecon(1,3) +// // plot RecObjects, made in event 1, PHOS module 3 +// root [2] a->DrawRecon(1,3,"PHOSRP","another PID") +// // plot RecObjets made in the event 1, PHOS module 3, +// // produced in the another reconstruction pass, +// // which produced PHOS RecParticles ("PHOSRP") with +// // title "another PID". +// root [3] a->InvariantMass() +// // Calculates "REAL" and "MIXED" invariant mass +// // distributions of kGAMMA and (kGAMMA+kNEUTRALEM) +// // and APPENDS this to the file "invmass.root" +// root [4] a->PositionResolution() +// // calculates two dimentional histos: energy of the primary +// // photon vs distance betwin incedence point and reconstructed +// // poisition. One can analyse the produced file position.root +// // with macro PhotonPosition.C +// root [5] a->EnergyResolution() +// // calculates two dimentional histos: energy of the primary +// // photon vs energy of the reconstructed particle. One can +// // analyse the produced file energy.root +// // with macro PhotonEnergy.C +// root [6] a->Contamination() +// // fills spectra of primary photons and several kinds of +// // reconstructed particles, so that analyzing them one can +// // estimate conatmination, efficiency of registration etc. +//*-- +//*-- Author: Dmitri Peressounko (SUBATECH & RRC Kurchatov Institute) ////////////////////////////////////////////////////////////////////////////// + // --- ROOT system --- #include "TFile.h" #include "TH1.h" -#include "TPad.h" +#include "TH2.h" #include "TH2.h" #include "TParticle.h" +#include "TDatabasePDG.h" #include "TClonesArray.h" -#include "TTree.h" #include "TMath.h" -#include "TCanvas.h" +#include "TROOT.h" // --- Standard library --- -#include -#include - // --- AliRoot header files --- -#include "AliRun.h" +#include "AliLog.h" +#include "AliStack.h" +#include "AliPHOSGeometry.h" #include "AliPHOSAnalyze.h" -#include "AliPHOSClusterizerv1.h" -#include "AliPHOSTrackSegmentMakerv1.h" -#include "AliPHOSPIDv1.h" -#include "AliPHOSReconstructioner.h" #include "AliPHOSDigit.h" +#include "AliPHOSSDigitizer.h" +#include "AliPHOSEmcRecPoint.h" +#include "AliPHOSCpvRecPoint.h" #include "AliPHOSTrackSegment.h" #include "AliPHOSRecParticle.h" +#include "AliPHOSLoader.h" + ClassImp(AliPHOSAnalyze) +//____________________________________________________________________________ +AliPHOSAnalyze::AliPHOSAnalyze(): + fCorrection(1.2), //Value calculated for default parameters of reconstruction + fEvt(0), + ffileName(), + fRunLoader(0) +{ + // default ctor (useless) +} //____________________________________________________________________________ - AliPHOSAnalyze::AliPHOSAnalyze() +AliPHOSAnalyze::AliPHOSAnalyze(Text_t * fileName): + fCorrection(1.05), //Value calculated for default parameters of reconstruction + fEvt(0), + ffileName(fileName), + fRunLoader(0) { - // ctor - - fRootFile = 0 ; + // ctor: analyze events from root file "name" + fRunLoader = AliRunLoader::Open(fileName,"AliPHOSAnalyze"); + if (fRunLoader == 0x0) + { + AliError(Form("Error Loading session")); + } } //____________________________________________________________________________ -AliPHOSAnalyze::AliPHOSAnalyze(Text_t * name) +AliPHOSAnalyze::AliPHOSAnalyze(const AliPHOSAnalyze & ana): + TObject(ana), + fCorrection(0.), + fEvt(0), + ffileName(), + fRunLoader(0) { - // ctor - - Bool_t ok = OpenRootFile(name) ; - if ( !ok ) { - cout << " AliPHOSAnalyze > Error opening " << name << endl ; - } - else { - gAlice = (AliRun*) fRootFile->Get("gAlice"); - fPHOS = (AliPHOSv0 *)gAlice->GetDetector("PHOS") ; - fGeom = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() ) ; - fEvt = -999 ; - } + // copy ctor + ( (AliPHOSAnalyze &)ana ).Copy(*this) ; } //____________________________________________________________________________ @@ -80,760 +134,1175 @@ AliPHOSAnalyze::~AliPHOSAnalyze() { // dtor - fRootFile->Close() ; - delete fRootFile ; - fRootFile = 0 ; +} +//____________________________________________________________________________ +void AliPHOSAnalyze::DrawRecon(Int_t Nevent,Int_t Nmod){ + //Draws pimary particles and reconstructed + //digits, RecPoints, RecPartices etc + //for event Nevent in the module Nmod. + + //========== Create ObjectLoader + if (fRunLoader == 0x0) + { + AliError(Form("Error Loading session")); + return; + } + + AliPHOSLoader* gime = dynamic_cast(fRunLoader->GetLoader("PHOSLoader")); + if ( gime == 0 ) + { + AliError(Form("Could not obtain the Loader object !")); + return ; + } + + + if(Nevent >= fRunLoader->GetNumberOfEvents() ) { + AliError(Form("There is no event %d only %d events available", Nevent, fRunLoader->GetNumberOfEvents() )) ; + return ; + } + AliPHOSGeometry * phosgeom = AliPHOSGeometry::GetInstance() ; + fRunLoader->GetEvent(Nevent); + + Int_t nx = phosgeom->GetNPhi() ; + Int_t nz = phosgeom->GetNZ() ; + const Float_t * cri= phosgeom->GetEMCAGeometry()->GetCrystalHalfSize() ; + Float_t x = nx*cri[0] ; + Float_t z = nz*cri[2] ; + Int_t nxCPV = (Int_t) (nx*phosgeom->GetPadSizePhi()/(2.*cri[0])) ; + Int_t nzCPV = (Int_t) (nz*phosgeom->GetPadSizeZ()/(2.*cri[2])) ; + + TH2F * emcDigits = (TH2F*) gROOT->FindObject("emcDigits") ; + if(emcDigits) + emcDigits->Delete() ; + emcDigits = new TH2F("emcDigits","EMC digits", nx,-x,x,nz,-z,z); + TH2F * emcSdigits =(TH2F*) gROOT->FindObject("emcSdigits") ; + if(emcSdigits) + emcSdigits->Delete() ; + emcSdigits = new TH2F("emcSdigits","EMC sdigits", nx,-x,x,nz,-z,z); + TH2F * emcRecPoints = (TH2F*)gROOT->FindObject("emcRecPoints") ; + if(emcRecPoints) + emcRecPoints->Delete() ; + emcRecPoints = new TH2F("emcRecPoints","EMC RecPoints",nx,-x,x,nz,-z,z); + TH2F * cpvSdigits =(TH2F*) gROOT->FindObject("cpvSdigits") ; + if(cpvSdigits) + cpvSdigits->Delete() ; + cpvSdigits = new TH2F("cpvSdigits","CPV sdigits", nx,-x,x,nz,-z,z); + TH2F * cpvDigits = (TH2F*)gROOT->FindObject("cpvDigits") ; + if(cpvDigits) + cpvDigits->Delete() ; + cpvDigits = new TH2F("cpvDigits","CPV digits", nxCPV,-x,x,nzCPV,-z,z) ; + TH2F * cpvRecPoints= (TH2F*)gROOT->FindObject("cpvRecPoints") ; + if(cpvRecPoints) + cpvRecPoints->Delete() ; + cpvRecPoints = new TH2F("cpvRecPoints","CPV RecPoints", nxCPV,-x,x,nzCPV,-z,z) ; + + TH2F * phot = (TH2F*)gROOT->FindObject("phot") ; + if(phot) + phot->Delete() ; + phot = new TH2F("phot","Primary Photon", nx,-x,x,nz,-z,z); + TH2F * recPhot = (TH2F*)gROOT->FindObject("recPhot") ; + if(recPhot) + recPhot->Delete() ; + recPhot = new TH2F("recPhot","RecParticles with primary Photon",nx,-x,x,nz,-z,z); + + //Get Vertex + Double_t vtx[3]={0.,0.,0.} ; +//DP: extract vertex either from Generator or from data - delete fPHOS ; - fPHOS = 0 ; + + //Plot Primary Particles + + if (fRunLoader->Stack() == 0x0) fRunLoader->LoadKinematics("READ"); + - delete fClu ; - fClu = 0 ; + const TParticle * primary ; + Int_t iPrimary ; + for ( iPrimary = 0 ; iPrimary < fRunLoader->Stack()->GetNprimary() ; iPrimary++) + { + primary = fRunLoader->Stack()->Particle(iPrimary); + + Int_t primaryType = primary->GetPdgCode(); +// if( (primaryType == 211)||(primaryType == -211)||(primaryType == 2212)||(primaryType == -2212) +// ||(primaryType == 11)||(primaryType == -11) ) { +// Int_t moduleNumber ; +// Double_t primX, primZ ; +// phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ; +// if(moduleNumber==Nmod) +// charg->Fill(primZ,primX,primary->Energy()) ; +// } + if( primaryType == 22 ) { + Int_t moduleNumber ; + Double_t primX, primZ ; + phosgeom->ImpactOnEmc(vtx,primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ; + if(moduleNumber==Nmod) + phot->Fill(primZ,primX,primary->Energy()) ; + } +// else{ +// if( primaryType == -2112 ) { +// Int_t moduleNumber ; +// Double_t primX, primZ ; +// phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ; +// if(moduleNumber==Nmod) +// nbar->Fill(primZ,primX,primary->Energy()) ; +// } +// } + } - delete fPID ; - fPID = 0 ; + + Int_t iSDigit ; + AliPHOSDigit * sdigit ; + const TClonesArray * sdigits = gime->SDigits() ; + Int_t nsdig[5] = {0,0,0,0,0} ; + if(sdigits){ + for(iSDigit = 0; iSDigit < sdigits->GetEntriesFast() ; iSDigit++) + { + sdigit = (AliPHOSDigit *) sdigits->At(iSDigit) ; + Int_t relid[4]; + phosgeom->AbsToRelNumbering(sdigit->GetId(), relid) ; + Float_t xd,zd ; + phosgeom->RelPosInModule(relid,xd,zd); + Float_t e = sdigit->GetEnergy() ; + nsdig[relid[0]-1]++ ; + if(relid[0]==Nmod){ + if(relid[1]==0) //EMC + emcSdigits->Fill(xd,zd,e) ; + if( relid[1]!=0 ) + cpvSdigits->Fill(xd,zd,e) ; + } + } + } + TString message ; + message = "Number of EMC + CPV SDigits per module: \n" ; + message += "%d %d %d %d %d\n"; + AliInfo(Form(message.Data(), nsdig[0], nsdig[1], nsdig[2], nsdig[3], nsdig[4] )) ; + + //Plot digits + Int_t iDigit ; + AliPHOSDigit * digit ; + const TClonesArray * digits = gime->Digits(); + if(digits) { + for(iDigit = 0; iDigit < digits->GetEntriesFast(); iDigit++) + { + digit = (AliPHOSDigit *) digits->At(iDigit) ; + Int_t relid[4]; + phosgeom->AbsToRelNumbering(digit->GetId(), relid) ; + Float_t xd,zd ; + phosgeom->RelPosInModule(relid,xd,zd) ; + Float_t e = digit->GetEnergy() ; + if(relid[0]==Nmod){ + if(relid[1]==0) //EMC + emcDigits->Fill(xd,zd,e) ; + if( relid[1]!=0 ) + cpvDigits->Fill(xd,zd,e) ; + } + } + } + + + //Plot RecPoints + Int_t irecp ; + TVector3 pos ; + TObjArray * emcrp = gime->EmcRecPoints() ; + if(emcrp) { + for(irecp = 0; irecp < emcrp->GetEntriesFast() ; irecp ++){ + AliPHOSEmcRecPoint * emc = (AliPHOSEmcRecPoint *) emcrp->At(irecp) ; + if(emc->GetPHOSMod()==Nmod){ + emc->GetLocalPosition(pos) ; + emcRecPoints->Fill(pos.X(),pos.Z(),emc->GetEnergy()); + } + } + } + + TObjArray * cpvrp = gime->CpvRecPoints() ; + if(cpvrp) { + for(irecp = 0; irecp < cpvrp->GetEntriesFast() ; irecp ++){ + AliPHOSRecPoint * cpv = (AliPHOSCpvRecPoint *) cpvrp->At(irecp) ; + if(cpv->GetPHOSMod()==Nmod){ + cpv->GetLocalPosition(pos) ; + cpvRecPoints->Fill(pos.X(),pos.Z(),cpv->GetEnergy()); + } + } + } + + //Plot RecParticles + AliPHOSRecParticle * recParticle ; + Int_t iRecParticle ; + TClonesArray * rp = gime->RecParticles() ; + TClonesArray * ts = gime->TrackSegments() ; + if(rp && ts && emcrp) { + for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast() ; iRecParticle++ ) + { + recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ; + Int_t moduleNumberRec ; + Double_t recX, recZ ; + phosgeom->ImpactOnEmc(vtx,recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ; + if(moduleNumberRec == Nmod){ + + Double_t minDistance = 5. ; + Int_t closestPrimary = -1 ; + + //extract list of primaries: it is stored at EMC RecPoints + Int_t emcIndex = ((AliPHOSTrackSegment *) ts->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ; + Int_t numberofprimaries ; + Int_t * listofprimaries = ((AliPHOSRecPoint*) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries) ; + Int_t index ; + const TParticle * primPart ; + Double_t distance = minDistance ; + + for ( index = 0 ; index < numberofprimaries ; index++){ + primPart = fRunLoader->Stack()->Particle(listofprimaries[index]) ; + Int_t moduleNumber ; + Double_t primX, primZ ; + phosgeom->ImpactOnEmc(vtx,primPart->Theta(), primPart->Phi(), moduleNumber, primX, primZ) ; + if(moduleNumberRec == moduleNumber) + distance = TMath::Sqrt((recX-primX)*(recX-primX)+(recZ-primZ)*(recZ-primZ) ) ; + if(minDistance > distance) + { + minDistance = distance ; + closestPrimary = listofprimaries[index] ; + } + } + + if(closestPrimary >=0 ){ + + Int_t primaryType = fRunLoader->Stack()->Particle(closestPrimary)->GetPdgCode() ; + + if(primaryType==22) + recPhot->Fill(recZ,recX,recParticle->Energy()) ; +// else +// if(primaryType==-2112) +// recNbar->Fill(recZ,recX,recParticle->Energy()) ; + } + } + } - delete fRec ; - fRec = 0 ; + } + + //Plot made histograms + emcSdigits->Draw("box") ; + emcDigits->SetLineColor(5) ; + emcDigits->Draw("boxsame") ; + emcRecPoints->SetLineColor(2) ; + emcRecPoints->Draw("boxsame") ; + cpvSdigits->SetLineColor(1) ; + cpvSdigits->Draw("boxsame") ; + +} +//____________________________________________________________________________ +void AliPHOSAnalyze::Ls(){ + //lists branches and titles of PHOS-related branches of TreeR, TreeD, TreeS + + if (fRunLoader == 0x0) + { + AliError(Form("Error Loading session")); + return; + } + + AliPHOSLoader* gime = dynamic_cast(fRunLoader->GetLoader("PHOSLoader")); + if ( gime == 0 ) + { + AliError(Form("Could not obtain the Loader object !")); + return ; + } - delete fTrs ; - fTrs = 0 ; -} + Int_t ibranch; + TObjArray * branches; + + if (gime->TreeS() == 0x0) + { + if (gime->LoadSDigits("READ")) + { + AliError(Form("Problems with loading summable digits")); + return; + } + } + branches = gime->TreeS()->GetListOfBranches() ; + + TString message ; + message = "TreeS:\n" ; + for(ibranch = 0;ibranch GetEntries();ibranch++){ + TBranch * branch=(TBranch *) branches->At(ibranch) ; + if(strstr(branch->GetName(),"PHOS") ){ + message += " " ; + message += branch->GetName() ; + message += " " ; + message += branch->GetTitle() ; + message += "\n" ; + } + } + if (gime->TreeD() == 0x0) + { + if (gime->LoadDigits("READ")) + { + AliError(Form("Problems with loading digits")); + return; + } + } + + branches = gime->TreeD()->GetListOfBranches() ; + + message += "TreeD:\n" ; + for(ibranch = 0;ibranch GetEntries();ibranch++){ + TBranch * branch=(TBranch *) branches->At(ibranch) ; + if(strstr(branch->GetName(),"PHOS") ) { + message += " "; + message += branch->GetName() ; + message += " " ; + message += branch->GetTitle() ; + message +="\n" ; + } + } + + if (gime->TreeR() == 0x0) + { + if (gime->LoadRecPoints("READ")) + { + AliError(Form("Problems with loading rec points")); + return; + } + } + branches = gime->TreeR()->GetListOfBranches() ; + + message += "TreeR: \n" ; + for(ibranch = 0;ibranch GetEntries();ibranch++){ + TBranch * branch=(TBranch *) branches->At(ibranch) ; + if(strstr(branch->GetName(),"PHOS") ) { + message += " " ; + message += branch->GetName() ; + message += " " ; + message += branch->GetTitle() ; + message += "\n" ; + } + } + AliInfo(Form(message.Data())) ; +} //____________________________________________________________________________ -void AliPHOSAnalyze::AnalyzeOneEvent(Int_t evt) + void AliPHOSAnalyze::InvariantMass() { - Bool_t ok = Init(evt) ; + // Calculates Real and Mixed invariant mass distributions + if (fRunLoader == 0x0) + { + AliError(Form("Error Loading session")); + return; + } - if ( ok ) { - //=========== Get the number of entries in the Digits array - - Int_t nId = fPHOS->Digits()->GetEntries(); - printf("AnalyzeOneEvent > Number of entries in the Digit array is %d \n",nId); - - //=========== Do the reconstruction - - cout << "AnalyzeOneEvent > Found " << nId << " digits in PHOS" << endl ; + AliPHOSLoader* gime = dynamic_cast(fRunLoader->GetLoader("PHOSLoader")); + if ( gime == 0 ) + { + AliError(Form("Could not obtain the Loader object !")); + return ; + } + + gime->LoadRecParticles("READ"); + + Int_t nMixedEvents = 4 ; //# of events used for calculation of 'mixed' distribution + + + //opening file + TFile * mfile = new TFile("invmass.root","update"); + + //========== Reading /Booking Histograms + TH2D * hRealEM = 0 ; + hRealEM = (TH2D*) mfile->Get("hRealEM") ; + if(hRealEM == 0) + hRealEM = new TH2D("hRealEM", "Real for EM particles", 250,0.,1.,40,0.,4.) ; + TH2D * hRealPhot = 0 ; + + hRealPhot = (TH2D*)mfile->Get("hRealPhot"); + if(hRealPhot == 0) + hRealPhot = new TH2D("hRealPhot", "Real for kPhoton particles", 250,0.,1.,40,0.,4.) ; + + TH2D * hMixedEM = 0 ; + hMixedEM = (TH2D*) mfile->Get("hMixedEM") ; + if(hMixedEM == 0) + hMixedEM = new TH2D("hMixedEM", "Mixed for EM particles", 250,0.,1.,40,0.,4.) ; + + TH2D * hMixedPhot = 0 ; + hMixedPhot = (TH2D*) mfile->Get("hMixedPhot") ; + if(hMixedPhot == 0) + hMixedPhot = new TH2D("hMixedPhot","Mixed for kPhoton particles",250,0.,1.,40,0.,4.) ; + + + //reading event and copyng it to TConesArray of all photons + + TClonesArray * allRecParticleList = new TClonesArray("AliPHOSRecParticle", 1000) ; + Int_t * nRecParticles = new Int_t[nMixedEvents] ; // to mark boundaries of each event in the total list + for(Int_t index = 0; index < nMixedEvents; index ++) + nRecParticles[index] = 0 ; + Int_t iRecPhot = 0 ; // number of EM particles in total list + + //scan over all events + Int_t event ; + + + if (fRunLoader->TreeE() == 0x0) fRunLoader->LoadHeader(); + + + Int_t maxevent = (Int_t)fRunLoader->TreeE()->GetEntries(); + // for(event = 0; event < gime->MaxEvent(); event++ ){ + + + + for(event = 0; event < maxevent; event++ ){ + fRunLoader->GetEvent(event); //will read only TreeR - fPHOS->Reconstruction(fRec); + //copy EM RecParticles to the "total" list + const AliPHOSRecParticle * recParticle ; + Int_t iRecParticle ; + TClonesArray * rp = gime->RecParticles() ; + if(!rp){ + AliError(Form("Can't find RecParticles")) ; + return ; + } + + for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast(); iRecParticle++ ) + { + recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ; + if((recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST)|| + (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW)) + new( (*allRecParticleList)[iRecPhot++] ) AliPHOSRecParticle(*recParticle) ; + } - // =========== End of reconstruction + Int_t mevent = event%nMixedEvents ; //event number in the "mixed" cicle + nRecParticles[mevent] = iRecPhot-1 ; - cout << "AnalyzeOneEvent > event # " << fEvt << " processed" << endl ; - } // ok - else - cout << "AnalyzeOneEvent > filed to process event # " << evt << endl ; + //check, if it is time to calculate invariant mass? + if((mevent == 0) && (event +1 == maxevent)){ + + // if((mevent == 0) && (event +1 == gime->MaxEvent())){ + + //calculate invariant mass: + Int_t irp1,irp2 ; + Int_t nCurEvent = 0 ; + + for(irp1 = 0; irp1 < allRecParticleList->GetEntries()-1; irp1++){ + AliPHOSRecParticle * rp1 = (AliPHOSRecParticle *)allRecParticleList->At(irp1) ; + + for(irp2 = irp1+1; irp2 < allRecParticleList->GetEntries(); irp2++){ + AliPHOSRecParticle * rp2 = (AliPHOSRecParticle *)allRecParticleList->At(irp2) ; + + Double_t invMass ; + invMass = (rp1->Energy()+rp2->Energy())*(rp1->Energy()+rp2->Energy())- + (rp1->Px()+rp2->Px())*(rp1->Px()+rp2->Px())- + (rp1->Py()+rp2->Py())*(rp1->Py()+rp2->Py())- + (rp1->Pz()+rp2->Pz())*(rp1->Pz()+rp2->Pz()) ; + + if(invMass> 0) + invMass = TMath::Sqrt(invMass); + + Double_t pt ; + pt = TMath::Sqrt((rp1->Px()+rp2->Px() )*( rp1->Px()+rp2->Px() ) + + (rp1->Py()+rp2->Py() )*( rp1->Py()+rp2->Py() ) ); + + if(irp1 > nRecParticles[nCurEvent]) + nCurEvent++; + + if(irp2 <= nRecParticles[nCurEvent]){ //'Real' event + hRealEM->Fill(invMass,pt); + if((rp1->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST)&& + (rp2->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST) ) + hRealPhot->Fill(invMass,pt); + } + else{ + hMixedEM->Fill(invMass,pt); + if((rp1->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST)&& + (rp2->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST) ) + hMixedPhot->Fill(invMass,pt); + } //real-mixed + + } //loop over second rp + }//loop over first rp + + //Make some cleanings + for(Int_t index = 0; index < nMixedEvents; index ++) + nRecParticles[index] = 0 ; + iRecPhot = 0 ; + allRecParticleList->Clear() ; + + } + } + delete allRecParticleList ; + + //writing output + mfile->cd(); + + hRealEM->Write(0,kOverwrite) ; + hRealPhot->Write(0,kOverwrite) ; + hMixedEM->Write(0,kOverwrite) ; + hMixedPhot->Write(0,kOverwrite) ; + + mfile->Write(); + mfile->Close(); + delete mfile ; + delete nRecParticles; } //____________________________________________________________________________ - void AliPHOSAnalyze::AnalyzeManyEvents(Int_t Nevents, Int_t module) // analyzes many events + void AliPHOSAnalyze::EnergyResolution() { + //fills two dimentional histo: energy of primary vs. energy of reconstructed - if ( fRootFile == 0 ) - cout << "AnalyzeManyEvents > " << "Root File not openned" << endl ; - else - { - //========== Get AliRun object from file - gAlice = (AliRun*) fRootFile->Get("gAlice") ; - //=========== Get the PHOS object and associated geometry from the file - fPHOS = (AliPHOSv0 *)gAlice->GetDetector("PHOS") ; - fGeom = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() ); - //========== Booking Histograms - cout << "AnalyzeManyEvents > " << "Booking Histograms" << endl ; - BookingHistograms(); - Int_t ievent; - Int_t relid[4] ; - AliPHOSDigit * digit ; - AliPHOSEmcRecPoint * emc ; - AliPHOSPpsdRecPoint * ppsd ; - // AliPHOSTrackSegment * tracksegment ; - AliPHOSRecParticle * recparticle; - for ( ievent=0; ievent " << "Starting Analyzing " << endl ; - //========== Create the Clusterizer - fClu = new AliPHOSClusterizerv1() ; - fClu->SetEmcEnergyThreshold(0.025) ; - fClu->SetEmcClusteringThreshold(0.50) ; - fClu->SetPpsdEnergyThreshold (0.0000002) ; - fClu->SetPpsdClusteringThreshold(0.0000001) ; - fClu->SetLocalMaxCut(0.03) ; - fClu->SetCalibrationParameters(0., 0.00000001) ; - //========== Creates the track segment maker - fTrs = new AliPHOSTrackSegmentMakerv1() ; - fTrs->UnsetUnfoldFlag() ; - //========== Creates the particle identifier - fPID = new AliPHOSPIDv1() ; - fPID->SetShowerProfileCuts(0.3, 1.8, 0.3, 1.8 ) ; - fPID->Print() ; - //========== Creates the Reconstructioner - fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ; - //========== Event Number - if ( ( log10(ievent+1) - (Int_t)(log10(ievent+1)) ) == 0. ) cout << "AnalyzeManyEvents > " << "Event is " << ievent << endl ; - //=========== Connects the various Tree's for evt - gAlice->GetEvent(ievent); - //=========== Gets the Digit TTree - gAlice->TreeD()->GetEvent(0) ; - //=========== Gets the number of entries in the Digits array - TIter nextdigit(fPHOS->Digits()) ; - while( ( digit = (AliPHOSDigit *)nextdigit() ) ) - { - fGeom->AbsToRelNumbering(digit->GetId(), relid) ; - if (fClu->IsInEmc(digit)) fhEmcDigit->Fill(fClu->Calibrate(digit->GetAmp())) ; - else - { - if (relid[1]<17) fhVetoDigit->Fill(fClu->Calibrate(digit->GetAmp())); - if (relid[1]>16) fhConvertorDigit->Fill(fClu->Calibrate(digit->GetAmp())); - } - } - //=========== Do the reconstruction - fPHOS->Reconstruction(fRec); - //=========== Cluster in module - TIter nextEmc(fPHOS->EmcClusters() ) ; - while((emc = (AliPHOSEmcRecPoint *)nextEmc())) - { - if ( emc->GetPHOSMod() == module ) - { - fhEmcCluster->Fill( emc->GetTotalEnergy() ); - TIter nextPpsd( fPHOS->PpsdClusters()) ; - while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsd())) - { - if ( ppsd->GetPHOSMod() == module ) - { - if (!ppsd->GetUp()) fhConvertorEmc->Fill(emc->GetTotalEnergy(),ppsd->GetTotalEnergy()) ; - } - } - } - } - //=========== Cluster in module PPSD Down - TIter nextPpsd(fPHOS->PpsdClusters() ) ; - while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsd())) - { - if ( ppsd->GetPHOSMod() == module ) - { - if (!ppsd->GetUp()) fhConvertorCluster->Fill(ppsd->GetTotalEnergy()) ; - if (ppsd->GetUp()) fhVetoCluster ->Fill(ppsd->GetTotalEnergy()) ; - } - } - //========== TRackSegments in the event - TIter nextRecParticle(fPHOS->RecParticles() ) ; - while((recparticle = (AliPHOSRecParticle *)nextRecParticle())) - { - if ( recparticle->GetPHOSTrackSegment()->GetPHOSMod() == module ) - { - cout << "Particle type is " << recparticle->GetType() << endl ; - switch(recparticle->GetType()) - { - case kGAMMA: - fhPhotonEnergy->Fill(recparticle->Energy() ) ; - //fhPhotonPositionX->Fill(recpart. ) ; - //fhPhotonPositionY->Fill(recpart. ) ; - cout << "PHOTON" << endl; - break; - case kELECTRON: - fhElectronEnergy->Fill(recparticle->Energy() ) ; - //fhElectronPositionX->Fill(recpart. ) ; - //fhElectronPositionY->Fill(recpart. ) ; - cout << "ELECTRON" << endl; - break; - case kNEUTRALHADRON: - fhNeutralHadronEnergy->Fill(recparticle->Energy() ) ; - //fhNeutralHadronPositionX->Fill(recpart. ) ; - //fhNeutralHadronPositionY->Fill(recpart. ) ; - cout << "NEUTRAl HADRON" << endl; - break ; - case kNEUTRALEM: - fhNeutralEMEnergy->Fill(recparticle->Energy() ) ; - //fhNeutralEMPositionX->Fill(recpart. ) ; - //fhNeutralEMPositionY->Fill(recpart. ) ; - //cout << "NEUTRAL EM" << endl; - break ; - case kCHARGEDHADRON: - fhChargedHadronEnergy->Fill(recparticle->Energy() ) ; - //fhChargedHadronPositionX->Fill(recpart. ) ; - //fhChargedHadronPositionY->Fill(recpart. ) ; - cout << "CHARGED HADRON" << endl; - break ; - case kGAMMAHADRON: - fhPhotonHadronEnergy->Fill(recparticle->Energy() ) ; - //fhPhotonHadronPositionX->Fill(recpart. ) ; - //fhPhotonHadronPositionY->Fill(recpart. ) ; - cout << "PHOTON HADRON" << endl; - break ; - } - } - } - // Deleting fClu, fTrs, fPID et fRec - fClu->Delete(); - fTrs->Delete(); - fPID->Delete(); - fRec->Delete(); - - } // endfor - SavingHistograms(); - } // endif -} // endfunction + TH2F * hAllEnergy = 0 ; //all reconstructed with primary photon + TH2F * hPhotEnergy= 0 ; //kGamma with primary photon + TH2F * hEMEnergy = 0 ; //electromagnetic with primary photon + //opening file and reading histograms if any + TFile * efile = new TFile("energy.root","update"); -//____________________________________________________________________________ -void AliPHOSAnalyze::BookingHistograms() -{ - if (fhEmcDigit ) - delete fhEmcDigit ; - if (fhVetoDigit ) - delete fhVetoDigit ; - if (fhConvertorDigit ) - delete fhConvertorDigit ; - if (fhEmcCluster ) - delete fhEmcCluster ; - if (fhVetoCluster ) - delete fhVetoCluster ; - if (fhConvertorCluster ) - delete fhConvertorCluster ; - if (fhConvertorEmc ) - delete fhConvertorEmc ; - - fhEmcDigit = new TH1F("hEmcDigit", "hEmcDigit", 1000, 0. , 25.); - fhVetoDigit = new TH1F("hVetoDigit", "hVetoDigit", 500, 0. , 3.e-5); - fhConvertorDigit = new TH1F("hConvertorDigit","hConvertorDigit", 500, 0. , 3.e-5); - fhEmcCluster = new TH1F("hEmcCluster", "hEmcCluster", 1000, 0. , 30.); - fhVetoCluster = new TH1F("hVetoCluster", "hVetoCluster", 500, 0. , 3.e-5); - fhConvertorCluster = new TH1F("hConvertorCluster","hConvertorCluster",500, 0. , 3.e-5); - fhConvertorEmc = new TH2F("hConvertorEmc", "hConvertorEmc", 200, 1. , 3., 200, 0., 3.e-5); - fhPhotonEnergy = new TH1F("hPhotonEnergy", "hPhotonEnergy", 1000, 0. , 30.); - fhElectronEnergy = new TH1F("hElectronEnergy","hElectronEnergy", 1000, 0. , 30.); - fhNeutralHadronEnergy = new TH1F("hNeutralHadronEnergy", "hNeutralHadronEnergy", 1000, 0. , 30.); - fhNeutralEMEnergy = new TH1F("hNeutralEMEnergy", "hNeutralEMEnergy", 1000, 0. , 30.); - fhChargedHadronEnergy = new TH1F("hChargedHadronEnergy", "hChargedHadronEnergy", 1000, 0. , 30.); - fhPhotonHadronEnergy = new TH1F("hPhotonHadronEnergy","hPhotonHadronEnergy",500,-80. , 80.); - fhPhotonPositionX = new TH1F("hPhotonPositionX","hPhotonPositionX", 500,-80. , 80.); - fhElectronPositionX = new TH1F("hElectronPositionX","hElectronPositionX",500,-80. , 80.); - fhNeutralHadronPositionX = new TH1F("hNeutralHadronPositionX","hNeutralHadronPositionX",500,-80. , 80.); - fhNeutralEMPositionX = new TH1F("hNeutralEMPositionX","hNeutralEMPositionX",500,-80. , 80.); - fhChargedHadronPositionX = new TH1F("hChargedHadronPositionX","hChargedHadronPositionX",500,-80. , 80.); - fhPhotonHadronPositionX = new TH1F("hPhotonHadronPositionX","hPhotonHadronPositionX",500,-80. , 80.); - fhPhotonPositionY = new TH1F("hPhotonPositionY","hPhotonPositionY", 500,-80. , 80.); - fhElectronPositionY = new TH1F("hElectronPositionY","hElectronPositionY",500,-80. , 80.); - fhNeutralHadronPositionY = new TH1F("hNeutralHadronPositionY","hNeutralHadronPositionY",500,-80. , 80.); - fhNeutralEMPositionY = new TH1F("hNeutralEMPositionY","hNeutralEMPositionY",500,-80. , 80.); - fhChargedHadronPositionY = new TH1F("hChargedHadronPositionY","hChargedHadronPositionY",500,-80. , 80.); - fhPhotonHadronPositionY = new TH1F("hPhotonHadronPositionY","hPhotonHadronPositionY",500,-80. , 80.); + hAllEnergy = (TH2F*)efile->Get("hAllEnergy") ; + if(hAllEnergy == 0) + hAllEnergy = new TH2F("hAllEnergy", "Energy of any RP with primary photon",100, 0., 5., 100, 0., 5.); + hPhotEnergy =(TH2F*) efile->Get("hPhotEnergy") ; + if(hPhotEnergy == 0) + hPhotEnergy = new TH2F("hPhotEnergy", "Energy of kGAMMA with primary photon",100, 0., 5., 100, 0., 5.); -} -//____________________________________________________________________________ -Bool_t AliPHOSAnalyze::Init(Int_t evt) -{ + hEMEnergy =(TH2F*) efile->Get("hEMEnergy"); + if(hEMEnergy == 0) + hEMEnergy = new TH2F("hEMEnergy", "Energy of EM with primary photon", 100, 0., 5., 100, 0., 5.); - Bool_t ok = kTRUE ; + + if (fRunLoader == 0x0) + { + AliError(Form("Error Loading session")); + return; + } + + AliPHOSLoader* gime = dynamic_cast(fRunLoader->GetLoader("PHOSLoader")); + if ( gime == 0 ) + { + AliError(Form("Could not obtain the Loader object !")); + return ; + } + + + AliPHOSGeometry * phosgeom = AliPHOSGeometry::GetInstance() ; + + Int_t ievent; + Int_t maxevent = (Int_t)fRunLoader->TreeE()->GetEntries(); + + fRunLoader->LoadKinematics("READ"); + gime->LoadTracks("READ"); - //========== Open galice root file + for ( ievent=0; ievent < maxevent ; ievent++){ - if ( fRootFile == 0 ) { - Text_t * name = new Text_t[80] ; - cout << "AnalyzeOneEvent > Enter file root file name : " ; - cin >> name ; - Bool_t ok = OpenRootFile(name) ; - if ( !ok ) - cout << " AliPHOSAnalyze > Error opening " << name << endl ; - else { - //========== Get AliRun object from file + //read the current event + fRunLoader->GetEvent(ievent) ; + + Double_t vtx[3]={0.,0.,0.} ; + + const AliPHOSRecParticle * recParticle ; + Int_t iRecParticle ; + TClonesArray * rp = gime->RecParticles() ; + if(!rp) { + AliError(Form("Event %d, Can't find RecParticles ", ievent)) ; + return ; + } + TClonesArray * ts = gime->TrackSegments() ; + if(!ts) { + AliError(Form("Event %d, Can't find TrackSegments", ievent)) ; + return ; + } + TObjArray * emcrp = gime->EmcRecPoints() ; + if(!emcrp){ + AliError(Form("Event %d, Can't find EmcRecPoints", ievent)) ; + return ; + } - gAlice = (AliRun*) fRootFile->Get("gAlice") ; + for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast() ;iRecParticle++ ){ + recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ; - //=========== Get the PHOS object and associated geometry from the file + //find the closest primary + Int_t moduleNumberRec ; + Double_t recX, recZ ; + phosgeom->ImpactOnEmc(vtx,recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ; - fPHOS = (AliPHOSv0 *)gAlice->GetDetector("PHOS") ; - fGeom = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() ); - } // else !ok - } // if fRootFile - - if ( ok ) { - - //========== Create the Clusterizer - - fClu = new AliPHOSClusterizerv1() ; - fClu->SetEmcEnergyThreshold(0.025) ; - fClu->SetEmcClusteringThreshold(0.50) ; - fClu->SetPpsdEnergyThreshold (0.0000002) ; - fClu->SetPpsdClusteringThreshold(0.0000001) ; - fClu->SetLocalMaxCut(0.03) ; - fClu->SetCalibrationParameters(0., 0.00000001) ; - cout << "AnalyzeOneEvent > using clusterizer " << fClu->GetName() << endl ; - fClu->PrintParameters() ; - - //========== Creates the track segment maker - - fTrs = new AliPHOSTrackSegmentMakerv1() ; - cout << "AnalyzeOneEvent > using tack segment maker " << fTrs->GetName() << endl ; - - //========== Creates the particle identifier - - fPID = new AliPHOSPIDv1() ; - cout << "AnalyzeOneEvent > using particle identifier " << fPID->GetName() << endl ; - - //========== Creates the Reconstructioner - - fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ; - - //=========== Connect the various Tree's for evt - - if ( evt == -999 ) { - cout << "AnalyzeOneEvent > Enter event number : " ; - cin >> evt ; - cout << evt << endl ; - } - fEvt = evt ; - gAlice->GetEvent(evt); - - //=========== Get the Digit TTree - - gAlice->TreeD()->GetEvent(0) ; - - } // ok - - return ok ; -} - + Double_t minDistance = 100. ; + Int_t closestPrimary = -1 ; + + //extract list of primaries: it is stored at EMC RecPoints + Int_t emcIndex = ((AliPHOSTrackSegment*) ts->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ; + Int_t numberofprimaries ; + Int_t * listofprimaries = ((AliPHOSEmcRecPoint*) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries) ; + + Int_t index ; + const TParticle * primary ; + Double_t distance = minDistance ; + Double_t dX, dZ; + Double_t dXmin = 0.; + Double_t dZmin = 0. ; + for ( index = 0 ; index < numberofprimaries ; index++){ + + primary = fRunLoader->Stack()->Particle(listofprimaries[index]) ; + + Int_t moduleNumber ; + Double_t primX, primZ ; + phosgeom->ImpactOnEmc(vtx,primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ; + if(moduleNumberRec == moduleNumber) { + dX = recX - primX; + dZ = recZ - primZ; + distance = TMath::Sqrt(dX*dX + dZ*dZ) ; + if(minDistance > distance) { + minDistance = distance ; + dXmin = dX; + dZmin = dZ; + closestPrimary = listofprimaries[index] ; + } + } + } -//____________________________________________________________________________ -void AliPHOSAnalyze::DisplayKineEvent(Int_t evt) -{ - if (evt == -999) - evt = fEvt ; - - Int_t module ; - cout << "DisplayKineEvent > which module (1-5, -1: all) ? " ; - cin >> module ; cout << module << endl ; - - Int_t testparticle ; - cout << " 22 : PHOTON " << endl - << " (-)11 : (POSITRON)ELECTRON " << endl - << " (-)2112 : (ANTI)NEUTRON " << endl - << " -999 : Everything else " << endl ; - cout << "DisplayKineEvent > enter PDG particle code to display " ; - cin >> testparticle ; cout << testparticle << endl ; - - Text_t histoname[80] ; - sprintf(histoname,"Event %d: Incident particles in module %d", evt, module) ; - - Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module - fGeom->EmcModuleCoverage(module, tm, tM, pm, pM, kDegre) ; - - Double_t theta, phi ; - fGeom->EmcXtalCoverage(theta, phi, kDegre) ; - - Int_t tdim = (Int_t)( (tM - tm) / theta ) ; - Int_t pdim = (Int_t)( (pM - pm) / phi ) ; - - tm -= theta ; - tM += theta ; - pm -= phi ; - pM += phi ; - - TH2F * histoparticle = new TH2F("histoparticle", histoname, - pdim, pm, pM, tdim, tm, tM) ; - histoparticle->SetStats(kFALSE) ; - - // Get pointers to Alice Particle TClonesArray - - TParticle * particle; - TClonesArray * particlearray = gAlice->Particles(); - - Text_t canvasname[80]; - sprintf(canvasname,"Particles incident in PHOS/EMC module # %d",module) ; - TCanvas * kinecanvas = new TCanvas("kinecanvas", canvasname, 650, 500) ; - - // get the KINE Tree - - TTree * kine = gAlice->TreeK() ; - Stat_t nParticles = kine->GetEntries() ; - cout << "DisplayKineEvent > events in kine " << nParticles << endl ; - - // loop over particles - - Double_t kRADDEG = 180. / TMath::Pi() ; - Int_t index1 ; - Int_t nparticlein = 0 ; - for (index1 = 0 ; index1 < nParticles ; index1++){ - Int_t nparticle = particlearray->GetEntriesFast() ; - Int_t index2 ; - for( index2 = 0 ; index2 < nparticle ; index2++) { - particle = (TParticle*)particlearray->UncheckedAt(index2) ; - Int_t particletype = particle->GetPdgCode() ; - if (testparticle == -999 || testparticle == particletype) { - Double_t phi = particle->Phi() ; - Double_t theta = particle->Theta() ; - Int_t mod ; - Double_t x, z ; - fGeom->ImpactOnEmc(theta, phi, mod, z, x) ; - if ( mod == module ) { - nparticlein++ ; - histoparticle->Fill(phi*kRADDEG, theta*kRADDEG, particle->Energy() ) ; - } - } + //if found primary, fill histograms + if(closestPrimary >=0 ){ + primary = fRunLoader->Stack()->Particle(closestPrimary) ; + if(primary->GetPdgCode() == 22){ + hAllEnergy->Fill(primary->Energy(), recParticle->Energy()) ; + if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST){ + hPhotEnergy->Fill(primary->Energy(), recParticle->Energy() ) ; + hEMEnergy->Fill(primary->Energy(), recParticle->Energy() ) ; + } + else + if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW) + hEMEnergy->Fill(primary->Energy(), recParticle->Energy() ) ; + } + } } } - kinecanvas->Draw() ; - histoparticle->Draw("color") ; - TPaveText * pavetext = new TPaveText(294, 100, 300, 101); - Text_t text[40] ; - sprintf(text, "Particles: %d ", nparticlein) ; - pavetext->AddText(text) ; - pavetext->Draw() ; - kinecanvas->Update(); + + //write filled histograms + efile->cd() ; + hAllEnergy->Write(0,kOverwrite) ; + hPhotEnergy->Write(0,kOverwrite) ; + hEMEnergy->Write(0,kOverwrite) ; + // efile->Write() ; + efile->Close() ; + delete efile ; } //____________________________________________________________________________ -void AliPHOSAnalyze::DisplayRecParticles() +void AliPHOSAnalyze::PositionResolution() { - if (fEvt == -999) { - cout << "DisplayRecParticles > Analyze an event first ... (y/n) " ; - Text_t answer[1] ; - cin >> answer ; cout << answer ; - if ( answer == "y" ) - AnalyzeOneEvent() ; - } - if (fEvt != -999) { + //fills two dimentional histo: energy vs. primary - reconstructed distance + + + + TH2F * hAllPosition = 0; // Position of any RP with primary photon + TH2F * hPhotPosition = 0; // Position of kGAMMA with primary photon + TH2F * hEMPosition = 0; // Position of EM with primary photon + + TH1F * hAllPositionX = 0; // X-Position Resolution of photons with photon primary + TH1F * hAllPositionZ = 0; // Z-Position Resolution of photons with photon primary + + + //opening file and reading histograms if any + TFile * pfile = new TFile("position.root","update"); + + hAllPosition = (TH2F*)pfile->Get("hAllPosition"); + if(hAllPosition == 0) + hAllPosition = new TH2F("hAllPosition", + "Position of any RP with primary photon",100, 0., 5., 100, 0., 5.); + hPhotPosition= (TH2F*)pfile->Get("hPhotPosition"); + if(hPhotPosition == 0) + hPhotPosition = new TH2F("hPhotPosition", + "Position of kGAMMA with primary photon",100, 0., 5., 100, 0., 5.); + hEMPosition= (TH2F*)pfile->Get("hEMPosition") ; + if(hEMPosition == 0) + hEMPosition = new TH2F("hEMPosition", + "Position of EM with primary photon", 100, 0., 5., 100, 0., 5.); + hAllPositionX = (TH1F*)pfile->Get("hAllPositionX") ; + if(hAllPositionX == 0) + hAllPositionX = new TH1F("hAllPositionX", + "Delta X of any RP with primary photon",100, -2., 2.); + hAllPositionZ =(TH1F*) pfile->Get("hAllPositionZ") ; + if(hAllPositionZ == 0) + hAllPositionZ = new TH1F("hAllPositionZ", + "Delta X of any RP with primary photon",100, -2., 2.); + + if (fRunLoader == 0x0) + { + AliError(Form("Error Loading session")); + return; + } + + AliPHOSLoader* gime = dynamic_cast(fRunLoader->GetLoader("PHOSLoader")); + if ( gime == 0 ) + { + AliError(Form("Could not obtain the Loader object !")); + return ; + } + + if (fRunLoader->TreeE() == 0x0) fRunLoader->LoadHeader(); + + AliPHOSGeometry * phosgeom = AliPHOSGeometry::GetInstance() ; + + Int_t ievent; + Int_t maxevent = (Int_t)fRunLoader->TreeE()->GetEntries() ; + for ( ievent=0; ievent < maxevent ; ievent++){ + + //read the current event + fRunLoader->GetEvent(ievent) ; + + //DP:Extract vertex position + Double_t vtx[3]={0.,0.,0.} ; + + TClonesArray * rp = gime->RecParticles() ; + if(!rp) { + AliError(Form("Event %d, Can't find RecParticles", ievent)) ; + return ; + } + TClonesArray * ts = gime->TrackSegments() ; + if(!ts) { + AliError(Form("Event %d, Can't find TrackSegments", ievent)) ; + return ; + } + TObjArray * emcrp = gime->EmcRecPoints() ; + if(!emcrp){ + AliError(Form("Event %d, Can't find EmcRecPoints", ievent)) ; + return ; + } + + + const AliPHOSRecParticle * recParticle ; + Int_t iRecParticle ; + for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast(); iRecParticle++ ){ + recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ; - Int_t module ; - cout << "DisplayRecParticles > which module (1-5, -1: all) ? " ; - cin >> module ; cout << module << endl ; - Text_t histoname[80] ; - sprintf(histoname,"Event %d: Reconstructed particles in module %d", fEvt, module) ; - Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module - fGeom->EmcModuleCoverage(module, tm, tM, pm, pM, kDegre) ; - Double_t theta, phi ; - fGeom->EmcXtalCoverage(theta, phi, kDegre) ; - Int_t tdim = (Int_t)( (tM - tm) / theta ) ; - Int_t pdim = (Int_t)( (pM - pm) / phi ) ; - tm -= theta ; - tM += theta ; - pm -= phi ; - TH2F * histoRparticle = new TH2F("histoRparticle", histoname, - pdim, pm, pM, tdim, tm, tM) ; - histoRparticle->SetStats(kFALSE) ; - Text_t canvasname[80] ; - sprintf(canvasname, "Reconstructed particles in PHOSmodule # %d", module) ; - TCanvas * rparticlecanvas = new TCanvas("RparticleCanvas", canvasname, 650, 500) ; - RecParticlesList * rpl = fPHOS->RecParticles() ; - Int_t nRecParticles = rpl->GetEntries() ; - Int_t nRecParticlesInModule = 0 ; - TIter nextRecPart(rpl) ; - AliPHOSRecParticle * rp ; - cout << "DisplayRecParticles > " << nRecParticles << " reconstructed particles " << endl ; - Double_t kRADDEG = 180. / TMath::Pi() ; - while ( (rp = (AliPHOSRecParticle *)nextRecPart() ) ) { - AliPHOSTrackSegment * ts = rp->GetPHOSTrackSegment() ; - if ( ts->GetPHOSMod() == module ) { - nRecParticlesInModule++ ; - Double_t theta = rp->Theta() * kRADDEG ; - Double_t phi = rp->Phi() * kRADDEG ; - Double_t energy = rp->Energy() ; - histoRparticle->Fill(phi, theta, energy) ; - } - } - histoRparticle->Draw("color") ; + //find the closest primary + Int_t moduleNumberRec ; + Double_t recX, recZ ; + phosgeom->ImpactOnEmc(vtx,recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ; + + Double_t minDistance = 100. ; + Int_t closestPrimary = -1 ; + + //extract list of primaries: it is stored at EMC RecPoints + Int_t emcIndex = ((AliPHOSTrackSegment*) ts->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ; + Int_t numberofprimaries ; + Int_t * listofprimaries = ((AliPHOSEmcRecPoint *) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries) ; - nextRecPart.Reset() ; - while ( (rp = (AliPHOSRecParticle *)nextRecPart() ) ) { - AliPHOSTrackSegment * ts = rp->GetPHOSTrackSegment() ; - if ( ts->GetPHOSMod() == module ) - rp->Draw("P") ; + Int_t index ; + const TParticle * primary ; + Double_t distance = minDistance ; + Double_t dX = 1000; // incredible number + Double_t dZ = 1000; // for the case if no primary will be found + Double_t dXmin = 0.; + Double_t dZmin = 0. ; + for ( index = 0 ; index < numberofprimaries ; index++){ + primary = fRunLoader->Stack()->Particle(listofprimaries[index]) ; + Int_t moduleNumber ; + Double_t primX, primZ ; + phosgeom->ImpactOnEmc(vtx,primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ; + if(moduleNumberRec == moduleNumber) { + dX = recX - primX; + dZ = recZ - primZ; + distance = TMath::Sqrt(dX*dX + dZ*dZ) ; + if(minDistance > distance) { + minDistance = distance ; + dXmin = dX; + dZmin = dZ; + closestPrimary = listofprimaries[index] ; + } + } + } + + //if found primary, fill histograms + if(closestPrimary >=0 ){ + primary = fRunLoader->Stack()->Particle(closestPrimary) ; + if(primary->GetPdgCode() == 22){ + hAllPosition->Fill(primary->Energy(), minDistance) ; + hAllPositionX->Fill(primary->Energy(), dX) ; + hAllPositionZ->Fill(primary->Energy(), dZ) ; + if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST){ + hPhotPosition->Fill(primary->Energy(), minDistance ) ; + hEMPosition->Fill(primary->Energy(), minDistance ) ; + } + else + if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW) + hEMPosition->Fill(primary->Energy(), minDistance ) ; + } } - - Text_t text[80] ; - sprintf(text, "reconstructed particles: %d", nRecParticlesInModule) ; - TPaveText * pavetext = new TPaveText(292, 100, 300, 101); - pavetext->AddText(text) ; - pavetext->Draw() ; - rparticlecanvas->Update() ; } -} + } + + //Write output histgrams + pfile->cd() ; + hAllPosition->Write(0,kOverwrite) ; + hAllPositionX->Write(0,kOverwrite) ; + hAllPositionZ->Write(0,kOverwrite) ; + hPhotPosition->Write(0,kOverwrite) ; + hEMPosition->Write(0,kOverwrite) ; + pfile->Write() ; + pfile->Close() ; + delete pfile ; + +} //____________________________________________________________________________ -void AliPHOSAnalyze::DisplayRecPoints() -{ - if (fEvt == -999) { - cout << "DisplayRecPoints > Analyze an event first ... (y/n) " ; - Text_t answer[1] ; - cin >> answer ; cout << answer ; - if ( answer == "y" ) - AnalyzeOneEvent() ; - } - if (fEvt != -999) { - - Int_t module ; - cout << "DisplayRecPoints > which module (1-5, -1: all) ? " ; - cin >> module ; cout << module << endl ; +void AliPHOSAnalyze::Contamination(){ +// fills spectra of primary photons and several kinds of +// reconstructed particles, so that analyzing them one can +// estimate conatmination, efficiency of registration etc. + + //define several general histograms + TH1F * hPrimary = 0; //spectrum (P_t distribution) of primary photons + TH1F * hAllRP = 0; //spectrum of all RecParticles in PHOS + TH1F * hPhot = 0; //spectrum of kGAMMA RecParticles + TH1F * hShape = 0; //spectrum of all EM RecParticles + TH1F * hVeto = 0; //spectrum of all neutral RecParticles + + //Now separate histograms in accoradance with primary + //primary - photon + TH1F * hPhotReg = 0; //Registeres as photon + TH1F * hPhotEM = 0; //Registered as EM + + //primary - n + TH1F * hNReg = 0; //Registeres as photon + TH1F * hNEM = 0; //Registered as EM + + //primary - nBar + TH1F * hNBarReg = 0; //Registeres as photon + TH1F * hNBarEM = 0; //Registered as EM + + //primary - charged hadron (pBar excluded) + TH1F * hChargedReg = 0; //Registeres as photon + TH1F * hChargedEM = 0; //Registered as EM + + //primary - pBar + TH1F * hPbarReg = 0; //Registeres as photon + TH1F * hPbarEM = 0; //Registered as EM + + + //Reading histograms from the file + TFile * cfile = new TFile("contamination.root","update") ; + + //read general histograms + hPrimary = (TH1F*) cfile->Get("hPrimary") ; + if(hPrimary == 0) + hPrimary= new TH1F("hPrimary", "Primary photon spectrum", 100, 0., 5.); + hAllRP = (TH1F*)cfile->Get("hAllRP") ; + if(hAllRP == 0) + hAllRP = new TH1F("hAllRP","All Reconstructed particles", 100, 0., 5.); + hPhot = (TH1F*)cfile->Get("hPhot") ; + if(hPhot == 0) + hPhot = new TH1F("hPhot","All kGAMMA RecParticles",100, 0., 5.); + hShape = (TH1F*) cfile->Get("hShape") ; + if(hShape == 0) + hShape = new TH1F("hShape","All particles with EM shower",100, 0., 5.); + hVeto= (TH1F*)cfile->Get("hVeto") ; + if(hVeto == 0) + hVeto = new TH1F("hVeto", "All uncharged particles", 100, 0., 5.); + + + //primary - photon + hPhotReg = (TH1F*)cfile->Get("hPhotReg"); + if(hPhotReg == 0) + hPhotReg = new TH1F("hPhotReg","Photon registered as photon",100, 0., 5.); + hPhotEM =(TH1F*)cfile->Get("hPhotEM"); + if(hPhotEM== 0) + hPhotEM = new TH1F("hPhotEM", "Photon registered as EM", 100, 0., 5.); + + //primary - n + hNReg = (TH1F*)cfile->Get("hNReg"); + if(hNReg== 0) + hNReg = new TH1F("hNReg", "N registered as photon", 100, 0., 5.); + hNEM = (TH1F*)cfile->Get("hNEM"); + if(hNEM== 0) + hNEM = new TH1F("hNEM", "N registered as EM", 100, 0., 5.); + + //primary - nBar + hNBarReg =(TH1F*)cfile->Get("hNBarReg"); + if(hNBarReg== 0) + hNBarReg = new TH1F("hNBarReg", "NBar registered as photon", 100, 0., 5.); + hNBarEM =(TH1F*)cfile->Get("hNBarEM"); + if(hNBarEM== 0) + hNBarEM = new TH1F("hNBarEM", "NBar registered as EM", 100, 0., 5.); + + //primary - charged hadron (pBar excluded) + hChargedReg = (TH1F*)cfile->Get("hChargedReg"); + if(hChargedReg== 0) + hChargedReg= new TH1F("hChargedReg", "Charged hadron registered as photon",100, 0., 5.); + hChargedEM = (TH1F*)cfile->Get("hChargedEM"); + if(hChargedEM== 0) + hChargedEM= new TH1F("hChargedEM","Charged registered as EM",100, 0., 5.); + + //primary - pBar + hPbarReg = (TH1F*)cfile->Get("hPbarReg"); + if(hPbarReg== 0) + hPbarReg= new TH1F("hPbarReg", "pBar registered as photon",100, 0., 5.); + hPbarEM = (TH1F*)cfile->Get("hPbarEM"); + if(hPbarEM== 0) + hPbarEM= new TH1F("hPbarEM","Pbar registered as EM",100, 0., 5.); + - Text_t canvasname[80]; - sprintf(canvasname,"Digits in PHOS/EMC module # %d",module) ; - TCanvas * modulecanvas = new TCanvas("module", canvasname, 650, 500) ; - modulecanvas->Draw() ; + //Now make some initializations - //=========== Creating 2d-histogram of the PHOS module - // a little bit junkie but is used to test Geom functinalities + Int_t counter[8][5] ; //# of registered particles + Int_t i1,i2 ; + for(i1 = 0; i1<8; i1++) + for(i2 = 0; i2<5; i2++) + counter[i1][i2] = 0 ; - Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module - - fGeom->EmcModuleCoverage(module, tm, tM, pm, pM); - // convert angles into coordinates local to the EMC module of interest - - Int_t emcModuleNumber ; - Double_t emcModulexm, emcModulezm ; // minimum local coordinate in a given EMCA module - Double_t emcModulexM, emcModulezM ; // maximum local coordinate in a given EMCA module - fGeom->ImpactOnEmc(tm, pm, emcModuleNumber, emcModulezm, emcModulexm) ; - fGeom->ImpactOnEmc(tM, pM, emcModuleNumber, emcModulezM, emcModulexM) ; - Int_t xdim = (Int_t)( ( emcModulexM - emcModulexm ) / fGeom->GetCrystalSize(0) ) ; - Int_t zdim = (Int_t)( ( emcModulezM - emcModulezm ) / fGeom->GetCrystalSize(2) ) ; - Float_t xmin = emcModulexm - fGeom->GetCrystalSize(0) ; - Float_t xMax = emcModulexM + fGeom->GetCrystalSize(0) ; - Float_t zmin = emcModulezm - fGeom->GetCrystalSize(2) ; - Float_t zMax = emcModulezM + fGeom->GetCrystalSize(2) ; - Text_t histoname[80]; - sprintf(histoname,"Event %d: Digits and RecPoints in module %d", fEvt, module) ; - TH2F * hModule = new TH2F("HistoReconstructed", histoname, - xdim, xmin, xMax, zdim, zmin, zMax) ; - hModule->SetMaximum(2.0); - hModule->SetMinimum(0.0); - hModule->SetStats(kFALSE); - - TIter next(fPHOS->Digits()) ; - Float_t energy, y, z; - Float_t etot=0.; - Int_t relid[4]; Int_t nDigits = 0 ; - AliPHOSDigit * digit ; - while((digit = (AliPHOSDigit *)next())) - { - fGeom->AbsToRelNumbering(digit->GetId(), relid) ; - if (relid[0] == module) - { - if (energy > 0.025 ){ - nDigits++ ; - energy = fClu->Calibrate(digit->GetAmp()) ; - etot += energy ; - fGeom->RelPosInModule(relid,y,z) ; - hModule->Fill(y, z, energy) ; - } - } - } - cout <<"DrawRecPoints > Found in module " - << module << " " << nDigits << " digits with total energy " << etot << endl ; - hModule->Draw("col2") ; - - //=========== Cluster in module - - TClonesArray * emcRP = fPHOS->EmcClusters() ; - etot = 0.; - Int_t totalnClusters = 0 ; - Int_t nClusters = 0 ; - TIter nextemc(emcRP) ; - AliPHOSEmcRecPoint * emc ; - while((emc = (AliPHOSEmcRecPoint *)nextemc())) - { - // Int_t numberofprimaries ; - // Int_t * primariesarray = new Int_t[10] ; - // emc->GetPrimaries(numberofprimaries, primariesarray) ; - totalnClusters++ ; - if ( emc->GetPHOSMod() == module ) - { - nClusters++ ; - energy = emc->GetTotalEnergy() ; - etot+= energy ; - emc->Draw("M") ; - } - } - cout << "DrawRecPoints > Found " << totalnClusters << " EMC Clusters in PHOS" << endl ; - cout << "DrawRecPoints > Found in module " << module << " " << nClusters << " EMC Clusters " << endl ; - cout << "DrawRecPoints > total energy " << etot << endl ; - - TPaveText * pavetext = new TPaveText(22, 80, 83, 90); - Text_t text[40] ; - sprintf(text, "digits: %d; clusters: %d", nDigits, nClusters) ; - pavetext->AddText(text) ; - pavetext->Draw() ; - modulecanvas->Update(); - - //=========== Cluster in module PPSD Down - - TClonesArray * ppsdRP = fPHOS->PpsdClusters() ; - etot = 0.; - TIter nextPpsd(ppsdRP) ; - AliPHOSPpsdRecPoint * ppsd ; - while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsd())) - { - totalnClusters++ ; - if ( ppsd->GetPHOSMod() == module ) - { - nClusters++ ; - energy = ppsd->GetEnergy() ; - etot+=energy ; - if (!ppsd->GetUp()) ppsd->Draw("P") ; - } - } - cout << "DrawRecPoints > Found " << totalnClusters << " Ppsd Down Clusters in PHOS" << endl ; - cout << "DrawRecPoints > Found in module " << module << " " << nClusters << " Ppsd Down Clusters " << endl ; - cout << "DrawRecPoints > total energy " << etot << endl ; - - //=========== Cluster in module PPSD Up - - ppsdRP = fPHOS->PpsdClusters() ; - etot = 0.; - TIter nextPpsdUp(ppsdRP) ; - while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsdUp())) - { - totalnClusters++ ; - if ( ppsd->GetPHOSMod() == module ) - { - nClusters++ ; - energy = ppsd->GetEnergy() ; - etot+=energy ; - if (ppsd->GetUp()) ppsd->Draw("P") ; - } - } - cout << "DrawRecPoints > Found " << totalnClusters << " Ppsd Up Clusters in PHOS" << endl ; - cout << "DrawRecPoints > Found in module " << module << " " << nClusters << " Ppsd Up Clusters " << endl ; - cout << "DrawRecPoints > total energy " << etot << endl ; + + + if (fRunLoader == 0x0) + { + AliError(Form("Error Loading session")); + return; + } + + AliPHOSLoader* gime = dynamic_cast(fRunLoader->GetLoader("PHOSLoader")); + if ( gime == 0 ) + { + AliError(Form("Could not obtain the Loader object !")); + return ; + } + + if (fRunLoader->TreeE() == 0x0) fRunLoader->LoadHeader(); + AliPHOSGeometry * phosgeom = AliPHOSGeometry::GetInstance() ; + + Int_t ievent; + Int_t maxevent = (Int_t)fRunLoader->TreeE()->GetEntries() ; + for ( ievent=0; ievent < maxevent ; ievent++){ - } // if !-999 -} + fRunLoader->GetEvent(ievent) ; + + //DP:Extract vertex position + Double_t vtx[3]={0.,0.,0.} ; -//____________________________________________________________________________ -void AliPHOSAnalyze::DisplayTrackSegments() -{ - if (fEvt == -999) { - cout << "DisplayTrackSegments > Analyze an event first ... (y/n) " ; - Text_t answer[1] ; - cin >> answer ; cout << answer ; - if ( answer == "y" ) - AnalyzeOneEvent() ; - } - if (fEvt != -999) { - - Int_t module ; - cout << "DisplayTrackSegments > which module (1-5, -1: all) ? " ; - cin >> module ; cout << module << endl ; - //=========== Creating 2d-histogram of the PHOS module - // a little bit junkie but is used to test Geom functinalities + TClonesArray * rp = gime->RecParticles() ; + if(!rp) { + AliError(Form("Event %d, Can't find RecParticles", ievent)) ; + return ; + } + TClonesArray * ts = gime->TrackSegments() ; + if(!ts) { + AliError(Form("Event %d, Can't find TrackSegments", ievent)) ; + return ; + } + TObjArray * emcrp = gime->EmcRecPoints() ; + if(!emcrp){ + AliError(Form("Event %d, Can't find EmcRecPoints", ievent)) ; + return ; + } + + + //=========== Make spectrum of the primary photons + const TParticle * primary ; + Int_t iPrimary ; + for( iPrimary = 0 ; iPrimary < fRunLoader->Stack()->GetNprimary() ; iPrimary++){ + primary = fRunLoader->Stack()->Particle(iPrimary) ; + Int_t primaryType = primary->GetPdgCode() ; + if( primaryType == 22 ) { + //check, if photons folls onto PHOS + Int_t moduleNumber ; + Double_t primX, primZ ; + phosgeom->ImpactOnEmc(vtx,primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ; + if(moduleNumber) + hPrimary->Fill(primary->Energy()) ; + + } - Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module + } + + //========== Now scan over RecParticles + const AliPHOSRecParticle * recParticle ; + Int_t iRecParticle ; + for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast(); iRecParticle++ ){ + recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ; + //fill histo spectrum of all RecParticles + hAllRP->Fill(CorrectedEnergy(recParticle->Energy())) ; - fGeom->EmcModuleCoverage(module, tm, tM, pm, pM); - // convert angles into coordinates local to the EMC module of interest + //==========find the closest primary + Int_t moduleNumberRec ; + Double_t recX, recZ ; + phosgeom->ImpactOnEmc(vtx,recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ; - Int_t emcModuleNumber ; - Double_t emcModulexm, emcModulezm ; // minimum local coordinate in a given EMCA module - Double_t emcModulexM, emcModulezM ; // maximum local coordinate in a given EMCA module - fGeom->ImpactOnEmc(tm, pm, emcModuleNumber, emcModulezm, emcModulexm) ; - fGeom->ImpactOnEmc(tM, pM, emcModuleNumber, emcModulezM, emcModulexM) ; - Int_t xdim = (Int_t)( ( emcModulexM - emcModulexm ) / fGeom->GetCrystalSize(0) ) ; - Int_t zdim = (Int_t)( ( emcModulezM - emcModulezm ) / fGeom->GetCrystalSize(2) ) ; - Float_t xmin = emcModulexm - fGeom->GetCrystalSize(0) ; - Float_t xMax = emcModulexM + fGeom->GetCrystalSize(0) ; - Float_t zmin = emcModulezm - fGeom->GetCrystalSize(2) ; - Float_t zMax = emcModulezM + fGeom->GetCrystalSize(2) ; - Text_t histoname[80]; - sprintf(histoname,"Event %d: Track Segments in module %d", fEvt, module) ; - TH2F * histotrack = new TH2F("histotrack", histoname, - xdim, xmin, xMax, zdim, zmin, zMax) ; - histotrack->SetStats(kFALSE); - Text_t canvasname[80]; - sprintf(canvasname,"Track segments in PHOS/EMC-PPSD module # %d", module) ; - TCanvas * trackcanvas = new TCanvas("TrackSegmentCanvas", canvasname, 650, 500) ; - histotrack->Draw() ; - - TrackSegmentsList * trsegl = fPHOS->TrackSegments() ; - AliPHOSTrackSegment * trseg ; - - Int_t nTrackSegments = trsegl->GetEntries() ; + Double_t minDistance = 100. ; + Int_t closestPrimary = -1 ; + + //extract list of primaries: it is stored at EMC RecPoints + Int_t emcIndex = ((AliPHOSTrackSegment *) ts->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ; + Int_t numberofprimaries ; + Int_t * listofprimaries = ((AliPHOSEmcRecPoint *) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries) ; Int_t index ; - Float_t etot = 0 ; - Int_t nTrackSegmentsInModule = 0 ; - for(index = 0; index < nTrackSegments ; index++){ - trseg = (AliPHOSTrackSegment * )trsegl->At(index) ; - etot+= trseg->GetEnergy() ; - if ( trseg->GetPHOSMod() == module ) { - nTrackSegmentsInModule++ ; - trseg->Draw("P"); - } - } - Text_t text[80] ; - sprintf(text, "track segments: %d", nTrackSegmentsInModule) ; - TPaveText * pavetext = new TPaveText(22, 80, 83, 90); - pavetext->AddText(text) ; - pavetext->Draw() ; - trackcanvas->Update() ; - cout << "DisplayTrackSegments > Found " << trsegl->GetEntries() << " Track segments with total energy "<< etot << endl ; - - } -} -//____________________________________________________________________________ -Bool_t AliPHOSAnalyze::OpenRootFile(Text_t * name) -{ - fRootFile = new TFile(name) ; - return fRootFile->IsOpen() ; -} -//____________________________________________________________________________ -void AliPHOSAnalyze::SavingHistograms() -{ - Text_t outputname[80] ;// = fRootFile->GetName(); - sprintf(outputname,"%s.analyzed",fRootFile->GetName()); - TFile output(outputname,"RECREATE"); - output.cd(); - if (fhEmcDigit ) - fhEmcDigit->Write() ; - if (fhVetoDigit ) - fhVetoDigit->Write() ; - if (fhConvertorDigit ) - fhConvertorDigit->Write() ; - if (fhEmcCluster ) - fhEmcCluster->Write() ; - if (fhVetoCluster ) - fhVetoCluster->Write() ; - if (fhConvertorCluster ) - fhConvertorCluster->Write() ; - if (fhConvertorEmc ) - fhConvertorEmc->Write() ; - if (fhPhotonEnergy) - fhPhotonEnergy->Write() ; - if (fhPhotonPositionX) - fhPhotonPositionX->Write() ; - if (fhPhotonPositionY) - fhPhotonPositionX->Write() ; - if (fhElectronEnergy) - fhElectronEnergy->Write() ; - if (fhElectronPositionX) - fhElectronPositionX->Write() ; - if (fhElectronPositionY) - fhElectronPositionX->Write() ; - if (fhNeutralHadronEnergy) - fhNeutralHadronEnergy->Write() ; - if (fhNeutralHadronPositionX) - fhNeutralHadronPositionX->Write() ; - if (fhNeutralHadronPositionY) - fhNeutralHadronPositionX->Write() ; - if (fhNeutralEMEnergy) - fhNeutralEMEnergy->Write() ; - if (fhNeutralEMPositionX) - fhNeutralEMPositionX->Write() ; - if (fhNeutralEMPositionY) - fhNeutralEMPositionX->Write() ; - if (fhChargedHadronEnergy) - fhChargedHadronEnergy->Write() ; - if (fhChargedHadronPositionX) - fhChargedHadronPositionX->Write() ; - if (fhChargedHadronPositionY) - fhChargedHadronPositionX->Write() ; - if (fhPhotonHadronEnergy) - fhPhotonHadronEnergy->Write() ; - if (fhPhotonHadronPositionX) - fhPhotonHadronPositionX->Write() ; - if (fhPhotonHadronPositionY) - fhPhotonHadronPositionX->Write() ; - - output.Write(); - output.Close(); + Double_t distance = minDistance ; + Double_t dX, dZ; + Double_t dXmin = 0.; + Double_t dZmin = 0. ; + for ( index = 0 ; index < numberofprimaries ; index++){ + primary = fRunLoader->Stack()->Particle(listofprimaries[index]) ; + Int_t moduleNumber ; + Double_t primX, primZ ; + phosgeom->ImpactOnEmc(vtx,primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ; + if(moduleNumberRec == moduleNumber) { + dX = recX - primX; + dZ = recZ - primZ; + distance = TMath::Sqrt(dX*dX + dZ*dZ) ; + if(minDistance > distance) { + minDistance = distance ; + dXmin = dX; + dZmin = dZ; + closestPrimary = listofprimaries[index] ; + } + } + } + + //===========define the "type" of closest primary + if(closestPrimary >=0 ){ + Int_t primaryCode = -1; + primary = fRunLoader->Stack()->Particle(closestPrimary) ; + Int_t primaryType = primary->GetPdgCode() ; + if(primaryType == 22) // photon ? + primaryCode = 0 ; + else + if(primaryType == 2112) // neutron + primaryCode = 1 ; + else + if(primaryType == -2112) // Anti neutron + primaryCode = 2 ; + else + if(primaryType == -2122) //Anti proton + primaryCode = 4 ; + else { + TParticle tempo(*primary) ; + if(tempo.GetPDG()->Charge()) + primaryCode = 3 ; + } + + //==========Now look at the type of RecParticle + Float_t energy = CorrectedEnergy(recParticle->Energy()) ; + if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST){ + hPhot->Fill(energy ) ; + switch(primaryCode){ + case 0: + hPhotReg->Fill(energy ) ; + break ; + case 1: + hNReg->Fill(energy ) ; + break ; + case 2: + hNBarReg->Fill(energy ) ; + break ; + case 3: + hChargedReg->Fill(energy ) ; + break ; + case 4: + hPbarReg->Fill(energy ) ; + break ; + default: + break ; + } + } + if((recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST)|| + (recParticle->GetType() == AliPHOSFastRecParticle::kCHARGEDEMFAST)|| + (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW)|| + (recParticle->GetType() == AliPHOSFastRecParticle::kCHARGEDEMSLOW) ){ //with EM shower + hShape->Fill(energy ) ; + switch(primaryCode){ + case 0: + hPhotEM->Fill(energy ) ; + break ; + case 1: + hNEM->Fill(energy ) ; + break ; + case 2: + hNBarEM->Fill(energy ) ; + break ; + case 3: + hChargedEM->Fill(energy ) ; + break ; + case 4: + hPbarEM->Fill(energy ) ; + break ; + default: + break ; + } + } + + if((recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST)|| + (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALHAFAST) || + (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW) || + (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALHASLOW) ) //nuetral + hVeto->Fill(energy ) ; + + //fill number of primaries identified as ... + if(primaryCode >= 0) // Primary code defined + counter[recParticle->GetType()][primaryCode]++ ; + + } + + } // no closest primary found + } + + + //=================== SaveHistograms + cfile->cd() ; + hPrimary->Write(0,kOverwrite); + hAllRP->Write(0,kOverwrite); + hPhot->Write(0,kOverwrite); + hShape->Write(0,kOverwrite); + hVeto->Write(0,kOverwrite); + hPhotReg->Write(0,kOverwrite); + hPhotEM->Write(0,kOverwrite); + hNReg ->Write(0,kOverwrite); + hNEM ->Write(0,kOverwrite); + hNBarReg ->Write(0,kOverwrite); + hNBarEM ->Write(0,kOverwrite); + hChargedReg ->Write(0,kOverwrite); + hChargedEM ->Write(0,kOverwrite); + hPbarReg ->Write(0,kOverwrite); + hPbarEM ->Write(0,kOverwrite); + + cfile->Write(0,kOverwrite); + cfile->Close(); + delete cfile ; + + + //print Final Table + maxevent = (Int_t)AliRunLoader::Instance()->TreeE()->GetEntries() ; + + TString message ; + message = "Resolutions: Analyzed %d event(s)\n" ; + + message += " Primary: Photon Neutron Antineutron Charged hadron AntiProton\n" ; + message += "--------------------------------------------------------------------------------\n" ; + message += " kGAMMA: " ; + message += "%d %d %d %d %d\n" ; + message += " kGAMMAHA: " ; + message += "%d %d %d %d %d\n" ; + message += " kNEUTRALEM: " ; + message += "%d %d %d %d %d\n" ; + message += " kNEUTRALHA: " ; + message += "%d %d %d %d %d\n" ; + message += " kABSURDEM: "; + message += "%d %d %d %d %d\n" ; + message += " kABSURDHA: " ; + message += "%d %d %d %d %d\n" ; + message += " kELECTRON: " ; + message += "%d %d %d %d %d\n" ; + message += " kCHARGEDHA: " ; + message += "%d %d %d %d %d\n" ; + + message += "--------------------------------------------------------------------------------" ; + + + Int_t totalInd = 0 ; + for(i1 = 0; i1<8; i1++) + for(i2 = 0; i2<5; i2++) + totalInd+=counter[i1][i2] ; + message += "Indentified particles: %d" ; + + AliInfo(Form(message.Data(), maxevent, + counter[2][0], counter[2][1], counter[2][2], counter[2][3], counter[2][4], + counter[3][0], counter[3][1], counter[3][2], counter[3][3], counter[3][4], + counter[0][0], counter[0][1], counter[0][2], counter[0][3], counter[0][4], + counter[1][0], counter[1][1], counter[1][2], counter[1][3], counter[1][4], + counter[4][0], counter[4][1], counter[4][2], counter[4][3], counter[4][4], + counter[5][0], counter[5][1], counter[5][2], counter[5][3], counter[5][4], + counter[6][0], counter[6][1], counter[6][2], counter[6][3], counter[6][4], + counter[7][0], counter[7][1], counter[7][2], counter[7][3], counter[7][4], + totalInd )) ; + }