Compilation on Windows/Cygwin
[u/mrichter/AliRoot.git] / PWG4 / AliAnaGamma.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.4  2007/10/29 13:48:42  gustavo
21  * Corrected coding violations
22  *
23  * Revision 1.2  2007/08/17 12:40:04  schutz
24  * New analysis classes by Gustavo Conesa
25  *
26  * Revision 1.1.2.1  2007/07/26 10:32:09  schutz
27  * new analysis classes in the the new analysis framework
28  *
29  *
30  */
31
32 //_________________________________________________________________________
33 // Base class for gamma and correlation analysis
34 // It is called by the task class AliAnalysisGammaTask and it connects the input (ESD/AOD/MonteCarlo)
35 // got with AliGammaReader (produces TClonesArrays of TParticles), with the analysis classes 
36 // AliAnaGammaDirect, AliAnaGammaCorrelation ....
37 //
38 //*-- Author: Gustavo Conesa (INFN-LNF)
39
40 // --- ROOT system ---
41
42 #include <TParticle.h>
43 #include <TH2.h>
44
45 //---- AliRoot system ---- 
46 #include "AliAnaGamma.h" 
47 #include "AliGammaReader.h" 
48 #include "AliAnaGammaDirect.h" 
49 #include "AliAnaGammaCorrelation.h" 
50 #include "AliAnaGammaSelection.h" 
51 #include "AliNeutralMesonSelection.h"
52 #include "AliAODCaloCluster.h"
53 #include "AliAODEvent.h"
54 #include "Riostream.h"
55 #include "AliLog.h"
56
57 ClassImp(AliAnaGamma)
58
59
60 //____________________________________________________________________________
61   AliAnaGamma::AliAnaGamma() : 
62     TObject(),
63     fOutputContainer(0x0), 
64     fAnaType(0),  fCalorimeter(0), fData(0x0), fKine(0x0), 
65     fReader(0x0), fGammaDirect(0x0), fGammaCorrelation(0x0), fGammaSelection(0x0),
66     fNeutralMesonSelection(0x0), fAODclusters(0x0), fNAODclusters(0)
67 {
68   //Default Ctor
69
70   //Initialize parameters, pointers and histograms
71   if(!fReader)
72     fReader = new AliGammaReader();
73   if(!fGammaDirect)
74     fGammaDirect = new AliAnaGammaDirect();
75   if(!fGammaCorrelation)
76     fGammaCorrelation = new AliAnaGammaCorrelation();
77   if(!fGammaSelection)
78     fGammaSelection = new AliAnaGammaSelection();
79   if(!fNeutralMesonSelection)
80     fNeutralMesonSelection = new AliNeutralMesonSelection();
81
82   fAODclusters = 0;
83
84   InitParameters();
85   
86 }
87
88 //____________________________________________________________________________
89 AliAnaGamma::AliAnaGamma(const AliAnaGamma & g) :   
90   TObject(),
91   fOutputContainer(g. fOutputContainer), 
92   fAnaType(g.fAnaType),  fCalorimeter(g.fCalorimeter), 
93   fData(g.fData), fKine(g.fKine),fReader(g.fReader),
94   fGammaDirect(g.fGammaDirect), fGammaCorrelation(g.fGammaCorrelation),
95   fGammaSelection(g.fGammaSelection),
96   fNeutralMesonSelection(g.fNeutralMesonSelection),  
97   fAODclusters(g. fAODclusters), fNAODclusters(g.fNAODclusters)
98 {
99   // cpy ctor
100   
101 }
102
103 //_________________________________________________________________________
104 AliAnaGamma & AliAnaGamma::operator = (const AliAnaGamma & source)
105 {
106   // assignment operator
107
108   if(this == &source)return *this;
109   ((TObject *)this)->operator=(source);
110
111   fOutputContainer = source.fOutputContainer ;
112   fAnaType = source.fAnaType;
113   fCalorimeter = source.fCalorimeter ;
114   fData = source.fData ; 
115   fKine = source.fKine ;
116   fReader = source.fReader ;
117   fGammaDirect = source.fGammaDirect ;
118   fGammaCorrelation = source.fGammaCorrelation ;
119   fGammaSelection = source.fGammaSelection ;
120   fNeutralMesonSelection = source.fNeutralMesonSelection ;
121
122   return *this;
123
124 }
125
126 //____________________________________________________________________________
127 AliAnaGamma::~AliAnaGamma() 
128 {
129   // Remove all pointers.
130
131   // Protection added in case of NULL pointers (MG)
132   if (fOutputContainer) {
133      fOutputContainer->Clear();
134      delete fOutputContainer ;
135   }   
136
137   if (fData) delete fData ; 
138   if (fKine) delete fKine ;
139   if (fReader) delete fReader ;
140   if (fGammaDirect) delete fGammaDirect ;
141   if (fGammaCorrelation) delete fGammaCorrelation ;
142   if (fGammaSelection) delete fGammaSelection ;
143   if (fNeutralMesonSelection) delete fNeutralMesonSelection ;
144
145 }
146
147 //________________________________________________________________________
148 void AliAnaGamma::Init()
149 {  
150
151   //Init container histograms and other common variables
152
153   //Histograms container
154   fOutputContainer = new TList ;
155   
156   //Fill container with appropriate histograms
157   
158   //Set selection  analysis histograms
159   TList * selectcontainer =  fGammaSelection->GetCreateOutputObjects(); 
160   for(Int_t i = 0; i < selectcontainer->GetEntries(); i++){
161     Bool_t  add = kTRUE ;
162     TString name = (selectcontainer->At(i))->GetName();   
163     if(!fReader->IsEMCALOn() && name.Contains("EMCAL")) add = kFALSE;
164     if(!fReader->IsPHOSOn() && name.Contains("PHOS"))   add = kFALSE;
165     if(!fReader->IsCTSOn() &&  !fGammaSelection->FillCTS() && name.Contains("CTS"))   add = kFALSE;
166     if(add) fOutputContainer->Add(selectcontainer->At(i)) ;
167   }  //Set selection  analysis histograms
168  
169   
170   //Set prompt photon analysis histograms
171   TList * promptcontainer =  fGammaDirect->GetCreateOutputObjects(); 
172   for(Int_t i = 0; i < promptcontainer->GetEntries(); i++)
173     fOutputContainer->Add(promptcontainer->At(i)) ;
174   
175   //Check if selected options are correct or set them when necessary
176   if(fReader->GetDataType() == AliGammaReader::kMCData){
177     fGammaDirect->SetMC();//Only useful with AliGammaMCDataReader, by default kFALSE
178     fGammaSelection->SetMC();//Only useful with AliGammaMCDataReader, by default kFALSE
179   }
180   
181   if(fAnaType == kCorrelation){
182     
183     //Check if selected options are correct
184     if (fGammaDirect->GetICMethod()==AliAnaGammaDirect::kSeveralIC)
185       AliFatal("Correlation not allowed with multiple isolation cuts, kCorrelation and kSeveralIC do not go together");
186     
187     if(fGammaCorrelation->GetCorrelationType() == AliAnaGammaCorrelation::kParton && 
188        fReader->GetDataType() != AliGammaReader::kMC)
189       AliFatal("kParton must be analyzed with data kMC");
190     
191     //Set the parameters for the neutral pair selection depending on the analysis, 
192     fNeutralMesonSelection->SetDeltaPhiCutRange(fGammaCorrelation->GetDeltaPhiMinCut(), 
193                                                 fGammaCorrelation->GetDeltaPhiMaxCut());  
194     
195
196     if(fGammaCorrelation->GetCorrelationType() == AliAnaGammaCorrelation::kHadron){
197       fNeutralMesonSelection->SetPhiPtSelection(AliNeutralMesonSelection::kSelectPhiMinPt);
198       fNeutralMesonSelection->SetMinPt(fGammaCorrelation->GetMinPtHadron());
199       
200     }
201     
202     if(fGammaCorrelation->GetCorrelationType() == AliAnaGammaCorrelation::kJetLeadCone){
203       fNeutralMesonSelection->SetPhiPtSelection(AliNeutralMesonSelection::kSelectPhiPtRatio);
204       fNeutralMesonSelection->SetRatioCutRange(fGammaCorrelation->GetRatioMinCut(), 
205                                                fGammaCorrelation->GetRatioMaxCut());
206     }
207
208     //Set the neutral mesosn selection histograms
209     TList * neutralmesoncontainer =  fNeutralMesonSelection->GetCreateOutputObjects();
210     if(fNeutralMesonSelection->AreNeutralMesonSelectionHistosKept()){
211       for(Int_t i = 0; i < neutralmesoncontainer->GetEntries(); i++)
212         fOutputContainer->Add(neutralmesoncontainer->At(i)) ;
213     }
214     
215     //Set correlation histograms
216     TList * correlationcontainer =  fGammaCorrelation->GetCreateOutputObjects();
217     for(Int_t i = 0; i < correlationcontainer->GetEntries(); i++)
218       fOutputContainer->Add(correlationcontainer->At(i)) ;
219     fGammaCorrelation->SetOutputContainer(fOutputContainer);
220     fGammaCorrelation->SetNeutralMesonSelection(fNeutralMesonSelection);
221  
222   }//kCorrelation  
223   
224 }
225
226 //____________________________________________________________________________
227 void AliAnaGamma::InitParameters()
228 {
229
230   //Init data members
231  
232   fAnaType = kPrompt;
233   fCalorimeter = "EMCAL";
234
235 }
236
237 //__________________________________________________________________
238 void AliAnaGamma::Print(const Option_t * opt) const
239 {
240
241   //Print some relevant parameters set for the analysis
242   if(! opt)
243     return;
244   
245   Info("Print", "%s %s", GetName(), GetTitle() ) ;
246   printf("Analysis type           =     %d\n", fAnaType) ;
247   printf("Calorimeter           =     %s\n", fCalorimeter.Data()) ;
248
249   switch(fAnaType)
250     {
251     case kPrompt:
252       {
253         fGammaDirect->Print("");
254       }// case kIsolationCut
255       break;
256    
257     case kCorrelation:
258       {
259         fGammaCorrelation->Print("");
260       }//  case kCorrelation
261       break;
262       
263     }//switch
264
265
266
267 //____________________________________________________________________________
268 Bool_t AliAnaGamma::ProcessEvent(){
269   //Process analysis for this event
270
271   if(!fOutputContainer)
272     AliFatal("Histograms not initialized");
273
274   //CreateTLists with arrays of TParticles. Filled with particles only relevant for the analysis.
275
276   TClonesArray * plCTS      = new TClonesArray("TParticle",1000); // All particles refitted in Central Tracking System (ITS+TPC)
277   TClonesArray * plEMCAL    = new TClonesArray("TParticle",1000);   // All particles measured in Jet Calorimeter (EMCAL)
278   TClonesArray * plPHOS     = new TClonesArray("TParticle",1000);  // All particles measured  Gamma calorimeter
279   TClonesArray * plPrimCTS      = new TClonesArray("TParticle",1000); // primary tracks
280   TClonesArray * plPrimEMCAL    = new TClonesArray("TParticle",1000);   // primary emcal clusters
281   TClonesArray * plPrimPHOS     = new TClonesArray("TParticle",1000);  // primary phos clusters
282   TClonesArray * plParton   = new TClonesArray("TParticle",1000);  // All partons
283   //Fill lists with photons, neutral particles and charged particles
284   //look for the highest energy photon in the event inside fCalorimeter
285     
286   //Fill particle lists 
287   if(fReader->GetDataType() == AliGammaReader::kData){
288     AliDebug(1,"Data analysis");
289     fReader->CreateParticleList(fData, NULL,plCTS,plEMCAL,plPHOS,NULL,NULL,NULL); 
290   }
291   else if( fReader->GetDataType()== AliGammaReader::kMC){
292     AliDebug(1,"Kinematics analysis");
293     fReader->CreateParticleList(fKine, NULL,plCTS,plEMCAL,plPHOS,plParton,NULL,NULL); 
294   }
295   else if(fReader->GetDataType() == AliGammaReader::kMCData) {
296    AliDebug(1,"Data + Kinematics analysis");
297    fReader->CreateParticleList(fData, fKine,plCTS,plEMCAL,plPHOS,plPrimCTS,plPrimEMCAL,plPrimPHOS); 
298   }
299   else
300     AliError("Option not implemented");
301
302   //Fill AOD with calorimeter
303   //Temporal solution, just for testing
304   //FillAODs(plPHOS,plEMCAL);  
305
306   //Select particles to do the final analysis.
307   if(fReader->IsEMCALOn()) fGammaSelection->Selection("EMCAL",plEMCAL,plPrimEMCAL);
308   if(fReader->IsPHOSOn())  fGammaSelection->Selection("PHOS",plPHOS,plPrimPHOS);
309   if(fReader->IsCTSOn() &&  fGammaSelection->FillCTS()) fGammaSelection->Selection("CTS",plCTS,plPrimCTS);
310
311
312   //Search highest energy prompt gamma in calorimeter
313   if(fCalorimeter == "PHOS")
314     MakeAnalysis(plPHOS, plEMCAL, plCTS, plParton, plPrimPHOS) ; 
315   else if (fCalorimeter == "EMCAL")
316     MakeAnalysis(plEMCAL, plPHOS, plCTS,plParton, plPrimEMCAL) ; 
317   else
318     AliFatal("Wrong calorimeter name");
319
320   plCTS->Clear() ;
321   plEMCAL->Clear() ;
322   plPHOS->Clear() ;
323   plParton->Clear() ;
324   plPrimCTS->Clear() ;
325   plPrimEMCAL->Clear() ;
326   plPrimPHOS->Clear() ;
327
328   delete plCTS ;
329   delete plPHOS ;
330   delete plEMCAL ;
331   delete plParton ;
332   delete plPrimCTS ;
333   delete plPrimPHOS ;
334   delete plPrimEMCAL ;
335
336   return kTRUE;
337
338 }
339
340 //____________________________________________________________________________
341 void AliAnaGamma::MakeAnalysis(TClonesArray * plCalo, TClonesArray * plNe, TClonesArray * plCTS, TClonesArray * plParton, TClonesArray * plPrimCalo)  {
342   
343   TParticle * pGamma = new TParticle ;
344   Bool_t isInCalo = kFALSE ;
345   
346   switch(fAnaType)
347     {
348       
349       //Only Prompt photon analysis
350     case kPrompt:
351       {   
352         AliDebug(1,"kPrompt analysis");
353         switch(fGammaDirect->GetICMethod())
354           {
355             
356           case AliAnaGammaDirect::kSeveralIC:
357             {
358               fGammaDirect->MakeSeveralICAnalysis(plCalo, plCTS);
359               AliDebug(1,"kSeveralIC analysis");
360             }
361             break;
362             
363           default :
364             {
365               fGammaDirect->GetPromptGamma(plCalo, plCTS,plPrimCalo, pGamma,isInCalo);
366               if(!isInCalo)
367                 AliDebug(1,"Prompt gamma not found");
368             }
369             break;
370           }//IC method
371       }// case kPrompt:
372       break;
373       
374       //Correlate prompt photon with something: parton, hadron, jet.
375     case kCorrelation:
376       {
377         AliDebug(1,"kCorrelation analysis");
378         //Find prompt photon    
379         fGammaDirect->GetPromptGamma(plCalo, plCTS,plPrimCalo, pGamma,isInCalo);
380
381         if(isInCalo){//If prompt photon found, do correlation
382           
383           switch(fGammaCorrelation->GetCorrelationType())
384             {
385             case AliAnaGammaCorrelation::kParton:
386               {
387                 AliDebug(1,"kParton correlation");
388                 fGammaCorrelation->MakeGammaCorrelation(pGamma, plParton, NULL);
389               }//  case kParton
390               break;
391       
392             case AliAnaGammaCorrelation::kHadron:
393               {
394                 AliDebug(1,"kHadron correlation");
395                 fGammaCorrelation->MakeGammaCorrelation(pGamma, plCTS, plNe);
396               }//  case kHadron
397               break;
398
399             case AliAnaGammaCorrelation::kJetLeadCone:
400               {         
401                 AliDebug(1,"kJetLeadCone correlation");
402                 fGammaCorrelation->MakeGammaCorrelation(pGamma, plCTS, plNe);
403               }//  case kJetLeadCone
404               break;
405               
406             case AliAnaGammaCorrelation::kJetFinder:
407               { 
408                 AliDebug(1,"kJetFinder correlation");
409                 printf("Analysis not implemented \n");
410               }//  case kJetFinder
411               break;
412             }// switch correlation
413         }// is in calo
414         else  AliDebug(2,"Prompt gamma not found");
415       }// case kCorrelation
416       break;
417
418     } //switch(fAnaType)
419
420   delete pGamma ; 
421   
422 }
423
424 //____________________________________________________
425 void AliAnaGamma::AddCluster(AliAODCaloCluster p)
426 {
427 // Add new jet to the list
428   cout<<"AOD list pointer "<<fAODclusters<<" nAODclusters "<<fNAODclusters<<endl;
429   cout<<"list entries "<<fAODclusters->GetEntries()<<endl;
430   new ((*fAODclusters)[fNAODclusters++]) AliAODCaloCluster(p);
431 }
432
433 //___________________________________________________
434 void AliAnaGamma::ConnectAOD(AliAODEvent* aod)
435 {
436 // Connect to the AOD
437   fAODclusters = aod->GetCaloClusters();
438 }
439
440 //____________________________________________________________________________
441 void AliAnaGamma::FillAODs(TClonesArray * plPHOS, TClonesArray * plEMCAL){
442   //Fill AOD caloClusters
443   //Temporal method, just for testing AOD creation
444   
445   //Fill AOD with PHOS clusters
446   Int_t nphos =  plPHOS->GetEntries() ;
447   cout<<"PHOS entries "<<nphos<<endl;
448   for(Int_t ipr = 0;ipr < nphos ; ipr ++ ){
449     TParticle * particle = dynamic_cast<TParticle *>(plPHOS->At(ipr)) ;
450     Float_t pos []= {0,0,0};
451     AliAODCaloCluster phos(ipr,0,0x0,particle->Energy(),pos,0x0,AliAODCluster::kPHOSNeutral,0);
452     AddCluster(phos);
453   }
454   
455   //Fill AOD with EMCAL clusters
456   cout<<"EMCAL entries "<<plEMCAL->GetEntries()<<endl;
457   for(Int_t ipr = 0;ipr < plEMCAL->GetEntries() ; ipr ++ ){
458     TParticle * particle = dynamic_cast<TParticle *>(plEMCAL->At(ipr)) ;
459     Float_t pos []= {0,0,0};
460     AliAODCaloCluster emcal(ipr+nphos,0,0x0,particle->Energy(),pos,0x0,AliAODCluster::kEMCALClusterv1,0);
461     AddCluster(emcal);
462   }
463 }