1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
17 /* History of cvs commits:
20 * Revision 1.9 2007/11/17 16:39:49 gustavo
21 * removed deleting of not owned data and deleting of histograms which are exported to the output file (MG)
23 * Revision 1.8 2007/10/29 13:48:42 gustavo
24 * Corrected coding violations
26 * Revision 1.6 2007/08/17 12:40:04 schutz
27 * New analysis classes by Gustavo Conesa
29 * Revision 1.4.4.4 2007/07/26 10:32:09 schutz
30 * new analysis classes in the the new analysis framework
35 //_________________________________________________________________________
36 // Class for the prompt gamma analysis, isolation cut
38 // Class created from old AliPHOSGammaJet
39 // (see AliRoot versions previous Release 4-09)
41 //*-- Author: Gustavo Conesa (LNF-INFN)
42 //////////////////////////////////////////////////////////////////////////////
45 // --- ROOT system ---
46 #include <TParticle.h>
49 #include "Riostream.h"
52 // --- AliRoot system ---
53 #include "AliAnaGammaDirect.h"
56 ClassImp(AliAnaGammaDirect)
58 //____________________________________________________________________________
59 AliAnaGammaDirect::AliAnaGammaDirect() :
62 fConeSize(0.),fPtThreshold(0.),fPtSumThreshold(0),
63 fICMethod(0), fAnaMC(0), fIsolatePi0(0),
64 fhNGamma(0),fhPhiGamma(0),fhEtaGamma(0), fhConeSumPt(0),
67 fNCones(0),fNPtThres(0), fConeSizes(), fPtThresholds(),
68 fhPtThresIsolated(), fhPtSumIsolated(), fntSeveralIC()
72 //Initialize parameters
77 //____________________________________________________________________________
78 AliAnaGammaDirect::AliAnaGammaDirect(const AliAnaGammaDirect & g) :
80 fMinGammaPt(g.fMinGammaPt),
81 fConeSize(g.fConeSize),
82 fPtThreshold(g.fPtThreshold),
83 fPtSumThreshold(g.fPtSumThreshold),
84 fICMethod(g.fICMethod),
86 fIsolatePi0(g.fIsolatePi0),
87 fhNGamma(g.fhNGamma),fhPhiGamma(g.fhPhiGamma),
88 fhEtaGamma(g.fhEtaGamma), fhConeSumPt(g.fhConeSumPt),
89 fntuplePrompt(g.fntuplePrompt),
91 fNCones(g.fNCones),fNPtThres(g.fNPtThres), fConeSizes(),fPtThresholds(),
92 fhPtThresIsolated(), fhPtSumIsolated()
97 for(Int_t i = 0; i < fNCones ; i++){
98 fConeSizes[i] = g.fConeSizes[i];
99 fhPtSumIsolated[i] = g.fhPtSumIsolated[i];
100 fntSeveralIC[i]= g.fntSeveralIC[i];
101 for(Int_t j = 0; j < fNPtThres ; j++)
102 fhPtThresIsolated[i][j] = g.fhPtThresIsolated[i][j];
105 for(Int_t i = 0; i < fNPtThres ; i++)
106 fPtThresholds[i]= g.fPtThresholds[i];
111 //_________________________________________________________________________
112 AliAnaGammaDirect & AliAnaGammaDirect::operator = (const AliAnaGammaDirect & source)
114 // assignment operator
116 if(&source == this) return *this;
118 fMinGammaPt = source.fMinGammaPt ;
119 fConeSize = source.fConeSize ;
120 fPtThreshold = source.fPtThreshold ;
121 fPtSumThreshold = source.fPtSumThreshold ;
122 fICMethod = source.fICMethod ;
123 fAnaMC = source.fAnaMC ;
124 fIsolatePi0 = source.fIsolatePi0 ;
126 fhNGamma = source.fhNGamma ;
127 fhPhiGamma = source.fhPhiGamma ;
128 fhEtaGamma = source.fhEtaGamma ;
129 fhConeSumPt = source.fhConeSumPt ;
131 fntuplePrompt = source.fntuplePrompt ;
134 fNCones = source.fNCones ;
135 fNPtThres = source.fNPtThres ;
137 for(Int_t i = 0; i < fNCones ; i++){
138 fConeSizes[i] = source.fConeSizes[i];
139 fhPtSumIsolated[i] = source.fhPtSumIsolated[i] ;
140 fntSeveralIC[i]= source.fntSeveralIC[i];
141 for(Int_t j = 0; j < fNPtThres ; j++)
142 fhPtThresIsolated[i][j] = source.fhPtThresIsolated[i][j] ;
145 for(Int_t i = 0; i < fNPtThres ; i++)
146 fPtThresholds[i]= source.fPtThresholds[i];
152 //____________________________________________________________________________
153 AliAnaGammaDirect::~AliAnaGammaDirect()
155 // Remove all pointers except analysis output pointers.
159 //________________________________________________________________________
160 TList * AliAnaGammaDirect::GetCreateOutputObjects()
163 // Create histograms to be saved in output file and
164 // store them in outputContainer
165 TList * outputContainer = new TList() ;
166 outputContainer->SetName("DirectGammaHistos") ;
168 //Histograms of highest gamma identified in Event
169 fhNGamma = new TH1F("NGamma","Number of #gamma over calorimeter",240,0,120);
170 fhNGamma->SetYTitle("N");
171 fhNGamma->SetXTitle("p_{T #gamma}(GeV/c)");
172 outputContainer->Add(fhNGamma) ;
174 fhPhiGamma = new TH2F
175 ("PhiGamma","#phi_{#gamma}",200,0,120,200,0,7);
176 fhPhiGamma->SetYTitle("#phi");
177 fhPhiGamma->SetXTitle("p_{T #gamma} (GeV/c)");
178 outputContainer->Add(fhPhiGamma) ;
180 fhEtaGamma = new TH2F
181 ("EtaGamma","#phi_{#gamma}",200,0,120,200,-0.8,0.8);
182 fhEtaGamma->SetYTitle("#eta");
183 fhEtaGamma->SetXTitle("p_{T #gamma} (GeV/c)");
184 outputContainer->Add(fhEtaGamma) ;
186 fhConeSumPt = new TH2F
187 ("ConePtSum","#Sigma p_{T} in cone ",200,0,120,100,0,100);
188 fhConeSumPt->SetYTitle("#Sigma p_{T}");
189 fhConeSumPt->SetXTitle("p_{T #gamma} (GeV/c)");
190 outputContainer->Add(fhConeSumPt) ;
193 fntuplePrompt = new TNtuple("ntuplePromptGamma", "Tree of prompt #gamma", "ptcluster:phicluster:etacluster:pdg:status:ptprimary:phiprimary:etaprimary:pdgprimary:statusprimary");
194 outputContainer->Add(fntuplePrompt) ;
196 if(fICMethod == kSeveralIC){
199 for(Int_t icone = 0; icone<fNCones; icone++){
200 sprintf(name,"PtSumIsolated_Cone_%d",icone);
201 sprintf(title,"Candidate cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
202 fhPtSumIsolated[icone] = new TH2F(name, title,240,0,120,120,0,10);
203 fhPtSumIsolated[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
204 fhPtSumIsolated[icone]->SetXTitle("p_{T} (GeV/c)");
205 outputContainer->Add(fhPtSumIsolated[icone]) ;
207 sprintf(name,"nt_Cone_%d",icone);
208 sprintf(title,"ntuple for cone size %d",icone);
209 fntSeveralIC[icone] = new TNtuple(name, title, "ptcand:phicand:etacand:pdg:status:ptsum:ncone0:ncone1:ncone2:ncone3:ncone4:ncone5");
210 outputContainer->Add(fntSeveralIC[icone]) ;
212 for(Int_t ipt = 0; ipt<fNPtThres;ipt++){
213 sprintf(name,"PtThresIsol_Cone_%d_Pt%d",icone,ipt);
214 sprintf(title,"Isolated candidate p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
215 fhPtThresIsolated[icone][ipt] = new TH1F(name, title,240,0,120);
216 fhPtThresIsolated[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
217 outputContainer->Add(fhPtThresIsolated[icone][ipt]) ;
224 return outputContainer ;
228 //____________________________________________________________________________
229 void AliAnaGammaDirect::GetPromptGamma(TClonesArray * plCalo, TClonesArray * plCTS, TClonesArray * plPrimCalo, TParticle *pGamma, Bool_t &found) const
231 //Search for the prompt photon in Calorimeter with pt > fMinGammaPt
236 Float_t coneptsum = 0 ;
238 for(Int_t ipr = 0;ipr < plCalo->GetEntries() ; ipr ++ ){
239 TParticle * particle = dynamic_cast<TParticle *>(plCalo->At(ipr)) ;
241 if((particle->Pt() > fMinGammaPt) && (particle->Pt() > pt) &&
242 (particle->GetPdgCode() == 22 || (fIsolatePi0 && particle->GetPdgCode() == 111))){
245 pGamma->SetMomentum(particle->Px(),particle->Py(),particle->Pz(),particle->Energy());
246 pGamma->SetStatusCode(particle->GetStatusCode());
247 pGamma->SetPdgCode(particle->GetPdgCode());
253 if( ( fICMethod == kPtIC || fICMethod == kSumPtIC ) && found)
255 Bool_t icPtThres = kFALSE;
256 Bool_t icPtSum = kFALSE;
257 MakeIsolationCut(plCTS,plCalo, pGamma, index,n,
258 icPtThres, icPtSum,coneptsum);
259 if(fICMethod == kPtIC) //Pt thres method
261 if(fICMethod == kSumPtIC) //Pt cone sum method
266 AliDebug(1,Form("Cluster with p_{T} larger than the pt treshold %f found in calorimeter ", fMinGammaPt)) ;
267 AliDebug(1,Form("Gamma: pt %f, phi %f, eta %f", pGamma->Pt(),pGamma->Phi(),pGamma->Eta())) ;
268 //Fill prompt gamma histograms
269 Float_t ptcluster = pGamma->Pt();
270 Float_t phicluster = pGamma->Phi();
271 Float_t etacluster = pGamma->Eta();
272 Int_t statuscluster = pGamma->GetStatusCode();
273 Int_t pdgcluster = pGamma->GetPdgCode();
275 fhNGamma ->Fill(ptcluster);
276 fhPhiGamma ->Fill(ptcluster,phicluster);
277 fhEtaGamma ->Fill(ptcluster,etacluster);
278 fhConeSumPt->Fill(ptcluster,coneptsum);
280 Float_t ptprimary = 0 ;
281 Float_t phiprimary = 0 ;
282 Float_t etaprimary = 0 ;
283 Int_t pdgprimary = 0 ;
284 Int_t statusprimary = 0 ;
287 TParticle * primary = dynamic_cast<TParticle *>(plPrimCalo->At(index)) ;
288 ptprimary = primary->Pt();
289 phiprimary = primary->Phi();
290 etaprimary = primary->Eta();
291 pdgprimary = TMath::Abs(primary->GetPdgCode()) ;
292 statusprimary = primary->GetStatusCode(); // = 2 means decay gamma!!!
294 AliDebug(1, Form("Identified prompt Gamma pT %2.2f; Primary, pdg %d, pT %2.2f",ptcluster,pdgprimary,ptprimary));
295 //printf("Identified prompt Gamma pT %2.2f; Primary, pdg %d, pT %2.2f \n",ptcluster,pdgprimary,ptprimary);
298 //Fill ntuple with cluster / MC data
300 fntuplePrompt->Fill(ptcluster,phicluster,etacluster,pdgcluster,statuscluster,ptprimary,phiprimary, etaprimary,pdgprimary,statusprimary);
303 AliDebug(1,Form("NO Cluster with pT larger than %f found in calorimeter ", fMinGammaPt)) ;
306 //____________________________________________________________________________
307 void AliAnaGammaDirect::InitParameters()
310 //Initialize the parameters of the analysis.
313 //Fill particle lists when PID is ok
316 fPtSumThreshold = 1.;
318 fICMethod = kNoIC; // 0 don't isolate, 1 pt thresh method, 2 cone pt sum method
320 fIsolatePi0 = kFALSE ;
321 //-----------kSeveralIC-----------------
324 fConeSizes[0] = 0.1; fConeSizes[1] = 0.2; fConeSizes[2] = 0.3; fConeSizes[3] = 0.4; fConeSizes[4] = 0.5;
325 fPtThresholds[0]=0.; fPtThresholds[1]=1.; fPtThresholds[2]=2.; fPtThresholds[3]=3.; fPtThresholds[4]=4.;fPtThresholds[5]=5.;
329 //__________________________________________________________________
330 void AliAnaGammaDirect::MakeIsolationCut(TClonesArray * plCTS,
332 TParticle * pCandidate,
333 Int_t index, Int_t & n,
334 Bool_t &icmpt, Bool_t &icms,
335 Float_t &coneptsum) const
337 //Search in cone around a candidate particle if it is isolated
338 Float_t phiC = pCandidate->Phi() ;
339 Float_t etaC = pCandidate->Eta() ;
341 Float_t eta = -100. ;
342 Float_t phi = -100. ;
344 TParticle * particle = new TParticle;
351 //Check charged particles in cone.
352 for(Int_t ipr = 0;ipr < plCTS->GetEntries() ; ipr ++ ){
353 particle = dynamic_cast<TParticle *>(plCTS->At(ipr)) ;
355 eta = particle->Eta();
356 phi = particle->Phi() ;
358 //Check if there is any particle inside cone with pt larger than fPtThreshold
359 rad = TMath::Sqrt((eta-etaC)*(eta-etaC)+ (phi-phiC)*(phi-phiC));
361 AliDebug(3,Form("charged in cone pt %f, phi %f, eta %f, R %f ",pt,phi,eta,rad));
363 if(pt > fPtThreshold ) n++;
365 }// charged particle loop
367 //Check neutral particles in cone.
368 for(Int_t ipr = 0;ipr < plNe->GetEntries() ; ipr ++ ){
369 if(ipr != index){//Do not count the candidate
370 particle = dynamic_cast<TParticle *>(plNe->At(ipr)) ;
372 eta = particle->Eta();
373 phi = particle->Phi() ;
375 //Check if there is any particle inside cone with pt larger than fPtThreshold
376 rad = TMath::Sqrt((eta-etaC)*(eta-etaC)+ (phi-phiC)*(phi-phiC));
378 AliDebug(3,Form("charged in cone pt %f, phi %f, eta %f, R %f ",pt,phi,eta,rad));
380 if(pt > fPtThreshold ) n++;
383 }// neutral particle loop
387 if(coneptsum <= fPtSumThreshold)
392 //__________________________________________________________________
393 void AliAnaGammaDirect::MakeSeveralICAnalysis(TClonesArray * plCalo, TClonesArray * plCTS)
395 //Isolation Cut Analysis for both methods and different pt cuts and cones
396 if (fICMethod != kSeveralIC)
397 AliFatal("Remember to set in config file: directGamma->SetICMethod(kSeveralIC)");
399 //Search maximum energy photon in the event
402 TParticle * pCandidate = new TParticle();
403 Bool_t found = kFALSE;
405 for(Int_t ipr = 0;ipr < plCalo->GetEntries() ; ipr ++ ){
406 TParticle * particle = dynamic_cast<TParticle *>(plCalo->At(ipr)) ;
408 if((particle->Pt() > fMinGammaPt) && (particle->Pt() > ptC) &&
409 (particle->GetPdgCode() == 22 || (fIsolatePi0 && particle->GetPdgCode() == 111))){
411 ptC = particle->Pt();
412 pCandidate = particle ;
417 //If there is a large cluster, larger than threshold, study isolation cut
421 fhPhiGamma->Fill(ptC,pCandidate->Phi());
422 fhEtaGamma->Fill(ptC,pCandidate->Eta());
424 Int_t ncone[10][10];//[fNCones][fNPtThres];
425 Bool_t icPtThres = kFALSE;
426 Bool_t icPtSum = kFALSE;
428 for(Int_t icone = 0; icone<fNCones; icone++){
429 fConeSize=fConeSizes[icone] ;
430 Float_t coneptsum = 0 ;
431 for(Int_t ipt = 0; ipt<fNPtThres;ipt++){
433 fPtThreshold=fPtThresholds[ipt] ;
434 MakeIsolationCut(plCTS,plCalo, pCandidate, index,
435 ncone[icone][ipt], icPtThres, icPtSum,coneptsum);
436 AliDebug(1,Form("Candidate pt %f, pt in cone %f, Isolated? ICPt %d, ICSum %d",
437 pCandidate->Pt(), coneptsum, icPtThres, icPtSum));
438 // if(ptC >15 && ptC < 25 && (icPtThres || icPtSum) && ipt ==0){
439 // 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",
440 // fConeSize, fPtThreshold, fPtSumThreshold, ptC, pCandidate->Eta(), pCandidate->Phi()*TMath::RadToDeg(), coneptsum, icPtThres, icPtSum);
441 // //cout<<"mother label "<<pCandidate->GetMother(0)<<endl;
444 fhPtThresIsolated[icone][ipt]->Fill(ptC);
446 fhPtSumIsolated[icone]->Fill(ptC,coneptsum) ;
448 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]);
450 }//found high energy gamma in the event
453 //__________________________________________________________________
454 void AliAnaGammaDirect::Print(const Option_t * opt) const
457 //Print some relevant parameters set for the analysis
461 Info("Print", "%s %s", GetName(), GetTitle() ) ;
463 printf("Min Gamma pT = %f\n", fMinGammaPt) ;
464 printf("IC method = %d\n", fICMethod) ;
465 printf("Cone Size = %f\n", fConeSize) ;
466 if(fICMethod == kPtIC) printf("pT threshold = %f\n", fPtThreshold) ;
467 if(fICMethod == kSumPtIC) printf("pT sum threshold = %f\n", fPtSumThreshold) ;
469 if(fICMethod == kSeveralIC){
470 printf("N Cone Sizes = %d\n", fNCones) ;
471 printf("N pT thresholds = %d\n", fNPtThres) ;
472 printf("Cone Sizes = \n") ;
473 for(Int_t i = 0; i < fNCones; i++)
474 printf(" %f;", fConeSizes[i]) ;
476 for(Int_t i = 0; i < fNPtThres; i++)
477 printf(" %f;", fPtThresholds[i]) ;