X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PHOS%2FAliPHOSAnalyze.cxx;h=69c72ca7c8edf8400a2bc7a08d0b73228d5359c7;hb=146aa588911dc12a6bae78254624f277975a3b85;hp=e4e48774044bac3e798a6fddedda8068c36dbaa5;hpb=1b5d8c54b89d252178996872bd70626e0fcbc611;p=u%2Fmrichter%2FAliRoot.git diff --git a/PHOS/AliPHOSAnalyze.cxx b/PHOS/AliPHOSAnalyze.cxx index e4e48774044..69c72ca7c8e 100644 --- a/PHOS/AliPHOSAnalyze.cxx +++ b/PHOS/AliPHOSAnalyze.cxx @@ -16,1670 +16,1293 @@ /* $Id$ */ //_________________________________________________________________________ -// Algorythm class to analyze PHOSv1 events: -// Construct histograms and displays them. -// Use the macro EditorBar.C for best access to the functionnalities +// 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: Y. Schutz (SUBATECH) & Gines Martinez (SUBATECH) +//*-- 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 "TStyle.h" +#include "TROOT.h" // --- Standard library --- -#include -#include - // --- AliRoot header files --- -#include "AliRun.h" -#include "AliPHOSv1.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 "AliPHOSIndexToObject.h" -#include "AliPHOSHit.h" -#include "AliPHOSCpvRecPoint.h" -#include "AliPHOSPpsdRecPoint.h" +#include "AliPHOSLoader.h" + ClassImp(AliPHOSAnalyze) //____________________________________________________________________________ - AliPHOSAnalyze::AliPHOSAnalyze() +AliPHOSAnalyze::AliPHOSAnalyze(): + fCorrection(1.2), //Value calculated for default parameters of reconstruction + fEvt(0), + ffileName(), + fRunLoader(0) { // default ctor (useless) - - fRootFile = 0 ; } //____________________________________________________________________________ -AliPHOSAnalyze::AliPHOSAnalyze(Text_t * name) +AliPHOSAnalyze::AliPHOSAnalyze(Text_t * fileName): + fCorrection(1.05), //Value calculated for default parameters of reconstruction + fEvt(0), + ffileName(fileName), + fRunLoader(0) { // ctor: analyze events from root file "name" - - Bool_t ok = OpenRootFile(name) ; - if ( !ok ) { - cout << " AliPHOSAnalyze > Error opening " << name << endl ; - } - else { - //========== Get AliRun object from file - gAlice = (AliRun*) fRootFile->Get("gAlice") ; - - //=========== Get the PHOS object and associated geometry from the file - fPHOS = (AliPHOSv1 *)gAlice->GetDetector("PHOS") ; - fGeom = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() ); - - //========== Initializes the Index to Object converter - fObjGetter = AliPHOSIndexToObject::GetInstance() ; // <--- To be redone - //========== Current event number - fEvt = -999 ; - - } - fDebugLevel = 0; - // fClu = 0 ; - // fPID = 0 ; - // fTrs = 0 ; - // fRec = 0 ; - ResetHistograms() ; + fRunLoader = AliRunLoader::Open(fileName,"AliPHOSAnalyze"); + if (fRunLoader == 0x0) + { + AliError(Form("Error Loading session")); + } } //____________________________________________________________________________ -AliPHOSAnalyze::AliPHOSAnalyze(const AliPHOSAnalyze & ana) +AliPHOSAnalyze::AliPHOSAnalyze(const AliPHOSAnalyze & ana): + TObject(ana), + fCorrection(0.), + fEvt(0), + ffileName(), + fRunLoader(0) { // copy ctor ( (AliPHOSAnalyze &)ana ).Copy(*this) ; } -//____________________________________________________________________________ -void AliPHOSAnalyze::Copy(TObject & obj) -{ - // copy an analysis into an other one - TObject::Copy(obj) ; - // I do nothing more because the copy is silly but the Code checkers requires one -} - //____________________________________________________________________________ AliPHOSAnalyze::~AliPHOSAnalyze() { // dtor - if(fRootFile->IsOpen()) fRootFile->Close() ; - if(fRootFile) {delete fRootFile ; fRootFile=0 ;} - if(fPHOS) {delete fPHOS ; fPHOS =0 ;} - // if(fClu) {delete fClu ; fClu =0 ;} - // if(fPID) {delete fPID ; fPID =0 ;} - // if(fRec) {delete fRec ; fRec =0 ;} - // if(fTrs) {delete fTrs ; fTrs =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. - -// TH2F * digitOccupancy = new TH2F("digitOccupancy","EMC digits", 64,-71.,71.,64,-71.,71.); -// TH2F * sdigitOccupancy = new TH2F("sdigitOccupancy","EMC sdigits", 64,-71.,71.,64,-71.,71.); -// TH2F * emcOccupancy = new TH2F("emcOccupancy","EMC RecPoints",64,-71.,71.,64,-71.,71.); -// TH2F * ppsdUp = new TH2F("ppsdUp","PPSD Up digits", 128,-71.,71.,128,-71.,71.) ; -// TH2F * ppsdUpCl = new TH2F("ppsdUpCl","PPSD Up RecPoints",128,-71.,71.,128,-71.,71.) ; -// TH2F * ppsdLow = new TH2F("ppsdLow","PPSD Low digits", 128,-71.,71.,128,-71.,71.) ; -// TH2F * ppsdLowCl = new TH2F("ppsdLowCl","PPSD Low RecPoints",128,-71.,71.,128,-71.,71.) ; -// TH2F * nbar = new TH2F("nbar","Primary nbar", 64,-71.,71.,64,-71.,71.); -// TH2F * phot = new TH2F("phot","Primary Photon", 64,-71.,71.,64,-71.,71.); -// TH2F * charg = new TH2F("charg","Primary charged",64,-71.,71.,64,-71.,71.); -// TH2F * recPhot = new TH2F("recPhot","RecParticles with primary Photon",64,-71.,71.,64,-71.,71.); -// TH2F * recNbar = new TH2F("recNbar","RecParticles with primary Nbar", 64,-71.,71.,64,-71.,71.); - -// //========== Create the Clusterizer -// // fClu = new AliPHOSClusterizerv1() ; + //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 + + + //Plot Primary Particles -// gAlice->GetEvent(Nevent); + if (fRunLoader->Stack() == 0x0) fRunLoader->LoadKinematics("READ"); -// TParticle * primary ; -// Int_t iPrimary ; -// for ( iPrimary = 0 ; iPrimary < gAlice->GetNtrack() ; iPrimary++) -// { -// primary = gAlice->Particle(iPrimary) ; -// Int_t primaryType = primary->GetPdgCode() ; -// if( (primaryType == 211)||(primaryType == -211)||(primaryType == 2212)||(primaryType == -2212) ) { + + 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 ; -// fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, 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 ; -// fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ; -// if(moduleNumber==Nmod) -// phot->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 ; -// fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ; +// phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ; // if(moduleNumber==Nmod) // nbar->Fill(primZ,primX,primary->Energy()) ; // } // } -// } - -// //Set TreeS here and get AliPHOSSdigitizer + } - -// gAlice->TreeS()->GetEvent(0) ; - -// Int_t iSDigit ; -// AliPHOSDigit * sdigit ; -// if(fPHOS->SDigits()){ -// for(iSDigit = 0; iSDigit < fPHOS->SDigits()->GetEntries(); iSDigit++) -// { -// sdigit = (AliPHOSDigit *) fPHOS->SDigits()->At(iSDigit) ; -// Int_t relid[4]; -// fGeom->AbsToRelNumbering(sdigit->GetId(), relid) ; -// Float_t x,z ; -// fGeom->RelPosInModule(relid,x,z) ; -// Float_t e = 1 ; //<--- fPHOS->Calibrate(sdigit->GetAmp()) ; -// if(relid[0]==Nmod){ -// if(relid[1]==0) //EMC -// sdigitOccupancy->Fill(x,z,e) ; -// if((relid[1]>0)&&(relid[1]<17)) -// ppsdUp->Fill(x,z,e) ; -// if(relid[1]>16) -// ppsdLow->Fill(x,z,e) ; -// } -// } -// } -// else{ -// cout << "No SDigits read " << endl ; -// } + 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) ; + } + } + } -// gAlice->TreeD()->GetEvent(0) ; -// if(fPHOS->Digits()){ -// Int_t iDigit ; -// AliPHOSDigit * digit ; -// for(iDigit = 0; iDigit < fPHOS->Digits()->GetEntries(); iDigit++) -// { -// digit = (AliPHOSDigit *) fPHOS->Digits()->At(iDigit) ; -// Int_t relid[4]; -// fGeom->AbsToRelNumbering(digit->GetId(), relid) ; -// Float_t x,z ; -// fGeom->RelPosInModule(relid,x,z) ; -// Float_t e = 1; //<--- fClu->Calibrate(digit->GetAmp()) ; -// if(relid[0]==Nmod){ -// if(relid[1]==0) //EMC -// digitOccupancy->Fill(x,z,e) ; -// if((relid[1]>0)&&(relid[1]<17)) -// ppsdUp->Fill(x,z,e) ; -// if(relid[1]>16) -// ppsdLow->Fill(x,z,e) ; -// } -// } -// } -// else{ -// cout << "No Digits read " << endl ; -// } - -// gAlice->TreeR()->GetEvent(0) ; + //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 * emcRecPoints = fPHOS->EmcRecPoints() ; -// TObjArray * ppsdRecPoints = fPHOS->PpsdRecPoints() ; -// TClonesArray * recParticleList = fPHOS->RecParticles() ; + 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()) ; + } + } + } + } -// Int_t irecp ; -// TVector3 pos ; + //Plot made histograms + emcSdigits->Draw("box") ; + emcDigits->SetLineColor(5) ; + emcDigits->Draw("boxsame") ; + emcRecPoints->SetLineColor(2) ; + emcRecPoints->Draw("boxsame") ; + cpvSdigits->SetLineColor(1) ; + cpvSdigits->Draw("boxsame") ; -// if(emcRecPoints ){ -// for(irecp = 0; irecp < emcRecPoints->GetEntries() ; irecp ++){ -// AliPHOSEmcRecPoint * emc= (AliPHOSEmcRecPoint*)emcRecPoints->At(irecp) ; -// if(emc->GetPHOSMod()==Nmod){ -// emc->GetLocalPosition(pos) ; -// emcOccupancy->Fill(pos.X(),pos.Z(),emc->GetEnergy()); -// } -// } -// } -// else{ -// cout << "No EMC rec points read " << endl ; -// } - -// if(ppsdRecPoints ){ -// for(irecp = 0; irecp < ppsdRecPoints->GetEntries() ; irecp ++){ -// AliPHOSPpsdRecPoint * ppsd= (AliPHOSPpsdRecPoint *)ppsdRecPoints->At(irecp) ; - -// ppsd->GetLocalPosition(pos) ; -// cout << "PPSD " << irecp << " " << ppsd->GetPHOSMod() << " " << pos.X() << " " << pos.Z() << endl ; - -// if(ppsd->GetPHOSMod()==Nmod){ -// ppsd->GetLocalPosition(pos) ; -// if(ppsd->GetUp()) -// ppsdUpCl->Fill(pos.X(),pos.Z(),ppsd->GetEnergy()); -// else -// ppsdLowCl->Fill(pos.X(),pos.Z(),ppsd->GetEnergy()); -// } -// } -// } -// else{ -// cout << "No PPSD/CPV rec points read " << endl ; -// } - -// AliPHOSRecParticle * recParticle ; -// Int_t iRecParticle ; -// if(recParticleList ){ -// for(iRecParticle = 0; iRecParticle < recParticleList->GetEntries() ;iRecParticle++ ) -// { -// recParticle = (AliPHOSRecParticle *) recParticleList->At(iRecParticle) ; - -// Int_t moduleNumberRec ; -// Double_t recX, recZ ; -// fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ; -// if(moduleNumberRec == Nmod){ - -// Double_t minDistance = 5. ; -// Int_t closestPrimary = -1 ; - -// Int_t numberofprimaries ; -// Int_t * listofprimaries = recParticle->GetPrimaries(numberofprimaries) ; -// Int_t index ; -// TParticle * primary ; -// Double_t distance = minDistance ; - -// for ( index = 0 ; index < numberofprimaries ; index++){ -// primary = gAlice->Particle(listofprimaries[index]) ; -// Int_t moduleNumber ; -// Double_t primX, primZ ; -// fGeom->ImpactOnEmc(primary->Theta(), primary->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 = gAlice->Particle(closestPrimary)->GetPdgCode() ; - -// if(primaryType==22) -// recPhot->Fill(recZ,recX,recParticle->Energy()) ; -// else -// if(primaryType==-2112) -// recNbar->Fill(recZ,recX,recParticle->Energy()) ; -// } -// } -// } -// } -// else{ -// cout << "Not Rec Prticles read " << endl ; -// } +} +//____________________________________________________________________________ +void AliPHOSAnalyze::Ls(){ + //lists branches and titles of PHOS-related branches of TreeR, TreeD, TreeS -// digitOccupancy->Draw("box") ; -// sdigitOccupancy->SetLineColor(5) ; -// sdigitOccupancy->Draw("box") ; -// emcOccupancy->SetLineColor(2) ; -// emcOccupancy->Draw("boxsame") ; -// ppsdUp->SetLineColor(3) ; -// ppsdUp->Draw("boxsame") ; -// ppsdLow->SetLineColor(4) ; -// ppsdLow->Draw("boxsame") ; -// phot->SetLineColor(8) ; -// phot->Draw("boxsame") ; -// nbar->SetLineColor(6) ; -// nbar->Draw("boxsame") ; + if (fRunLoader == 0x0) + { + AliError(Form("Error Loading session")); + return; + } -} -// //____________________________________________________________________________ -// void AliPHOSAnalyze::Reconstruct(Int_t nevents,Int_t firstEvent ) -// { - -// // Performs reconstruction of EMC and CPV (GPS2, IHEP or MIXT) -// // for events from FirstEvent to Nevents - -// Int_t ievent ; -// for ( ievent=firstEvent; ievent Starting Reconstructing " << endl ; -// //========== Create the Clusterizer -// fClu = new AliPHOSClusterizerv1() ; - -// //========== 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 ) ; - -// //========== Creates the Reconstructioner -// fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ; -// if (fDebugLevel != 0) fRec -> SetDebugReconstruction(kTRUE); -// } - -// if (fDebugLevel != 0 || -// (ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0) -// cout << "======= Analyze ======> Event " << ievent+1 << endl ; - - -// gAlice->GetEvent(ievent) ; -// gAlice->SetEvent(ievent) ; - -// if(gAlice->TreeS() == 0) gAlice->MakeTree("S"); -// fPHOS->MakeBranch("S") ; - -// fPHOS->Hits2SDigits() ; - -// if(gAlice->TreeD() == 0) gAlice->MakeTree("D"); -// fPHOS->MakeBranch("D") ; + AliPHOSLoader* gime = dynamic_cast(fRunLoader->GetLoader("PHOSLoader")); + if ( gime == 0 ) + { + AliError(Form("Could not obtain the Loader object !")); + return ; + } -// fPHOS->SDigits2Digits() ; -// if(gAlice->TreeR() == 0) gAlice->MakeTree("R"); - -// fPHOS->Reconstruction(fRec); - -// gAlice->TreeS()->Fill() ; -// gAlice->TreeS()->Write(0,TObject::kOverwrite); - -// gAlice->TreeD()->Fill() ; -// gAlice->TreeD()->Write(0,TObject::kOverwrite); - -// } + Int_t ibranch; + TObjArray * branches; -// if(fClu) {delete fClu ; fClu =0 ;} -// if(fPID) {delete fPID ; fPID =0 ;} -// if(fRec) {delete fRec ; fRec =0 ;} -// if(fTrs) {delete fTrs ; fTrs =0 ;} + 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() ; -// } - -//------------------------------------------------------------------------------------- -void AliPHOSAnalyze::ReadAndPrintCPV(Int_t EvFirst, Int_t EvLast) -{ -// // -// // Read and print generated and reconstructed hits in CPV -// // for events from EvFirst to Nevent. -// // If only EvFirst is defined, print only this one event. -// // Author: Yuri Kharlov -// // 12 October 2000 -// // - -// if (EvFirst!=0 && EvLast==0) EvLast=EvFirst; -// for ( Int_t ievent=EvFirst; ievent<=EvLast; ievent++) { - -// //========== Event Number> -// cout << endl << "==== ReadAndPrintCPV ====> Event is " << ievent+1 << endl ; - -// //=========== Connects the various Tree's for evt -// Int_t ntracks = gAlice->GetEvent(ievent); - -// //========== Creating branches =================================== -// AliPHOSRecPoint::RecPointsList ** emcRecPoints = fPHOS->EmcRecPoints() ; -// gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP" , emcRecPoints ) ; - -// AliPHOSRecPoint::RecPointsList ** cpvRecPoints = fPHOS->PpsdRecPoints() ; -// gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", cpvRecPoints ) ; - -// // Read and print CPV hits - -// AliPHOSCPVModule cpvModule; -// TClonesArray *cpvHits; -// Int_t nCPVhits; -// AliPHOSCPVHit *cpvHit; -// TLorentzVector p; -// Float_t xgen, zgen; -// Int_t ipart; -// Int_t nGenHits = 0; -// for (Int_t itrack=0; itrackResetHits(); -// gAlice->TreeH()->GetEvent(itrack); -// Int_t iModule = 0 ; -// for (iModule=0; iModule < fGeom->GetNCPVModules(); iModule++) { -// cpvModule = fPHOS->GetCPVModule(iModule); -// cpvHits = cpvModule.Hits(); -// nCPVhits = cpvHits->GetEntriesFast(); -// for (Int_t ihit=0; ihitUncheckedAt(ihit); -// p = cpvHit->GetMomentum(); -// xgen = cpvHit->X(); -// zgen = cpvHit->Y(); -// ipart = cpvHit->GetIpart(); -// printf("CPV hit in module %d: ",iModule+1); -// printf(" p = (%f, %f, %f, %f) GeV,\n", -// p.Px(),p.Py(),p.Pz(),p.Energy()); -// printf(" (X,Z) = (%8.4f, %8.4f) cm, ipart = %d\n", -// xgen,zgen,ipart); -// } -// } -// } - -// // Read and print CPV reconstructed points - -// //=========== Gets the Reconstruction TTree -// gAlice->TreeR()->GetEvent(0) ; -// printf("Recpoints: %d\n",(*fPHOS->CpvRecPoints())->GetEntries()); -// TIter nextRP(*fPHOS->CpvRecPoints() ) ; -// AliPHOSCpvRecPoint *cpvRecPoint ; -// Int_t nRecPoints = 0; -// while( ( cpvRecPoint = (AliPHOSCpvRecPoint *)nextRP() ) ) { -// nRecPoints++; -// TVector3 locpos; -// cpvRecPoint->GetLocalPosition(locpos); -// Int_t phosModule = cpvRecPoint->GetPHOSMod(); -// printf("CPV recpoint in module %d: (X,Z) = (%f,%f) cm\n", -// phosModule,locpos.X(),locpos.Z()); -// } -// printf("This event has %d generated hits and %d reconstructed points\n", -// nGenHits,nRecPoints); -// } + 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::AnalyzeCPV(Int_t Nevents) + void AliPHOSAnalyze::InvariantMass() { -// // -// // Analyzes CPV characteristics -// // Author: Yuri Kharlov -// // 9 October 2000 -// // - -// // Book histograms - -// TH1F *hDx = new TH1F("hDx" ,"CPV x-resolution@reconstruction",100,-5. , 5.); -// TH1F *hDz = new TH1F("hDz" ,"CPV z-resolution@reconstruction",100,-5. , 5.); -// TH1F *hDr = new TH1F("hDr" ,"CPV r-resolution@reconstruction",100, 0. , 5.); -// TH1S *hNrp = new TH1S("hNrp" ,"CPV rec.point multiplicity", 21,-0.5,20.5); -// TH1S *hNrpX = new TH1S("hNrpX","CPV rec.point Phi-length" , 21,-0.5,20.5); -// TH1S *hNrpZ = new TH1S("hNrpZ","CPV rec.point Z-length" , 21,-0.5,20.5); - -// cout << "Start CPV Analysis"<< endl ; -// for ( Int_t ievent=0; ievent -// // if ( (ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0) -// cout << endl << "==== AnalyzeCPV ====> Event is " << ievent+1 << endl ; - -// //=========== Connects the various Tree's for evt -// Int_t ntracks = gAlice->GetEvent(ievent); - -// //========== Creating branches =================================== -// AliPHOSRecPoint::RecPointsList ** emcRecPoints = fPHOS->EmcRecPoints() ; -// gAlice->TreeR()->SetBranchAddress( "PHOSEmcRP" , emcRecPoints ) ; - -// AliPHOSRecPoint::RecPointsList ** cpvRecPoints = fPHOS->PpsdRecPoints() ; -// gAlice->TreeR()->SetBranchAddress( "PHOSPpsdRP", cpvRecPoints ) ; - -// // Create and fill arrays of hits for each CPV module - -// Int_t nOfModules = fGeom->GetNModules(); -// TClonesArray **hitsPerModule = new TClonesArray *[nOfModules]; -// Int_t iModule = 0; -// for (iModule=0; iModule < nOfModules; iModule++) -// hitsPerModule[iModule] = new TClonesArray("AliPHOSCPVHit",100); - -// AliPHOSCPVModule cpvModule; -// TClonesArray *cpvHits; -// Int_t nCPVhits; -// AliPHOSCPVHit *cpvHit; -// TLorentzVector p; -// Float_t xzgen[2]; -// Int_t ipart; - -// // First go through all primary tracks and fill the arrays -// // of hits per each CPV module - -// for (Int_t itrack=0; itrackResetHits(); -// gAlice->TreeH()->GetEvent(itrack); -// for (Int_t iModule=0; iModule < nOfModules; iModule++) { -// cpvModule = fPHOS->GetCPVModule(iModule); -// cpvHits = cpvModule.Hits(); -// nCPVhits = cpvHits->GetEntriesFast(); -// for (Int_t ihit=0; ihitUncheckedAt(ihit); -// p = cpvHit->GetMomentum(); -// xzgen[0] = cpvHit->X(); -// xzgen[1] = cpvHit->Y(); -// ipart = cpvHit->GetIpart(); -// TClonesArray &lhits = *(TClonesArray *)hitsPerModule[iModule]; -// new(lhits[hitsPerModule[iModule]->GetEntriesFast()]) AliPHOSCPVHit(*cpvHit); -// } -// cpvModule.Clear(); -// } -// } -// for (iModule=0; iModule < nOfModules; iModule++) { -// Int_t nsum = hitsPerModule[iModule]->GetEntriesFast(); -// printf("Module %d has %d hits\n",iModule,nsum); -// } - -// // Then go through reconstructed points and for each find -// // the closeset hit -// // The distance from the rec.point to the closest hit -// // gives the coordinate resolution of the CPV - -// // Get the Reconstruction Tree -// gAlice->TreeR()->GetEvent(0) ; -// TIter nextRP(*fPHOS->PpsdRecPoints() ) ; -// AliPHOSCpvRecPoint *cpvRecPoint ; -// Float_t xgen, zgen; -// while( ( cpvRecPoint = (AliPHOSCpvRecPoint *)nextRP() ) ) { -// TVector3 locpos; -// cpvRecPoint->GetLocalPosition(locpos); -// Int_t phosModule = cpvRecPoint->GetPHOSMod(); -// Int_t rpMult = cpvRecPoint->GetDigitsMultiplicity(); -// Int_t rpMultX, rpMultZ; -// cpvRecPoint->GetClusterLengths(rpMultX,rpMultZ); -// Float_t xrec = locpos.X(); -// Float_t zrec = locpos.Z(); -// Float_t dxmin = 1.e+10; -// Float_t dzmin = 1.e+10; -// Float_t r2min = 1.e+10; -// Float_t r2; - -// cpvHits = hitsPerModule[phosModule-1]; -// Int_t nCPVhits = cpvHits->GetEntriesFast(); -// for (Int_t ihit=0; ihitUncheckedAt(ihit); -// xgen = cpvHit->X(); -// zgen = cpvHit->Y(); -// r2 = TMath::Power((xgen-xrec),2) + TMath::Power((zgen-zrec),2); -// if ( r2 < r2min ) { -// r2min = r2; -// dxmin = xgen - xrec; -// dzmin = zgen - zrec; -// } -// } -// hDx ->Fill(dxmin); -// hDz ->Fill(dzmin); -// hDr ->Fill(TMath::Sqrt(r2min)); -// hNrp ->Fill(rpMult); -// hNrpX->Fill(rpMultX); -// hNrpZ->Fill(rpMultZ); -// } -// delete [] hitsPerModule; -// } -// // Save histograms - -// Text_t outputname[80] ; -// sprintf(outputname,"%s.analyzed",fRootFile->GetName()); -// TFile output(outputname,"RECREATE"); -// output.cd(); - -// hDx ->Write() ; -// hDz ->Write() ; -// hDr ->Write() ; -// hNrp ->Write() ; -// hNrpX->Write() ; -// hNrpZ->Write() ; - -// // Plot histograms - -// TCanvas *cpvCanvas = new TCanvas("CPV","CPV analysis",20,20,800,400); -// gStyle->SetOptStat(111111); -// gStyle->SetOptFit(1); -// gStyle->SetOptDate(1); -// cpvCanvas->Divide(3,2); - -// cpvCanvas->cd(1); -// gPad->SetFillColor(10); -// hNrp->SetFillColor(16); -// hNrp->Draw(); - -// cpvCanvas->cd(2); -// gPad->SetFillColor(10); -// hNrpX->SetFillColor(16); -// hNrpX->Draw(); - -// cpvCanvas->cd(3); -// gPad->SetFillColor(10); -// hNrpZ->SetFillColor(16); -// hNrpZ->Draw(); - -// cpvCanvas->cd(4); -// gPad->SetFillColor(10); -// hDx->SetFillColor(16); -// hDx->Fit("gaus"); -// hDx->Draw(); - -// cpvCanvas->cd(5); -// gPad->SetFillColor(10); -// hDz->SetFillColor(16); -// hDz->Fit("gaus"); -// hDz->Draw(); - -// cpvCanvas->cd(6); -// gPad->SetFillColor(10); -// hDr->SetFillColor(16); -// hDr->Draw(); - -// cpvCanvas->Print("CPV.ps"); + // Calculates Real and Mixed invariant mass distributions + 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 ; + } + + 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 -//____________________________________________________________________________ - void AliPHOSAnalyze::InvariantMass(Int_t Nevents ) -{ -// // Calculates Real and Mixed invariant mass distributions + 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 ; + -// const Int_t knMixedEvents = 4 ; //# of events used for calculation of 'mixed' distribution -// Int_t mixedLoops = (Int_t )TMath::Ceil(Nevents/knMixedEvents) ; + if (fRunLoader->TreeE() == 0x0) fRunLoader->LoadHeader(); -// //========== Booking Histograms -// TH2D * hRealEM = new TH2D("hRealEM", "Real for EM particles", 250,0.,1.,40,0.,4.) ; -// TH2D * hRealPhot = new TH2D("hRealPhot", "Real for kPhoton particles", 250,0.,1.,40,0.,4.) ; -// TH2D * hMixedEM = new TH2D("hMixedEM", "Mixed for EM particles", 250,0.,1.,40,0.,4.) ; -// TH2D * hMixedPhot= new TH2D("hMixedPhot","Mixed for kPhoton particles",250,0.,1.,40,0.,4.) ; -// Int_t ievent; -// Int_t eventInMixedLoop ; + Int_t maxevent = (Int_t)fRunLoader->TreeE()->GetEntries(); + // for(event = 0; event < gime->MaxEvent(); event++ ){ -// Int_t nRecParticles[4];//knMixedEvents] ; -// AliPHOSRecParticle::RecParticlesList * allRecParticleList = new TClonesArray("AliPHOSRecParticle", knMixedEvents*1000) ; -// for(eventInMixedLoop = 0; eventInMixedLoop < mixedLoops; eventInMixedLoop++ ){ -// Int_t iRecPhot = 0 ; + for(event = 0; event < maxevent; event++ ){ + fRunLoader->GetEvent(event); //will read only TreeR -// for ( ievent=0; ievent < knMixedEvents; ievent++){ - -// Int_t absEventNumber = eventInMixedLoop*knMixedEvents + ievent ; + //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) ; + } + + Int_t mevent = event%nMixedEvents ; //event number in the "mixed" cicle + nRecParticles[mevent] = iRecPhot-1 ; + + //check, if it is time to calculate invariant mass? + if((mevent == 0) && (event +1 == maxevent)){ -// //=========== Connects the various Tree's for evt -// gAlice->GetEvent(absEventNumber); - -// //========== Creating branches =================================== -// fPHOS->SetTreeAddress() ; + // if((mevent == 0) && (event +1 == gime->MaxEvent())){ -// gAlice->TreeD()->GetEvent(0) ; -// gAlice->TreeR()->GetEvent(0) ; + //calculate invariant mass: + Int_t irp1,irp2 ; + Int_t nCurEvent = 0 ; -// TClonesArray * recParticleList = fPHOS->RecParticles() ; + 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 - -// AliPHOSRecParticle * recParticle ; -// Int_t iRecParticle ; -// for(iRecParticle = 0; iRecParticle < recParticleList->GetEntries() ;iRecParticle++ ) -// { -// recParticle = (AliPHOSRecParticle *) recParticleList->At(iRecParticle) ; -// if((recParticle->GetType() == AliPHOSFastRecParticle::kGAMMA)|| -// (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEM)){ -// new( (*allRecParticleList)[iRecPhot] ) AliPHOSRecParticle(*recParticle) ; -// iRecPhot++; -// } -// } + //Make some cleanings + for(Int_t index = 0; index < nMixedEvents; index ++) + nRecParticles[index] = 0 ; + iRecPhot = 0 ; + allRecParticleList->Clear() ; -// nRecParticles[ievent] = iRecPhot-1 ; -// } - -// //Now 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::kGAMMA)&&(rp2->GetType() == AliPHOSFastRecParticle::kGAMMA)) -// hRealPhot->Fill(invMass,pt); -// } -// else{ -// hMixedEM->Fill(invMass,pt); -// if((rp1->GetType() == AliPHOSFastRecParticle::kGAMMA)&&(rp2->GetType() == AliPHOSFastRecParticle::kGAMMA)) -// hMixedPhot->Fill(invMass,pt); -// } //real-mixed - -// } //loop over second rp -// }//loop over first rp -// allRecParticleList->Delete() ; -// } //Loop over events - -// delete allRecParticleList ; + } + } + delete allRecParticleList ; -// //writing output -// TFile output("invmass.root","RECREATE"); -// output.cd(); + //writing output + mfile->cd(); -// hRealEM->Write() ; -// hRealPhot->Write() ; -// hMixedEM->Write() ; -// hMixedPhot->Write() ; + hRealEM->Write(0,kOverwrite) ; + hRealPhot->Write(0,kOverwrite) ; + hMixedEM->Write(0,kOverwrite) ; + hMixedPhot->Write(0,kOverwrite) ; -// output.Write(); -// output.Close(); + mfile->Write(); + mfile->Close(); + delete mfile ; + delete nRecParticles; } //____________________________________________________________________________ - void AliPHOSAnalyze::ReadAndPrintEMC(Int_t EvFirst, Int_t EvLast) + void AliPHOSAnalyze::EnergyResolution() { -// // -// // Read and print generated and reconstructed hits in EMC -// // for events from EvFirst to Nevent. -// // If only EvFirst is defined, print only this one event. -// // Author: Yuri Kharlov -// // 24 November 2000 -// // - -// if (EvFirst!=0 && EvLast==0) EvLast=EvFirst; -// Int_t ievent; -// for (ievent=EvFirst; ievent<=EvLast; ievent++) { - -// //========== Event Number> -// cout << endl << "==== ReadAndPrintEMC ====> Event is " << ievent+1 << endl ; + //fills two dimentional histo: energy of primary vs. energy of reconstructed -// //=========== Connects the various Tree's for evt -// Int_t ntracks = gAlice->GetEvent(ievent); -// fPHOS->SetTreeAddress() ; - -// gAlice->TreeD()->GetEvent(0) ; -// gAlice->TreeR()->GetEvent(0) ; + TH2F * hAllEnergy = 0 ; //all reconstructed with primary photon + TH2F * hPhotEnergy= 0 ; //kGamma with primary photon + TH2F * hEMEnergy = 0 ; //electromagnetic with primary photon -// // Loop over reconstructed particles - -// TClonesArray ** recParticleList = fPHOS->RecParticles() ; -// AliPHOSRecParticle * recParticle ; -// Int_t iRecParticle ; -// Int_t *primList; -// Int_t nPrimary; -// for(iRecParticle = 0; iRecParticle < (*recParticleList)->GetEntries() ;iRecParticle++ ) { -// recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(iRecParticle) ; -// Float_t recE = recParticle->Energy(); -// primList = recParticle->GetPrimaries(nPrimary); -// Int_t moduleNumberRec ; -// Double_t recX, recZ ; -// fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ; -// printf("Rec point: module %d, (X,Z) = (%8.4f,%8.4f) cm, E = %.3f GeV, primary = %d\n", -// moduleNumberRec,recX,recZ,recE,*primList); -// } - -// // Read and print EMC hits from EMCn branches - -// AliPHOSCPVModule emcModule; -// TClonesArray *emcHits; -// Int_t nEMChits; -// AliPHOSCPVHit *emcHit; -// TLorentzVector p; -// Float_t xgen, zgen; -// Int_t ipart, primary; -// Int_t nGenHits = 0; -// for (Int_t itrack=0; itrackResetHits(); -// gAlice->TreeH()->GetEvent(itrack); -// Int_t iModule = 0 ; -// for (iModule=0; iModule < fGeom->GetNModules(); iModule++) { -// emcModule = fPHOS->GetEMCModule(iModule); -// emcHits = emcModule.Hits(); -// nEMChits = emcHits->GetEntriesFast(); -// for (Int_t ihit=0; ihitUncheckedAt(ihit); -// p = emcHit->GetMomentum(); -// xgen = emcHit->X(); -// zgen = emcHit->Y(); -// ipart = emcHit->GetIpart(); -// primary= emcHit->GetTrack(); -// printf("EMC hit A: module %d, ",iModule+1); -// printf(" p = (%f .4, %f .4, %f .4, %f .4) GeV,\n", -// p.Px(),p.Py(),p.Pz(),p.Energy()); -// printf(" (X,Z) = (%8.4f, %8.4f) cm, ipart = %d, primary = %d\n", -// xgen,zgen,ipart,primary); -// } -// } -// } - -// // // Read and print EMC hits from PHOS branch - -// // for (Int_t itrack=0; itrackResetHits(); -// // gAlice->TreeH()->GetEvent(itrack); -// // TClonesArray *hits = fPHOS->Hits(); -// // AliPHOSHit *hit ; -// // Int_t ihit; -// // for ( ihit = 0 ; ihit < hits->GetEntries() ; ihit++ ) { -// // hit = (AliPHOSHit*)hits->At(ihit) ; -// // Float_t hitXYZ[3]; -// // hitXYZ[0] = hit->X(); -// // hitXYZ[1] = hit->Y(); -// // hitXYZ[2] = hit->Z(); -// // ipart = hit->GetPid(); -// // primary = hit->GetPrimary(); -// // Int_t absId = hit->GetId(); -// // Int_t relId[4]; -// // fGeom->AbsToRelNumbering(absId, relId) ; -// // Int_t module = relId[0]; -// // if (relId[1]==0 && !(hitXYZ[0]==0 && hitXYZ[2]==0)) -// // printf("EMC hit B: module %d, (X,Z) = (%8.4f, %8.4f) cm, ipart = %d, primary = %d\n", -// // module,hitXYZ[0],hitXYZ[2],ipart,primary); -// // } -// // } - -// } -} + //opening file and reading histograms if any + TFile * efile = new TFile("energy.root","update"); -//____________________________________________________________________________ - void AliPHOSAnalyze::AnalyzeEMC(Int_t Nevents) -{ -// // -// // Read generated and reconstructed hits in EMC for Nevents events. -// // Plots the coordinate and energy resolution histograms. -// // Coordinate resolution is a difference between the reconstructed -// // coordinate and the exact coordinate on the face of the PHOS -// // Author: Yuri Kharlov -// // 27 November 2000 -// // - -// // Book histograms - -// TH1F *hDx1 = new TH1F("hDx1" ,"EMC x-resolution", 100,-5. , 5.); -// TH1F *hDz1 = new TH1F("hDz1" ,"EMC z-resolution", 100,-5. , 5.); -// TH1F *hDE1 = new TH1F("hDE1" ,"EMC E-resolution", 100,-2. , 2.); - -// TH2F *hDx2 = new TH2F("hDx2" ,"EMC x-resolution", 100, 0., 10., 100,-5. , 5.); -// TH2F *hDz2 = new TH2F("hDz2" ,"EMC z-resolution", 100, 0., 10., 100,-5. , 5.); -// TH2F *hDE2 = new TH2F("hDE2" ,"EMC E-resolution", 100, 0., 10., 100, 0. , 5.); - -// cout << "Start EMC Analysis"<< endl ; -// for (Int_t ievent=0; ievent -// if ( (ievent+1) % (Int_t)TMath::Power( 10, (Int_t)TMath::Log10(ievent+1) ) == 0) -// cout << "==== AnalyzeEMC ====> Event is " << ievent+1 << endl ; - -// //=========== Connects the various Tree's for evt -// Int_t ntracks = gAlice->GetEvent(ievent); + 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.); -// fPHOS->SetTreeAddress() ; - -// gAlice->TreeD()->GetEvent(0) ; -// gAlice->TreeR()->GetEvent(0) ; + hPhotEnergy =(TH2F*) efile->Get("hPhotEnergy") ; + if(hPhotEnergy == 0) + hPhotEnergy = new TH2F("hPhotEnergy", "Energy of kGAMMA with primary photon",100, 0., 5., 100, 0., 5.); -// // Create and fill arrays of hits for each EMC module - -// Int_t nOfModules = fGeom->GetNModules(); -// TClonesArray **hitsPerModule = new TClonesArray *[nOfModules]; -// Int_t iModule; -// for (iModule=0; iModule < nOfModules; iModule++) -// hitsPerModule[iModule] = new TClonesArray("AliPHOSCPVHit",100); - -// AliPHOSCPVModule emcModule; -// TClonesArray *emcHits; -// Int_t nEMChits; -// AliPHOSCPVHit *emcHit; - -// // First go through all primary tracks and fill the arrays -// // of hits per each EMC module - -// for (Int_t itrack=0; itrackResetHits(); -// gAlice->TreeH()->GetEvent(itrack); -// for (Int_t iModule=0; iModule < nOfModules; iModule++) { -// emcModule = fPHOS->GetEMCModule(iModule); -// emcHits = emcModule.Hits(); -// nEMChits = emcHits->GetEntriesFast(); -// for (Int_t ihit=0; ihitUncheckedAt(ihit); -// TClonesArray &lhits = *(TClonesArray *)hitsPerModule[iModule]; -// new(lhits[hitsPerModule[iModule]->GetEntriesFast()]) AliPHOSCPVHit(*emcHit); -// } -// emcModule.Clear(); -// } -// } + hEMEnergy =(TH2F*) efile->Get("hEMEnergy"); + if(hEMEnergy == 0) + hEMEnergy = new TH2F("hEMEnergy", "Energy of EM with primary photon", 100, 0., 5., 100, 0., 5.); -// // Loop over reconstructed particles - -// TClonesArray ** recParticleList = fPHOS->RecParticles() ; -// AliPHOSRecParticle * recParticle ; -// Int_t nEMCrecs = (*recParticleList)->GetEntries(); -// if (nEMCrecs == 1) { -// recParticle = (AliPHOSRecParticle *) (*recParticleList)->At(0) ; -// Float_t recE = recParticle->Energy(); -// Int_t phosModule; -// Double_t recX, recZ ; -// fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), phosModule, recX, recZ) ; - -// // for this rec.point take the hit list in the same PHOS module - -// emcHits = hitsPerModule[phosModule-1]; -// Int_t nEMChits = emcHits->GetEntriesFast(); -// if (nEMChits == 1) { -// Float_t genX, genZ, genE; -// for (Int_t ihit=0; ihitUncheckedAt(ihit); -// genX = emcHit->X(); -// genZ = emcHit->Y(); -// genE = emcHit->GetMomentum().E(); -// } -// Float_t dx = recX - genX; -// Float_t dz = recZ - genZ; -// Float_t de = recE - genE; -// hDx1 ->Fill(dx); -// hDz1 ->Fill(dz); -// hDE1 ->Fill(de); -// hDx2 ->Fill(genE,dx); -// hDz2 ->Fill(genE,dz); -// hDE2 ->Fill(genE,recE); -// } -// } -// delete [] hitsPerModule; -// } -// // Save histograms - -// Text_t outputname[80] ; -// sprintf(outputname,"%s.analyzed",fRootFile->GetName()); -// TFile output(outputname,"RECREATE"); -// output.cd(); - -// hDx1 ->Write() ; -// hDz1 ->Write() ; -// hDE1 ->Write() ; -// hDx2 ->Write() ; -// hDz2 ->Write() ; -// hDE2 ->Write() ; - -// // Plot histograms - -// TCanvas *emcCanvas = new TCanvas("EMC","EMC analysis",20,20,700,300); -// gStyle->SetOptStat(111111); -// gStyle->SetOptFit(1); -// gStyle->SetOptDate(1); -// emcCanvas->Divide(3,1); - -// emcCanvas->cd(1); -// gPad->SetFillColor(10); -// hDx1->SetFillColor(16); -// hDx1->Draw(); - -// emcCanvas->cd(2); -// gPad->SetFillColor(10); -// hDz1->SetFillColor(16); -// hDz1->Draw(); - -// emcCanvas->cd(3); -// gPad->SetFillColor(10); -// hDE1->SetFillColor(16); -// hDE1->Draw(); - -// emcCanvas->Print("EMC.ps"); -} + 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 ; + } -//____________________________________________________________________________ - void AliPHOSAnalyze::AnalyzeResolutions(Int_t Nevents ) -{ -// // analyzes Nevents events and calculate Energy and Position resolution as well as -// // probaility of correct indentifiing of the incident particle - -// //========== Booking Histograms -// cout << "AnalyzeResolutions > " << "Booking Histograms" << endl ; -// BookResolutionHistograms(); - -// Int_t counter[9][5] ; -// Int_t i1,i2,totalInd = 0 ; -// for(i1 = 0; i1<9; i1++) -// for(i2 = 0; i2<5; i2++) -// counter[i1][i2] = 0 ; + + AliPHOSGeometry * phosgeom = AliPHOSGeometry::GetInstance() ; + + Int_t ievent; + Int_t maxevent = (Int_t)fRunLoader->TreeE()->GetEntries(); + + fRunLoader->LoadKinematics("READ"); + gime->LoadTracks("READ"); -// Int_t totalPrimary = 0 ; -// Int_t totalRecPart = 0 ; -// Int_t totalRPwithPrim = 0 ; -// Int_t ievent; - -// cout << "Start Analysing"<< endl ; -// for ( ievent=0; ieventGetEvent(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")) ; + return ; + } -// //========== Event Number> -// // if ( ( log10((Float_t)(ievent+1)) - (Int_t)(log10((Float_t)(ievent+1))) ) == 0. ) -// cout << "AnalyzeResolutions > " << "Event is " << ievent << endl ; + for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast() ;iRecParticle++ ){ + recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ; -// //=========== Connects the various Tree's for evt -// gAlice->GetEvent(ievent); - -// //=========== Gets the Kine TTree -// gAlice->TreeK()->GetEvent(0) ; + //find the closest primary + Int_t moduleNumberRec ; + Double_t recX, recZ ; + phosgeom->ImpactOnEmc(vtx,recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ; -// //=========== Gets the list of Primari Particles - -// TParticle * primary ; -// Int_t iPrimary ; -// for ( iPrimary = 0 ; iPrimary < gAlice->GetNtrack() ; iPrimary++) -// { -// primary = gAlice->Particle(iPrimary) ; -// Int_t primaryType = primary->GetPdgCode() ; -// if( primaryType == 22 ) { -// Int_t moduleNumber ; -// Double_t primX, primZ ; -// fGeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ; -// if(moduleNumber){ -// fhPrimary->Fill(primary->Energy()) ; -// if(primary->Energy() > 0.3) -// totalPrimary++ ; -// } -// } -// } + Double_t minDistance = 100. ; + Int_t closestPrimary = -1 ; -// fPHOS->SetTreeAddress() ; + //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) ; -// gAlice->TreeD()->GetEvent(0) ; -// gAlice->TreeR()->GetEvent(0) ; + 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] ; + } + } + } + + //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() ) ; + } + } + } + } + + //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::PositionResolution() +{ + //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) ; -// TClonesArray * recParticleList = fPHOS->RecParticles() ; + //find the closest primary + Int_t moduleNumberRec ; + Double_t recX, recZ ; + phosgeom->ImpactOnEmc(vtx,recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ; -// AliPHOSRecParticle * recParticle ; -// Int_t iRecParticle ; -// for(iRecParticle = 0; iRecParticle < recParticleList->GetEntries() ;iRecParticle++ ) -// { -// recParticle = (AliPHOSRecParticle *) recParticleList->At(iRecParticle) ; -// fhAllRP->Fill(CorrectEnergy(recParticle->Energy())) ; - -// Int_t moduleNumberRec ; -// Double_t recX, recZ ; -// fGeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ; - -// Double_t minDistance = 100. ; -// Int_t closestPrimary = -1 ; - -// Int_t numberofprimaries ; -// Int_t * listofprimaries = recParticle->GetPrimaries(numberofprimaries) ; -// Int_t index ; -// 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 = gAlice->Particle(listofprimaries[index]) ; -// Int_t moduleNumber ; -// Double_t primX, primZ ; -// fGeom->ImpactOnEmc(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] ; -// } -// } -// } -// totalRecPart++ ; - -// if(closestPrimary >=0 ){ -// totalRPwithPrim++; - -// Int_t primaryType = gAlice->Particle(closestPrimary)->GetPdgCode() ; -// // TParticlePDG* pDGparticle = gAlice->ParticleAt(closestPrimary)->GetPDG(); -// // Double_t charge = PDGparticle->Charge() ; -// // if(charge) -// // cout <<"Primary " <At(closestPrimary))->Energy() << endl ; -// Int_t primaryCode ; -// switch(primaryType) -// { -// case 22: -// primaryCode = 0; //Photon -// fhAllEnergy ->Fill(gAlice->Particle(closestPrimary)->Energy(), recParticle->Energy()) ; -// fhAllPosition ->Fill(gAlice->Particle(closestPrimary)->Energy(), minDistance) ; -// fhAllPositionX->Fill(dXmin); -// fhAllPositionZ->Fill(dZmin); -// break; -// case 11 : -// primaryCode = 1; //Electron -// break; -// case -11 : -// primaryCode = 1; //positron -// break; -// case 321 : -// primaryCode = 4; //K+ -// break; -// case -321 : -// primaryCode = 4; //K- -// break; -// case 310 : -// primaryCode = 4; //K0s -// break; -// case 130 : -// primaryCode = 4; //K0l -// break; -// case 211 : -// primaryCode = 2; //K0l -// break; -// case -211 : -// primaryCode = 2; //K0l -// break; -// case 2212 : -// primaryCode = 2; //K0l -// break; -// case -2212 : -// primaryCode = 2; //K0l -// break; -// default: -// primaryCode = 3; //ELSE -// break; -// } - -// switch(recParticle->GetType()) -// { -// case AliPHOSFastRecParticle::kGAMMA: -// if(primaryType == 22){ -// fhPhotEnergy->Fill(gAlice->Particle(closestPrimary)->Energy(), recParticle->Energy() ) ; -// fhEMEnergy->Fill(gAlice->Particle(closestPrimary)->Energy(), recParticle->Energy() ) ; -// fhPPSDEnergy->Fill(gAlice->Particle(closestPrimary)->Energy(), recParticle->Energy() ) ; - -// fhPhotPosition->Fill(gAlice->Particle(closestPrimary)->Energy(),minDistance) ; -// fhEMPosition->Fill(gAlice->Particle(closestPrimary)->Energy(),minDistance) ; -// fhPPSDPosition->Fill(gAlice->Particle(closestPrimary)->Energy(),minDistance) ; - -// fhPhotReg->Fill(CorrectEnergy(recParticle->Energy()) ) ; -// fhPhotEM->Fill(CorrectEnergy(recParticle->Energy()) ) ; -// fhPhotPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ; - -// fhPhotPhot->Fill(CorrectEnergy(recParticle->Energy()) ) ; -// } -// if(primaryType == 2112){ //neutron -// fhNReg->Fill(CorrectEnergy(recParticle->Energy()) ) ; -// fhNEM->Fill(CorrectEnergy(recParticle->Energy()) ) ; -// fhNPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ; -// } - -// if(primaryType == -2112){ //neutron ~ -// fhNBarReg->Fill(CorrectEnergy(recParticle->Energy()) ) ; -// fhNBarEM->Fill(CorrectEnergy(recParticle->Energy()) ) ; -// fhNBarPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ; - -// } -// if(primaryCode == 2){ -// fhChargedReg->Fill(CorrectEnergy(recParticle->Energy()) ) ; -// fhChargedEM->Fill(CorrectEnergy(recParticle->Energy()) ) ; -// fhChargedPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ; -// } - -// fhAllReg->Fill(CorrectEnergy(recParticle->Energy()) ) ; -// fhAllEM->Fill(CorrectEnergy(recParticle->Energy()) ) ; -// fhAllPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ; -// fhShape->Fill(CorrectEnergy(recParticle->Energy()) ) ; -// fhVeto->Fill(CorrectEnergy(recParticle->Energy()) ) ; -// fhPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ; -// counter[0][primaryCode]++; -// break; -// case AliPHOSFastRecParticle::kELECTRON: -// if(primaryType == 22){ -// fhPhotElec->Fill(CorrectEnergy(recParticle->Energy()) ) ; -// fhEMEnergy->Fill(gAlice->Particle(closestPrimary)->Energy(), recParticle->Energy() ) ; -// fhEMPosition->Fill(gAlice->Particle(closestPrimary)->Energy(),minDistance) ; -// fhPhotEM->Fill(CorrectEnergy(recParticle->Energy()) ) ; -// fhPhotPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ; -// } -// if(primaryType == 2112){ //neutron -// fhNEM->Fill(CorrectEnergy(recParticle->Energy()) ) ; -// fhNPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ; -// } - -// if(primaryType == -2112){ //neutron ~ -// fhNBarEM->Fill(CorrectEnergy(recParticle->Energy()) ) ; -// fhNBarPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ; - -// } -// if(primaryCode == 2){ -// fhChargedEM->Fill(CorrectEnergy(recParticle->Energy()) ) ; -// fhChargedPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ; -// } - -// fhAllEM->Fill(CorrectEnergy(recParticle->Energy()) ) ; -// fhAllPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ; -// fhShape->Fill(CorrectEnergy(recParticle->Energy()) ) ; -// fhPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ; -// counter[1][primaryCode]++; -// break; -// case AliPHOSFastRecParticle::kNEUTRALHA: -// if(primaryType == 22) -// fhPhotNeuH->Fill(CorrectEnergy(recParticle->Energy()) ) ; - -// fhVeto->Fill(CorrectEnergy(recParticle->Energy()) ) ; -// counter[2][primaryCode]++; -// break ; -// case AliPHOSFastRecParticle::kNEUTRALEM: -// if(primaryType == 22){ -// fhEMEnergy->Fill(gAlice->Particle(closestPrimary)->Energy(),recParticle->Energy() ) ; -// fhEMPosition->Fill(gAlice->Particle(closestPrimary)->Energy(),minDistance ) ; - -// fhPhotNuEM->Fill(CorrectEnergy(recParticle->Energy()) ) ; -// fhPhotEM->Fill(CorrectEnergy(recParticle->Energy()) ) ; -// } -// if(primaryType == 2112) //neutron -// fhNEM->Fill(CorrectEnergy(recParticle->Energy()) ) ; - -// if(primaryType == -2112) //neutron ~ -// fhNBarEM->Fill(CorrectEnergy(recParticle->Energy()) ) ; - -// if(primaryCode == 2) -// fhChargedEM->Fill(CorrectEnergy(recParticle->Energy()) ) ; - -// fhAllEM->Fill(CorrectEnergy(recParticle->Energy()) ) ; -// fhShape->Fill(CorrectEnergy(recParticle->Energy()) ) ; -// fhVeto->Fill(CorrectEnergy(recParticle->Energy()) ) ; - -// counter[3][primaryCode]++; -// break ; -// case AliPHOSFastRecParticle::kCHARGEDHA: -// if(primaryType == 22) //photon -// fhPhotChHa->Fill(CorrectEnergy(recParticle->Energy()) ) ; - -// counter[4][primaryCode]++ ; -// break ; -// case AliPHOSFastRecParticle::kGAMMAHA: -// if(primaryType == 22){ //photon -// fhPhotGaHa->Fill(CorrectEnergy(recParticle->Energy()) ) ; -// fhPPSDEnergy->Fill(gAlice->Particle(closestPrimary)->Energy(), recParticle->Energy() ) ; -// fhPPSDPosition->Fill(gAlice->Particle(closestPrimary)->Energy(),minDistance) ; -// fhPhotPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ; -// } -// if(primaryType == 2112){ //neutron -// fhNPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ; -// } - -// if(primaryType == -2112){ //neutron ~ -// fhNBarPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ; -// } -// if(primaryCode == 2){ -// fhChargedPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ; -// } - -// fhAllPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ; -// fhVeto->Fill(CorrectEnergy(recParticle->Energy()) ) ; -// fhPPSD->Fill(CorrectEnergy(recParticle->Energy()) ) ; -// counter[5][primaryCode]++ ; -// break ; -// case AliPHOSFastRecParticle::kABSURDEM: -// counter[6][primaryCode]++ ; -// fhShape->Fill(CorrectEnergy(recParticle->Energy()) ) ; -// break; -// case AliPHOSFastRecParticle::kABSURDHA: -// counter[7][primaryCode]++ ; -// break; -// default: -// counter[8][primaryCode]++ ; -// break; -// } -// } -// } -// } // endfor -// SaveHistograms(); -// cout << "Resolutions: Analyzed " << Nevents << " event(s)" << endl ; -// cout << "Resolutions: Total primary " << totalPrimary << endl ; -// cout << "Resoluitons: Total reconstracted " << totalRecPart << endl ; -// cout << "TotalReconstructed with Primarie " << totalRPwithPrim << endl ; -// cout << " Primary: Photon Electron Ch. Hadr. Neutr. Hadr Kaons" << endl ; -// cout << " Detected as photon " << counter[0][0] << " " << counter[0][1] << " " << counter[0][2] << " " <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 = 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] ; + } + } + } -} // endfunction + //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 ) ; + } + } + } + } + + //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::BookingHistograms() -{ - // Books the histograms where the results of the analysis are stored (to be changed) - - delete fhEmcDigit ; - delete fhVetoDigit ; - delete fhConvertorDigit ; - delete fhEmcCluster ; - delete fhVetoCluster ; - delete fhConvertorCluster ; - delete fhConvertorEmc ; +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.); - 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); -} -//____________________________________________________________________________ -void AliPHOSAnalyze::BookResolutionHistograms() -{ - // Books the histograms where the results of the Resolution analysis are stored - -// if(fhAllEnergy) -// delete fhAllEnergy ; -// if(fhPhotEnergy) -// delete fhPhotEnergy ; -// if(fhEMEnergy) -// delete fhEMEnergy ; -// if(fhPPSDEnergy) -// delete fhPPSDEnergy ; - - - fhAllEnergy = new TH2F("hAllEnergy", "Energy of any RP with primary photon",100, 0., 5., 100, 0., 5.); - fhPhotEnergy = new TH2F("hPhotEnergy", "Energy of kGAMMA with primary photon",100, 0., 5., 100, 0., 5.); - fhEMEnergy = new TH2F("hEMEnergy", "Energy of EM with primary photon", 100, 0., 5., 100, 0., 5.); - fhPPSDEnergy = new TH2F("hPPSDEnergy", "Energy of PPSD with primary photon", 100, 0., 5., 100, 0., 5.); - -// if(fhAllPosition) -// delete fhAllPosition ; -// if(fhPhotPosition) -// delete fhPhotPosition ; -// if(fhEMPosition) -// delete fhEMPosition ; -// if(fhPPSDPosition) -// delete fhPPSDPosition ; - - - fhAllPosition = new TH2F("hAllPosition", "Position of any RP with primary photon",100, 0., 5., 100, 0., 5.); - fhPhotPosition = new TH2F("hPhotPosition", "Position of kGAMMA with primary photon",100, 0., 5., 100, 0., 5.); - fhEMPosition = new TH2F("hEMPosition", "Position of EM with primary photon", 100, 0., 5., 100, 0., 5.); - fhPPSDPosition = new TH2F("hPPSDPosition", "Position of PPSD with primary photon", 100, 0., 5., 100, 0., 5.); - - fhAllPositionX = new TH1F("hAllPositionX", "#Delta X of any RP with primary photon",100, -2., 2.); - fhAllPositionZ = new TH1F("hAllPositionZ", "#Delta X of any RP with primary photon",100, -2., 2.); - -// if(fhAllReg) -// delete fhAllReg ; -// if(fhPhotReg) -// delete fhPhotReg ; -// if(fhNReg) -// delete fhNReg ; -// if(fhNBarReg) -// delete fhNBarReg ; -// if(fhChargedReg) -// delete fhChargedReg ; + //Now make some initializations + + 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 ; + + + + if (fRunLoader == 0x0) + { + AliError(Form("Error Loading session")); + return; + } - fhAllReg = new TH1F("hAllReg", "All primaries registered as photon", 100, 0., 5.); - fhPhotReg = new TH1F("hPhotReg", "Photon registered as photon", 100, 0., 5.); - fhNReg = new TH1F("hNReg", "N registered as photon", 100, 0., 5.); - fhNBarReg = new TH1F("hNBarReg", "NBar registered as photon", 100, 0., 5.); - fhChargedReg= new TH1F("hChargedReg", "Charged hadron registered as photon",100, 0., 5.); + AliPHOSLoader* gime = dynamic_cast(fRunLoader->GetLoader("PHOSLoader")); + if ( gime == 0 ) + { + AliError(Form("Could not obtain the Loader object !")); + return ; + } -// if(fhAllEM) -// delete fhAllEM ; -// if(fhPhotEM) -// delete fhPhotEM ; -// if(fhNEM) -// delete fhNEM ; -// if(fhNBarEM) -// delete fhNBarEM ; -// if(fhChargedEM) -// delete fhChargedEM ; + if (fRunLoader->TreeE() == 0x0) fRunLoader->LoadHeader(); + AliPHOSGeometry * phosgeom = AliPHOSGeometry::GetInstance() ; - fhAllEM = new TH1F("hAllEM", "All primary registered as EM",100, 0., 5.); - fhPhotEM = new TH1F("hPhotEM", "Photon registered as EM", 100, 0., 5.); - fhNEM = new TH1F("hNEM", "N registered as EM", 100, 0., 5.); - fhNBarEM = new TH1F("hNBarEM", "NBar registered as EM", 100, 0., 5.); - fhChargedEM= new TH1F("hChargedEM","Charged registered as EM",100, 0., 5.); - -// if(fhAllPPSD) -// delete fhAllPPSD ; -// if(fhPhotPPSD) -// delete fhPhotPPSD ; -// if(fhNPPSD) -// delete fhNPPSD ; -// if(fhNBarPPSD) -// delete fhNBarPPSD ; -// if(fhChargedPPSD) -// delete fhChargedPPSD ; + Int_t ievent; + Int_t maxevent = (Int_t)fRunLoader->TreeE()->GetEntries() ; + for ( ievent=0; ievent < maxevent ; ievent++){ + + 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 ; + } + + + //=========== 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()) ; + + } + + } + + //========== 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())) ; + + //==========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) ; + Int_t index ; + 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 + } - fhAllPPSD = new TH1F("hAllPPSD", "All primary registered as PPSD",100, 0., 5.); - fhPhotPPSD = new TH1F("hPhotPPSD", "Photon registered as PPSD", 100, 0., 5.); - fhNPPSD = new TH1F("hNPPSD", "N registered as PPSD", 100, 0., 5.); - fhNBarPPSD = new TH1F("hNBarPPSD", "NBar registered as PPSD", 100, 0., 5.); - fhChargedPPSD= new TH1F("hChargedPPSD","Charged registered as PPSD",100, 0., 5.); -// if(fhPrimary) -// delete fhPrimary ; - fhPrimary= new TH1F("hPrimary", "hPrimary", 100, 0., 5.); - -// if(fhAllRP) -// delete fhAllRP ; -// if(fhVeto) -// delete fhVeto ; -// if(fhShape) -// delete fhShape ; -// if(fhPPSD) -// delete fhPPSD ; - - fhAllRP = new TH1F("hAllRP","All Reconstructed particles", 100, 0., 5.); - fhVeto = new TH1F("hVeto", "All uncharged particles", 100, 0., 5.); - fhShape = new TH1F("hShape","All particles with EM shaower",100, 0., 5.); - fhPPSD = new TH1F("hPPSD", "All PPSD photon particles", 100, 0., 5.); - - -// if(fhPhotPhot) -// delete fhPhotPhot ; -// if(fhPhotElec) -// delete fhPhotElec ; -// if(fhPhotNeuH) -// delete fhPhotNeuH ; -// if(fhPhotNuEM) -// delete fhPhotNuEM ; -// if(fhPhotChHa) -// delete fhPhotChHa ; -// if(fhPhotGaHa) -// delete fhPhotGaHa ; - - fhPhotPhot = new TH1F("hPhotPhot","hPhotPhot", 100, 0., 5.); //Photon registered as photon - fhPhotElec = new TH1F("hPhotElec","hPhotElec", 100, 0., 5.); //Photon registered as Electron - fhPhotNeuH = new TH1F("hPhotNeuH","hPhotNeuH", 100, 0., 5.); //Photon registered as Neutral Hadron - fhPhotNuEM = new TH1F("hPhotNuEM","hPhotNuEM", 100, 0., 5.); //Photon registered as Neutral EM - fhPhotChHa = new TH1F("hPhotChHa","hPhotChHa", 100, 0., 5.); //Photon registered as Charged Hadron - fhPhotGaHa = new TH1F("hPhotGaHa","hPhotGaHa", 100, 0., 5.); //Photon registered as Gamma-Hadron -} - -//____________________________________________________________________________ -Bool_t AliPHOSAnalyze::OpenRootFile(Text_t * name) -{ - // Open the root file named "name" + //=================== 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); - fRootFile = new TFile(name, "update") ; - return fRootFile->IsOpen() ; -} - -//____________________________________________________________________________ -void AliPHOSAnalyze::SaveHistograms() -{ - // Saves the histograms in a root file named "name.analyzed" - - Text_t outputname[80] ; - sprintf(outputname,"%s.analyzed",fRootFile->GetName()); - TFile output(outputname,"RECREATE"); - output.cd(); - - if (fhAllEnergy) - fhAllEnergy->Write() ; - if (fhPhotEnergy) - fhPhotEnergy->Write() ; - if(fhEMEnergy) - fhEMEnergy->Write() ; - if(fhPPSDEnergy) - fhPPSDEnergy->Write() ; - if(fhAllPosition) - fhAllPosition->Write() ; - if(fhAllPositionX) - fhAllPositionX->Write() ; - if(fhAllPositionZ) - fhAllPositionZ->Write() ; - if(fhPhotPosition) - fhPhotPosition->Write() ; - if(fhEMPosition) - fhEMPosition->Write() ; - if(fhPPSDPosition) - fhPPSDPosition->Write() ; - if (fhAllReg) - fhAllReg->Write() ; - if (fhPhotReg) - fhPhotReg->Write() ; - if(fhNReg) - fhNReg->Write() ; - if(fhNBarReg) - fhNBarReg->Write() ; - if(fhChargedReg) - fhChargedReg->Write() ; - if (fhAllEM) - fhAllEM->Write() ; - if (fhPhotEM) - fhPhotEM->Write() ; - if(fhNEM) - fhNEM->Write() ; - if(fhNBarEM) - fhNBarEM->Write() ; - if(fhChargedEM) - fhChargedEM->Write() ; - if (fhAllPPSD) - fhAllPPSD->Write() ; - if (fhPhotPPSD) - fhPhotPPSD->Write() ; - if(fhNPPSD) - fhNPPSD->Write() ; - if(fhNBarPPSD) - fhNBarPPSD->Write() ; - if(fhChargedPPSD) - fhChargedPPSD->Write() ; - if(fhPrimary) - fhPrimary->Write() ; - if(fhAllRP) - fhAllRP->Write() ; - if(fhVeto) - fhVeto->Write() ; - if(fhShape) - fhShape->Write() ; - if(fhPPSD) - fhPPSD->Write() ; - if(fhPhotPhot) - fhPhotPhot->Write() ; - if(fhPhotElec) - fhPhotElec->Write() ; - if(fhPhotNeuH) - fhPhotNeuH->Write() ; - if(fhPhotNuEM) - fhPhotNuEM->Write() ; - if(fhPhotNuEM) - fhPhotNuEM->Write() ; - if(fhPhotChHa) - fhPhotChHa->Write() ; - if(fhPhotGaHa) - fhPhotGaHa->Write() ; - if(fhEnergyCorrelations) - fhEnergyCorrelations->Write() ; + cfile->Write(0,kOverwrite); + cfile->Close(); + delete cfile ; + - output.Write(); - output.Close(); -} -//____________________________________________________________________________ -Float_t AliPHOSAnalyze::CorrectEnergy(Float_t ERecPart) -{ - return ERecPart/0.8783 ; -} + //print Final Table + maxevent = (Int_t)AliRunLoader::Instance()->TreeE()->GetEntries() ; -//____________________________________________________________________________ -void AliPHOSAnalyze::ResetHistograms() -{ - fhEnergyCorrelations = 0 ; //Energy correlations between Eloss in Convertor and PPSD(2) - - fhEmcDigit = 0 ; // Histo of digit energies in the Emc - fhVetoDigit = 0 ; // Histo of digit energies in the Veto - fhConvertorDigit = 0 ; // Histo of digit energies in the Convertor - fhEmcCluster = 0 ; // Histo of Cluster energies in Emc - fhVetoCluster = 0 ; // Histo of Cluster energies in Veto - fhConvertorCluster = 0 ; // Histo of Cluster energies in Convertor - fhConvertorEmc = 0 ; // 2d Convertor versus Emc energies - - fhAllEnergy = 0 ; - fhPhotEnergy = 0 ; // Total spectrum of detected photons - fhEMEnergy = 0 ; // Spectrum of detected electrons with electron primary - fhPPSDEnergy = 0 ; - fhAllPosition = 0 ; - fhAllPositionX = 0 ; - fhAllPositionZ = 0 ; - fhPhotPosition = 0 ; - fhEMPosition = 0 ; - fhPPSDPosition = 0 ; - - fhPhotReg = 0 ; - fhAllReg = 0 ; - fhNReg = 0 ; - fhNBarReg = 0 ; - fhChargedReg = 0 ; - fhPhotEM = 0 ; - fhAllEM = 0 ; - fhNEM = 0 ; - fhNBarEM = 0 ; - fhChargedEM = 0 ; - fhPhotPPSD = 0 ; - fhAllPPSD = 0 ; - fhNPPSD = 0 ; - fhNBarPPSD = 0 ; - fhChargedPPSD = 0 ; - - fhPrimary = 0 ; - - fhPhotPhot = 0 ; - fhPhotElec = 0 ; - fhPhotNeuH = 0 ; - fhPhotNuEM = 0 ; - fhPhotChHa = 0 ; - fhPhotGaHa = 0 ; + 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 )) ; }