New Gamma package
[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  *
21  */
22
23 //_________________________________________________________________________
24 // Class for the analysis of gamma correlations (gamma-jet, 
25 // gamma-hadron and isolation cut.
26 // This class contains 3 main methods: one to fill lists of particles (ESDs) comming 
27 //  from the CTS (ITS+TPC) and the calorimeters;  the other one tags a candidate 
28 //  cluster as isolated;  the last one search in the 
29 //  corresponing calorimeter for the highest energy cluster, identified it as 
30 //  prompt photon;
31 //
32 //  Class created from old AliPHOSGammaJet 
33 //  (see AliRoot versions previous Release 4-09)
34 //
35 //*-- Author: Gustavo Conesa (LNF-INFN) 
36 //////////////////////////////////////////////////////////////////////////////
37
38
39 // --- ROOT system ---
40
41 #include <TFile.h>
42 #include <TParticle.h>
43 #include <TH2.h>
44
45 #include "AliAnaGammaDirect.h" 
46 #include "AliESD.h"
47 #include "AliESDtrack.h"
48 #include "AliESDCaloCluster.h"
49 #include "Riostream.h"
50 #include "AliLog.h"
51
52 ClassImp(AliAnaGammaDirect)
53
54 //____________________________________________________________________________
55   AliAnaGammaDirect::AliAnaGammaDirect(const char *name) : 
56     AliAnalysisTask(name,""), fChain(0), fESD(0),
57     fOutputContainer(new TObjArray(100)), 
58     fPrintInfo(0), fMinGammaPt(0.),
59     fCalorimeter(""), fEMCALPID(0),fPHOSPID(0),
60     fConeSize(0.),fPtThreshold(0.),fPtSumThreshold(0), 
61     fNCones(0),fNPtThres(0),fMakeICMethod(0)
62 {
63   //Ctor        
64   TList * list = gDirectory->GetListOfKeys() ; 
65   TIter next(list) ; 
66   TH2F * h = 0 ;
67   Int_t index ; 
68   for (index = 0 ; index < list->GetSize()-1 ; index++) { 
69     //-1 to avoid GammaJet Task
70     h = dynamic_cast<TH2F*>(gDirectory->Get(list->At(index)->GetName())) ; 
71     fOutputContainer->Add(h) ; 
72   }
73   
74   for(Int_t i = 0; i < 10 ; i++){
75     fConeSizes[i]=0;
76     fPtThresholds[i]=0;
77   }
78   
79   // Input slot #0 works with an Ntuple
80   DefineInput(0, TChain::Class());
81   // Output slot #0 writes into a TH1 container
82   DefineOutput(0,  TObjArray::Class()) ; 
83   
84 }
85
86
87 //____________________________________________________________________________
88 AliAnaGammaDirect::AliAnaGammaDirect(const AliAnaGammaDirect & g) : 
89   AliAnalysisTask(g), fChain(g.fChain), fESD(g.fESD),
90   fOutputContainer(g. fOutputContainer),  fPrintInfo(g.fPrintInfo),
91   fMinGammaPt(g.fMinGammaPt), fCalorimeter(g.fCalorimeter),
92   fEMCALPID(g.fEMCALPID),fPHOSPID(g.fPHOSPID),
93   fConeSize(g.fConeSize),
94   fPtThreshold(g.fPtThreshold),
95   fPtSumThreshold(g.fPtSumThreshold), 
96   fNCones(g.fNCones),fNPtThres(g.fNPtThres),
97   fMakeICMethod(g.fMakeICMethod)
98 {
99   // cpy ctor
100   SetName (g.GetName()) ; 
101   SetTitle(g.GetTitle()) ; 
102
103   for(Int_t i = 0; i < 10 ; i++){
104     fConeSizes[i]=  g.fConeSizes[i];
105     fPtThresholds[i]=   g.fPtThresholds[i];
106   }
107 }
108
109 //____________________________________________________________________________
110 AliAnaGammaDirect::~AliAnaGammaDirect() 
111 {
112   // Remove all pointers
113   fOutputContainer->Clear() ; 
114   delete fOutputContainer ;
115   delete fhNGamma    ;  
116   delete fhPhiGamma  ; 
117   delete fhEtaGamma   ;  
118   delete fhPtCandidate ;
119   delete [] fhPtThresIsolated ;
120   delete [] fhPtSumIsolated ;
121
122 }
123
124 //____________________________________________________________________________
125 void AliAnaGammaDirect::CreateParticleList(TClonesArray * pl, 
126                                         TClonesArray * plCTS, 
127                                         TClonesArray * plEMCAL,  
128                                         TClonesArray * plPHOS){
129   
130   //Create a list of particles from the ESD. These particles have been measured 
131   //by the Central Tracking system (TPC+ITS), PHOS and EMCAL 
132   
133   Int_t index = pl->GetEntries() ; 
134   Int_t npar  = 0 ;
135   Float_t *pid = new Float_t[AliPID::kSPECIESN];  
136   AliDebug(3,"Fill particle lists");
137   if(fPrintInfo)
138     AliInfo(Form("fCalorimeter %s",fCalorimeter.Data()));
139
140   Double_t v[3] ; //vertex ;
141   fESD->GetVertex()->GetXYZ(v) ; 
142
143   //########### PHOS ##############
144   if(fCalorimeter == "PHOS"){
145
146     Int_t begphos = fESD->GetFirstPHOSCluster();  
147     Int_t endphos = fESD->GetFirstPHOSCluster() + 
148       fESD->GetNumberOfPHOSClusters() ;  
149     Int_t indexNePHOS = plPHOS->GetEntries() ;
150     AliDebug(3,Form("First PHOS particle %d, last particle %d", begphos,endphos));
151
152     if(fCalorimeter == "PHOS"){
153       for (npar =  begphos; npar <  endphos; npar++) {//////////////PHOS track loop
154         AliESDCaloCluster * clus = fESD->GetCaloCluster(npar) ; // retrieve track from esd
155
156         //Create a TParticle to fill the particle list
157         
158         Float_t en = clus->GetClusterEnergy() ;
159         Float_t *p = new Float_t();
160         clus->GetGlobalPosition(p) ;
161         TVector3 pos(p[0],p[1],p[2]) ; 
162         Double_t phi  = pos.Phi();
163         Double_t theta= pos.Theta();
164         Double_t px = en*TMath::Cos(phi)*TMath::Sin(theta);;
165         Double_t py = en*TMath::Sin(phi)*TMath::Sin(theta);
166         Double_t pz = en*TMath::Cos(theta);
167
168         TParticle * particle = new TParticle() ;
169         particle->SetMomentum(px,py,pz,en) ;
170         AliDebug(4,Form("PHOS clusters: pt %f, phi %f, eta %f", particle->Pt(),particle->Phi(),particle->Eta()));
171         
172         //Select only photons
173         
174         pid=clus->GetPid();
175         //cout<<"pid "<<pid[AliPID::kPhoton]<<endl ;
176         if( !fPHOSPID)
177           new((*plPHOS)[indexNePHOS++])   TParticle(*particle) ;
178         else if( pid[AliPID::kPhoton] > 0.75)
179           new((*plPHOS)[indexNePHOS++])   TParticle(*particle) ;
180       }
181     }
182   }
183
184   //########### CTS (TPC+ITS) #####################
185   Int_t begtpc   = 0 ;  
186   Int_t endtpc   = fESD->GetNumberOfTracks() ;
187   Int_t indexCh  = plCTS->GetEntries() ;
188   AliDebug(3,Form("First CTS particle %d, last particle %d", begtpc,endtpc));
189   
190   for (npar =  begtpc; npar <  endtpc; npar++) {////////////// track loop
191     AliESDtrack * track = fESD->GetTrack(npar) ; // retrieve track from esd
192     //We want tracks fitted in the detectors:
193     ULong_t status=AliESDtrack::kTPCrefit;
194     status|=AliESDtrack::kITSrefit;
195    
196     //We want tracks whose PID bit is set:
197     //     ULong_t status =AliESDtrack::kITSpid;
198     //     status|=AliESDtrack::kTPCpid;
199   
200     if ( (track->GetStatus() & status) == status) {//Check if the bits we want are set
201       // Do something with the tracks which were successfully
202       // re-fitted 
203       Double_t en = 0; //track ->GetTPCsignal() ;
204       Double_t mom[3];
205       track->GetPxPyPz(mom) ;
206       Double_t px = mom[0];
207       Double_t py = mom[1];
208       Double_t pz = mom[2]; //Check with TPC people if this is correct.
209       Int_t pdg = 11; //Give any charged PDG code, in this case electron.
210       //I just want to tag the particle as charged
211        TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1, 
212                                                  px, py, pz, en, v[0], v[1], v[2], 0);
213   
214       //TParticle * particle = new TParticle() ;
215       //particle->SetMomentum(px,py,pz,en) ;
216  
217       new((*plCTS)[indexCh++])       TParticle(*particle) ;    
218       new((*pl)[index++])           TParticle(*particle) ;
219     }
220   }
221   
222   //################ EMCAL ##############
223   
224   Int_t begem = fESD->GetFirstEMCALCluster();  
225   Int_t endem = fESD->GetFirstEMCALCluster() + 
226     fESD->GetNumberOfEMCALClusters() ;  
227   Int_t indexNe  = plEMCAL->GetEntries() ; 
228   
229   AliDebug(3,Form("First EMCAL particle %d, last particle %d",begem,endem));
230     
231     for (npar =  begem; npar <  endem; npar++) {//////////////EMCAL track loop
232       AliESDCaloCluster * clus = fESD->GetCaloCluster(npar) ; // retrieve track from esd
233       Int_t clustertype= clus->GetClusterType();
234       if(clustertype == AliESDCaloCluster::kClusterv1){
235         Float_t en = clus->GetClusterEnergy() ;
236         Float_t *p = new Float_t();
237         clus->GetGlobalPosition(p) ;
238         TVector3 pos(p[0],p[1],p[2]) ;
239         Double_t phi  = pos.Phi();
240         Double_t theta= pos.Theta();
241         Double_t px = en*TMath::Cos(phi)*TMath::Sin(theta);;
242         Double_t py = en*TMath::Sin(phi)*TMath::Sin(theta);
243         Double_t pz = en*TMath::Cos(theta);
244
245         pid=clus->GetPid();
246         if(fCalorimeter == "EMCAL")
247           {
248             TParticle * particle = new TParticle() ;
249             particle->SetMomentum(px,py,pz,en) ;
250             AliDebug(4,Form("EMCAL clusters: pt %f, phi %f, eta %f", particle->Pt(),particle->Phi(),particle->Eta()));
251             if(!fEMCALPID) //Only identified particles
252               new((*plEMCAL)[indexNe++])       TParticle(*particle) ;
253             else if(pid[AliPID::kPhoton] > 0.75)
254               new((*plEMCAL)[indexNe++])       TParticle(*particle) ;       
255           }
256         else
257           {
258             Int_t pdg = 0;
259             if(fEMCALPID) 
260               {
261                 if( pid[AliPID::kPhoton] > 0.75) //This has to be fixen.
262                   pdg = 22;
263                 else if( pid[AliPID::kPi0] > 0.75)
264                   pdg = 111;
265               }
266             else
267               pdg = 22; //No PID, assume all photons
268             
269             TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1, 
270                                                  px, py, pz, en, v[0], v[1], v[2], 0);
271             AliDebug(4,Form("EMCAL clusters: pt %f, phi %f, eta %f", particle->Pt(),particle->Phi(),particle->Eta()));
272             
273             new((*plEMCAL)[indexNe++])       TParticle(*particle) ; 
274             new((*pl)[index++])           TParticle(*particle) ;
275           }
276       }
277     }
278
279     AliDebug(3,"Particle lists filled");
280     
281 }
282
283
284
285 //____________________________________________________________________________
286 void AliAnaGammaDirect::Exec(Option_t *) 
287 {
288   
289   // Processing of one event
290     
291   //Get ESDs
292   Long64_t entry = fChain->GetReadEntry() ;
293   
294   if (!fESD) {
295     AliError("fESD is not connected to the input!") ; 
296     return ; 
297   }
298   
299   if (fPrintInfo) 
300      AliInfo(Form("%s ----> Processing event # %lld",  (dynamic_cast<TChain *>(fChain))->GetFile()->GetName(), entry)) ; 
301
302   //CreateTLists with arrays of TParticles. Filled with particles only relevant for the analysis.
303
304   TClonesArray * particleList = new TClonesArray("TParticle",1000); // All particles refitted in CTS and detected in EMCAL (jet)
305   TClonesArray * plCTS         = new TClonesArray("TParticle",1000); // All particles refitted in Central Tracking System (ITS+TPC)
306   TClonesArray * plNe         = new TClonesArray("TParticle",1000);   // All particles measured in Jet Calorimeter (EMCAL)
307   TClonesArray * plPHOS     = new TClonesArray("TParticle",1000);  // All particles measured in PHOS as Gamma calorimeter
308   TClonesArray * plEMCAL   = new TClonesArray("TParticle",1000);  // All particles measured in EMCAL as Gamma calorimeter
309
310   TParticle *pGamma = new TParticle(); //It will contain the kinematics of the found prompt gamma
311   TParticle *pLeading = new TParticle(); //It will contain the kinematics of the found leading particle
312   
313   //Fill lists with photons, neutral particles and charged particles
314   //look for the highest energy photon in the event inside fCalorimeter
315
316   AliDebug(2, "Fill particle lists, get prompt gamma");
317
318   //Fill particle lists 
319   CreateParticleList(particleList, plCTS,plEMCAL,plPHOS); 
320
321   if(fCalorimeter == "PHOS")
322     plNe = plPHOS;
323   if(fCalorimeter == "EMCAL")
324     plNe = plEMCAL;
325
326
327   //_______________Analysis 1__________________________
328   //Look for the prompt photon in the selected calorimeter
329   if(fMakeICMethod < 3){
330     
331     Bool_t iIsInPHOSorEMCAL = kFALSE ; //To check if Gamma was in any calorimeter
332     //Search highest energy prompt gamma in calorimeter
333     GetPromptGamma(plNe,  plCTS, pGamma, iIsInPHOSorEMCAL) ; 
334     
335     AliDebug(1, Form("Is Gamma in %s? %d",fCalorimeter.Data(),iIsInPHOSorEMCAL));
336     
337     //If there is any photon candidate in fCalorimeter
338     if(iIsInPHOSorEMCAL){
339       if (fPrintInfo)
340         AliInfo(Form("Prompt Gamma: pt %f, phi %f, eta %f", pGamma->Pt(),pGamma->Phi(),pGamma->Eta())) ;
341       
342     }//Gamma in Calo
343     
344   }//Analysis 1
345
346
347   //_______________Analysis 2__________________________
348   //Look for the prompt photon in the selected calorimeter
349   //Isolation Cut Analysis for both methods and different pt cuts and cones
350   if(fMakeICMethod == 3){
351     
352     for(Int_t ipr = 0; ipr < plNe->GetEntries() ; ipr ++ ){
353       TParticle * pCandidate = dynamic_cast<TParticle *>(plNe->At(ipr)) ;
354       
355       if(pCandidate->Pt() > fMinGammaPt){
356         
357         Bool_t  icPtThres   = kFALSE;
358         Bool_t  icPtSum     = kFALSE;
359         
360         Float_t ptC             = pCandidate->Pt() ;
361         
362         fhPtCandidate->Fill(ptC);
363         
364         for(Int_t icone = 0; icone<fNCones; icone++){
365           fConeSize  = fConeSizes[icone] ;
366           Float_t coneptsum = 0 ;
367           for(Int_t ipt = 0; ipt<fNPtThres;ipt++){ 
368             fPtThreshold   =  fPtThresholds[ipt] ;
369             MakeIsolationCut(plCTS,plNe, pCandidate, ipr, icPtThres, icPtSum,coneptsum);
370             AliDebug(4,Form("Candidate pt %f, pt in cone %f, Isolated? ICPt %d, ICSum %d",
371                             pCandidate->Pt(), coneptsum, icPtThres, icPtSum));
372             
373             fhPtThresIsolated[icone][ipt]->Fill(ptC); 
374           }//pt thresh loop
375           fhPtSumIsolated[icone]->Fill(ptC,coneptsum) ;
376         }//cone size loop
377       }//min pt candidate
378     }//candidate loop
379   }//Analysis 2
380   
381   AliDebug(2, "End of analysis, delete pointers");
382   
383   particleList->Delete() ; 
384   plCTS->Delete() ;
385   plNe->Delete() ;
386   plPHOS->Delete() ;
387   pLeading->Delete();
388   pGamma->Delete();
389
390   delete plNe ;
391   delete plCTS ;
392   delete particleList ;
393   //  delete pLeading;
394   //  delete pGamma;
395
396   PostData(0, fOutputContainer);
397 }    
398
399
400 //____________________________________________________________________________
401 void AliAnaGammaDirect::GetPromptGamma(TClonesArray * pl, TClonesArray * plCTS, TParticle *pGamma, Bool_t &Is) const 
402 {
403   //Search for the prompt photon in Calorimeter with pt > fMinGammaPt
404
405   Double_t pt = 0;
406   Int_t index = -1; 
407   for(Int_t ipr = 0;ipr < pl->GetEntries() ; ipr ++ ){
408     TParticle * particle = dynamic_cast<TParticle *>(pl->At(ipr)) ;
409
410     if((particle->Pt() > fMinGammaPt) && (particle->Pt() > pt)){
411       index = ipr ;
412       pt = particle->Pt();
413       pGamma->SetMomentum(particle->Px(),particle->Py(),particle->Pz(),particle->Energy());
414       AliDebug(4,Form("Cluster in calo: pt %f, phi %f, eta %f", pGamma->Pt(),pGamma->Phi(),pGamma->Eta())) ;
415       Is  = kTRUE;
416     }
417   }
418
419   //Do Isolation?
420   if(fMakeICMethod && Is)
421     {
422       Float_t coneptsum = 0 ;
423       Bool_t  icPtThres   = kFALSE;
424       Bool_t  icPtSum     = kFALSE;
425       MakeIsolationCut(plCTS,pl, pGamma, index, 
426                        icPtThres, icPtSum,coneptsum);
427       if(fMakeICMethod == 1) //Pt thres method
428         Is = icPtThres ;
429       if(fMakeICMethod == 2) //Pt cone sum method
430         Is = icPtSum ;
431     }
432   
433   if(Is){
434     AliDebug(3,Form("Cluster with p_{T} larger than %f found in calorimeter ", fMinGammaPt)) ;
435     AliDebug(3,Form("Gamma: pt %f, phi %f, eta %f", pGamma->Pt(),pGamma->Phi(),pGamma->Eta())) ;
436     //Fill prompt gamma histograms
437     fhNGamma  ->Fill(pGamma->Pt());
438     fhPhiGamma->Fill( pGamma->Pt(),pGamma->Phi());
439     fhEtaGamma->Fill(pGamma->Pt(),pGamma->Eta());
440   }
441   else
442     AliDebug(1,Form("NO Cluster with pT larger than %f found in calorimeter ", fMinGammaPt)) ;
443 }
444
445   //____________________________________________________________________________
446 void AliAnaGammaDirect::Init(const Option_t * )
447 {
448   // Initialisation of branch container 
449  
450   AliDebug(2,Form("*** Initialization of %s", GetName())) ; 
451   
452   // Get input data
453   fChain = dynamic_cast<TChain *>(GetInputData(0)) ;
454   if (!fChain) {
455     AliError(Form("Input 0 for %s not found\n", GetName()));
456     return ;
457   }
458   
459   if (!fESD) {
460     // One should first check if the branch address was taken by some other task
461     char ** address = (char **)GetBranchAddress(0, "ESD") ;
462     if (address) 
463       fESD = (AliESD *)(*address) ; 
464     if (!fESD) 
465       fChain->SetBranchAddress("ESD", &fESD) ;  
466   }
467   // The output objects will be written to 
468   TDirectory * cdir = gDirectory ; 
469   // Open a file for output #0
470   char outputName[1024] ; 
471   sprintf(outputName, "%s.root", GetName() ) ; 
472   OpenFile(0, outputName , "RECREATE") ; 
473   if (cdir) 
474     cdir->cd() ; 
475
476  //  //Initialize the parameters of the analysis.
477   fCalorimeter="PHOS";
478   fPrintInfo           = kTRUE;
479   fMinGammaPt  = 10. ;
480
481   //Fill particle lists when PID is ok
482   fEMCALPID = kFALSE;
483   fPHOSPID = kFALSE;
484
485   fConeSize             = 0.2 ; 
486   fPtThreshold         = 2.0; 
487   fPtSumThreshold  = 1.; 
488
489   fNCones           = 4 ; 
490   fNPtThres         = 4 ; 
491   fConeSizes[0] = 0.1; fConeSizes[0] = 0.2; fConeSizes[2] = 0.3; fConeSizes[3] = 0.4;
492   fPtThresholds[0]=1.; fPtThresholds[0]=2.; fPtThresholds[0]=3.; fPtThresholds[0]=4.;
493
494   fMakeICMethod = 1; // 0 don't isolate, 1 pt thresh method, 2 cone pt sum method, 3 make isolation study
495
496   //Initialization of histograms 
497   MakeHistos() ;
498 }
499
500 //__________________________________________________________________
501 void  AliAnaGammaDirect::MakeIsolationCut(TClonesArray * plCTS, 
502                                            TClonesArray * plNe, 
503                                            TParticle * pCandidate, 
504                                            Int_t index, 
505                                            Bool_t  &icmpt,  Bool_t  &icms, 
506                                            Float_t &coneptsum) const 
507 {  
508   //Search in cone around a candidate particle if it is isolated 
509   Float_t phiC  = pCandidate->Phi() ;
510   Float_t etaC = pCandidate->Eta() ;
511   Float_t pt     = -100. ;
512   Float_t eta   = -100.  ;
513   Float_t phi    = -100.  ;
514   Float_t rad   = -100 ;
515   Int_t    n        = 0 ;
516   TParticle * particle  = new TParticle();
517
518   coneptsum = 0; 
519   icmpt = kFALSE;
520   icms  = kFALSE;
521
522   //Check charged particles in cone.
523   for(Int_t ipr = 0;ipr < plCTS->GetEntries() ; ipr ++ ){
524     particle = dynamic_cast<TParticle *>(plCTS->At(ipr)) ;
525     pt    = particle->Pt();
526     eta  = particle->Eta();
527     phi  = particle->Phi() ;
528     
529     //Check if there is any particle inside cone with pt larger than  fPtThreshold
530     rad = TMath::Sqrt(TMath::Power((eta-etaC),2) +
531                       TMath::Power((phi-phiC),2));
532     if(rad<fConeSize){
533       AliDebug(3,Form("charged in cone pt %f, phi %f, eta %f, R %f ",pt,phi,eta,rad));
534       coneptsum+=pt;
535       if(pt > fPtThreshold ) n++;
536     }
537   }// charged particle loop
538   
539   //Check neutral particles in cone.
540   for(Int_t ipr = 0;ipr < plNe->GetEntries() ; ipr ++ ){
541     if(ipr != index){//Do not count the candidate
542       particle = dynamic_cast<TParticle *>(plNe->At(ipr)) ;
543       pt    = particle->Pt();
544       eta  = particle->Eta();
545       phi  = particle->Phi() ;
546       
547       //Check if there is any particle inside cone with pt larger than  fPtThreshold
548       rad = TMath::Sqrt(TMath::Power((eta-etaC),2) +
549                         TMath::Power((phi-phiC),2));
550       if(rad<fConeSize){
551         AliDebug(3,Form("charged in cone pt %f, phi %f, eta %f, R %f ",pt,phi,eta,rad));
552         coneptsum+=pt;
553         if(pt > fPtThreshold ) n++;
554       }
555     }
556   }// neutral particle loop
557   
558   if(n == 0) 
559     icmpt =  kTRUE ;
560   if(coneptsum < fPtSumThreshold)
561     icms  =  kTRUE ;
562
563 }
564
565 //___________________________________________________________________
566 void AliAnaGammaDirect::MakeHistos()
567 {
568   // Create histograms to be saved in output file and 
569   // stores them in fOutputContainer
570   
571   fOutputContainer = new TObjArray(10000) ;
572
573   //Histograms of highest gamma identified in Event
574   fhNGamma  = new TH1F("NGamma","Number of #gamma over PHOS",240,0,120); 
575   fhNGamma->SetYTitle("N");
576   fhNGamma->SetXTitle("p_{T #gamma}(GeV/c)");
577   fOutputContainer->Add(fhNGamma) ; 
578   
579   fhPhiGamma  = new TH2F
580     ("PhiGamma","#phi_{#gamma}",200,0,120,200,0,7); 
581   fhPhiGamma->SetYTitle("#phi");
582   fhPhiGamma->SetXTitle("p_{T #gamma} (GeV/c)");
583   fOutputContainer->Add(fhPhiGamma) ; 
584   
585   fhEtaGamma  = new TH2F
586     ("EtaGamma","#phi_{#gamma}",200,0,120,200,-0.8,0.8); 
587   fhEtaGamma->SetYTitle("#eta");
588   fhEtaGamma->SetXTitle("p_{T #gamma} (GeV/c)");
589   fOutputContainer->Add(fhEtaGamma) ;
590
591   if( fMakeICMethod== 3 ){
592
593     //Isolation cut histograms
594     fhPtCandidate  = new TH1F
595       ("PtCandidate","p_{T} of candidate particles for isolation",240,0,120); 
596     fhPtCandidate->SetXTitle("p_{T} (GeV/c)");
597     fOutputContainer->Add(fhPtCandidate) ;
598     
599     char name[128];
600     char title[128];
601     for(Int_t icone = 0; icone<fNCones; icone++){
602       sprintf(name,"PtSumIsolated_Cone_%d",icone);
603       sprintf(title,"Candidate cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
604       fhPtSumIsolated[icone]  = new TH2F(name, title,240,0,120,120,0,10);
605       fhPtSumIsolated[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
606       fhPtSumIsolated[icone]->SetXTitle("p_{T} (GeV/c)");
607       fOutputContainer->Add(fhPtSumIsolated[icone]) ; 
608       
609       for(Int_t ipt = 0; ipt<fNPtThres;ipt++){ 
610         sprintf(name,"PtThresIsol_Cone_%d_Pt%d",icone,ipt);
611         sprintf(title,"Isolated candidate p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
612         fhPtThresIsolated[icone][ipt]  = new TH1F(name, title,240,0,120);
613         fhPtThresIsolated[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
614         fOutputContainer->Add(fhPtThresIsolated[icone][ipt]) ; 
615       }//icone loop
616   }//ipt loop
617   }
618 }
619
620 void AliAnaGammaDirect::Print(const Option_t * opt) const
621 {
622
623   //Print some relevant parameters set for the analysis
624   if(! opt)
625     return;
626
627   Info("Print", "%s %s", GetName(), GetTitle() ) ;
628   printf("Cone Size               =     %f\n", fConeSize) ; 
629   printf("pT threshold           =     %f\n", fPtThreshold) ;
630   printf("pT sum threshold   =     %f\n", fPtSumThreshold) ; 
631   printf("Min Gamma pT      =     %f\n", fMinGammaPt) ; 
632   printf("Calorimeter            =     %s\n", fCalorimeter.Data()) ; 
633
634
635 void AliAnaGammaDirect::Terminate(Option_t *)
636 {
637    // The Terminate() function is the last function to be called during
638    // a query. It always runs on the client, it can be used to present
639    // the results graphically or save the results to file.
640     
641
642 }