Compilation on Windows/Cygwin
[u/mrichter/AliRoot.git] / PWG4 / AliAnaGammaDirect.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
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  **************************************************************************/
15 /* $Id$ */
16
17 /* History of cvs commits:
18  *
19  * $Log$
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)
22  *
23  * Revision 1.8  2007/10/29 13:48:42  gustavo
24  * Corrected coding violations
25  *
26  * Revision 1.6  2007/08/17 12:40:04  schutz
27  * New analysis classes by Gustavo Conesa
28  *
29  * Revision 1.4.4.4  2007/07/26 10:32:09  schutz
30  * new analysis classes in the the new analysis framework
31  *
32  *
33  */
34
35 //_________________________________________________________________________
36 // Class for the prompt gamma analysis, isolation cut
37 //
38 //  Class created from old AliPHOSGammaJet 
39 //  (see AliRoot versions previous Release 4-09)
40 //
41 //*-- Author: Gustavo Conesa (LNF-INFN) 
42 //////////////////////////////////////////////////////////////////////////////
43   
44   
45 // --- ROOT system --- 
46 #include <TParticle.h>
47 #include <TH2.h>
48 #include <TList.h>
49 #include "Riostream.h"
50 #include "TROOT.h"
51
52 // --- AliRoot system --- 
53 #include "AliAnaGammaDirect.h" 
54 #include "AliLog.h"
55
56 ClassImp(AliAnaGammaDirect)
57   
58 //____________________________________________________________________________
59   AliAnaGammaDirect::AliAnaGammaDirect() : 
60     TObject(),
61     fMinGammaPt(0.),
62     fConeSize(0.),fPtThreshold(0.),fPtSumThreshold(0), 
63     fICMethod(0), fAnaMC(0), fIsolatePi0(0),
64     fhNGamma(0),fhPhiGamma(0),fhEtaGamma(0), fhConeSumPt(0), 
65     fntuplePrompt(0),
66     //kSeveralIC
67     fNCones(0),fNPtThres(0), fConeSizes(),  fPtThresholds(), 
68     fhPtThresIsolated(), fhPtSumIsolated(), fntSeveralIC()
69 {
70   //default ctor
71   
72   //Initialize parameters
73   InitParameters();
74
75 }
76
77 //____________________________________________________________________________
78 AliAnaGammaDirect::AliAnaGammaDirect(const AliAnaGammaDirect & g) : 
79   TObject(g),
80   fMinGammaPt(g.fMinGammaPt), 
81   fConeSize(g.fConeSize),
82   fPtThreshold(g.fPtThreshold),
83   fPtSumThreshold(g.fPtSumThreshold), 
84   fICMethod(g.fICMethod),
85   fAnaMC( g.fAnaMC), 
86   fIsolatePi0(g.fIsolatePi0),
87   fhNGamma(g.fhNGamma),fhPhiGamma(g.fhPhiGamma),
88   fhEtaGamma(g.fhEtaGamma), fhConeSumPt(g.fhConeSumPt),    
89   fntuplePrompt(g.fntuplePrompt),
90   //kSeveralIC
91   fNCones(g.fNCones),fNPtThres(g.fNPtThres), fConeSizes(),fPtThresholds(), 
92   fhPtThresIsolated(), fhPtSumIsolated()
93 {
94   // cpy ctor
95   
96   //kSeveralIC
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]; 
103   }
104   
105   for(Int_t i = 0; i < fNPtThres ; i++)
106     fPtThresholds[i]=   g.fPtThresholds[i];
107
108
109 }
110
111 //_________________________________________________________________________
112 AliAnaGammaDirect & AliAnaGammaDirect::operator = (const AliAnaGammaDirect & source)
113 {
114   // assignment operator
115   
116   if(&source == this) return *this;
117
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 ;
125
126   fhNGamma = source.fhNGamma ; 
127   fhPhiGamma = source.fhPhiGamma ;
128   fhEtaGamma = source.fhEtaGamma ;
129   fhConeSumPt = source.fhConeSumPt ;
130
131   fntuplePrompt = source.fntuplePrompt ;
132
133   //kSeveralIC
134   fNCones = source.fNCones ;
135   fNPtThres = source.fNPtThres ; 
136    
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] ;
143   }
144   
145   for(Int_t i = 0; i < fNPtThres ; i++)
146     fPtThresholds[i]=   source.fPtThresholds[i];
147
148   return *this;
149   
150 }
151
152 //____________________________________________________________________________
153 AliAnaGammaDirect::~AliAnaGammaDirect() 
154 {
155   // Remove all pointers except analysis output pointers.
156
157 }
158
159 //________________________________________________________________________
160 TList *  AliAnaGammaDirect::GetCreateOutputObjects()
161 {  
162
163   // Create histograms to be saved in output file and 
164   // store them in outputContainer
165   TList * outputContainer = new TList() ; 
166   outputContainer->SetName("DirectGammaHistos") ; 
167   
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) ; 
173   
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) ; 
179   
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) ;
185
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) ;
191   
192   //NTUPLE
193   fntuplePrompt = new TNtuple("ntuplePromptGamma", "Tree of prompt #gamma", "ptcluster:phicluster:etacluster:pdg:status:ptprimary:phiprimary:etaprimary:pdgprimary:statusprimary");
194   outputContainer->Add(fntuplePrompt) ;
195
196   if(fICMethod == kSeveralIC){
197     char name[128];
198     char title[128];
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]) ; 
206       
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]) ; 
211
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]) ; 
218       }//icone loop
219     }//ipt loop
220   }
221
222 //  gROOT->cd();
223
224   return outputContainer ;
225
226 }
227
228 //____________________________________________________________________________
229 void AliAnaGammaDirect::GetPromptGamma(TClonesArray * plCalo, TClonesArray * plCTS, TClonesArray * plPrimCalo, TParticle *pGamma, Bool_t &found) const 
230 {
231   //Search for the prompt photon in Calorimeter with pt > fMinGammaPt
232
233   Double_t pt = 0;
234   Int_t n = 0;
235   Int_t index = -1; 
236   Float_t coneptsum = 0 ;
237
238   for(Int_t ipr = 0;ipr < plCalo->GetEntries() ; ipr ++ ){
239     TParticle * particle = dynamic_cast<TParticle *>(plCalo->At(ipr)) ;
240
241     if((particle->Pt() > fMinGammaPt) && (particle->Pt() > pt) && 
242        (particle->GetPdgCode() == 22 || (fIsolatePi0 && particle->GetPdgCode() == 111))){
243       index = ipr ;
244       pt = particle->Pt();
245       pGamma->SetMomentum(particle->Px(),particle->Py(),particle->Pz(),particle->Energy());
246       pGamma->SetStatusCode(particle->GetStatusCode());
247       pGamma->SetPdgCode(particle->GetPdgCode());
248       found  = kTRUE;
249     }
250   }
251  
252   //Do Isolation?
253   if( ( fICMethod == kPtIC  ||  fICMethod == kSumPtIC )  && found)
254     {
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
260         found = icPtThres ;
261       if(fICMethod == kSumPtIC) //Pt cone sum method
262         found = icPtSum ;
263     }
264   
265   if(found){
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();
274
275     fhNGamma   ->Fill(ptcluster);
276     fhPhiGamma ->Fill(ptcluster,phicluster);
277     fhEtaGamma ->Fill(ptcluster,etacluster);
278     fhConeSumPt->Fill(ptcluster,coneptsum);
279
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 ;
285
286     if(fAnaMC){
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!!!
293  
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);
296     }
297
298     //Fill ntuple with cluster / MC data
299 //    gROOT->cd();
300     fntuplePrompt->Fill(ptcluster,phicluster,etacluster,pdgcluster,statuscluster,ptprimary,phiprimary, etaprimary,pdgprimary,statusprimary);
301   }
302   else
303     AliDebug(1,Form("NO Cluster with pT larger than %f found in calorimeter ", fMinGammaPt)) ;
304 }
305
306   //____________________________________________________________________________
307 void AliAnaGammaDirect::InitParameters()
308 {
309  
310   //Initialize the parameters of the analysis.
311   fMinGammaPt  = 5. ;
312
313   //Fill particle lists when PID is ok
314   fConeSize             = 0.4 ; 
315   fPtThreshold         = 1.; 
316   fPtSumThreshold  = 1.; 
317
318   fICMethod = kNoIC; // 0 don't isolate, 1 pt thresh method, 2 cone pt sum method
319   fAnaMC = kFALSE ;
320   fIsolatePi0 = kFALSE ;
321  //-----------kSeveralIC-----------------
322   fNCones           = 5 ; 
323   fNPtThres         = 6 ; 
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.;
326
327 }
328
329 //__________________________________________________________________
330 void  AliAnaGammaDirect::MakeIsolationCut(TClonesArray * plCTS, 
331                                           TClonesArray * plNe, 
332                                           TParticle * pCandidate, 
333                                           Int_t index, Int_t & n,
334                                           Bool_t  &icmpt,  Bool_t  &icms, 
335                                           Float_t &coneptsum) const 
336 {  
337   //Search in cone around a candidate particle if it is isolated 
338   Float_t phiC  = pCandidate->Phi() ;
339   Float_t etaC = pCandidate->Eta() ;
340   Float_t pt     = -100. ;
341   Float_t eta   = -100.  ;
342   Float_t phi    = -100.  ;
343   Float_t rad   = -100 ;
344   TParticle * particle  = new TParticle;
345
346   n = 0 ;
347   coneptsum = 0.; 
348   icmpt = kFALSE;
349   icms  = kFALSE;
350
351   //Check charged particles in cone.
352   for(Int_t ipr = 0;ipr < plCTS->GetEntries() ; ipr ++ ){
353     particle = dynamic_cast<TParticle *>(plCTS->At(ipr)) ;
354     pt    = particle->Pt();
355     eta  = particle->Eta();
356     phi  = particle->Phi() ;
357     
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));
360     if(rad < fConeSize){
361       AliDebug(3,Form("charged in cone pt %f, phi %f, eta %f, R %f ",pt,phi,eta,rad));
362       coneptsum+=pt;
363       if(pt > fPtThreshold ) n++;
364     }
365   }// charged particle loop
366   
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)) ;
371       pt    = particle->Pt();
372       eta  = particle->Eta();
373       phi  = particle->Phi() ;
374       
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));
377       if(rad < fConeSize){
378         AliDebug(3,Form("charged in cone pt %f, phi %f, eta %f, R %f ",pt,phi,eta,rad));
379         coneptsum+=pt;
380         if(pt > fPtThreshold ) n++;
381       }
382     }
383   }// neutral particle loop
384   
385   if(n == 0) 
386     icmpt =  kTRUE ;
387   if(coneptsum <= fPtSumThreshold)
388     icms  =  kTRUE ;
389
390 }
391
392 //__________________________________________________________________
393 void  AliAnaGammaDirect::MakeSeveralICAnalysis(TClonesArray * plCalo, TClonesArray * plCTS) 
394 {
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)");
398
399   //Search maximum energy photon in the event
400   Double_t ptC = 0;
401   Int_t index = -1; 
402   TParticle * pCandidate = new TParticle();
403   Bool_t found = kFALSE;
404   
405   for(Int_t ipr = 0;ipr < plCalo->GetEntries() ; ipr ++ ){
406     TParticle * particle = dynamic_cast<TParticle *>(plCalo->At(ipr)) ;
407     
408     if((particle->Pt() > fMinGammaPt) && (particle->Pt() > ptC) && 
409        (particle->GetPdgCode() == 22 ||  (fIsolatePi0 && particle->GetPdgCode() == 111))){
410       index = ipr ;
411       ptC = particle->Pt();
412       pCandidate = particle ;
413       found  = kTRUE;
414     }
415   }
416   
417   //If there is a large cluster, larger than threshold, study isolation cut
418   if(found){
419     
420     fhNGamma->Fill(ptC);
421     fhPhiGamma->Fill(ptC,pCandidate->Phi());
422     fhEtaGamma->Fill(ptC,pCandidate->Eta());
423     
424     Int_t ncone[10][10];//[fNCones][fNPtThres];
425     Bool_t  icPtThres   = kFALSE;
426     Bool_t  icPtSum     = kFALSE;
427     
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++){
432         ncone[icone][ipt]=0;
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;
442 //      }
443         
444         fhPtThresIsolated[icone][ipt]->Fill(ptC); 
445       }//pt thresh loop
446       fhPtSumIsolated[icone]->Fill(ptC,coneptsum) ;
447 //      gROOT->cd();
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]);
449     }//cone size loop
450   }//found high energy gamma in the event
451 }
452
453 //__________________________________________________________________
454 void AliAnaGammaDirect::Print(const Option_t * opt) const
455 {
456   
457   //Print some relevant parameters set for the analysis
458   if(! opt)
459     return;
460   
461   Info("Print", "%s %s", GetName(), GetTitle() ) ;
462   
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) ;
468    
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]) ;
475     printf("    \n") ;
476     for(Int_t i = 0; i < fNPtThres; i++)
477       printf("   %f;",  fPtThresholds[i]) ;
478   }
479
480   printf("    \n") ;
481   
482