/* History of cvs commits:
*
* $Log$
+ * Revision 1.9 2007/11/17 16:39:49 gustavo
+ * removed deleting of not owned data and deleting of histograms which are exported to the output file (MG)
+ *
+ * Revision 1.8 2007/10/29 13:48:42 gustavo
+ * Corrected coding violations
+ *
+ * Revision 1.6 2007/08/17 12:40:04 schutz
+ * New analysis classes by Gustavo Conesa
+ *
+ * Revision 1.4.4.4 2007/07/26 10:32:09 schutz
+ * new analysis classes in the the new analysis framework
+ *
*
*/
//_________________________________________________________________________
-// Class for the analysis of gamma correlations (gamma-jet,
-// gamma-hadron and isolation cut.
-// This class contains 3 main methods: one to fill lists of particles (ESDs) comming
-// from the CTS (ITS+TPC) and the calorimeters; the other one tags a candidate
-// cluster as isolated; the last one search in the
-// corresponing calorimeter for the highest energy cluster, identified it as
-// prompt photon;
+// Class for the prompt gamma analysis, isolation cut
//
// Class created from old AliPHOSGammaJet
// (see AliRoot versions previous Release 4-09)
//
//*-- Author: Gustavo Conesa (LNF-INFN)
//////////////////////////////////////////////////////////////////////////////
-
-
-// --- ROOT system ---
-
-#include <TFile.h>
+
+
+// --- ROOT system ---
#include <TParticle.h>
#include <TH2.h>
+#include <TList.h>
+#include "Riostream.h"
+#include "TROOT.h"
+// --- AliRoot system ---
#include "AliAnaGammaDirect.h"
-#include "AliESD.h"
-#include "AliESDtrack.h"
-#include "AliESDCaloCluster.h"
-#include "Riostream.h"
#include "AliLog.h"
ClassImp(AliAnaGammaDirect)
-
+
//____________________________________________________________________________
- AliAnaGammaDirect::AliAnaGammaDirect(const char *name) :
- AliAnalysisTask(name,""), fChain(0), fESD(0),
- fOutputContainer(new TObjArray(100)),
- fPrintInfo(0), fMinGammaPt(0.),
- fCalorimeter(""), fEMCALPID(0),fPHOSPID(0),
+ AliAnaGammaDirect::AliAnaGammaDirect() :
+ TObject(),
+ fMinGammaPt(0.),
fConeSize(0.),fPtThreshold(0.),fPtSumThreshold(0),
- fNCones(0),fNPtThres(0),fMakeICMethod(0)
+ fICMethod(0), fAnaMC(0), fIsolatePi0(0),
+ fhNGamma(0),fhPhiGamma(0),fhEtaGamma(0), fhConeSumPt(0),
+ fntuplePrompt(0),
+ //kSeveralIC
+ fNCones(0),fNPtThres(0), fConeSizes(), fPtThresholds(),
+ fhPtThresIsolated(), fhPtSumIsolated(), fntSeveralIC()
{
- //Ctor
- TList * list = gDirectory->GetListOfKeys() ;
- TIter next(list) ;
- TH2F * h = 0 ;
- Int_t index ;
- for (index = 0 ; index < list->GetSize()-1 ; index++) {
- //-1 to avoid GammaJet Task
- h = dynamic_cast<TH2F*>(gDirectory->Get(list->At(index)->GetName())) ;
- fOutputContainer->Add(h) ;
- }
-
- for(Int_t i = 0; i < 10 ; i++){
- fConeSizes[i]=0;
- fPtThresholds[i]=0;
- }
-
- // Input slot #0 works with an Ntuple
- DefineInput(0, TChain::Class());
- // Output slot #0 writes into a TH1 container
- DefineOutput(0, TObjArray::Class()) ;
+ //default ctor
-}
+ //Initialize parameters
+ InitParameters();
+}
//____________________________________________________________________________
AliAnaGammaDirect::AliAnaGammaDirect(const AliAnaGammaDirect & g) :
- AliAnalysisTask(g), fChain(g.fChain), fESD(g.fESD),
- fOutputContainer(g. fOutputContainer), fPrintInfo(g.fPrintInfo),
- fMinGammaPt(g.fMinGammaPt), fCalorimeter(g.fCalorimeter),
- fEMCALPID(g.fEMCALPID),fPHOSPID(g.fPHOSPID),
+ TObject(g),
+ fMinGammaPt(g.fMinGammaPt),
fConeSize(g.fConeSize),
fPtThreshold(g.fPtThreshold),
fPtSumThreshold(g.fPtSumThreshold),
- fNCones(g.fNCones),fNPtThres(g.fNPtThres),
- fMakeICMethod(g.fMakeICMethod)
+ fICMethod(g.fICMethod),
+ fAnaMC( g.fAnaMC),
+ fIsolatePi0(g.fIsolatePi0),
+ fhNGamma(g.fhNGamma),fhPhiGamma(g.fhPhiGamma),
+ fhEtaGamma(g.fhEtaGamma), fhConeSumPt(g.fhConeSumPt),
+ fntuplePrompt(g.fntuplePrompt),
+ //kSeveralIC
+ fNCones(g.fNCones),fNPtThres(g.fNPtThres), fConeSizes(),fPtThresholds(),
+ fhPtThresIsolated(), fhPtSumIsolated()
{
// cpy ctor
- SetName (g.GetName()) ;
- SetTitle(g.GetTitle()) ;
-
- for(Int_t i = 0; i < 10 ; i++){
- fConeSizes[i]= g.fConeSizes[i];
- fPtThresholds[i]= g.fPtThresholds[i];
+
+ //kSeveralIC
+ for(Int_t i = 0; i < fNCones ; i++){
+ fConeSizes[i] = g.fConeSizes[i];
+ fhPtSumIsolated[i] = g.fhPtSumIsolated[i];
+ fntSeveralIC[i]= g.fntSeveralIC[i];
+ for(Int_t j = 0; j < fNPtThres ; j++)
+ fhPtThresIsolated[i][j] = g.fhPtThresIsolated[i][j];
}
-}
+
+ for(Int_t i = 0; i < fNPtThres ; i++)
+ fPtThresholds[i]= g.fPtThresholds[i];
-//____________________________________________________________________________
-AliAnaGammaDirect::~AliAnaGammaDirect()
-{
- // Remove all pointers
- fOutputContainer->Clear() ;
- delete fOutputContainer ;
- delete fhNGamma ;
- delete fhPhiGamma ;
- delete fhEtaGamma ;
- delete fhPtCandidate ;
- delete [] fhPtThresIsolated ;
- delete [] fhPtSumIsolated ;
}
-//____________________________________________________________________________
-void AliAnaGammaDirect::CreateParticleList(TClonesArray * pl,
- TClonesArray * plCTS,
- TClonesArray * plEMCAL,
- TClonesArray * plPHOS){
-
- //Create a list of particles from the ESD. These particles have been measured
- //by the Central Tracking system (TPC+ITS), PHOS and EMCAL
-
- Int_t index = pl->GetEntries() ;
- Int_t npar = 0 ;
- Float_t *pid = new Float_t[AliPID::kSPECIESN];
- AliDebug(3,"Fill particle lists");
- if(fPrintInfo)
- AliInfo(Form("fCalorimeter %s",fCalorimeter.Data()));
-
- Double_t v[3] ; //vertex ;
- fESD->GetVertex()->GetXYZ(v) ;
-
- //########### PHOS ##############
- if(fCalorimeter == "PHOS"){
-
- Int_t begphos = fESD->GetFirstPHOSCluster();
- Int_t endphos = fESD->GetFirstPHOSCluster() +
- fESD->GetNumberOfPHOSClusters() ;
- Int_t indexNePHOS = plPHOS->GetEntries() ;
- AliDebug(3,Form("First PHOS particle %d, last particle %d", begphos,endphos));
-
- if(fCalorimeter == "PHOS"){
- for (npar = begphos; npar < endphos; npar++) {//////////////PHOS track loop
- AliESDCaloCluster * clus = fESD->GetCaloCluster(npar) ; // retrieve track from esd
-
- //Create a TParticle to fill the particle list
-
- Float_t en = clus->GetClusterEnergy() ;
- Float_t *p = new Float_t();
- clus->GetGlobalPosition(p) ;
- TVector3 pos(p[0],p[1],p[2]) ;
- Double_t phi = pos.Phi();
- Double_t theta= pos.Theta();
- Double_t px = en*TMath::Cos(phi)*TMath::Sin(theta);;
- Double_t py = en*TMath::Sin(phi)*TMath::Sin(theta);
- Double_t pz = en*TMath::Cos(theta);
-
- TParticle * particle = new TParticle() ;
- particle->SetMomentum(px,py,pz,en) ;
- AliDebug(4,Form("PHOS clusters: pt %f, phi %f, eta %f", particle->Pt(),particle->Phi(),particle->Eta()));
-
- //Select only photons
-
- pid=clus->GetPid();
- //cout<<"pid "<<pid[AliPID::kPhoton]<<endl ;
- if( !fPHOSPID)
- new((*plPHOS)[indexNePHOS++]) TParticle(*particle) ;
- else if( pid[AliPID::kPhoton] > 0.75)
- new((*plPHOS)[indexNePHOS++]) TParticle(*particle) ;
- }
- }
- }
-
- //########### CTS (TPC+ITS) #####################
- Int_t begtpc = 0 ;
- Int_t endtpc = fESD->GetNumberOfTracks() ;
- Int_t indexCh = plCTS->GetEntries() ;
- AliDebug(3,Form("First CTS particle %d, last particle %d", begtpc,endtpc));
+//_________________________________________________________________________
+AliAnaGammaDirect & AliAnaGammaDirect::operator = (const AliAnaGammaDirect & source)
+{
+ // assignment operator
- for (npar = begtpc; npar < endtpc; npar++) {////////////// track loop
- AliESDtrack * track = fESD->GetTrack(npar) ; // retrieve track from esd
- //We want tracks fitted in the detectors:
- ULong_t status=AliESDtrack::kTPCrefit;
- status|=AliESDtrack::kITSrefit;
+ if(&source == this) return *this;
+
+ fMinGammaPt = source.fMinGammaPt ;
+ fConeSize = source.fConeSize ;
+ fPtThreshold = source.fPtThreshold ;
+ fPtSumThreshold = source.fPtSumThreshold ;
+ fICMethod = source.fICMethod ;
+ fAnaMC = source.fAnaMC ;
+ fIsolatePi0 = source.fIsolatePi0 ;
+
+ fhNGamma = source.fhNGamma ;
+ fhPhiGamma = source.fhPhiGamma ;
+ fhEtaGamma = source.fhEtaGamma ;
+ fhConeSumPt = source.fhConeSumPt ;
+
+ fntuplePrompt = source.fntuplePrompt ;
+
+ //kSeveralIC
+ fNCones = source.fNCones ;
+ fNPtThres = source.fNPtThres ;
- //We want tracks whose PID bit is set:
- // ULong_t status =AliESDtrack::kITSpid;
- // status|=AliESDtrack::kTPCpid;
-
- if ( (track->GetStatus() & status) == status) {//Check if the bits we want are set
- // Do something with the tracks which were successfully
- // re-fitted
- Double_t en = 0; //track ->GetTPCsignal() ;
- Double_t mom[3];
- track->GetPxPyPz(mom) ;
- Double_t px = mom[0];
- Double_t py = mom[1];
- Double_t pz = mom[2]; //Check with TPC people if this is correct.
- Int_t pdg = 11; //Give any charged PDG code, in this case electron.
- //I just want to tag the particle as charged
- TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1,
- px, py, pz, en, v[0], v[1], v[2], 0);
-
- //TParticle * particle = new TParticle() ;
- //particle->SetMomentum(px,py,pz,en) ;
-
- new((*plCTS)[indexCh++]) TParticle(*particle) ;
- new((*pl)[index++]) TParticle(*particle) ;
- }
+ for(Int_t i = 0; i < fNCones ; i++){
+ fConeSizes[i] = source.fConeSizes[i];
+ fhPtSumIsolated[i] = source.fhPtSumIsolated[i] ;
+ fntSeveralIC[i]= source.fntSeveralIC[i];
+ for(Int_t j = 0; j < fNPtThres ; j++)
+ fhPtThresIsolated[i][j] = source.fhPtThresIsolated[i][j] ;
}
- //################ EMCAL ##############
-
- Int_t begem = fESD->GetFirstEMCALCluster();
- Int_t endem = fESD->GetFirstEMCALCluster() +
- fESD->GetNumberOfEMCALClusters() ;
- Int_t indexNe = plEMCAL->GetEntries() ;
-
- AliDebug(3,Form("First EMCAL particle %d, last particle %d",begem,endem));
-
- for (npar = begem; npar < endem; npar++) {//////////////EMCAL track loop
- AliESDCaloCluster * clus = fESD->GetCaloCluster(npar) ; // retrieve track from esd
- Int_t clustertype= clus->GetClusterType();
- if(clustertype == AliESDCaloCluster::kClusterv1){
- Float_t en = clus->GetClusterEnergy() ;
- Float_t *p = new Float_t();
- clus->GetGlobalPosition(p) ;
- TVector3 pos(p[0],p[1],p[2]) ;
- Double_t phi = pos.Phi();
- Double_t theta= pos.Theta();
- Double_t px = en*TMath::Cos(phi)*TMath::Sin(theta);;
- Double_t py = en*TMath::Sin(phi)*TMath::Sin(theta);
- Double_t pz = en*TMath::Cos(theta);
-
- pid=clus->GetPid();
- if(fCalorimeter == "EMCAL")
- {
- TParticle * particle = new TParticle() ;
- particle->SetMomentum(px,py,pz,en) ;
- AliDebug(4,Form("EMCAL clusters: pt %f, phi %f, eta %f", particle->Pt(),particle->Phi(),particle->Eta()));
- if(!fEMCALPID) //Only identified particles
- new((*plEMCAL)[indexNe++]) TParticle(*particle) ;
- else if(pid[AliPID::kPhoton] > 0.75)
- new((*plEMCAL)[indexNe++]) TParticle(*particle) ;
- }
- else
- {
- Int_t pdg = 0;
- if(fEMCALPID)
- {
- if( pid[AliPID::kPhoton] > 0.75) //This has to be fixen.
- pdg = 22;
- else if( pid[AliPID::kPi0] > 0.75)
- pdg = 111;
- }
- else
- pdg = 22; //No PID, assume all photons
-
- TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1,
- px, py, pz, en, v[0], v[1], v[2], 0);
- AliDebug(4,Form("EMCAL clusters: pt %f, phi %f, eta %f", particle->Pt(),particle->Phi(),particle->Eta()));
-
- new((*plEMCAL)[indexNe++]) TParticle(*particle) ;
- new((*pl)[index++]) TParticle(*particle) ;
- }
- }
- }
+ for(Int_t i = 0; i < fNPtThres ; i++)
+ fPtThresholds[i]= source.fPtThresholds[i];
- AliDebug(3,"Particle lists filled");
-
+ return *this;
+
}
-
-
//____________________________________________________________________________
-void AliAnaGammaDirect::Exec(Option_t *)
+AliAnaGammaDirect::~AliAnaGammaDirect()
{
-
- // Processing of one event
-
- //Get ESDs
- Long64_t entry = fChain->GetReadEntry() ;
-
- if (!fESD) {
- AliError("fESD is not connected to the input!") ;
- return ;
- }
-
- if (fPrintInfo)
- AliInfo(Form("%s ----> Processing event # %lld", (dynamic_cast<TChain *>(fChain))->GetFile()->GetName(), entry)) ;
+ // Remove all pointers except analysis output pointers.
- //CreateTLists with arrays of TParticles. Filled with particles only relevant for the analysis.
+}
- TClonesArray * particleList = new TClonesArray("TParticle",1000); // All particles refitted in CTS and detected in EMCAL (jet)
- TClonesArray * plCTS = new TClonesArray("TParticle",1000); // All particles refitted in Central Tracking System (ITS+TPC)
- TClonesArray * plNe = new TClonesArray("TParticle",1000); // All particles measured in Jet Calorimeter (EMCAL)
- TClonesArray * plPHOS = new TClonesArray("TParticle",1000); // All particles measured in PHOS as Gamma calorimeter
- TClonesArray * plEMCAL = new TClonesArray("TParticle",1000); // All particles measured in EMCAL as Gamma calorimeter
+//________________________________________________________________________
+TList * AliAnaGammaDirect::GetCreateOutputObjects()
+{
- TParticle *pGamma = new TParticle(); //It will contain the kinematics of the found prompt gamma
- TParticle *pLeading = new TParticle(); //It will contain the kinematics of the found leading particle
+ // Create histograms to be saved in output file and
+ // store them in outputContainer
+ TList * outputContainer = new TList() ;
+ outputContainer->SetName("DirectGammaHistos") ;
- //Fill lists with photons, neutral particles and charged particles
- //look for the highest energy photon in the event inside fCalorimeter
-
- AliDebug(2, "Fill particle lists, get prompt gamma");
-
- //Fill particle lists
- CreateParticleList(particleList, plCTS,plEMCAL,plPHOS);
-
- if(fCalorimeter == "PHOS")
- plNe = plPHOS;
- if(fCalorimeter == "EMCAL")
- plNe = plEMCAL;
+ //Histograms of highest gamma identified in Event
+ fhNGamma = new TH1F("NGamma","Number of #gamma over calorimeter",240,0,120);
+ fhNGamma->SetYTitle("N");
+ fhNGamma->SetXTitle("p_{T #gamma}(GeV/c)");
+ outputContainer->Add(fhNGamma) ;
+
+ fhPhiGamma = new TH2F
+ ("PhiGamma","#phi_{#gamma}",200,0,120,200,0,7);
+ fhPhiGamma->SetYTitle("#phi");
+ fhPhiGamma->SetXTitle("p_{T #gamma} (GeV/c)");
+ outputContainer->Add(fhPhiGamma) ;
+
+ fhEtaGamma = new TH2F
+ ("EtaGamma","#phi_{#gamma}",200,0,120,200,-0.8,0.8);
+ fhEtaGamma->SetYTitle("#eta");
+ fhEtaGamma->SetXTitle("p_{T #gamma} (GeV/c)");
+ outputContainer->Add(fhEtaGamma) ;
+ fhConeSumPt = new TH2F
+ ("ConePtSum","#Sigma p_{T} in cone ",200,0,120,100,0,100);
+ fhConeSumPt->SetYTitle("#Sigma p_{T}");
+ fhConeSumPt->SetXTitle("p_{T #gamma} (GeV/c)");
+ outputContainer->Add(fhConeSumPt) ;
+
+ //NTUPLE
+ fntuplePrompt = new TNtuple("ntuplePromptGamma", "Tree of prompt #gamma", "ptcluster:phicluster:etacluster:pdg:status:ptprimary:phiprimary:etaprimary:pdgprimary:statusprimary");
+ outputContainer->Add(fntuplePrompt) ;
- //_______________Analysis 1__________________________
- //Look for the prompt photon in the selected calorimeter
- if(fMakeICMethod < 3){
-
- Bool_t iIsInPHOSorEMCAL = kFALSE ; //To check if Gamma was in any calorimeter
- //Search highest energy prompt gamma in calorimeter
- GetPromptGamma(plNe, plCTS, pGamma, iIsInPHOSorEMCAL) ;
-
- AliDebug(1, Form("Is Gamma in %s? %d",fCalorimeter.Data(),iIsInPHOSorEMCAL));
-
- //If there is any photon candidate in fCalorimeter
- if(iIsInPHOSorEMCAL){
- if (fPrintInfo)
- AliInfo(Form("Prompt Gamma: pt %f, phi %f, eta %f", pGamma->Pt(),pGamma->Phi(),pGamma->Eta())) ;
+ if(fICMethod == kSeveralIC){
+ char name[128];
+ char title[128];
+ for(Int_t icone = 0; icone<fNCones; icone++){
+ sprintf(name,"PtSumIsolated_Cone_%d",icone);
+ sprintf(title,"Candidate cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
+ fhPtSumIsolated[icone] = new TH2F(name, title,240,0,120,120,0,10);
+ fhPtSumIsolated[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
+ fhPtSumIsolated[icone]->SetXTitle("p_{T} (GeV/c)");
+ outputContainer->Add(fhPtSumIsolated[icone]) ;
- }//Gamma in Calo
-
- }//Analysis 1
-
+ sprintf(name,"nt_Cone_%d",icone);
+ sprintf(title,"ntuple for cone size %d",icone);
+ fntSeveralIC[icone] = new TNtuple(name, title, "ptcand:phicand:etacand:pdg:status:ptsum:ncone0:ncone1:ncone2:ncone3:ncone4:ncone5");
+ outputContainer->Add(fntSeveralIC[icone]) ;
- //_______________Analysis 2__________________________
- //Look for the prompt photon in the selected calorimeter
- //Isolation Cut Analysis for both methods and different pt cuts and cones
- if(fMakeICMethod == 3){
-
- for(Int_t ipr = 0; ipr < plNe->GetEntries() ; ipr ++ ){
- TParticle * pCandidate = dynamic_cast<TParticle *>(plNe->At(ipr)) ;
-
- if(pCandidate->Pt() > fMinGammaPt){
-
- Bool_t icPtThres = kFALSE;
- Bool_t icPtSum = kFALSE;
-
- Float_t ptC = pCandidate->Pt() ;
-
- fhPtCandidate->Fill(ptC);
-
- for(Int_t icone = 0; icone<fNCones; icone++){
- fConeSize = fConeSizes[icone] ;
- Float_t coneptsum = 0 ;
- for(Int_t ipt = 0; ipt<fNPtThres;ipt++){
- fPtThreshold = fPtThresholds[ipt] ;
- MakeIsolationCut(plCTS,plNe, pCandidate, ipr, icPtThres, icPtSum,coneptsum);
- AliDebug(4,Form("Candidate pt %f, pt in cone %f, Isolated? ICPt %d, ICSum %d",
- pCandidate->Pt(), coneptsum, icPtThres, icPtSum));
-
- fhPtThresIsolated[icone][ipt]->Fill(ptC);
- }//pt thresh loop
- fhPtSumIsolated[icone]->Fill(ptC,coneptsum) ;
- }//cone size loop
- }//min pt candidate
- }//candidate loop
- }//Analysis 2
-
- AliDebug(2, "End of analysis, delete pointers");
-
- particleList->Delete() ;
- plCTS->Delete() ;
- plNe->Delete() ;
- plPHOS->Delete() ;
- pLeading->Delete();
- pGamma->Delete();
+ for(Int_t ipt = 0; ipt<fNPtThres;ipt++){
+ sprintf(name,"PtThresIsol_Cone_%d_Pt%d",icone,ipt);
+ sprintf(title,"Isolated candidate p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
+ fhPtThresIsolated[icone][ipt] = new TH1F(name, title,240,0,120);
+ fhPtThresIsolated[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
+ outputContainer->Add(fhPtThresIsolated[icone][ipt]) ;
+ }//icone loop
+ }//ipt loop
+ }
- delete plNe ;
- delete plCTS ;
- delete particleList ;
- // delete pLeading;
- // delete pGamma;
+// gROOT->cd();
- PostData(0, fOutputContainer);
-}
+ return outputContainer ;
+}
//____________________________________________________________________________
-void AliAnaGammaDirect::GetPromptGamma(TClonesArray * pl, TClonesArray * plCTS, TParticle *pGamma, Bool_t &Is) const
+void AliAnaGammaDirect::GetPromptGamma(TClonesArray * plCalo, TClonesArray * plCTS, TClonesArray * plPrimCalo, TParticle *pGamma, Bool_t &found) const
{
//Search for the prompt photon in Calorimeter with pt > fMinGammaPt
Double_t pt = 0;
+ Int_t n = 0;
Int_t index = -1;
- for(Int_t ipr = 0;ipr < pl->GetEntries() ; ipr ++ ){
- TParticle * particle = dynamic_cast<TParticle *>(pl->At(ipr)) ;
+ Float_t coneptsum = 0 ;
+
+ for(Int_t ipr = 0;ipr < plCalo->GetEntries() ; ipr ++ ){
+ TParticle * particle = dynamic_cast<TParticle *>(plCalo->At(ipr)) ;
- if((particle->Pt() > fMinGammaPt) && (particle->Pt() > pt)){
+ if((particle->Pt() > fMinGammaPt) && (particle->Pt() > pt) &&
+ (particle->GetPdgCode() == 22 || (fIsolatePi0 && particle->GetPdgCode() == 111))){
index = ipr ;
pt = particle->Pt();
pGamma->SetMomentum(particle->Px(),particle->Py(),particle->Pz(),particle->Energy());
- AliDebug(4,Form("Cluster in calo: pt %f, phi %f, eta %f", pGamma->Pt(),pGamma->Phi(),pGamma->Eta())) ;
- Is = kTRUE;
+ pGamma->SetStatusCode(particle->GetStatusCode());
+ pGamma->SetPdgCode(particle->GetPdgCode());
+ found = kTRUE;
}
}
-
+
//Do Isolation?
- if(fMakeICMethod && Is)
+ if( ( fICMethod == kPtIC || fICMethod == kSumPtIC ) && found)
{
- Float_t coneptsum = 0 ;
Bool_t icPtThres = kFALSE;
Bool_t icPtSum = kFALSE;
- MakeIsolationCut(plCTS,pl, pGamma, index,
+ MakeIsolationCut(plCTS,plCalo, pGamma, index,n,
icPtThres, icPtSum,coneptsum);
- if(fMakeICMethod == 1) //Pt thres method
- Is = icPtThres ;
- if(fMakeICMethod == 2) //Pt cone sum method
- Is = icPtSum ;
+ if(fICMethod == kPtIC) //Pt thres method
+ found = icPtThres ;
+ if(fICMethod == kSumPtIC) //Pt cone sum method
+ found = icPtSum ;
}
- if(Is){
- AliDebug(3,Form("Cluster with p_{T} larger than %f found in calorimeter ", fMinGammaPt)) ;
- AliDebug(3,Form("Gamma: pt %f, phi %f, eta %f", pGamma->Pt(),pGamma->Phi(),pGamma->Eta())) ;
- //Fill prompt gamma histograms
- fhNGamma ->Fill(pGamma->Pt());
- fhPhiGamma->Fill( pGamma->Pt(),pGamma->Phi());
- fhEtaGamma->Fill(pGamma->Pt(),pGamma->Eta());
+ if(found){
+ AliDebug(1,Form("Cluster with p_{T} larger than the pt treshold %f found in calorimeter ", fMinGammaPt)) ;
+ AliDebug(1,Form("Gamma: pt %f, phi %f, eta %f", pGamma->Pt(),pGamma->Phi(),pGamma->Eta())) ;
+ //Fill prompt gamma histograms
+ Float_t ptcluster = pGamma->Pt();
+ Float_t phicluster = pGamma->Phi();
+ Float_t etacluster = pGamma->Eta();
+ Int_t statuscluster = pGamma->GetStatusCode();
+ Int_t pdgcluster = pGamma->GetPdgCode();
+
+ fhNGamma ->Fill(ptcluster);
+ fhPhiGamma ->Fill(ptcluster,phicluster);
+ fhEtaGamma ->Fill(ptcluster,etacluster);
+ fhConeSumPt->Fill(ptcluster,coneptsum);
+
+ Float_t ptprimary = 0 ;
+ Float_t phiprimary = 0 ;
+ Float_t etaprimary = 0 ;
+ Int_t pdgprimary = 0 ;
+ Int_t statusprimary = 0 ;
+
+ if(fAnaMC){
+ TParticle * primary = dynamic_cast<TParticle *>(plPrimCalo->At(index)) ;
+ ptprimary = primary->Pt();
+ phiprimary = primary->Phi();
+ etaprimary = primary->Eta();
+ pdgprimary = TMath::Abs(primary->GetPdgCode()) ;
+ statusprimary = primary->GetStatusCode(); // = 2 means decay gamma!!!
+
+ AliDebug(1, Form("Identified prompt Gamma pT %2.2f; Primary, pdg %d, pT %2.2f",ptcluster,pdgprimary,ptprimary));
+ //printf("Identified prompt Gamma pT %2.2f; Primary, pdg %d, pT %2.2f \n",ptcluster,pdgprimary,ptprimary);
+ }
+
+ //Fill ntuple with cluster / MC data
+// gROOT->cd();
+ fntuplePrompt->Fill(ptcluster,phicluster,etacluster,pdgcluster,statuscluster,ptprimary,phiprimary, etaprimary,pdgprimary,statusprimary);
}
else
AliDebug(1,Form("NO Cluster with pT larger than %f found in calorimeter ", fMinGammaPt)) ;
}
//____________________________________________________________________________
-void AliAnaGammaDirect::Init(const Option_t * )
+void AliAnaGammaDirect::InitParameters()
{
- // Initialisation of branch container
- AliDebug(2,Form("*** Initialization of %s", GetName())) ;
-
- // Get input data
- fChain = dynamic_cast<TChain *>(GetInputData(0)) ;
- if (!fChain) {
- AliError(Form("Input 0 for %s not found\n", GetName()));
- return ;
- }
-
- if (!fESD) {
- // One should first check if the branch address was taken by some other task
- char ** address = (char **)GetBranchAddress(0, "ESD") ;
- if (address)
- fESD = (AliESD *)(*address) ;
- if (!fESD)
- fChain->SetBranchAddress("ESD", &fESD) ;
- }
- // The output objects will be written to
- TDirectory * cdir = gDirectory ;
- // Open a file for output #0
- char outputName[1024] ;
- sprintf(outputName, "%s.root", GetName() ) ;
- OpenFile(0, outputName , "RECREATE") ;
- if (cdir)
- cdir->cd() ;
-
- // //Initialize the parameters of the analysis.
- fCalorimeter="PHOS";
- fPrintInfo = kTRUE;
- fMinGammaPt = 10. ;
+ //Initialize the parameters of the analysis.
+ fMinGammaPt = 5. ;
//Fill particle lists when PID is ok
- fEMCALPID = kFALSE;
- fPHOSPID = kFALSE;
-
- fConeSize = 0.2 ;
- fPtThreshold = 2.0;
+ fConeSize = 0.4 ;
+ fPtThreshold = 1.;
fPtSumThreshold = 1.;
- fNCones = 4 ;
- fNPtThres = 4 ;
- fConeSizes[0] = 0.1; fConeSizes[0] = 0.2; fConeSizes[2] = 0.3; fConeSizes[3] = 0.4;
- fPtThresholds[0]=1.; fPtThresholds[0]=2.; fPtThresholds[0]=3.; fPtThresholds[0]=4.;
+ fICMethod = kNoIC; // 0 don't isolate, 1 pt thresh method, 2 cone pt sum method
+ fAnaMC = kFALSE ;
+ fIsolatePi0 = kFALSE ;
+ //-----------kSeveralIC-----------------
+ fNCones = 5 ;
+ fNPtThres = 6 ;
+ fConeSizes[0] = 0.1; fConeSizes[1] = 0.2; fConeSizes[2] = 0.3; fConeSizes[3] = 0.4; fConeSizes[4] = 0.5;
+ fPtThresholds[0]=0.; fPtThresholds[1]=1.; fPtThresholds[2]=2.; fPtThresholds[3]=3.; fPtThresholds[4]=4.;fPtThresholds[5]=5.;
- fMakeICMethod = 1; // 0 don't isolate, 1 pt thresh method, 2 cone pt sum method, 3 make isolation study
-
- //Initialization of histograms
- MakeHistos() ;
}
//__________________________________________________________________
void AliAnaGammaDirect::MakeIsolationCut(TClonesArray * plCTS,
- TClonesArray * plNe,
- TParticle * pCandidate,
- Int_t index,
- Bool_t &icmpt, Bool_t &icms,
- Float_t &coneptsum) const
+ TClonesArray * plNe,
+ TParticle * pCandidate,
+ Int_t index, Int_t & n,
+ Bool_t &icmpt, Bool_t &icms,
+ Float_t &coneptsum) const
{
//Search in cone around a candidate particle if it is isolated
Float_t phiC = pCandidate->Phi() ;
Float_t eta = -100. ;
Float_t phi = -100. ;
Float_t rad = -100 ;
- Int_t n = 0 ;
- TParticle * particle = new TParticle();
+ TParticle * particle = new TParticle;
- coneptsum = 0;
+ n = 0 ;
+ coneptsum = 0.;
icmpt = kFALSE;
icms = kFALSE;
phi = particle->Phi() ;
//Check if there is any particle inside cone with pt larger than fPtThreshold
- rad = TMath::Sqrt(TMath::Power((eta-etaC),2) +
- TMath::Power((phi-phiC),2));
- if(rad<fConeSize){
+ rad = TMath::Sqrt((eta-etaC)*(eta-etaC)+ (phi-phiC)*(phi-phiC));
+ if(rad < fConeSize){
AliDebug(3,Form("charged in cone pt %f, phi %f, eta %f, R %f ",pt,phi,eta,rad));
coneptsum+=pt;
if(pt > fPtThreshold ) n++;
phi = particle->Phi() ;
//Check if there is any particle inside cone with pt larger than fPtThreshold
- rad = TMath::Sqrt(TMath::Power((eta-etaC),2) +
- TMath::Power((phi-phiC),2));
- if(rad<fConeSize){
+ rad = TMath::Sqrt((eta-etaC)*(eta-etaC)+ (phi-phiC)*(phi-phiC));
+ if(rad < fConeSize){
AliDebug(3,Form("charged in cone pt %f, phi %f, eta %f, R %f ",pt,phi,eta,rad));
coneptsum+=pt;
if(pt > fPtThreshold ) n++;
if(n == 0)
icmpt = kTRUE ;
- if(coneptsum < fPtSumThreshold)
+ if(coneptsum <= fPtSumThreshold)
icms = kTRUE ;
}
-//___________________________________________________________________
-void AliAnaGammaDirect::MakeHistos()
+//__________________________________________________________________
+void AliAnaGammaDirect::MakeSeveralICAnalysis(TClonesArray * plCalo, TClonesArray * plCTS)
{
- // Create histograms to be saved in output file and
- // stores them in fOutputContainer
-
- fOutputContainer = new TObjArray(10000) ;
+ //Isolation Cut Analysis for both methods and different pt cuts and cones
+ if (fICMethod != kSeveralIC)
+ AliFatal("Remember to set in config file: directGamma->SetICMethod(kSeveralIC)");
- //Histograms of highest gamma identified in Event
- fhNGamma = new TH1F("NGamma","Number of #gamma over PHOS",240,0,120);
- fhNGamma->SetYTitle("N");
- fhNGamma->SetXTitle("p_{T #gamma}(GeV/c)");
- fOutputContainer->Add(fhNGamma) ;
+ //Search maximum energy photon in the event
+ Double_t ptC = 0;
+ Int_t index = -1;
+ TParticle * pCandidate = new TParticle();
+ Bool_t found = kFALSE;
- fhPhiGamma = new TH2F
- ("PhiGamma","#phi_{#gamma}",200,0,120,200,0,7);
- fhPhiGamma->SetYTitle("#phi");
- fhPhiGamma->SetXTitle("p_{T #gamma} (GeV/c)");
- fOutputContainer->Add(fhPhiGamma) ;
+ for(Int_t ipr = 0;ipr < plCalo->GetEntries() ; ipr ++ ){
+ TParticle * particle = dynamic_cast<TParticle *>(plCalo->At(ipr)) ;
+
+ if((particle->Pt() > fMinGammaPt) && (particle->Pt() > ptC) &&
+ (particle->GetPdgCode() == 22 || (fIsolatePi0 && particle->GetPdgCode() == 111))){
+ index = ipr ;
+ ptC = particle->Pt();
+ pCandidate = particle ;
+ found = kTRUE;
+ }
+ }
- fhEtaGamma = new TH2F
- ("EtaGamma","#phi_{#gamma}",200,0,120,200,-0.8,0.8);
- fhEtaGamma->SetYTitle("#eta");
- fhEtaGamma->SetXTitle("p_{T #gamma} (GeV/c)");
- fOutputContainer->Add(fhEtaGamma) ;
-
- if( fMakeICMethod== 3 ){
-
- //Isolation cut histograms
- fhPtCandidate = new TH1F
- ("PtCandidate","p_{T} of candidate particles for isolation",240,0,120);
- fhPtCandidate->SetXTitle("p_{T} (GeV/c)");
- fOutputContainer->Add(fhPtCandidate) ;
+ //If there is a large cluster, larger than threshold, study isolation cut
+ if(found){
+
+ fhNGamma->Fill(ptC);
+ fhPhiGamma->Fill(ptC,pCandidate->Phi());
+ fhEtaGamma->Fill(ptC,pCandidate->Eta());
+
+ Int_t ncone[10][10];//[fNCones][fNPtThres];
+ Bool_t icPtThres = kFALSE;
+ Bool_t icPtSum = kFALSE;
- char name[128];
- char title[128];
for(Int_t icone = 0; icone<fNCones; icone++){
- sprintf(name,"PtSumIsolated_Cone_%d",icone);
- sprintf(title,"Candidate cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
- fhPtSumIsolated[icone] = new TH2F(name, title,240,0,120,120,0,10);
- fhPtSumIsolated[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
- fhPtSumIsolated[icone]->SetXTitle("p_{T} (GeV/c)");
- fOutputContainer->Add(fhPtSumIsolated[icone]) ;
-
- for(Int_t ipt = 0; ipt<fNPtThres;ipt++){
- sprintf(name,"PtThresIsol_Cone_%d_Pt%d",icone,ipt);
- sprintf(title,"Isolated candidate p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
- fhPtThresIsolated[icone][ipt] = new TH1F(name, title,240,0,120);
- fhPtThresIsolated[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
- fOutputContainer->Add(fhPtThresIsolated[icone][ipt]) ;
- }//icone loop
- }//ipt loop
- }
+ fConeSize=fConeSizes[icone] ;
+ Float_t coneptsum = 0 ;
+ for(Int_t ipt = 0; ipt<fNPtThres;ipt++){
+ ncone[icone][ipt]=0;
+ fPtThreshold=fPtThresholds[ipt] ;
+ MakeIsolationCut(plCTS,plCalo, pCandidate, index,
+ ncone[icone][ipt], icPtThres, icPtSum,coneptsum);
+ AliDebug(1,Form("Candidate pt %f, pt in cone %f, Isolated? ICPt %d, ICSum %d",
+ pCandidate->Pt(), coneptsum, icPtThres, icPtSum));
+// if(ptC >15 && ptC < 25 && (icPtThres || icPtSum) && ipt ==0){
+// printf("R %0.1f, ptthres %1.1f, ptsum %1.1f, Candidate pt %2.2f, eta %2.2f, phi %2.2f, pt in cone %2.2f, Isolated? ICPt %d, ICSum %d\n",
+// fConeSize, fPtThreshold, fPtSumThreshold, ptC, pCandidate->Eta(), pCandidate->Phi()*TMath::RadToDeg(), coneptsum, icPtThres, icPtSum);
+// //cout<<"mother label "<<pCandidate->GetMother(0)<<endl;
+// }
+
+ fhPtThresIsolated[icone][ipt]->Fill(ptC);
+ }//pt thresh loop
+ fhPtSumIsolated[icone]->Fill(ptC,coneptsum) ;
+// gROOT->cd();
+ fntSeveralIC[icone]->Fill(ptC,pCandidate->Phi(),pCandidate->Eta(), pCandidate->GetPdgCode(), pCandidate->GetStatusCode(), coneptsum,ncone[icone][0],ncone[icone][1],ncone[icone][2],ncone[icone][3],ncone[icone][4],ncone[icone][5]);
+ }//cone size loop
+ }//found high energy gamma in the event
}
+//__________________________________________________________________
void AliAnaGammaDirect::Print(const Option_t * opt) const
{
-
+
//Print some relevant parameters set for the analysis
if(! opt)
return;
-
+
Info("Print", "%s %s", GetName(), GetTitle() ) ;
+
+ printf("Min Gamma pT = %f\n", fMinGammaPt) ;
+ printf("IC method = %d\n", fICMethod) ;
printf("Cone Size = %f\n", fConeSize) ;
- printf("pT threshold = %f\n", fPtThreshold) ;
- printf("pT sum threshold = %f\n", fPtSumThreshold) ;
- printf("Min Gamma pT = %f\n", fMinGammaPt) ;
- printf("Calorimeter = %s\n", fCalorimeter.Data()) ;
-}
-
-void AliAnaGammaDirect::Terminate(Option_t *)
-{
- // The Terminate() function is the last function to be called during
- // a query. It always runs on the client, it can be used to present
- // the results graphically or save the results to file.
-
+ if(fICMethod == kPtIC) printf("pT threshold = %f\n", fPtThreshold) ;
+ if(fICMethod == kSumPtIC) printf("pT sum threshold = %f\n", fPtSumThreshold) ;
+
+ if(fICMethod == kSeveralIC){
+ printf("N Cone Sizes = %d\n", fNCones) ;
+ printf("N pT thresholds = %d\n", fNPtThres) ;
+ printf("Cone Sizes = \n") ;
+ for(Int_t i = 0; i < fNCones; i++)
+ printf(" %f;", fConeSizes[i]) ;
+ printf(" \n") ;
+ for(Int_t i = 0; i < fNPtThres; i++)
+ printf(" %f;", fPtThresholds[i]) ;
+ }
-}
+ printf(" \n") ;
+
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff