/************************************************************************** * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * * * * Author: The ALICE Off-line Project. * * Contributors are mentioned in the code where appropriate. * * * * Permission to use, copy, modify and distribute this software and its * * documentation strictly for non-commercial purposes is hereby granted * * without fee, provided that the above copyright notice appears in all * * copies and that both the copyright notice and this permission notice * * appear in the supporting documentation. The authors make no claims * * about the suitability of this software for any purpose. It is * * provided "as is" without express or implied warranty. * **************************************************************************/ /* $Id$ */ //_________________________________________________________________________ // 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 "TH2.h" #include "TH2.h" #include "TParticle.h" #include "TDatabasePDG.h" #include "TClonesArray.h" #include "TMath.h" #include "TROOT.h" // --- Standard library --- // --- AliRoot header files --- #include "AliLog.h" #include "AliStack.h" #include "AliPHOSGeometry.h" #include "AliPHOSAnalyze.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(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" fRunLoader = AliRunLoader::Open(fileName,"AliPHOSAnalyze"); if (fRunLoader == 0x0) { AliError(Form("Error Loading session")); } } //____________________________________________________________________________ AliPHOSAnalyze::AliPHOSAnalyze(const AliPHOSAnalyze & ana): TObject(ana), fCorrection(0.), fEvt(0), ffileName(), fRunLoader(0) { // copy ctor ( (AliPHOSAnalyze &)ana ).Copy(*this) ; } //____________________________________________________________________________ AliPHOSAnalyze::~AliPHOSAnalyze() { // dtor } //____________________________________________________________________________ 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 //Plot Primary Particles if (fRunLoader->Stack() == 0x0) fRunLoader->LoadKinematics("READ"); 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()) ; // } // } } 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()) ; } } } } //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 ; } 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::InvariantMass() { // 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 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 //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)){ // 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::EnergyResolution() { //fills two dimentional histo: energy of primary vs. energy of reconstructed 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"); 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.); hEMEnergy =(TH2F*) efile->Get("hEMEnergy"); if(hEMEnergy == 0) hEMEnergy = new TH2F("hEMEnergy", "Energy of EM with primary photon", 100, 0., 5., 100, 0., 5.); 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"); for ( ievent=0; ievent < maxevent ; ievent++){ //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 ; } for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast() ;iRecParticle++ ){ recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ; //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 ; 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) ; //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 ; 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 ) ; } } } } //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::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.); //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; } 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++){ 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 } //=================== 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 )) ; }