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.6 2007/08/17 12:40:04 schutz
21 * New analysis classes by Gustavo Conesa
23 * Revision 1.4.4.4 2007/07/26 10:32:09 schutz
24 * new analysis classes in the the new analysis framework
29 //_________________________________________________________________________
30 // Class for the prompt gamma analysis, isolation cut
32 // Class created from old AliPHOSGammaJet
33 // (see AliRoot versions previous Release 4-09)
35 //*-- Author: Gustavo Conesa (LNF-INFN)
36 //////////////////////////////////////////////////////////////////////////////
39 // --- ROOT system ---
40 #include <TParticle.h>
43 #include "Riostream.h"
46 // --- AliRoot system ---
47 #include "AliAnaGammaDirect.h"
50 ClassImp(AliAnaGammaDirect)
52 //____________________________________________________________________________
53 AliAnaGammaDirect::AliAnaGammaDirect() :
56 fConeSize(0.),fPtThreshold(0.),fPtSumThreshold(0),
57 fICMethod(0), fAnaMC(0), fIsolatePi0(0),
58 fhNGamma(0),fhPhiGamma(0),fhEtaGamma(0), fhConeSumPt(0),
61 fNCones(0),fNPtThres(0), fConeSizes(), fPtThresholds(),
62 fhPtThresIsolated(), fhPtSumIsolated(), fntSeveralIC()
66 //Initialize parameters
71 //____________________________________________________________________________
72 AliAnaGammaDirect::AliAnaGammaDirect(const AliAnaGammaDirect & g) :
74 fMinGammaPt(g.fMinGammaPt),
75 fConeSize(g.fConeSize),
76 fPtThreshold(g.fPtThreshold),
77 fPtSumThreshold(g.fPtSumThreshold),
78 fICMethod(g.fICMethod),
80 fIsolatePi0(g.fIsolatePi0),
81 fhNGamma(g.fhNGamma),fhPhiGamma(g.fhPhiGamma),
82 fhEtaGamma(g.fhEtaGamma), fhConeSumPt(g.fhConeSumPt),
83 fntuplePrompt(g.fntuplePrompt),
85 fNCones(g.fNCones),fNPtThres(g.fNPtThres), fConeSizes(),fPtThresholds(),
86 fhPtThresIsolated(), fhPtSumIsolated()
91 for(Int_t i = 0; i < fNCones ; i++){
92 fConeSizes[i] = g.fConeSizes[i];
93 fhPtSumIsolated[i] = g.fhPtSumIsolated[i];
94 fntSeveralIC[i]= g.fntSeveralIC[i];
95 for(Int_t j = 0; j < fNPtThres ; j++)
96 fhPtThresIsolated[i][j] = g.fhPtThresIsolated[i][j];
99 for(Int_t i = 0; i < fNPtThres ; i++)
100 fPtThresholds[i]= g.fPtThresholds[i];
105 //_________________________________________________________________________
106 AliAnaGammaDirect & AliAnaGammaDirect::operator = (const AliAnaGammaDirect & source)
108 // assignment operator
110 if(&source == this) return *this;
112 fMinGammaPt = source.fMinGammaPt ;
113 fConeSize = source.fConeSize ;
114 fPtThreshold = source.fPtThreshold ;
115 fPtSumThreshold = source.fPtSumThreshold ;
116 fICMethod = source.fICMethod ;
117 fAnaMC = source.fAnaMC ;
118 fIsolatePi0 = source.fIsolatePi0 ;
120 fhNGamma = source.fhNGamma ;
121 fhPhiGamma = source.fhPhiGamma ;
122 fhEtaGamma = source.fhEtaGamma ;
123 fhConeSumPt = source.fhConeSumPt ;
125 fntuplePrompt = source.fntuplePrompt ;
128 fNCones = source.fNCones ;
129 fNPtThres = source.fNPtThres ;
131 for(Int_t i = 0; i < fNCones ; i++){
132 fConeSizes[i] = source.fConeSizes[i];
133 fhPtSumIsolated[i] = source.fhPtSumIsolated[i] ;
134 fntSeveralIC[i]= source.fntSeveralIC[i];
135 for(Int_t j = 0; j < fNPtThres ; j++)
136 fhPtThresIsolated[i][j] = source.fhPtThresIsolated[i][j] ;
139 for(Int_t i = 0; i < fNPtThres ; i++)
140 fPtThresholds[i]= source.fPtThresholds[i];
146 //____________________________________________________________________________
147 AliAnaGammaDirect::~AliAnaGammaDirect()
149 // Remove all pointers
155 delete fntuplePrompt ;
158 delete [] fhPtThresIsolated ;
159 delete [] fhPtSumIsolated ;
160 delete [] fntSeveralIC ;
164 //________________________________________________________________________
165 TList * AliAnaGammaDirect::GetCreateOutputObjects()
168 // Create histograms to be saved in output file and
169 // store them in outputContainer
170 TList * outputContainer = new TList() ;
171 outputContainer->SetName("DirectGammaHistos") ;
173 //Histograms of highest gamma identified in Event
174 fhNGamma = new TH1F("NGamma","Number of #gamma over calorimeter",240,0,120);
175 fhNGamma->SetYTitle("N");
176 fhNGamma->SetXTitle("p_{T #gamma}(GeV/c)");
177 outputContainer->Add(fhNGamma) ;
179 fhPhiGamma = new TH2F
180 ("PhiGamma","#phi_{#gamma}",200,0,120,200,0,7);
181 fhPhiGamma->SetYTitle("#phi");
182 fhPhiGamma->SetXTitle("p_{T #gamma} (GeV/c)");
183 outputContainer->Add(fhPhiGamma) ;
185 fhEtaGamma = new TH2F
186 ("EtaGamma","#phi_{#gamma}",200,0,120,200,-0.8,0.8);
187 fhEtaGamma->SetYTitle("#eta");
188 fhEtaGamma->SetXTitle("p_{T #gamma} (GeV/c)");
189 outputContainer->Add(fhEtaGamma) ;
191 fhConeSumPt = new TH2F
192 ("ConePtSum","#Sigma p_{T} in cone ",200,0,120,100,0,100);
193 fhConeSumPt->SetYTitle("#Sigma p_{T}");
194 fhConeSumPt->SetXTitle("p_{T #gamma} (GeV/c)");
195 outputContainer->Add(fhConeSumPt) ;
198 fntuplePrompt = new TNtuple("ntuplePromptGamma", "Tree of prompt #gamma", "ptcluster:phicluster:etacluster:pdg:status:ptprimary:phiprimary:etaprimary:pdgprimary:statusprimary");
199 outputContainer->Add(fntuplePrompt) ;
201 if(fICMethod == kSeveralIC){
204 for(Int_t icone = 0; icone<fNCones; icone++){
205 sprintf(name,"PtSumIsolated_Cone_%d",icone);
206 sprintf(title,"Candidate cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
207 fhPtSumIsolated[icone] = new TH2F(name, title,240,0,120,120,0,10);
208 fhPtSumIsolated[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
209 fhPtSumIsolated[icone]->SetXTitle("p_{T} (GeV/c)");
210 outputContainer->Add(fhPtSumIsolated[icone]) ;
212 sprintf(name,"nt_Cone_%d",icone);
213 sprintf(title,"ntuple for cone size %d",icone);
214 fntSeveralIC[icone] = new TNtuple(name, title, "ptcand:phicand:etacand:ptsum:type:ncone0:ncone1:ncone2:ncone3:ncone4:ncone5");
215 outputContainer->Add(fntSeveralIC[icone]) ;
217 for(Int_t ipt = 0; ipt<fNPtThres;ipt++){
218 sprintf(name,"PtThresIsol_Cone_%d_Pt%d",icone,ipt);
219 sprintf(title,"Isolated candidate p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
220 fhPtThresIsolated[icone][ipt] = new TH1F(name, title,240,0,120);
221 fhPtThresIsolated[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
222 outputContainer->Add(fhPtThresIsolated[icone][ipt]) ;
229 return outputContainer ;
233 //____________________________________________________________________________
234 void AliAnaGammaDirect::GetPromptGamma(TClonesArray * plCalo, TClonesArray * plCTS, TClonesArray * plPrimCalo, TParticle *pGamma, Bool_t &found) const
236 //Search for the prompt photon in Calorimeter with pt > fMinGammaPt
241 Float_t coneptsum = 0 ;
243 for(Int_t ipr = 0;ipr < plCalo->GetEntries() ; ipr ++ ){
244 TParticle * particle = dynamic_cast<TParticle *>(plCalo->At(ipr)) ;
246 if((particle->Pt() > fMinGammaPt) && (particle->Pt() > pt) &&
247 (particle->GetPdgCode() == 22 || (fIsolatePi0 && particle->GetPdgCode() == 111))){
250 pGamma->SetMomentum(particle->Px(),particle->Py(),particle->Pz(),particle->Energy());
251 pGamma->SetStatusCode(particle->GetStatusCode());
252 pGamma->SetPdgCode(particle->GetPdgCode());
258 if( ( fICMethod == kPtIC || fICMethod == kSumPtIC ) && found)
260 Bool_t icPtThres = kFALSE;
261 Bool_t icPtSum = kFALSE;
262 MakeIsolationCut(plCTS,plCalo, pGamma, index,n,
263 icPtThres, icPtSum,coneptsum);
264 if(fICMethod == kPtIC) //Pt thres method
266 if(fICMethod == kSumPtIC) //Pt cone sum method
271 AliDebug(1,Form("Cluster with p_{T} larger than the pt treshold %f found in calorimeter ", fMinGammaPt)) ;
272 AliDebug(1,Form("Gamma: pt %f, phi %f, eta %f", pGamma->Pt(),pGamma->Phi(),pGamma->Eta())) ;
273 //Fill prompt gamma histograms
274 Float_t ptcluster = pGamma->Pt();
275 Float_t phicluster = pGamma->Phi();
276 Float_t etacluster = pGamma->Eta();
277 Int_t statuscluster = pGamma->GetStatusCode();
278 Int_t pdgcluster = pGamma->GetPdgCode();
280 fhNGamma ->Fill(ptcluster);
281 fhPhiGamma ->Fill(ptcluster,phicluster);
282 fhEtaGamma ->Fill(ptcluster,etacluster);
283 fhConeSumPt->Fill(ptcluster,coneptsum);
285 Float_t ptprimary = 0 ;
286 Float_t phiprimary = 0 ;
287 Float_t etaprimary = 0 ;
288 Int_t pdgprimary = 0 ;
289 Int_t statusprimary = 0 ;
292 TParticle * primary = dynamic_cast<TParticle *>(plPrimCalo->At(index)) ;
293 ptprimary = primary->Pt();
294 phiprimary = primary->Phi();
295 etaprimary = primary->Eta();
296 pdgprimary = TMath::Abs(primary->GetPdgCode()) ;
297 statusprimary = primary->GetStatusCode(); // = 2 means decay gamma!!!
299 AliDebug(1, Form("Identified prompt Gamma pT %2.2f; Primary, pdg %d, pT %2.2f",ptcluster,pdgprimary,ptprimary));
300 //printf("Identified prompt Gamma pT %2.2f; Primary, pdg %d, pT %2.2f \n",ptcluster,pdgprimary,ptprimary);
303 //Fill ntuple with cluster / MC data
305 fntuplePrompt->Fill(ptcluster,phicluster,etacluster,pdgcluster,statuscluster,ptprimary,phiprimary, etaprimary,pdgprimary,statusprimary);
308 AliDebug(1,Form("NO Cluster with pT larger than %f found in calorimeter ", fMinGammaPt)) ;
311 //____________________________________________________________________________
312 void AliAnaGammaDirect::InitParameters()
315 //Initialize the parameters of the analysis.
318 //Fill particle lists when PID is ok
321 fPtSumThreshold = 1.;
323 fICMethod = kNoIC; // 0 don't isolate, 1 pt thresh method, 2 cone pt sum method
325 fIsolatePi0 = kFALSE ;
326 //-----------kSeveralIC-----------------
329 fConeSizes[0] = 0.1; fConeSizes[1] = 0.2; fConeSizes[2] = 0.3; fConeSizes[3] = 0.4; fConeSizes[4] = 0.5;
330 fPtThresholds[0]=0.; fPtThresholds[1]=1.; fPtThresholds[2]=2.; fPtThresholds[3]=3.; fPtThresholds[4]=4.;fPtThresholds[5]=5.;
334 //__________________________________________________________________
335 void AliAnaGammaDirect::MakeIsolationCut(TClonesArray * plCTS,
337 TParticle * pCandidate,
338 Int_t index, Int_t & n,
339 Bool_t &icmpt, Bool_t &icms,
340 Float_t &coneptsum) const
342 //Search in cone around a candidate particle if it is isolated
343 Float_t phiC = pCandidate->Phi() ;
344 Float_t etaC = pCandidate->Eta() ;
346 Float_t eta = -100. ;
347 Float_t phi = -100. ;
349 TParticle * particle = new TParticle;
356 //Check charged particles in cone.
357 for(Int_t ipr = 0;ipr < plCTS->GetEntries() ; ipr ++ ){
358 particle = dynamic_cast<TParticle *>(plCTS->At(ipr)) ;
360 eta = particle->Eta();
361 phi = particle->Phi() ;
363 //Check if there is any particle inside cone with pt larger than fPtThreshold
364 rad = TMath::Sqrt((eta-etaC)*(eta-etaC)+ (phi-phiC)*(phi-phiC));
366 AliDebug(3,Form("charged in cone pt %f, phi %f, eta %f, R %f ",pt,phi,eta,rad));
368 if(pt > fPtThreshold ) n++;
370 }// charged particle loop
372 //Check neutral particles in cone.
373 for(Int_t ipr = 0;ipr < plNe->GetEntries() ; ipr ++ ){
374 if(ipr != index){//Do not count the candidate
375 particle = dynamic_cast<TParticle *>(plNe->At(ipr)) ;
377 eta = particle->Eta();
378 phi = particle->Phi() ;
380 //Check if there is any particle inside cone with pt larger than fPtThreshold
381 rad = TMath::Sqrt((eta-etaC)*(eta-etaC)+ (phi-phiC)*(phi-phiC));
383 AliDebug(3,Form("charged in cone pt %f, phi %f, eta %f, R %f ",pt,phi,eta,rad));
385 if(pt > fPtThreshold ) n++;
388 }// neutral particle loop
392 if(coneptsum <= fPtSumThreshold)
397 //__________________________________________________________________
398 void AliAnaGammaDirect::MakeSeveralICAnalysis(TClonesArray * plCalo, TClonesArray * plCTS)
400 //Isolation Cut Analysis for both methods and different pt cuts and cones
401 if (fICMethod != kSeveralIC)
402 AliFatal("Remember to set in config file: directGamma->SetICMethod(kSeveralIC)");
405 //Search maximum energy photon in the event
408 TParticle * pCandidate = new TParticle();
409 Bool_t found = kFALSE;
411 for(Int_t ipr = 0;ipr < plCalo->GetEntries() ; ipr ++ ){
412 TParticle * particle = dynamic_cast<TParticle *>(plCalo->At(ipr)) ;
414 if((particle->Pt() > fMinGammaPt) && (particle->Pt() > ptC) &&
415 (particle->GetPdgCode() == 22 || (fIsolatePi0 && particle->GetPdgCode() == 111))){
417 ptC = particle->Pt();
418 pCandidate = particle ;
423 //If there is a large cluster, larger than threshold, study isolation cut
427 fhPhiGamma->Fill(ptC,pCandidate->Phi());
428 fhEtaGamma->Fill(ptC,pCandidate->Eta());
430 Int_t ncone[10][10];//[fNCones][fNPtThres];
431 Bool_t icPtThres = kFALSE;
432 Bool_t icPtSum = kFALSE;
434 for(Int_t icone = 0; icone<fNCones; icone++){
435 fConeSize=fConeSizes[icone] ;
436 Float_t coneptsum = 0 ;
437 for(Int_t ipt = 0; ipt<fNPtThres;ipt++){
439 fPtThreshold=fPtThresholds[ipt] ;
440 MakeIsolationCut(plCTS,plCalo, pCandidate, index,
441 ncone[icone][ipt], icPtThres, icPtSum,coneptsum);
442 AliDebug(1,Form("Candidate pt %f, pt in cone %f, Isolated? ICPt %d, ICSum %d",
443 pCandidate->Pt(), coneptsum, icPtThres, icPtSum));
444 // if(ptC >15 && ptC < 25 && (icPtThres || icPtSum) && ipt ==0){
445 // 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",
446 // fConeSize, fPtThreshold, fPtSumThreshold, ptC, pCandidate->Eta(), pCandidate->Phi()*TMath::RadToDeg(), coneptsum, icPtThres, icPtSum);
447 // //cout<<"mother label "<<pCandidate->GetMother(0)<<endl;
450 fhPtThresIsolated[icone][ipt]->Fill(ptC);
452 fhPtSumIsolated[icone]->Fill(ptC,coneptsum) ;
454 fntSeveralIC[icone]->Fill(ptC,pCandidate->Phi(),pCandidate->Eta(), coneptsum,type,ncone[icone][0],ncone[icone][1],ncone[icone][2],ncone[icone][3],ncone[icone][4],ncone[icone][5]);
456 }//found high energy gamma in the event
459 //__________________________________________________________________
460 void AliAnaGammaDirect::Print(const Option_t * opt) const
463 //Print some relevant parameters set for the analysis
467 Info("Print", "%s %s", GetName(), GetTitle() ) ;
469 printf("Min Gamma pT = %f\n", fMinGammaPt) ;
470 printf("IC method = %d\n", fICMethod) ;
471 printf("Cone Size = %f\n", fConeSize) ;
472 if(fICMethod == kPtIC) printf("pT threshold = %f\n", fPtThreshold) ;
473 if(fICMethod == kSumPtIC) printf("pT sum threshold = %f\n", fPtSumThreshold) ;
475 if(fICMethod == kSeveralIC){
476 printf("N Cone Sizes = %d\n", fNCones) ;
477 printf("N pT thresholds = %d\n", fNPtThres) ;
478 printf("Cone Sizes = \n") ;
479 for(Int_t i = 0; i < fNCones; i++)
480 printf(" %f;", fConeSizes[i]) ;
482 for(Int_t i = 0; i < fNPtThres; i++)
483 printf(" %f;", fPtThresholds[i]) ;