2110dd28874ec0eb32955074033e62209a10227f
[u/mrichter/AliRoot.git] / PWG4 / macros / AddTaskPi0EMCALPbPb.C
1
2 AliAnalysisTaskParticleCorrelation *AddTaskPi0EMCALPbPb(
3                                                     const TString data          = "AOD",
4                                                     const TString calorimeter   = "EMCAL", 
5                                                     const Bool_t  printSettings = kFALSE,
6                                                     const Bool_t  simulation    = kFALSE,
7                                                     const Bool_t  outputAOD     = kFALSE, 
8                                                     const TString outputfile    = "",
9                                                     const Int_t   year          = 2011,
10                                                     const TString col           = "PbPb",
11                                                     const TString trig          = "",
12                                                     const TString clustersArray = "" 
13                                                     )
14 {
15   // Creates a PartCorr task, configures it and adds it to the analysis manager.
16     
17   // Get the pointer to the existing analysis manager via the static access method.
18   
19   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
20   if (!mgr) {
21     ::Error("AddTask", "No analysis manager to connect to.");
22     return NULL;
23   }  
24   
25   // Check the analysis type using the event handlers connected to the analysis manager.
26   
27   if (!mgr->GetInputEventHandler()) {
28     ::Error("AddTask", "This task requires an input event handler");
29     return NULL;
30   }
31   Bool_t kInputDataType = "AOD";
32   if(!data.Contains("delta"))
33     kInputDataType = mgr->GetInputEventHandler()->GetDataType(); // can be "ESD" or "AOD"
34   
35   Bool_t kUseKinematics = kFALSE; 
36   if(simulation) { 
37     kUseKinematics = (mgr->GetMCtruthEventHandler())?kTRUE:kFALSE; 
38     if (!kUseKinematics && data=="AOD" && kInputDataType != "ESD") kUseKinematics = kTRUE; //AOD primary should be available ... 
39   } 
40   
41   //cout<<"********* ACCESS KINE? "<<kUseKinematics<<endl;
42   
43   // #### Configure analysis ####
44     
45   AliAnaPartCorrMaker * maker = new AliAnaPartCorrMaker();
46   
47   // General frame setting and configuration
48   maker->SetReader   (ConfigureReader(data,kInputDataType,calorimeter,clustersArray,col,outputAOD,kUseKinematics,printSettings)   ); 
49   maker->SetCaloUtils(ConfigureCaloUtils(year,col,clustersArray, simulation, printSettings)); 
50   
51   // Analysis tasks setting and configuration
52   Int_t n = 0;//Analysis number, order is important
53   maker->AddAnalysis(ConfigurePhotonAnalysis(data,calorimeter,year,col,clustersArray,simulation,"",printSettings), n++); // Photon cluster selection
54   maker->AddAnalysis(ConfigurePi0Analysis   (data,calorimeter,year,col,clustersArray,simulation,"",printSettings), n++); // Pi0 invariant mass analysis
55   
56   maker->SetAnaDebug(-1)  ;
57   maker->SwitchOnHistogramsMaker()  ;
58   if(data.Contains("delta")) maker->SwitchOffAODsMaker() ;
59   else                       maker->SwitchOnAODsMaker()  ;
60         
61   if(printSettings) maker->Print("");
62   
63   //printf("<< End Configuration of %d analysis for calorimeter %s >>\n",n, calorimeter.Data());
64   
65   // Create task
66   
67   AliAnalysisTaskParticleCorrelation * task = new AliAnalysisTaskParticleCorrelation (Form("PartCorr%s_Trig%s_Cl%s",calorimeter.Data(),trig.Data(),clustersArray.Data()));
68   task->SetConfigFileName(""); //Don't configure the analysis via configuration file.
69   //task->SetDebugLevel(-1);
70   task->SetBranches("ESD:AliESDRun.,AliESDHeader"); //just a trick to get Constantin's analysis to work
71   task->SetAnalysisMaker(maker);
72   mgr->AddTask(task);
73   
74   //Create containers
75   
76   if(outputfile.Length()==0)outputfile = AliAnalysisManager::GetCommonFileName(); 
77   TString name(Form("%s_Trig%s_Cl%s",calorimeter.Data(),trig.Data(),clustersArray.Data()));
78   AliAnalysisDataContainer *cout_pc   = mgr->CreateContainer(name.Data(), TList::Class(), 
79                                                              AliAnalysisManager::kOutputContainer, 
80                                                              Form("%s:%s",outputfile.Data(),name.Data()));
81         
82   AliAnalysisDataContainer *cout_cuts = mgr->CreateContainer(Form("%sCuts",name.Data()), TList::Class(), 
83                                                              AliAnalysisManager::kParamContainer, 
84                                                              Form("%s:%sCuts",outputfile.Data(),name.Data()));
85
86   // Create ONLY the output containers for the data produced by the task.
87   // Get and connect other common input/output containers via the manager as below
88   //==============================================================================
89   mgr->ConnectInput  (task, 0, mgr->GetCommonInputContainer());
90   // AOD output slot will be used in a different way in future
91   if(!data.Contains("delta")   && outputAOD) mgr->ConnectOutput (task, 0, mgr->GetCommonOutputContainer());
92   mgr->ConnectOutput (task, 1, cout_pc);
93   mgr->ConnectOutput (task, 2, cout_cuts);
94     
95   return task;
96 }
97
98 //____________________________________
99 AliCaloTrackReader * ConfigureReader(TString kData,TString kInputDataType, TString kCalorimeter, 
100                                      TString kClusterArray, TString kCollisions,
101                                      Bool_t kOutputAOD, Bool_t kUseKinematics, Bool_t kPrint)
102 {
103   
104   AliCaloTrackReader * reader = 0;
105   if     (kData.Contains("AOD"))   reader = new AliCaloTrackAODReader();
106   else if(kData=="ESD")            reader = new AliCaloTrackESDReader();
107   else if(kData=="MC" && 
108           kInputDataType == "ESD") reader = new AliCaloTrackMCReader();
109   
110   reader->SetDebug(-1);//10 for lots of messages
111   
112   //Delta AOD?
113   //reader->SetDeltaAODFileName("");
114   if(kOutputAOD) reader->SwitchOnWriteDeltaAOD()  ;
115   
116   // MC settings
117   if(kUseKinematics){
118     if(kInputDataType == "ESD"){
119       reader->SwitchOnStack();          
120       reader->SwitchOffAODMCParticles(); 
121     }
122     else if(kInputDataType == "AOD"){
123       reader->SwitchOffStack();          
124       reader->SwitchOnAODMCParticles(); 
125     }
126   }  
127   
128   //------------------------
129   // Detector input filling
130   //------------------------
131   
132   //Min cluster/track E
133   reader->SetEMCALPtMin(0.5); 
134   reader->SetEMCALPtMax(1000); 
135   reader->SetPHOSPtMin(0.3);
136   reader->SetPHOSPtMax(1000);
137   reader->SetCTSPtMin(0.1);
138   reader->SetCTSPtMax(1000);
139   
140   reader->SwitchOffFiducialCut();
141   
142   // Tracks
143   reader->SwitchOffCTS();
144   
145   //gROOT->LoadMacro("$ALICE_ROOT/PWG4/macros/CreateTrackCutsPWG4.C"); 
146   //AliESDtrackCuts * esdTrackCuts = CreateTrackCutsPWG4(10011004); 
147   //reader->SetTrackCuts(esdTrackCuts);
148   
149   // Calorimeter
150   
151   reader->SetEMCALClusterListName(kClusterArray);
152 //  if(kClusterArray == "") {
153 //    //printf("**************** Normal analysis **************** \n");
154 //    reader->SwitchOnClusterRecalculation(); // Bad map removal
155 //  }
156 //  else {
157 //    printf("**************** Input for analysis is Clusterizer %s **************** \n", kClusterArray.Data());
158 //    reader->SwitchOffClusterRecalculation();
159 //  }  
160   
161   if(kCalorimeter == "EMCAL") {
162     reader->SwitchOnEMCALCells();  
163     reader->SwitchOnEMCAL();
164   }
165   if(kCalorimeter == "PHOS") { 
166     reader->SwitchOnPHOSCells();  
167     reader->SwitchOnPHOS();
168   }
169   
170   // for case data="deltaAOD", no need to fill the EMCAL/PHOS cluster lists
171   if(kData.Contains("delta")){
172     reader->SwitchOffEMCAL();
173     reader->SwitchOffPHOS();
174     reader->SwitchOffEMCALCells(); 
175     reader->SwitchOffPHOSCells(); 
176   }
177   
178   //-----------------
179   // Event selection
180   //-----------------
181   
182   //if(!kUseKinematics) reader->SetFiredTriggerClassName("CEMC7EGA-B-NOPF-CENTNOTRD"); // L1 Gamma
183   
184   if     (kCollisions=="pp"  ) {
185     reader->SwitchOffEventSelection(); 
186     reader->SwitchOffV0ANDSelection() ;        // and besides v0 AND
187     reader->SwitchOffPrimaryVertexSelection(); // and besides primary vertex
188     reader->SetZvertexCut(50.);                // Open cut
189   }
190   else if(kCollisions=="PbPb") {
191     reader->SwitchOffEventSelection();         // remove pileup by default
192     reader->SwitchOffV0ANDSelection() ;        // and besides v0 AND
193     reader->SwitchOffPrimaryVertexSelection(); // and besides primary vertex
194     reader->SetZvertexCut(10.);                // Centrality defined in this range.
195     
196     // Centrality
197     reader->SetCentralityClass("V0M");
198     reader->SetCentralityOpt(10);  // 10 centrality bins
199     reader->SetCentralityBin(-1,-1); // Accept all events, if not select range
200     
201     // Event plane (only used in AliAnaPi0 for the moment)
202     reader->SetEventPlaneMethod("Q");
203   }
204   
205   if(kPrint) reader->Print("");
206   
207   return reader;
208   
209 }
210
211 //_______________________________________
212 AliCalorimeterUtils* ConfigureCaloUtils(Int_t kYears, TString kCollisions, TString kClusterArray, 
213                                         Bool_t kSimulation, Bool_t kPrint)
214 {
215   
216   AliCalorimeterUtils *cu = new AliCalorimeterUtils;
217   cu->SetDebug(-1);
218   
219   if(kYears==2010) cu->SetEMCALGeometryName("EMCAL_FIRSTYEARV1");
220   else             cu->SetEMCALGeometryName("EMCAL_COMPLETEV1");
221   
222   // Remove clusters close to borders, at least max energy cell is 1 cell away 
223   cu->SetNumberOfCellsFromEMCALBorder(1);
224   cu->SetNumberOfCellsFromPHOSBorder(2);
225   
226   //if(kClusterArray == "") 
227   //  cu->SwitchOffRecalculateClusterTrackMatching(); // Done in clusterization
228   //else            
229   //  cu->SwitchOnRecalculateClusterTrackMatching();
230   
231   //EMCAL only settings
232   AliEMCALRecoUtils * recou = cu->GetEMCALRecoUtils();
233   
234   
235   if(kCollisions=="pp"  ) { // Do only for pp for the moment
236     cu->SwitchOnCorrectClusterLinearity();
237     if(!kSimulation) recou->SetNonLinearityFunction(AliEMCALRecoUtils::kBeamTestCorrected);
238     else             recou->SetNonLinearityFunction(AliEMCALRecoUtils::kPi0MC);
239   }
240   
241   //recou->SwitchOnRejectExoticCell();
242   //if(kClusterArray == "") recou->SwitchOnRejectExoticCluster();
243   //else                    recou->SwitchOffRejectExoticCluster();
244   
245   if(kPrint) cu->Print("");
246   
247   return cu;
248   
249 }
250
251 //_____________________________________
252 AliAnaPhoton* ConfigurePhotonAnalysis(TString kData, TString kCalorimeter, Int_t kYears, TString kCollisions, TString kClusterArray, 
253                                       Bool_t kSimulation, TString kTrig = "", Bool_t kPrint)
254 {
255   
256   printf("Photon kData %s, kCalorimeter %s, kYears %d, kCollisions %s, kClusterArray %s, kSimulation %d, kTrig %s, kPrint %d \n",
257          kData.Data(), kCalorimeter.Data(), kYears, kCollisions.Data(), kClusterArray.Data(), kSimulation, kTrig.Data(), kPrint);
258   AliAnaPhoton *anaphoton = new AliAnaPhoton();
259   anaphoton->SetDebug(-1); //10 for lots of messages
260   
261   // cluster selection cuts
262   
263   anaphoton->SwitchOffFiducialCut();
264   
265   anaphoton->SetCalorimeter(kCalorimeter);
266   
267   if(kCalorimeter == "PHOS"){
268     anaphoton->SetNCellCut(2);// At least 2 cells
269     anaphoton->SetMinPt(0.3);
270     anaphoton->SetMinDistanceToBadChannel(2, 4, 5);
271   }
272   else {//EMCAL
273     anaphoton->SetNCellCut(1);// At least 2 cells
274     anaphoton->SetMinPt(0.5); // avoid mip peak at E = 260 MeV
275     anaphoton->SetMaxPt(1000); 
276     anaphoton->SetTimeCut(-1000,1000);// Time window of [400-900] ns if time recalibration is off, 
277     // restrict to less than 100 ns when time calibration is on 
278     anaphoton->SetMinDistanceToBadChannel(1, 2, 3); // For filtered AODs, new releases.
279   }
280   
281   anaphoton->SwitchOnTrackMatchRejection() ;
282   
283   //PID cuts (shower shape)
284   
285   AliCaloPID* caloPID = anaphoton->GetCaloPID();
286   anaphoton->SwitchOnCaloPID(); // if nothing else specified bayesian
287   anaphoton->SwitchOnCaloPIDRecalculation(); // off, get bayesian weights, on use simple cut
288   //caloPID->SetLambda0CutMax(0.30);
289   //caloPID->SetLambda0CutMin(0.10);
290   //anaphoton->SwitchOffFillShowerShapeHistograms();  // Filled before photon shower shape selection
291   
292   // Input / output delta AOD settings
293   
294   if(!kData.Contains("delta")) {
295     anaphoton->SetOutputAODName(Form("Photon%s_Trig%s_Cl%s",kCalorimeter.Data(), kTrig.Data(),kClusterArray.Data()));
296     anaphoton->SetOutputAODClassName("AliAODPWG4ParticleCorrelation");
297     //anaphoton->SetOutputAODClassName("AliAODPWG4Particle"); // use if no correlation done
298   }
299   else anaphoton->SetInputAODName(Form("Photon%s_Trig%s_Cl%s",kCalorimeter.Data(), kTrig.Data(),kClusterArray.Data()));
300   
301   //Set Histograms name tag, bins and ranges
302   
303   anaphoton->AddToHistogramsName("AnaPhoton_");
304   SetHistoRangeAndNBins(anaphoton, kCalorimeter, kYears, kCollisions, kSimulation); // see method below
305   
306   // Number of particle type MC histograms
307   //anaphoton->FillNOriginHistograms(8);
308   //anaphoton->FillNPrimaryHistograms(4);
309   
310   if(kPrint) anaphoton->Print("");
311   
312   return anaphoton;
313   
314 }
315
316
317
318 //_______________________________
319 AliAnaPi0* ConfigurePi0Analysis(TString kData, TString kCalorimeter, Int_t kYears, TString kCollisions, TString kClusterArray, 
320                                 Bool_t kSimulation, TString kTrig = "", Bool_t kPrint)
321 {
322   
323   AliAnaPi0 *anapi0 = new AliAnaPi0();
324   
325   anapi0->SetDebug(-1);//10 for lots of messages
326   if(kPrint) anapi0->Print("");
327   
328   // Input delta AOD settings
329   anapi0->SetInputAODName(Form("Photon%s_Trig%s_Cl%s",kCalorimeter.Data(), kTrig.Data(),kClusterArray.Data()));
330   
331   // Calorimeter settings
332   anapi0->SetCalorimeter(kCalorimeter);
333   if(kCalorimeter=="PHOS") anapi0->SetNumberOfModules(3); //PHOS first year
334   else {                   
335     if(kYears==2010) anapi0->SetNumberOfModules(4); //EMCAL first year
336     else             anapi0->SetNumberOfModules(10);
337   }
338   
339   //settings for pp collision mixing
340   anapi0->SwitchOnOwnMix(); //Off when mixing done with general mixing frame
341   
342   // Cuts 
343   if(kCalorimeter=="EMCAL") anapi0->SetPairTimeCut(70);
344   
345   if     (kCollisions=="pp"  ) {
346     printf("****** Configure Pi0 for pp analysis\n");
347     anapi0->SetNCentrBin(1);
348     anapi0->SetNZvertBin(10);
349     anapi0->SetNRPBin(1);
350     anapi0->SetNMaxEvMix(100);    
351   }
352   else if(kCollisions=="PbPb") {
353     printf("****** Configure Pi0 for PbPb analysis\n");
354     anapi0->SetNCentrBin(5);
355     anapi0->SetNZvertBin(3);
356     anapi0->SetNRPBin(1);
357     anapi0->SetNMaxEvMix(5);
358   }
359   
360   anapi0->SwitchOffMultipleCutAnalysis();
361   anapi0->SwitchOffSMCombinations();
362
363   //Set Histograms name tag, bins and ranges
364   
365   anapi0->AddToHistogramsName("AnaPi0_");
366   SetHistoRangeAndNBins(anapi0, kCalorimeter, kYears, kCollisions, kSimulation); // see method below
367   
368   return anapi0;
369   
370 }
371
372 //________________________________________________________
373 void SetHistoRangeAndNBins (AliAnaPartCorrBaseClass* ana, TString kCalorimeter, 
374                             Int_t kYears, TString kCollisions, Bool_t kSimulation)
375 {
376   // Set common bins for all analysis and MC histograms filling
377   
378   if(kSimulation) ana->SwitchOnDataMC() ;//Access MC stack and fill more histograms, AOD MC not implemented yet.
379   else            ana->SwitchOffDataMC() ;
380   
381   ana->SetHistoPtRangeAndNBins(0, 100, 250) ; // Energy and pt histograms
382   
383   if(kCalorimeter=="EMCAL"){
384     if(kYears==2010){
385       ana->SetHistoPhiRangeAndNBins(78*TMath::DegToRad(), 122*TMath::DegToRad(), 78) ;
386       //ana->SetHistoXRangeAndNBins(-230,90,120); // QA
387       //ana->SetHistoYRangeAndNBins(370,450,40);  // QA
388     }
389     else {           
390       ana->SetHistoPhiRangeAndNBins(78*TMath::DegToRad(), 182*TMath::DegToRad(), 108) ;
391       //ana->SetHistoXRangeAndNBins(-600,90,200); // QA
392       //ana->SetHistoYRangeAndNBins(100,450,100); // QA
393     }
394     
395     ana->SetHistoEtaRangeAndNBins(-0.72, 0.72, 144) ;
396   }
397   else{
398     ana->SetHistoPhiRangeAndNBins(260*TMath::DegToRad(), 320*TMath::DegToRad(), 60) ;
399     ana->SetHistoEtaRangeAndNBins(-0.13, 0.13, 130) ;
400     
401   }
402   
403   //ana->SetHistoShowerShapeRangeAndNBins(0, 3, 300);
404   
405   // Invariant mass analysis
406   ana->SetHistoMassRangeAndNBins(0., 1., 200) ;
407   ana->SetHistoAsymmetryRangeAndNBins(0., 1. , 100) ;
408   
409   // check if time calibration is on
410   //ana->SetHistoTimeRangeAndNBins(-1000.,1000,1000);
411   //ana->SetHistoDiffTimeRangeAndNBins(-200, 200, 800);
412   
413   // QA, electron, charged
414   //ana->SetHistoPOverERangeAndNBins(0,10.,100);
415   //ana->SetHistodEdxRangeAndNBins(0.,200.,200);
416   
417   // QA
418   //ana->SetHistoFinePtRangeAndNBins(0, 10, 200) ; // bining for fhAmpId
419   //ana->SetHistodRRangeAndNBins(0.,TMath::Pi(),150);
420   //ana->SetHistoRatioRangeAndNBins(0.,2.,100);
421   //ana->SetHistoVertexDistRangeAndNBins(0.,500.,500);
422   //ana->SetHistoNClusterCellRangeAndNBins(0,500,500);
423   //ana->SetHistoZRangeAndNBins(-400,400,200);
424   //ana->SetHistoRRangeAndNBins(400,450,25);
425   //ana->SetHistoV0SignalRangeAndNBins(0,5000,500);
426   //ana->SetHistoV0MultiplicityRangeAndNBins(0,5000,500);
427   //ana->SetHistoTrackMultiplicityRangeAndNBins(0,5000,500);
428   
429 }
430