AliPhysics  7f1bdba (7f1bdba)
AliAnalysisTaskCaloConv.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: Dmitri Peressounko (RRC KI) *
5  * *
6  * Permission to use, copy, modify and distribute this software and its *
7  * documentation strictly for non-commercial purposes is hereby granted *
8  * without fee, provided that the above copyright notice appears in all *
9  * copies and that both the copyright notice and this permission notice *
10  * appear in the supporting documentation. The authors make no claims *
11  * about the suitability of this software for any purpose. It is *
12  * provided "as is" without express or implied warranty. *
13  **************************************************************************/
14 
16 //---------------------------------------------
17 // Class used to do analysis on conversion + calorimeter pairs
18 // A lot of code cut-and-pasted from AliV0Reader and GammaConversion classes
19 //
20 //---------------------------------------------
22 
23 // root
24 #include "TChain.h"
25 #include "TH3.h"
26 #include "TH2.h"
27 #include "TDirectory.h"
28 #include "TLorentzVector.h"
29 #include "TRandom.h"
30 
31 // analysis
32 #include "AliAnalysisManager.h"
33 #include "AliESDInputHandler.h"
34 #include "AliAnalysisTaskCaloConv.h"
35 #include "AliLog.h"
36 #include "AliESDEvent.h"
37 #include "AliESDpid.h"
38 #include "AliESDtrackCuts.h"
39 #include "AliESDtrackCuts.h"
40 #include "AliCFContainer.h" // for CF
41 #include "AliESDCaloCluster.h"
42 #include "AliPHOSGeoUtils.h"
43 #include "AliEMCALGeometry.h"
44 #include "AliCFContainer.h"
45 #include "AliMCEventHandler.h"
46 #include "AliMCEvent.h"
47 #include "AliESDv0.h"
48 #include "AliKFParticle.h"
49 #include "AliKFVertex.h"
50 
51 class Riostream;
52 class TFile;
53 
54 ClassImp(AliAnalysisTaskCaloConv)
55 
56 
57 AliAnalysisTaskCaloConv::AliAnalysisTaskCaloConv():
58 AliAnalysisTaskSE(),
59  fMCEvent(NULL),
60  fESDEvent(NULL),
61  fESDpid(NULL),
62  fESDtrackCuts(NULL),
63  fOutputContainer(NULL),
64  fCFOutputContainer(NULL),
65  fConvCFCont(0x0),
66  fPHOSCFCont(0x0),
67  fEMCALCFCont(0x0),
68  fPi0CFCont(0x0),
69  fCentr(0.),
70  fTriggerCINT1B(kFALSE),
71  fToUseCF(kFALSE),
72  fMinOpeningAngleGhostCut(0.),
73  fPHOSgeom(0x0),
74  fEMCALgeom(0x0),
75  fPi0Thresh1(0.5),
76  fPi0Thresh2(1.),
77  fBadDistCutPHOS(3.3),
78  fBadDistCutEMCAL(6.),
79  fGammaV0s(),
80  fGammaPHOS(),
81  fGammaEMCAL(),
82  fConvEvent(NULL) ,
83  fPHOSEvent(NULL),
84  fEMCALEvent(NULL),
85  fnSigmaAboveElectronLine(5.),
86  fnSigmaBelowElectronLine(-3.),
87  fnSigmaAbovePionLine(0.),
88  fpnSigmaAbovePionLine(1.),
89  fprobCut(0.),
90  fmaxR(180.),
91  fmaxZ(240.),
92  fetaCut(0.9),
93  fptCut(0.02),
94  fchi2CutConversion(30.)
95 {
96  // Default constructor
97  Int_t nBin=10 ;
98  for(Int_t i=0;i<nBin;i++){
99  fPHOSEvents[i]=0 ;
100  fEMCALEvents[i]=0;
101  fConvEvents[i]=0;
102  }
103  char key[55] ;
104  for(Int_t i=0; i<6; i++){
105  snprintf(key,55,"PHOS_BadMap_mod%d",i) ;
106  fPHOSBadMap[i]=new TH2I(key,"Bad Modules map",64,0.,64.,56,0.,56.) ;
107  }
108  for(Int_t i=0; i<10; i++){
109  snprintf(key,55,"EMCAL_BadMap_mod%d",i) ;
110  fEMCALBadMap[i] = new TH2I(key,"Bad Modules map",24,0.,24.,48,0.,48.) ;
111  }
112 }
113 AliAnalysisTaskCaloConv::AliAnalysisTaskCaloConv(const char* name):
114  AliAnalysisTaskSE(name),
115  fMCEvent(NULL),
116  fESDEvent(NULL),
117  fESDpid(NULL),
118  fESDtrackCuts(NULL),
119  fOutputContainer(NULL),
120  fCFOutputContainer(NULL),
121  fConvCFCont(0x0),
122  fPHOSCFCont(0x0),
123  fEMCALCFCont(0x0),
124  fPi0CFCont(0x0),
125  fCentr(0.),
126  fTriggerCINT1B(kFALSE),
127  fToUseCF(kFALSE),
128  fMinOpeningAngleGhostCut(0.),
129  fPHOSgeom(0x0),
130  fEMCALgeom(0x0),
131  fPi0Thresh1(0.5),
132  fPi0Thresh2(1.),
133  fBadDistCutPHOS(3.3),
134  fBadDistCutEMCAL(6.),
135  fGammaV0s(),
136  fGammaPHOS(),
137  fGammaEMCAL(),
138  fConvEvent(NULL) ,
139  fPHOSEvent(NULL),
140  fEMCALEvent(NULL),
141  fnSigmaAboveElectronLine(5.),
142  fnSigmaBelowElectronLine(-3.),
143  fnSigmaAbovePionLine(0.),
144  fpnSigmaAbovePionLine(1.),
145  fprobCut(0.),
146  fmaxR(180.),
147  fmaxZ(240.),
148  fetaCut(0.9),
149  fptCut(0.02),
150  fchi2CutConversion(30.)
151 {
152  // Common I/O in slot 0
153  DefineInput (0, TChain::Class());
154  DefineOutput(0, TTree::Class());
155 
156  // Your private output
157  DefineOutput(1, TList::Class());
158  DefineOutput(2, TList::Class()); // for CF
159 
160  Int_t nBin=10 ;
161  for(Int_t i=0;i<nBin;i++){
162  fPHOSEvents[i]=0 ;
163  fEMCALEvents[i]=0;
164  fConvEvents[i]=0;
165  }
166  char key[55] ;
167  for(Int_t i=0; i<6; i++){
168  snprintf(key,55,"PHOS_BadMap_mod%d",i) ;
169  fPHOSBadMap[i]=new TH2I(key,"Bad Modules map",64,0.,64.,56,0.,56.) ;
170  }
171  for(Int_t i=0; i<10; i++){
172  snprintf(key,55,"EMCAL_BadMap_mod%d",i) ;
173  fEMCALBadMap[i] = new TH2I(key,"Bad Modules map",24,0.,24.,48,0.,48.) ;
174  }
175 // fESDpid = new AliESDpid;
176 }
177 //_____________________________________________________
178 AliAnalysisTaskCaloConv::~AliAnalysisTaskCaloConv()
179 {
180  // Remove all pointers
181 
182  if (AliAnalysisManager::GetAnalysisManager()->GetAnalysisType() !=
183  AliAnalysisManager::kProofAnalysis) {
184 
185  if(fOutputContainer){
186  fOutputContainer->Clear() ;
187  delete fOutputContainer ;
188  }
189  if(fCFOutputContainer){
190  fCFOutputContainer->Clear() ;
191  delete fCFOutputContainer ;
192  }
193 
194  if(fPHOSgeom){
195  delete fPHOSgeom ;
196  fPHOSgeom=0x0 ;
197  }
198 
199  if(fEMCALgeom){
200  delete fEMCALgeom ;
201  fEMCALgeom=0x0;
202  }
203 
204  for(Int_t ivtx=0; ivtx<10; ivtx++){
205  if(fPHOSEvents[ivtx]){
206  delete fPHOSEvents[ivtx] ;
207  fPHOSEvents[ivtx]=0x0 ;
208  }
209  if(fEMCALEvents[ivtx]){
210  delete fEMCALEvents[ivtx] ;
211  fEMCALEvents[ivtx]=0x0 ;
212  }
213  if(fConvEvents[ivtx]){
214  delete fConvEvents[ivtx] ;
215  fConvEvents[ivtx]=0x0 ;
216  }
217  }
218  for(Int_t i=0; i<6; i++)
219  if(fPHOSBadMap[i]){
220  delete fPHOSBadMap[i] ;
221  fPHOSBadMap[i]=0 ;
222  }
223  for(Int_t i=0; i<10; i++)
224  if(fEMCALBadMap[i]){
225  delete fEMCALBadMap[i];
226  fEMCALBadMap[i]=0 ;
227  }
228  }
229 }
230 //_____________________________________________________
231 void AliAnalysisTaskCaloConv::Init()
232 {
233  // Initialization
234  // AliLog::SetGlobalLogLevel(AliLog::kError);
235 }
236 //_____________________________________________________
237 void AliAnalysisTaskCaloConv::UserExec(Option_t */*option*/)
238 {
239  // Execute analysis for current event
240  // First select conversion and calorimeter photons
241  // then construct inv. mass distributions
242 
243  fMCEvent = MCEvent();
244 
245  AliESDInputHandler *esdHandler=dynamic_cast<AliESDInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
246  if(!fESDpid){
247  if( esdHandler && esdHandler->GetESDpid()){
248  fESDpid=new AliESDpid(*(esdHandler->GetESDpid())) ;
249  }
250  else {
251  fESDpid=new AliESDpid;
252  Double_t alephParameters[5];
253  if(static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())->MCEvent()){// simulation
254  alephParameters[0] = 2.15898e+00/50.;
255  alephParameters[1] = 1.75295e+01;
256  alephParameters[2] = 3.40030e-09;
257  alephParameters[3] = 1.96178e+00;
258  alephParameters[4] = 3.91720e+00;
259  fESDpid->GetTOFResponse().SetTimeResolution(80.);
260  }
261  else{// data
262  alephParameters[0] = 0.0283086;
263  alephParameters[1] = 2.63394e+01;
264  alephParameters[2] = 5.04114e-11;
265  alephParameters[3] = 2.12543e+00;
266  alephParameters[4] = 4.88663e+00;
267  fESDpid->GetTOFResponse().SetTimeResolution(130.);
268  fESDpid->GetTPCResponse().SetMip(47.9);
269  }
270 
271  fESDpid->GetTPCResponse().SetBetheBlochParameters(
272  alephParameters[0],alephParameters[1],alephParameters[2],
273  alephParameters[3],alephParameters[4]);
274  fESDpid->GetTPCResponse().SetSigma(3.79301e-03, 2.21280e+04);
275  }
276  }
277 
278  if(!fESDtrackCuts){
279 // fESDtrackCuts= AliESDtrackCuts::GetStandardITSTPCTrackCuts2010();
280  fESDtrackCuts = new AliESDtrackCuts;
281 
282  // TPC
283  fESDtrackCuts->SetMinNClustersTPC(70);
284  fESDtrackCuts->SetMaxChi2PerClusterTPC(4);
285  fESDtrackCuts->SetAcceptKinkDaughters(kFALSE);
286  fESDtrackCuts->SetRequireTPCRefit(kTRUE);
287  // ITS
288  fESDtrackCuts->SetRequireITSRefit(kTRUE);
289  fESDtrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
290  AliESDtrackCuts::kAny);
291 // if(selPrimaries) {
292  // 7*(0.0026+0.0050/pt^1.01)
293  fESDtrackCuts->SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");
294 // }
295  fESDtrackCuts->SetMaxDCAToVertexZ(2);
296  fESDtrackCuts->SetDCAToVertex2D(kFALSE);
297  fESDtrackCuts->SetRequireSigmaToVertex(kFALSE);
298 
299 
300 
301 /*
302  fESDtrackCuts = new AliESDtrackCuts("AliESDtrackCuts","AliESDtrackCuts");
303 
304  fESDtrackCuts->SetAcceptKinkDaughters(kFALSE);
305  fESDtrackCuts->SetMinNClustersTPC(70);
306  fESDtrackCuts->SetMaxChi2PerClusterTPC(4);
307  fESDtrackCuts->SetRequireTPCRefit(kTRUE);
308  fESDtrackCuts->SetRequireITSRefit(kTRUE);
309 // fESDtrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny); //TEMPORARY <-> REMOVE
310 */
311  }
312 
313 
314  fESDEvent=(AliESDEvent*)InputEvent();
315  FillHistogram("hEventsTrig",0.5) ;
316  Bool_t isSelected = esdHandler && ((esdHandler->IsEventSelected()& AliVEvent::kMB) == AliVEvent::kMB);
317  if(!isSelected){
318  printf("Not selected !!!!! \n") ;
319  PostData(1, fOutputContainer);
320  return ;
321  }
322  FillHistogram("hEventsTrig",1.5) ;
323 
324  //Take Only events with proper trigger
325  //No trigger in MC data => no check
326 // if(!fMCEvent && !fESDEvent->IsTriggerClassFired("CINT1B-ABCE-NOPF-ALL")){
327  //for LHC10e
328  if(!fMCEvent && !fESDEvent->IsTriggerClassFired("CINT1-B-NOPF-ALLNOTRD")){
329  PostData(1, fOutputContainer);
330  return ;
331  }
332  FillHistogram("hEventsTrig",2.5) ;
333 
334  //checks if we have a prim vertex
335  if(fESDEvent->GetPrimaryVertex()->GetNContributors()<=0) {
336  PostData(1, fOutputContainer);
337  return ;
338  }
339  FillHistogram("hEventsTrig",3.5) ;
340 
341  if(TMath::Abs(fESDEvent->GetPrimaryVertex()->GetZ())>10.){
342  PostData(1, fOutputContainer);
343  return ;
344  }
345  FillHistogram("hEventsTrig",4.5) ;
346 
347 
348  //Calculate charged multiplicity
349  Int_t trackCounter = 0;
350  for (Int_t i=0;i<fESDEvent->GetNumberOfTracks();++i) {
351  AliESDtrack *track = new AliESDtrack(*fESDEvent->GetTrack(i)) ;
352  if(fESDtrackCuts->AcceptTrack(track) && TMath::Abs(track->Eta())< 0.9)
353  trackCounter++;
354  delete track;
355  }
356  fCentr=trackCounter+0.5 ;
357  FillHistogram("hMult",fCentr) ;
358 
359  //Init geometry if not done yet
360  InitGeometry();
361 
362  //Select conversion and calorimeter photons
363  //Conversion photons should go first since they are used in calibration in SelectCALOPhotons()
364  SelectConvPhotons() ;
365  SelectPHOSPhotons() ;
366  SelectEMCALPhotons() ;
367  //Fill MC histograms if MC is present
368  ProcessMC();
369  FillRealMixed() ;
370 
371  PostData(1, fOutputContainer);
372  if(fToUseCF)
373  PostData(2, fCFOutputContainer); // for CF
374 
375 
376 }
377 //____________________________________________________________
378 void AliAnalysisTaskCaloConv::ConnectInputData(Option_t *option){
379  // see header file for documentation
380 
381  AliAnalysisTaskSE::ConnectInputData(option);
382 
383 }
384 //____________________________________________________________
385 void AliAnalysisTaskCaloConv::UserCreateOutputObjects()
386 {
387  //UserCreateOutputObjects
388  if(fDebug)gDirectory->Print() ;
389  // Create the output container
390  if(fOutputContainer != NULL){
391  delete fOutputContainer;
392  fOutputContainer = NULL;
393  }
394  fOutputContainer = new TList();
395  fOutputContainer->SetOwner(kTRUE);
396 
397  if(fCFOutputContainer != NULL){
398  delete fCFOutputContainer;
399  fCFOutputContainer = NULL;
400  }
401  //===========Correction Framework ======================
402  if(fToUseCF){
403  fCFOutputContainer = new TList();
404  fCFOutputContainer->SetOwner(kTRUE);
405 
406  //bins: pt,eta,mass
407  Int_t iBin[3]={500,40,100};
408  fConvCFCont = new AliCFContainer("ConvContainer","container for converted photons", 23,3,iBin);
409  fConvCFCont->SetBinLimits(0,0.,50.);
410  fConvCFCont->SetBinLimits(1,-2.,2.) ;
411  fConvCFCont->SetBinLimits(2,0.,1.);
412  fCFOutputContainer->Add(fConvCFCont) ;
413 
414  fPHOSCFCont = new AliCFContainer("PHOSContainer","container for PHOS photons", 10,2,iBin);
415  fPHOSCFCont->SetBinLimits(0,0.,50.);
416  fPHOSCFCont->SetBinLimits(1,-2.,2.) ;
417  fCFOutputContainer->Add(fPHOSCFCont) ;
418 
419  fEMCALCFCont = new AliCFContainer("EMCALContainer","container for EMCAL photons", 10,2,iBin);
420  fEMCALCFCont->SetBinLimits(0,0.,50.);
421  fEMCALCFCont->SetBinLimits(1,-2.,2.) ;
422  fCFOutputContainer->Add(fEMCALCFCont) ;
423 
424  fPi0CFCont = new AliCFContainer("Pi0Container","container for EMCAL photons", 10,2,iBin);
425  fPi0CFCont->SetBinLimits(0,0.,50.);
426  fPi0CFCont->SetBinLimits(1,-2.,2.) ;
427  fCFOutputContainer->Add(fPi0CFCont) ;
428 
429  }
430  //========================================================
431 
432  //Adding the histograms to the output container
433  Int_t firstRun= 125000 ;
434  Int_t lastRun = 135000 ;
435  Int_t nRuns =lastRun-firstRun+1 ;
436 
437  //Run QA histigrams
438  fOutputContainer->Add(new TH2F("hRunTrigger","Triggers fired",nRuns,float(firstRun),float(lastRun),2,0.,2.)) ;
439  fOutputContainer->Add(new TH1F("hRunEvents","Events per run",nRuns,float(firstRun),float(lastRun))) ;
440  fOutputContainer->Add(new TH1F("hRunConvs","Conversion photons per run",nRuns,float(firstRun),float(lastRun))) ;
441  fOutputContainer->Add(new TH1F("hRunPHOS","PHOS photons per run",nRuns,float(firstRun),float(lastRun))) ;
442  fOutputContainer->Add(new TH1F("hRunEMCAL","EMCAL photons per run",nRuns,float(firstRun),float(lastRun))) ;
443  fOutputContainer->Add(new TH1F("hVtxBin","Vtx distribution",10,0.,10.)) ;
444  fOutputContainer->Add(new TH1F("hEvents","Events processed",1,0.,1.)) ;
445  fOutputContainer->Add(new TH1F("hEventsTrig","Events processed",10,0.,10.)) ;
446 
447  fOutputContainer->Add(new TH2F("hQA_PHOS_mod1_soft","number of clusters per cell",64,0.,64.,56,0.,56.)) ;
448  fOutputContainer->Add(new TH2F("hQA_PHOS_mod2_soft","number of clusters per cell",64,0.,64.,56,0.,56.)) ;
449  fOutputContainer->Add(new TH2F("hQA_PHOS_mod3_soft","number of clusters per cell",64,0.,64.,56,0.,56.)) ;
450  fOutputContainer->Add(new TH2F("hQA_PHOS_mod4_soft","number of clusters per cell",64,0.,64.,56,0.,56.)) ;
451  fOutputContainer->Add(new TH2F("hQA_PHOS_mod5_soft","number of clusters per cell",64,0.,64.,56,0.,56.)) ;
452  fOutputContainer->Add(new TH2F("hQA_PHOS_mod1_hard","number of clusters per cell",64,0.,64.,56,0.,56.)) ;
453  fOutputContainer->Add(new TH2F("hQA_PHOS_mod2_hard","number of clusters per cell",64,0.,64.,56,0.,56.)) ;
454  fOutputContainer->Add(new TH2F("hQA_PHOS_mod3_hard","number of clusters per cell",64,0.,64.,56,0.,56.)) ;
455  fOutputContainer->Add(new TH2F("hQA_PHOS_mod4_hard","number of clusters per cell",64,0.,64.,56,0.,56.)) ;
456  fOutputContainer->Add(new TH2F("hQA_PHOS_mod5_hard","number of clusters per cell",64,0.,64.,56,0.,56.)) ;
457 
458  fOutputContainer->Add(new TH2F("hQA_EMCAL_SM0_soft","number of clusters per cell",24,0.,24,48,0.,48.)) ;
459  fOutputContainer->Add(new TH2F("hQA_EMCAL_SM1_soft","number of clusters per cell",24,0.,24,48,0.,48.)) ;
460  fOutputContainer->Add(new TH2F("hQA_EMCAL_SM2_soft","number of clusters per cell",24,0.,24,48,0.,48.)) ;
461  fOutputContainer->Add(new TH2F("hQA_EMCAL_SM3_soft","number of clusters per cell",24,0.,24,48,0.,48.)) ;
462  fOutputContainer->Add(new TH2F("hQA_EMCAL_SM4_soft","number of clusters per cell",24,0.,24,48,0.,48.)) ;
463  fOutputContainer->Add(new TH2F("hQA_EMCAL_SM5_soft","number of clusters per cell",24,0.,24,48,0.,48.)) ;
464  fOutputContainer->Add(new TH2F("hQA_EMCAL_SM6_soft","number of clusters per cell",24,0.,24,48,0.,48.)) ;
465  fOutputContainer->Add(new TH2F("hQA_EMCAL_SM7_soft","number of clusters per cell",24,0.,24,48,0.,48.)) ;
466  fOutputContainer->Add(new TH2F("hQA_EMCAL_SM8_soft","number of clusters per cell",24,0.,24,48,0.,48.)) ;
467  fOutputContainer->Add(new TH2F("hQA_EMCAL_SM9_soft","number of clusters per cell",24,0.,24,48,0.,48.)) ;
468  fOutputContainer->Add(new TH2F("hQA_EMCAL_SM0_hard","number of clusters per cell",24,0.,24,48,0.,48.)) ;
469  fOutputContainer->Add(new TH2F("hQA_EMCAL_SM1_hard","number of clusters per cell",24,0.,24,48,0.,48.)) ;
470  fOutputContainer->Add(new TH2F("hQA_EMCAL_SM2_hard","number of clusters per cell",24,0.,24,48,0.,48.)) ;
471  fOutputContainer->Add(new TH2F("hQA_EMCAL_SM3_hard","number of clusters per cell",24,0.,24,48,0.,48.)) ;
472  fOutputContainer->Add(new TH2F("hQA_EMCAL_SM4_hard","number of clusters per cell",24,0.,24,48,0.,48.)) ;
473  fOutputContainer->Add(new TH2F("hQA_EMCAL_SM5_hard","number of clusters per cell",24,0.,24,48,0.,48.)) ;
474  fOutputContainer->Add(new TH2F("hQA_EMCAL_SM6_hard","number of clusters per cell",24,0.,24,48,0.,48.)) ;
475  fOutputContainer->Add(new TH2F("hQA_EMCAL_SM7_hard","number of clusters per cell",24,0.,24,48,0.,48.)) ;
476  fOutputContainer->Add(new TH2F("hQA_EMCAL_SM8_hard","number of clusters per cell",24,0.,24,48,0.,48.)) ;
477  fOutputContainer->Add(new TH2F("hQA_EMCAL_SM9_hard","number of clusters per cell",24,0.,24,48,0.,48.)) ;
478 
479  fOutputContainer->Add(new TH2F("hQA_ConvPhiEta","Number of V0s phi eta",100,0.,TMath::TwoPi(),40,-1.5,1.5)) ;
480 
481  fOutputContainer->Add(new TH2F("hdEdx","dEdx of acceptaed electrons",1000,0.,10.,150,0.,150.)) ;
482 
483  fOutputContainer->Add(new TH1F("hMult","Multiplicity",200,0.,200.)) ;
484 
485  Int_t npt=200 ;
486  Double_t ptmax=20. ;
487  //Calibration of PHOS
488  fOutputContainer->Add(new TH3F("PHOS_mod1_th1","Inv.Mass distr. per channel",64,0.,64,56,0.,56,100,0.,0.5)) ;
489  fOutputContainer->Add(new TH3F("PHOS_mod2_th1","Inv.Mass distr. per channel",64,0.,64,56,0.,56,100,0.,0.5)) ;
490  fOutputContainer->Add(new TH3F("PHOS_mod3_th1","Inv.Mass distr. per channel",64,0.,64,56,0.,56,100,0.,0.5)) ;
491  fOutputContainer->Add(new TH3F("PHOS_mod4_th1","Inv.Mass distr. per channel",64,0.,64,56,0.,56,100,0.,0.5)) ;
492  fOutputContainer->Add(new TH3F("PHOS_mod5_th1","Inv.Mass distr. per channel",64,0.,64,56,0.,56,100,0.,0.5)) ;
493 
494  fOutputContainer->Add(new TH3F("PHOS_mod1_th2","Inv.Mass distr. per channel",64,0.,64,56,0.,56,100,0.,0.5)) ;
495  fOutputContainer->Add(new TH3F("PHOS_mod2_th2","Inv.Mass distr. per channel",64,0.,64,56,0.,56,100,0.,0.5)) ;
496  fOutputContainer->Add(new TH3F("PHOS_mod3_th2","Inv.Mass distr. per channel",64,0.,64,56,0.,56,100,0.,0.5)) ;
497  fOutputContainer->Add(new TH3F("PHOS_mod4_th2","Inv.Mass distr. per channel",64,0.,64,56,0.,56,100,0.,0.5)) ;
498  fOutputContainer->Add(new TH3F("PHOS_mod5_th2","Inv.Mass distr. per channel",64,0.,64,56,0.,56,100,0.,0.5)) ;
499 
500  //Pi0 histograms
501  //Vary Conversion cuts
502  fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_OnFly","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
503  fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_Offline","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
504  fOutputContainer->Add(new TH3F("PHOS_Re_mvsPt_OnFly_mult","Mass vs pt",400,0.,1.,npt,0.,ptmax,150,0.,150.)) ;
505  fOutputContainer->Add(new TH3F("PHOS_Re_mvsPt_Offline_mult","Mass vs pt",400,0.,1.,npt,0.,ptmax,150,0.,150.)) ;
506  fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_On_dEdx","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
507  fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_Off_dEdx","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
508 // fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_On_Prob","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
509 // fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_Off_Prob","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
510  fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_On_R120","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
511  fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_Off_R120","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
512  fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_On_Z","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
513  fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_Off_Z","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
514  fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_On_chi","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
515  fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_Off_chi","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
516  fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_On_Eta","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
517  fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_Off_Eta","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
518  fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_On_Wcut","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
519  fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_On_Wcut_Neu","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
520  fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_Off_Wcut","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
521  fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_Off_Wcut_Neu","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
522  fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_On_ArmQt","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
523  fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_Off_ArmQt","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
524 
525  fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_OnFly","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
526  fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_Offline","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
527  fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_On_dEdx","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
528  fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_Off_dEdx","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
529 // fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_On_Prob","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
530 // fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_Off_Prob","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
531  fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_On_R120","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
532  fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_Off_R120","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
533  fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_On_Z","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
534  fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_Off_Z","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
535  fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_On_chi","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
536  fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_Off_chi","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
537  fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_On_Eta","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
538  fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_Off_Eta","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
539  fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_On_Wcut","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
540  fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_Off_Wcut","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
541  fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_On_ArmQt","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
542  fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_Off_ArmQt","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
543 
544  fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_OnFly","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
545  fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_Offline","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
546  fOutputContainer->Add(new TH3F("EMCAL_Re_mvsPt_OnFly_mult","Mass vs pt",400,0.,1.,npt,0.,ptmax,30,0.,60.)) ;
547  fOutputContainer->Add(new TH3F("EMCAL_Re_mvsPt_Offline_mult","Mass vs pt",400,0.,1.,npt,0.,ptmax,30,0.,60.)) ;
548  fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_On_ArmQt","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
549  fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_Off_ArmQt","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
550  fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_On_dEdx","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
551  fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_Off_dEdx","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
552 // fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_On_Prob","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
553 // fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_Off_Prob","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
554  fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_On_R120","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
555  fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_Off_R120","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
556  fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_On_Z","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
557  fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_Off_Z","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
558  fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_On_chi","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
559  fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_Off_chi","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
560  fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_On_Eta","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
561  fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_Off_Eta","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
562  fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_On_Wcut","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
563  fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_Off_Wcut","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
564 
565  fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_OnFly","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
566  fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_Offline","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
567  fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_On_ArmQt","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
568  fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_Off_ArmQt","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
569  fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_On_dEdx","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
570  fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_Off_dEdx","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
571 // fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_On_Prob","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
572 // fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_Off_Prob","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
573  fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_On_R120","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
574  fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_Off_R120","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
575  fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_On_Z","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
576  fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_Off_Z","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
577  fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_On_chi","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
578  fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_Off_chi","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
579  fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_On_Eta","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
580  fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_Off_Eta","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
581  fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_On_Wcut","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
582  fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_Off_Wcut","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
583 
584  //PHOS PID variations
585  fOutputContainer->Add(new TH3F("PHOS_Re_mvsPt_alpha","Mass vs pt vs PHOS E",400,0.,1.,npt,0.,ptmax,10,0.,1.)) ;
586  fOutputContainer->Add(new TH3F("PHOS_Re_mvsPt_E","Mass vs pt vs PHOS E",400,0.,1.,npt,0.,ptmax,100,0.,10.)) ;
587  fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_all","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
588  fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_all_dist","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
589  fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_Disp","Mass vs pt, disp cut",400,0.,1.,npt,0.,ptmax)) ;
590  fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_TOF","Mass vs pt, TOF cut",400,0.,1.,npt,0.,ptmax)) ;
591  fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_Neutral","Mass vs pt, Neutral cut",400,0.,1.,npt,0.,ptmax)) ;
592  fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_DispNeutral","Mass vs pt, Disp & neutral cut",400,0.,1.,npt,0.,ptmax)) ;
593 
594  fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_all","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
595  fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_all_dist","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
596  fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_Disp","Mass vs pt, disp cut",400,0.,1.,npt,0.,ptmax)) ;
597  fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_TOF","Mass vs pt, TOF cut",400,0.,1.,npt,0.,ptmax)) ;
598  fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_Neutral","Mass vs pt, Neutral cut",400,0.,1.,npt,0.,ptmax)) ;
599  fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_DispNeutral","Mass vs pt, Disp & neutral cut",400,0.,1.,npt,0.,ptmax)) ;
600 
601  char key[155] ;
602  for(Int_t mod=1; mod<=5;mod++){
603  snprintf(key,155,"PHOS_Re_mvsPt_mod%d_single",mod) ;
604  fOutputContainer->Add(new TH2F(key,"Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
605  snprintf(key,155,"PHOS_Re_mvsPt_mod%d_all",mod) ;
606  fOutputContainer->Add(new TH2F(key,"Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
607  }
608 
609  //Single photon spectrum
610  //Conversion
611  fOutputContainer->Add(new TH1F("Single_conv_OnFly","Single photon spectrum",npt,0.,ptmax)) ;
612  fOutputContainer->Add(new TH1F("Single_conv_Offline","Single photon spectrum",npt,0.,ptmax)) ;
613  fOutputContainer->Add(new TH1F("Single_conv_On_ArmQt","Single photon spectrum",npt,0.,ptmax)) ;
614  fOutputContainer->Add(new TH1F("Single_conv_Off_ArmQt","Single photon spectrum",npt,0.,ptmax)) ;
615  fOutputContainer->Add(new TH1F("Single_conv_On_dEdx","Single photon spectrum",npt,0.,ptmax)) ;
616  fOutputContainer->Add(new TH1F("Single_conv_Off_dEdx","Single photon spectrum",npt,0.,ptmax)) ;
617 // fOutputContainer->Add(new TH1F("Single_conv_On_Prob","Single photon spectrum",npt,0.,ptmax)) ;
618 // fOutputContainer->Add(new TH1F("Single_conv_Off_Prob","Single photon spectrum",npt,0.,ptmax)) ;
619  fOutputContainer->Add(new TH1F("Single_conv_On_R120","Single photon spectrum",npt,0.,ptmax)) ;
620  fOutputContainer->Add(new TH1F("Single_conv_Off_R120","Single photon spectrum",npt,0.,ptmax)) ;
621  fOutputContainer->Add(new TH1F("Single_conv_On_Z","Single photon spectrum",npt,0.,ptmax)) ;
622  fOutputContainer->Add(new TH1F("Single_conv_Off_Z","Single photon spectrum",npt,0.,ptmax)) ;
623  fOutputContainer->Add(new TH1F("Single_conv_On_chi","Single photon spectrum",npt,0.,ptmax)) ;
624  fOutputContainer->Add(new TH1F("Single_conv_Off_chi","Single photon spectrum",npt,0.,ptmax)) ;
625  fOutputContainer->Add(new TH1F("Single_conv_On_Eta","Single photon spectrum",npt,0.,ptmax)) ;
626  fOutputContainer->Add(new TH1F("Single_conv_Off_Eta","Single photon spectrum",npt,0.,ptmax)) ;
627  fOutputContainer->Add(new TH1F("Single_conv_On_Wcut","Single photon spectrum",npt,0.,ptmax)) ;
628  fOutputContainer->Add(new TH1F("Single_conv_Off_Wcut","Single photon spectrum",npt,0.,ptmax)) ;
629 
630  //PHOS
631  fOutputContainer->Add(new TH2F("PHOS_single_all_mult","Single photon spectrum",npt,0.,ptmax,150,0.,150.)) ;
632  fOutputContainer->Add(new TH2F("PHOS_single_disp_mult","Single photon spectrum",npt,0.,ptmax,150,0.,150.)) ;
633  fOutputContainer->Add(new TH2F("PHOS_single_neu_mult","Single photon spectrum",npt,0.,ptmax,150,0.,150.)) ;
634 
635  for(Int_t mod=1; mod<=5;mod++){
636  snprintf(key,155,"PHOS_single_mod%d_all",mod) ;
637  fOutputContainer->Add(new TH1F(key,"Single photon spectrum",npt,0.,ptmax)) ;
638  snprintf(key,155,"PHOS_single_mod%d_disp",mod) ;
639  fOutputContainer->Add(new TH1F(key,"Single photon spectrum",npt,0.,ptmax)) ;
640  snprintf(key,155,"PHOS_single_mod%d_neutral",mod) ;
641  fOutputContainer->Add(new TH1F(key,"Single photon spectrum",npt,0.,ptmax)) ;
642  snprintf(key,155,"PHOS_single_mod%d_dispneutral",mod) ;
643  fOutputContainer->Add(new TH1F(key,"Single photon spectrum",npt,0.,ptmax)) ;
644  snprintf(key,155,"PHOS_single_mod%d_dist1",mod) ;
645  fOutputContainer->Add(new TH1F(key,"Single photon spectrum",npt,0.,ptmax)) ;
646  snprintf(key,155,"PHOS_single_mod%d_dist2",mod) ;
647  fOutputContainer->Add(new TH1F(key,"Single photon spectrum",npt,0.,ptmax)) ;
648  }
649 
650  for(Int_t mod=0; mod<4;mod++){
651  snprintf(key,155,"EMCAL_mod%d_th1",mod) ;
652  fOutputContainer->Add(new TH3F(key,"Inv.Mass distr. per channel",24,0.,24,48,0.,48,100,0.,0.5)) ;
653  snprintf(key,155,"EMCAL_mod%d_th2",mod) ;
654  fOutputContainer->Add(new TH3F(key,"Inv.Mass distr. per channel",24,0.,24,48,0.,48,100,0.,0.5)) ;
655 
656  snprintf(key,155,"EMCAL_Re_mvsPt_mod%d_single",mod) ;
657  fOutputContainer->Add(new TH2F(key,"Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
658 
659  snprintf(key,155,"EMCAL_Re_mvsPt_mod%d_all",mod) ;
660  fOutputContainer->Add(new TH2F(key,"Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
661  snprintf(key,155,"EMCAL_Re_mvsPt_mod%d_Disp",mod) ;
662  fOutputContainer->Add(new TH2F(key,"Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
663  snprintf(key,155,"EMCAL_Re_mvsPt_mod%d_TOF",mod) ;
664  fOutputContainer->Add(new TH2F(key,"Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
665  snprintf(key,155,"EMCAL_Re_mvsPt_mod%d_Neutral",mod) ;
666  fOutputContainer->Add(new TH2F(key,"Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
667  snprintf(key,155,"EMCAL_Re_mvsPt_mod%d_DispNeutral",mod) ;
668  fOutputContainer->Add(new TH2F(key,"Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
669 
670  snprintf(key,155,"EMCAL_Mi_mvsPt_mod%d_all",mod) ;
671  fOutputContainer->Add(new TH2F(key,"Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
672  snprintf(key,155,"EMCAL_Mi_mvsPt_mod%d_Disp",mod) ;
673  fOutputContainer->Add(new TH2F(key,"Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
674  snprintf(key,155,"EMCAL_Mi_mvsPt_mod%d_TOF",mod) ;
675  fOutputContainer->Add(new TH2F(key,"Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
676  snprintf(key,155,"EMCAL_Mi_mvsPt_mod%d_Neutral",mod) ;
677  fOutputContainer->Add(new TH2F(key,"Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
678  snprintf(key,155,"EMCAL_Mi_mvsPt_mod%d_DispNeutral",mod) ;
679  fOutputContainer->Add(new TH2F(key,"Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
680  }
681 
682  //MC info
683  fOutputContainer->Add(new TH2F("hMC_CaloConv_pi0_unitEta","Primary #pi^{0}",npt,0.,ptmax,150,0.,150.)) ;
684  fOutputContainer->Add(new TH2F("hMC_CaloConv_eta_unitEta","Primary #pi^{0}",npt,0.,ptmax,150,0.,150.)) ;
685  fOutputContainer->Add(new TH1F("hMC_CaloConv_allpi0","Primary #pi^{0}",npt,0.,ptmax)) ;
686  fOutputContainer->Add(new TH1F("hMC_CaloConv_alleta","Primary #pi^{0}",npt,0.,ptmax)) ;
687  fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0PHOSacc","#pi^{0} decayed in PHOS acc",npt,0.,ptmax)) ;
688  fOutputContainer->Add(new TH1F("hMC_CaloConv_etaPHOSacc","#pi^{0} decayed in PHOS acc",npt,0.,ptmax)) ;
689  fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0EMCALacc","#pi^{0} decayed in EMCAL acc",npt,0.,ptmax)) ;
690  fOutputContainer->Add(new TH1F("hMC_CaloConv_etaEMCALacc","#pi^{0} decayed in EMCAL acc",npt,0.,ptmax)) ;
691  fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_PHOS_conv","#pi^{0} decayed in PHOS acc asnd conv. photon",npt,0.,ptmax)) ;
692  fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_PHOS_conv","#pi^{0} decayed in PHOS acc asnd conv. photon",npt,0.,ptmax)) ;
693  fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_EMCAL_conv","#pi^{0} decayed in EMCAL acc asnd conv. photon",npt,0.,ptmax)) ;
694  fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_EMCAL_conv","#pi^{0} decayed in EMCAL acc asnd conv. photon",npt,0.,ptmax)) ;
695  fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_bothphot_conv","#pi^{0} both photons converted",npt,0.,ptmax)) ;
696  fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_bothphot_conv","#pi^{0} both photons converted",npt,0.,ptmax)) ;
697  fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_convPhotInCalo","#pi^{0} photon in calo converted",npt,0.,ptmax)) ;
698  fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_convPhotInCalo","#pi^{0} photon in calo converted",npt,0.,ptmax)) ;
699  fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0onfly_PHOSacc","#pi^{0} photon converted and V0 found",npt,0.,ptmax)) ;
700  fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_v0onfly_PHOSacc","#pi^{0} photon converted and V0 found",npt,0.,ptmax)) ;
701  fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0offline_PHOSacc","#pi^{0} photon converted and V0 found",npt,0.,ptmax)) ;
702  fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_v0offline_PHOSacc","#pi^{0} photon converted and V0 found",npt,0.,ptmax)) ;
703  fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0onfly_EMCALacc","#pi^{0} photon converted and V0 found",npt,0.,ptmax)) ;
704  fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_v0onfly_EMCALacc","#pi^{0} photon converted and V0 found",npt,0.,ptmax)) ;
705  fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0offline_EMCALacc","#pi^{0} photon converted and V0 found",npt,0.,ptmax)) ;
706  fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_v0offline_EMCALacc","#pi^{0} photon converted and V0 found",npt,0.,ptmax)) ;
707  fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0onfly_PHOSclu","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
708  fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_v0onfly_PHOSclu","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
709  fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0onfly_PHOSclu_pid","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
710  fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_v0onfly_PHOSclu_pid","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
711  fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0offline_PHOSclu","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
712  fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_v0offline_PHOSclu","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
713  fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0offline_PHOSclu_pid","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
714  fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_v0offline_PHOSclu_pid","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
715  fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0onfly_PHOSclu_good","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
716  fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_v0onfly_PHOSclu_good","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
717  fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0offline_PHOSclu_good","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
718  fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_v0offline_PHOSclu_good","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
719  fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0onfly_PHOSclu_mod1","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
720  fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_v0onfly_PHOSclu_mod1","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
721  fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0onfly_PHOSclu_mod2","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
722  fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_v0onfly_PHOSclu_mod2","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
723  fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0onfly_PHOSclu_mod3","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
724  fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_v0onfly_PHOSclu_mod3","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
725  fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0onfly_PHOSclu_mod4","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
726  fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_v0onfly_PHOSclu_mod4","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
727  fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0onfly_PHOSclu_mod5","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
728  fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_v0onfly_PHOSclu_mod5","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
729  fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0offline_PHOSclu_mod1","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
730  fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_v0offline_PHOSclu_mod1","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
731  fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0offline_PHOSclu_mod2","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
732  fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_v0offline_PHOSclu_mod2","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
733  fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0offline_PHOSclu_mod3","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
734  fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_v0offline_PHOSclu_mod3","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
735  fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0offline_PHOSclu_mod4","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
736  fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_v0offline_PHOSclu_mod4","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
737  fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0offline_PHOSclu_mod5","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
738  fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_v0offline_PHOSclu_mod5","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
739 
740  fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0on_PHOSclu_ptRec","#pi^{0} V0 and cluster in PHOS found(rec pt)",npt,0.,ptmax)) ;
741  fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_v0on_PHOSclu_ptRec","#pi^{0} V0 and cluster in PHOS found(rec pt)",npt,0.,ptmax)) ;
742  fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0off_PHOSclu_ptRec","#pi^{0} V0 and cluster in PHOS found(rec pt)",npt,0.,ptmax)) ;
743  fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_v0off_PHOSclu_ptRec","#pi^{0} V0 and cluster in PHOS found(rec pt)",npt,0.,ptmax)) ;
744  fOutputContainer->Add(new TH2F("hMC_CaloConv_pi0_v0on_PHOSclu_mvsPt","m vs pt for rec pi0s",400,0.,1.,npt,0.,ptmax)) ;
745  fOutputContainer->Add(new TH2F("hMC_CaloConv_eta_v0on_PHOSclu_mvsPt","m vs pt for rec pi0s",400,0.,1.,npt,0.,ptmax)) ;
746  fOutputContainer->Add(new TH2F("hMC_CaloConv_pi0_v0off_PHOSclu_mvsPt","m vs pt for rec pi0s",400,0.,1.,npt,0.,ptmax)) ;
747  fOutputContainer->Add(new TH2F("hMC_CaloConv_eta_v0off_PHOSclu_mvsPt","m vs pt for rec pi0s",400,0.,1.,npt,0.,ptmax)) ;
748 
749  fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0off_EMCALclu_ptRec","#pi^{0} V0 and cluster in EMCAL found (rec pt)",npt,0.,ptmax)) ;
750  fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0on_EMCALclu_ptRec","#pi^{0} V0 and cluster in EMCAL found (rec pt)",npt,0.,ptmax)) ;
751  fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0onfly_EMCALclu","#pi^{0} V0 and cluster in EMCAL found",npt,0.,ptmax)) ;
752  fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0onfly_EMCALclu_pid","#pi^{0} V0 and cluster in EMCAL found",npt,0.,ptmax)) ;
753  fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0offline_EMCALclu","#pi^{0} V0 and cluster in EMCAL found",npt,0.,ptmax)) ;
754  fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0offline_EMCALclu_pid","#pi^{0} V0 and cluster in EMCAL found",npt,0.,ptmax)) ;
755  fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0onfly_EMCALclu_good","#pi^{0} V0 and cluster in EMCAL found",npt,0.,ptmax)) ;
756  fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0offline_EMCALclu_good","#pi^{0} V0 and cluster in EMCAL found",npt,0.,ptmax)) ;
757  fOutputContainer->Add(new TH2F("hMC_CaloConv_pi0_v0on_EMCALclu_mvsPt","m vs pt for rec pi0s",400,0.,1.,npt,0.,ptmax)) ;
758  fOutputContainer->Add(new TH2F("hMC_CaloConv_pi0_v0off_EMCALclu_mvsPt","m vs pt for rec pi0s",400,0.,1.,npt,0.,ptmax)) ;
759 
760  fOutputContainer->Add(new TH3F("hMC_Resid_PHOS_Phot_mvsPt","m vs pt for rec pi0s",400,0.,1.,npt,0.,ptmax,10,0.,1.)) ;
761  fOutputContainer->Add(new TH3F("hMC_Resid_PHOS_Pi0_mvsPt","m vs pt for rec pi0s",400,0.,1.,npt,0.,ptmax,10,0.,1.)) ;
762  fOutputContainer->Add(new TH3F("hMC_Resid_PHOS_eta_mvsPt","m vs pt for rec pi0s",400,0.,1.,npt,0.,ptmax,10,0.,1.)) ;
763  fOutputContainer->Add(new TH3F("hMC_Resid_PHOS_K_mvsPt","m vs pt for rec pi0s",400,0.,1.,npt,0.,ptmax,10,0.,1.)) ;
764  fOutputContainer->Add(new TH3F("hMC_Resid_PHOS_pi_mvsPt","m vs pt for rec pi0s",400,0.,1.,npt,0.,ptmax,10,0.,1.)) ;
765  fOutputContainer->Add(new TH3F("hMC_Resid_PHOS_pbar_mvsPt","m vs pt for rec pi0s",400,0.,1.,npt,0.,ptmax,10,0.,1.)) ;
766  fOutputContainer->Add(new TH3F("hMC_Resid_PHOS_other_mvsPt","m vs pt for rec pi0s",400,0.,1.,npt,0.,ptmax,10,0.,1.)) ;
767  fOutputContainer->Add(new TH3F("hMC_Resid_EMCAL_Phot_mvsPt","m vs pt for rec pi0s",400,0.,1.,npt,0.,ptmax,10,0.,1.)) ;
768  fOutputContainer->Add(new TH3F("hMC_Resid_EMCAL_Pi0_mvsPt","m vs pt for rec pi0s",400,0.,1.,npt,0.,ptmax,10,0.,1.)) ;
769  fOutputContainer->Add(new TH3F("hMC_Resid_EMCAL_eta_mvsPt","m vs pt for rec pi0s",400,0.,1.,npt,0.,ptmax,10,0.,1.)) ;
770  fOutputContainer->Add(new TH3F("hMC_Resid_EMCAL_K_mvsPt","m vs pt for rec pi0s",400,0.,1.,npt,0.,ptmax,10,0.,1.)) ;
771  fOutputContainer->Add(new TH3F("hMC_Resid_EMCAL_pi_mvsPt","m vs pt for rec pi0s",400,0.,1.,npt,0.,ptmax,10,0.,1.)) ;
772  fOutputContainer->Add(new TH3F("hMC_Resid_EMCAL_pbar_mvsPt","m vs pt for rec pi0s",400,0.,1.,npt,0.,ptmax,10,0.,1.)) ;
773  fOutputContainer->Add(new TH3F("hMC_Resid_EMCAL_other_mvsPt","m vs pt for rec pi0s",400,0.,1.,npt,0.,ptmax,10,0.,1.)) ;
774 
775 
776  fOutputContainer->Add(new TH1F("hMC_CaloConv_phot","Primary photons",npt,0.,ptmax)) ;
777  fOutputContainer->Add(new TH1F("hMC_CaloConv_gammaPHOSacc","Photons in PHOS acc",npt,0.,ptmax)) ;
778  fOutputContainer->Add(new TH1F("hMC_CaloConv_gammaEMCALacc","Photons in EMCAL acc",npt,0.,ptmax)) ;
779  fOutputContainer->Add(new TH1F("hMC_CaloConv_gamma_conv","Converted photons",npt,0.,ptmax)) ;
780  fOutputContainer->Add(new TH1F("hMC_CaloConv_gamma_v0","Converted photons with V0",npt,0.,ptmax)) ;
781  fOutputContainer->Add(new TH2F("hMC_CaloConv_gamma_v0_devsE","Converted photons with V0",200,-1.,1.,npt,0.,ptmax)) ;
782  fOutputContainer->Add(new TH1F("hMC_CaloConv_gamma_PHOSclu","Photons with cluster in PHOS",npt,0.,ptmax)) ;
783  fOutputContainer->Add(new TH1F("hMC_CaloConv_gamma_PHOSclu_dist1","Photons with cluster in PHOS",npt,0.,ptmax)) ;
784  fOutputContainer->Add(new TH1F("hMC_CaloConv_gamma_PHOSclu_dist2","Photons with cluster in PHOS",npt,0.,ptmax)) ;
785  fOutputContainer->Add(new TH1F("hMC_CaloConv_gamma_EMCALclu","Photons with cluster in EMCAL",npt,0.,ptmax)) ;
786  fOutputContainer->Add(new TH1F("hMC_CaloConv_gamma_EMCALclu_dist1","Photons with cluster in EMCAL",npt,0.,ptmax)) ;
787  fOutputContainer->Add(new TH1F("hMC_CaloConv_gamma_EMCALclu_dist2","Photons with cluster in EMCAL",npt,0.,ptmax)) ;
788  fOutputContainer->Add(new TH1F("hMC_CaloConv_gamma_PHOSclu_recE","Photons with cluster in PHOS",npt,0.,ptmax)) ;
789  fOutputContainer->Add(new TH1F("hMC_CaloConv_gamma_EMCALclu_recE","Photons with cluster in EMCAL",npt,0.,ptmax)) ;
790  fOutputContainer->Add(new TH2F("hMC_CaloConv_gamma_PHOSclu_devsE","Photons with cluster in PHOS",200,-1.,1.,npt,0.,ptmax)) ;
791  fOutputContainer->Add(new TH2F("hMC_CaloConv_gamma_EMCALclu_devsE","Photons with cluster in EMCAL",200,-1.,1.,npt,0.,ptmax)) ;
792 
793 
794  //Non-linearity test
795  char keym[55] ;
796  for(Int_t iw=0;iw<10;iw++){ //resolution
797  for(Int_t in=0;in<10;in++){
798  snprintf(keym,55,"hMC_nonlinearity_w%d_n%d",iw,in) ;
799  fOutputContainer->Add(new TH2F(keym,"m vs pt, nonlinearity test" ,200,0.,0.5,npt,0.,ptmax)) ;
800  snprintf(keym,55,"hMC_nonlinearity_ConvPHOS_w%d_n%d",iw,in) ;
801  fOutputContainer->Add(new TH2F(keym,"m vs pt, nonlinearity test" ,200,0.,0.5,npt,0.,ptmax)) ;
802  snprintf(keym,55,"hMC_nonlinearity_EMCAL_w%d_n%d",iw,in) ;
803  fOutputContainer->Add(new TH2F(keym,"m vs pt, nonlinearity test" ,200,0.,0.5,npt,0.,ptmax)) ;
804  snprintf(keym,55,"hMC_CaloConv_pi0_v0onfly_PHOSclu_w%d_n%d",iw,in) ;
805  fOutputContainer->Add(new TH1F(keym,"#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
806  snprintf(keym,55,"hMC_CaloConv_pi0_v0onfly_ConvPHOSclu_w%d_n%d",iw,in) ;
807  fOutputContainer->Add(new TH1F(keym,"#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
808  snprintf(keym,55,"hMC_CaloConv_pi0_v0onfly_EMCALclu_w%d_n%d",iw,in) ;
809  fOutputContainer->Add(new TH1F(keym,"#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
810  }
811  }
812 
813 
814  fOutputContainer->Add(new TH3F("All_chi2_eta_pt","MC chi2 vs eta vs phi",100,0.,100.,200,-2.,2.,npt,0.,ptmax)) ;
815  fOutputContainer->Add(new TH2F("All_w_vs_m","MC w vs m",300,0.,TMath::Pi(),400,0.,1.)) ;
816  fOutputContainer->Add(new TH3F("MC_V0_pt_eta_phi","MC pt vs eta vs phi",npt,0.,ptmax,200,-2.,2.,200,0.,TMath::TwoPi())) ;
817  fOutputContainer->Add(new TH3F("MC_V0_m_eta_pt","MC m vs eta vs phi",400,0.,1.,200,-2.,2.,npt,0.,ptmax)) ;
818  fOutputContainer->Add(new TH3F("MC_V0_chi2_eta_pt","MC chi2 vs eta vs phi",100,0.,100.,200,-2.,2.,npt,0.,ptmax)) ;
819  fOutputContainer->Add(new TH2F("MC_V0_w_vs_m","MC w vs m",300,0.,TMath::Pi(),400,0.,1.)) ;
820 
821  fOutputContainer->SetName(GetName());
822 }
823 //______________________________________________________________________
824  void AliAnalysisTaskCaloConv::InitGeometry()
825 {
826  //If not done yet, create Geometry for PHOS and EMCAL
827  //and read misalignment matrixes from ESD/AOD (AOD not implemented yet)
828  //
829 
830  if(fPHOSgeom && fEMCALgeom){ //already initialized
831  return ;
832  }
833 
834  AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent()) ;
835  if(!esd )
836  AliFatal("Can not read geometry matrixes from ESD/AOD: NO ESD") ;
837  if(!fPHOSgeom){//reading PHOS matrixes
838  fPHOSgeom = new AliPHOSGeoUtils("IHEP","");
839  for(Int_t mod=0; mod<5; mod++){
840  if(esd){
841  const TGeoHMatrix* m=esd->GetPHOSMatrix(mod) ;
842  if(m)
843  fPHOSgeom->SetMisalMatrix(m, mod) ;
844  }
845  }
846  }
847  if(!fEMCALgeom){
848  fEMCALgeom = AliEMCALGeometry::GetInstance("EMCAL_FIRSTYEARV1");
849  for(Int_t mod=0; mod < (fEMCALgeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){ //<---Gustavo, could you check???
850  if(esd){
851  const TGeoHMatrix* m=esd->GetEMCALMatrix(mod) ;
852  if(m)
853  fEMCALgeom->SetMisalMatrix(m, mod) ;
854  }
855  }
856  }
857 }
858 //________________________________________________________________
859 void AliAnalysisTaskCaloConv::SelectPHOSPhotons(){
860  //SelectPHOSPhotons
861  // Loop over all CaloClusters
862  if(fPHOSEvent)
863  fPHOSEvent->Clear() ;
864  else
865  fPHOSEvent = new TClonesArray("TLorentzVector",10) ;
866  Int_t inPHOS = 0;
867  TLorentzVector pi0 ;
868 
869  //vertex
870  Double_t vtx[3];
871  vtx[0] = fESDEvent->GetPrimaryVertex()->GetX();
872  vtx[1] = fESDEvent->GetPrimaryVertex()->GetY();
873  vtx[2] = fESDEvent->GetPrimaryVertex()->GetZ();
874  for (Int_t i=0; i<fESDEvent->GetNumberOfCaloClusters(); i++) {
875  AliESDCaloCluster * clu = fESDEvent->GetCaloCluster(i);
876  if(!clu->IsPHOS())
877  continue ;
878  TLorentzVector p ;
879  clu ->GetMomentum(p ,vtx);
880  if(p.Energy()<0.25){
881  continue ;
882  }
883  if(clu->GetNCells()<=2){
884  continue ;
885  }
886 
887  Bool_t isNeutral = kTRUE ;
888  Bool_t isDispOK = kTRUE ;
889  Bool_t isTOFOK = kTRUE ;
890  Int_t iMod,iX,iZ ;
891 
892  isNeutral = clu->GetEmcCpvDistance()>5. ; //To be improved
893  isDispOK = kFALSE ;
894  Double_t l0=clu->GetM02(),l1=clu->GetM20() ;
895  if(l1>= 0 && l0>= 0 && l1 < 0.1 && l0 < 0.1) isDispOK=kFALSE ;
896  if(l1>= 0 && l0 > 0.5 && l1 < 0.1 && l0 < 1.5) isDispOK=kTRUE ;
897  if(l1>= 0 && l0 > 2.0 && l1 < 0.1 && l0 < 2.7) isDispOK=kFALSE ;
898  if(l1>= 0 && l0 > 2.7 && l1 < 0.1 && l0 < 4.0) isDispOK=kFALSE ;
899  if(l1 > 0.1 && l1 < 0.7 && l0 > 0.7 && l0 < 2.1) isDispOK=kTRUE ;
900  if(l1 > 0.1 && l1 < 0.3 && l0 > 3.0 && l0 < 5.0) isDispOK=kFALSE ;
901  if(l1 > 0.3 && l1 < 0.7 && l0 > 2.5 && l0 < 4.0) isDispOK=kFALSE ;
902  if(l1 > 0.7 && l1 < 1.3 && l0 > 1.0 && l0 < 1.6) isDispOK=kTRUE ;
903  if(l1 > 0.7 && l1 < 1.3 && l0 > 1.6 && l0 < 3.5) isDispOK=kTRUE ;
904  if(l1 > 1.3 && l1 < 3.5 && l0 > 1.3 && l0 < 3.5) isDispOK=kTRUE ;
905 
906  Float_t xyz[3] = {0,0,0};
907  clu->GetPosition(xyz); //Global position in ALICE system
908  TVector3 global(xyz) ;
909  Int_t relid[4] ;
910  if(!fPHOSgeom->GlobalPos2RelId(global,relid)){
911  printf("PHOS_beyond: x=%f, y=%f, z=%f \n",xyz[0],xyz[1],xyz[2]) ;
912  continue ;
913  }
914  iMod=relid[0] ;
915  iX=relid[2];
916  iZ=relid[3] ;
917  if(!IsGoodChannel("PHOS",iMod,iX,iZ))
918  continue ;
919 
920  Bool_t closeToBad=(clu->GetDistanceToBadChannel()>fBadDistCutPHOS) ;
921 
922  p.SetBit(kCaloPIDdisp,isDispOK) ;
923  p.SetBit(kCaloPIDtof,isTOFOK) ;
924  p.SetBit(kCaloPIDneutral,isNeutral) ;
925  p.SetBit(BIT(17+iMod),kTRUE) ;
926  p.SetBit(kCaloDistBad,closeToBad) ;
927  new((*fPHOSEvent)[inPHOS]) TLorentzVector(p) ;
928  fGammaPHOS[inPHOS] = i ;
929  inPHOS++ ;
930 
931 
932  //Single photon spectra
933  Double_t pt= p.Pt() ;
934  TString skey="PHOS_single_"; skey+="mod" ; skey+=iMod ; skey+="_all" ;
935  FillHistogram(skey,pt) ;
936  FillHistogram("PHOS_single_all_mult",pt,fCentr) ;
937  if(isDispOK){
938  skey="PHOS_single_"; skey+="mod" ; skey+=iMod ; skey+="_disp" ;
939  FillHistogram("PHOS_single_disp_mult",pt,fCentr) ;
940  FillHistogram(skey,pt) ;
941  }
942  if(isNeutral){
943  skey="PHOS_single_"; skey+="mod" ; skey+=iMod ; skey+="_neutral" ;
944  FillHistogram(skey,pt) ;
945  FillHistogram("PHOS_single_neu_mult",pt,fCentr) ;
946  }
947  if(isNeutral && isDispOK){
948  skey="PHOS_single_"; skey+="mod" ; skey+=iMod ; skey+="_dispneutral" ;
949  FillHistogram(skey,pt) ;
950  }
951  //Distance to bad channel
952  if(clu->GetDistanceToBadChannel()>fBadDistCutPHOS){
953  skey="PHOS_single_"; skey+="mod" ; skey+=iMod ; skey+="_dist1" ;
954  FillHistogram(skey,pt) ;
955  }
956  if(clu->GetDistanceToBadChannel()>2.*fBadDistCutPHOS){
957  skey="PHOS_single_"; skey+="mod" ; skey+=iMod ; skey+="_dist2" ;
958  FillHistogram(skey,pt) ;
959  }
960 
961  //Fill QA
962  if(clu->E()>0.5){
963  skey="hQA_PHOS_mod"; skey+=iMod; skey+="_soft" ;
964  FillHistogram(skey,iX-0.5, iZ-0.5,1.) ;
965  if(clu->E()>1.5){
966  skey="hQA_PHOS_mod"; skey+=iMod; skey+="_hard" ;
967  FillHistogram(skey,iX-0.5, iZ-0.5,1.) ;
968  }
969  }
970 
971  //Fill histogams for calibration
972  if(clu->E()<fPi0Thresh1 ) continue;
973  for(Int_t iconv=0;iconv<fConvEvent->GetEntriesFast();iconv++){
974  TLorentzVector *gammaConv=static_cast<TLorentzVector*>(fConvEvent->At(iconv)) ;
975  pi0=*gammaConv+p;
976  skey="PHOS_"; skey+="mod" ; skey+=iMod ; skey+="_th1" ;
977  FillHistogram(skey,iX-0.5, iZ-0.5,pi0.M());
978  if(isNeutral){
979  skey="PHOS_"; skey+="mod"; skey+=iMod ; skey+="_th2" ;
980  FillHistogram(skey,iX-0.5, iZ-0.5,pi0.M());
981  }
982  }
983  }
984 }
985 //____________________________________________________________
986 void AliAnalysisTaskCaloConv::SelectEMCALPhotons(){
987  //SelectEMCALPhotons
988  // Loop over all CaloClusters
989  if(fEMCALEvent)
990  fEMCALEvent->Clear() ;
991  else
992  fEMCALEvent = new TClonesArray("TLorentzVector",10) ;
993  Int_t inEMCAL = 0 ; //, inEMCALRecal=0;
994  TLorentzVector pi0 ;
995 
996  //vertex
997  Double_t vtx[3]={0.,0.,0.};
998  vtx[0] = fESDEvent->GetPrimaryVertex()->GetX();
999  vtx[1] = fESDEvent->GetPrimaryVertex()->GetY();
1000  vtx[2] = fESDEvent->GetPrimaryVertex()->GetZ();
1001  Int_t nEMCAL=0 ; //For QA
1002  for (Int_t i=0; i<fESDEvent->GetNumberOfCaloClusters(); i++) {
1003  AliESDCaloCluster * clu = fESDEvent->GetCaloCluster(i);
1004  if(clu->IsPHOS())
1005  continue ;
1006  TLorentzVector p ;
1007  TLorentzVector pRecal ;
1008  clu ->GetMomentum(p ,vtx);
1009  Bool_t isNeutral = kTRUE ;
1010  Bool_t isDispOK = kTRUE ;
1011  Bool_t isTOFOK = kTRUE ;
1012  Int_t iMod,iX,iZ ;
1013 
1014  if(clu->E()>0.1 ){
1015  nEMCAL++ ;
1016  isNeutral = clu->GetEmcCpvDistance()>10. ; //To be improved
1017  if(clu->GetTOF()>550.e-9 && clu->GetTOF()<750.e-9)
1018  isTOFOK=kTRUE ;
1019  else
1020  isTOFOK=kFALSE ;
1021  Float_t phi = p.Phi();
1022  if(phi < 0) phi+=TMath::TwoPi();
1023  Int_t absId = -1;
1024  fEMCALgeom->GetAbsCellIdFromEtaPhi(p.Eta(),phi, absId);
1025  iMod=fEMCALgeom->GetSuperModuleNumber(absId) ;
1026  Int_t imod = -1, iTower = -1, iIphi = -1, iIeta = -1, iphi = -1, ieta = -1;
1027  fEMCALgeom->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
1028  fEMCALgeom->GetCellPhiEtaIndexInSModule(imod,iTower, iIphi, iIeta,iphi,ieta);
1029  if(imod<0 || imod>5){
1030  printf("EMCAL: Beyond the geometry!\n") ;
1031  printf("phi=%f, eta=%f, absId=%d, SM=%d \n",p.Eta(),phi, absId, imod) ;
1032  continue ;
1033  }
1034 
1035  iX=iphi+1 ;
1036  iZ=ieta+1 ;
1037  if(!IsGoodChannel("EMCAL",iMod,iX,iZ))
1038  continue ;
1039  p.SetBit(kCaloPIDdisp,isDispOK) ;
1040  p.SetBit(kCaloPIDtof,isTOFOK) ;
1041  p.SetBit(kCaloPIDneutral,isNeutral) ;
1042  p.SetBit(BIT(17+imod),kTRUE) ;
1043  new((*fEMCALEvent)[inEMCAL]) TLorentzVector(p) ;
1044  fGammaEMCAL[inEMCAL] = i ;
1045  inEMCAL++ ;
1046 
1047 
1048  //Fill QA histograms
1049  if(clu->E()>0.5 && iMod>=0){ //Sometimes modules is negative not found??
1050  TString skey="hQA_EMCAL_SM";skey+=iMod ; skey+="_soft" ;
1051  FillHistogram(skey,iX-0.5, iZ-0.5,1.) ;
1052  if(clu->E()>1.5){
1053  skey="hQA_EMCAL_SM";skey+=iMod ; skey+="_hard" ;
1054  FillHistogram(skey,iX-0.5, iZ-0.5,1.) ;
1055  }
1056  }
1057 
1058  //Fill histograms for recalibration
1059  if(clu->E()<fPi0Thresh1) continue ;
1060  for(Int_t iconv=0;iconv<fConvEvent->GetEntriesFast();iconv++){
1061  TLorentzVector *gammaConv=static_cast<TLorentzVector*>(fConvEvent->At(iconv)) ;
1062  pi0=*gammaConv+p;
1063  TString skey="EMCAL_"; skey+="mod" ; skey+=iMod ; skey+="_th1" ;
1064  FillHistogram(skey,iX-0.5, iZ-0.5,pi0.M());
1065  if(clu->E()>fPi0Thresh2){
1066  skey="EMCAL_"; skey+="mod"; skey+=iMod ; skey+="_th2" ;
1067  FillHistogram(skey,iX-0.5, iZ-0.5,pi0.M());
1068  }
1069  }
1070  }
1071  }
1072 }
1073 //______________________________________________________________________
1074 void AliAnalysisTaskCaloConv::SelectConvPhotons(){
1075  //Fill list of conversion photons
1076  //that is scan v0s and select photon-like
1077 
1078  //set some constants
1079  const Double_t cutSigmaMass=0.0001; //Constraint on photon mass
1080  const Bool_t useImprovedVertex=kTRUE ; //Use verted with converted photon?
1081 // const Double_t zrSlope = TMath::Tan(2*TMath::ATan(TMath::Exp(-fetaCut)));
1082  const Double_t zrSlope12 = TMath::Tan(2*TMath::ATan(TMath::Exp(-1.2)));
1083  const Double_t zrSlope09 = TMath::Tan(2*TMath::ATan(TMath::Exp(-0.9)));
1084  const Double_t zOffset = 7.;
1085 
1086  if(!fConvEvent)
1087  fConvEvent = new TClonesArray("TLorentzVector",10) ;
1088  else
1089  fConvEvent->Clear() ;
1090 
1091  //No primary vertex in event => scip
1092  if(fESDEvent->GetPrimaryVertex()->GetNContributors()<=0) {
1093  return;
1094  }
1095 
1096  Int_t inConv=0 ;
1097  for(Int_t iv0=0; iv0<fESDEvent->GetNumberOfV0s();iv0++){
1098  AliESDv0 * v0 = fESDEvent->GetV0(iv0) ;
1099 
1100  AliESDtrack * pos = fESDEvent->GetTrack(v0->GetPindex()) ;
1101  AliESDtrack * neg = fESDEvent->GetTrack(v0->GetNindex()) ;
1102  const AliExternalTrackParam * paramPos = v0->GetParamP() ;
1103  const AliExternalTrackParam * paramNeg = v0->GetParamN() ;
1104  if(pos->GetSign() <0){//change tracks
1105  pos=neg ;
1106  neg=fESDEvent->GetTrack(v0->GetPindex()) ;
1107  paramPos=paramNeg ;
1108  paramNeg=v0->GetParamP() ;
1109  }
1110  AliKFParticle negKF(*paramNeg,11);
1111  AliKFParticle posKF(*paramPos,-11);
1112  AliKFParticle photKF(negKF,posKF) ;
1113  photKF.SetMassConstraint(0,cutSigmaMass);
1114 
1115  if(useImprovedVertex){
1116  AliKFVertex primaryVertexImproved(*(fESDEvent->GetPrimaryVertex()));
1117  //if Vtx do created
1118  if(primaryVertexImproved.GetNContributors()>1){
1119  primaryVertexImproved+=photKF;
1120  photKF.SetProductionVertex(primaryVertexImproved);
1121  }
1122  }
1123  Double_t m=0., width=0. ;
1124  photKF.GetMass(m,width);
1125 
1126  TLorentzVector photLV;
1127 // photLV.SetXYZM(negKF.Px()+posKF.Px(),negKF.Py()+posKF.Py(),negKF.Pz()+negKF.Pz(),0.) ;
1128 // photLV.SetXYZT(photKF.GetPx(),photKF.GetPy(),photKF.GetPz(),photKF.GetE()) ;
1129  photLV.SetXYZM(photKF.GetPx(),photKF.GetPy(),photKF.GetPz(),0.) ; //Produces slightly better pi0 width
1130 
1131  //Parameters for correction function
1132  Double_t a[3]={photLV.Pt(),photLV.Eta(),m} ;
1133  if(fToUseCF)
1134  fConvCFCont->Fill(a,0) ;
1135 
1136  //select V0 finder
1137  Bool_t isOnFly=kTRUE ;
1138  //select V0Finder
1139  if (v0->GetOnFlyStatus()){
1140  if(fToUseCF)
1141  fConvCFCont->Fill(a,1) ;
1142  }
1143  else{
1144  isOnFly=kFALSE ;
1145  if(fToUseCF)
1146  fConvCFCont->Fill(a,2) ;
1147  }
1148 
1149  //Number of TPC clusters
1150  if(neg->GetNcls(1) <2 || pos->GetNcls(1) <2){
1151  continue ;
1152  }
1153 
1154  //remove like sign pairs
1155  if(pos->GetSign() == neg->GetSign()){
1156  continue ;
1157  }
1158  if(fToUseCF){
1159  if(isOnFly)
1160  fConvCFCont->Fill(a,3) ;
1161  else
1162  fConvCFCont->Fill(a,4) ;
1163  }
1164 
1165  if( !(pos->GetStatus() & AliESDtrack::kTPCrefit) ||
1166  !(neg->GetStatus() & AliESDtrack::kTPCrefit) ){
1167  continue;
1168  }
1169  if(fToUseCF){
1170  if(isOnFly)
1171  fConvCFCont->Fill(a,5) ;
1172  else
1173  fConvCFCont->Fill(a,6) ;
1174  }
1175 
1176  if( neg->GetKinkIndex(0) > 0 ||
1177  pos->GetKinkIndex(0) > 0) {
1178  continue ;
1179  }
1180 
1181  //First rough PID
1182  if( fESDpid->NumberOfSigmasTPC(pos,AliPID::kElectron)<fnSigmaBelowElectronLine ||
1183  fESDpid->NumberOfSigmasTPC(pos,AliPID::kElectron)>fnSigmaAboveElectronLine ||
1184  fESDpid->NumberOfSigmasTPC(neg,AliPID::kElectron)<fnSigmaBelowElectronLine ||
1185  fESDpid->NumberOfSigmasTPC(neg,AliPID::kElectron)>fnSigmaAboveElectronLine ){
1186  continue ;
1187  }
1188  const Double_t minPnSigmaAbovePionLine = 1. ;
1189  const Double_t maxPnSigmaAbovePionLine = 3. ;
1190  const Double_t nSigmaAbovePionLine = 0 ;
1191  if(pos->P()>minPnSigmaAbovePionLine && pos->P()<maxPnSigmaAbovePionLine ){
1192  if(fESDpid->NumberOfSigmasTPC(pos,AliPID::kPion)<nSigmaAbovePionLine){
1193  continue ;
1194  }
1195  }
1196  if(neg->P()>minPnSigmaAbovePionLine && neg->P()<maxPnSigmaAbovePionLine){
1197  if(fESDpid->NumberOfSigmasTPC(neg,AliPID::kPion)<nSigmaAbovePionLine){
1198  continue ;
1199  }
1200  }
1201  //Strict dEdx
1202  Bool_t isdEdx=kTRUE;
1203  if(pos->P()>minPnSigmaAbovePionLine && pos->P()<maxPnSigmaAbovePionLine ){
1204  if(fESDpid->NumberOfSigmasTPC(pos,AliPID::kPion)<2.){
1205  isdEdx=kFALSE;
1206  }
1207  }
1208  if(neg->P()>minPnSigmaAbovePionLine && neg->P()<maxPnSigmaAbovePionLine){
1209  if(fESDpid->NumberOfSigmasTPC(neg,AliPID::kPion)<2.){
1210  isdEdx=kFALSE;
1211  }
1212  }
1213 
1214 
1215  //Kaon rejection
1216  const Double_t minPKaonRejection=1.5 ;
1217  const Double_t sigmaAroundLine=1. ;
1218  if(neg->P()<minPKaonRejection ){
1219  if(TMath::Abs(fESDpid->NumberOfSigmasTPC(neg,AliPID::kKaon))<sigmaAroundLine){
1220  isdEdx=kFALSE;
1221  }
1222  }
1223  if(pos->P()<minPKaonRejection ){
1224  if(TMath::Abs(fESDpid->NumberOfSigmasTPC(pos,AliPID::kKaon))<sigmaAroundLine){
1225  isdEdx=kFALSE;
1226  }
1227  }
1228 
1229  //Proton rejection
1230  const Double_t minPProtonRejection=2. ;
1231  if(neg->P()<minPProtonRejection){
1232  if(TMath::Abs(fESDpid->NumberOfSigmasTPC(neg,AliPID::kProton))<sigmaAroundLine){
1233  isdEdx=kFALSE;
1234  }
1235  }
1236  if(pos->P()<minPProtonRejection ){
1237  if(TMath::Abs(fESDpid->NumberOfSigmasTPC(pos,AliPID::kProton))<sigmaAroundLine){
1238  isdEdx=kFALSE;
1239  }
1240  }
1241 
1242  const Double_t minPPionRejection=0.5 ;
1243  if(neg->P()<minPPionRejection ){
1244  if(TMath::Abs(fESDpid->NumberOfSigmasTPC(neg,AliPID::kPion))<sigmaAroundLine){
1245  isdEdx=kFALSE;
1246  }
1247  }
1248  if(pos->P()<minPPionRejection ){
1249  if( TMath::Abs(fESDpid->NumberOfSigmasTPC(pos,AliPID::kPion))<sigmaAroundLine){
1250  isdEdx=kFALSE;
1251  }
1252  }
1253 
1254 
1255  if(isdEdx){
1256  FillHistogram("hdEdx",paramPos->GetP(),pos->GetTPCsignal()) ;
1257  FillHistogram("hdEdx",paramNeg->GetP(),neg->GetTPCsignal()) ;
1258  }
1259 
1260 
1261  //Check the pid probability
1262  Bool_t isProb=kTRUE ;
1263  Double_t posProbArray[10];
1264  Double_t negProbArray[10];
1265  neg->GetTPCpid(negProbArray);
1266  pos->GetTPCpid(posProbArray);
1267  if(negProbArray[AliPID::kElectron]<fprobCut || posProbArray[AliPID::kElectron]<fprobCut){
1268  isProb=kFALSE ;
1269  }
1270  if(fToUseCF){
1271  if(isOnFly)
1272  fConvCFCont->Fill(a,9) ;
1273  else
1274  fConvCFCont->Fill(a,10) ;
1275  }
1276 
1277  Double_t v0x=0.,v0y=0.,v0z=0.;
1278  v0->GetXYZ(v0x,v0y,v0z) ;
1279  Double_t r=TMath::Sqrt(v0x*v0x + v0y*v0y) ;
1280  //Remove Dalitz
1281  const Double_t rMin=2.8 ;
1282  if(r<rMin)
1283  continue ;
1284  if(r>fmaxR){ // cuts on distance from collision point
1285  continue;
1286  }
1287  Bool_t isStrictR=kFALSE ;
1288  if(r<120.)
1289  isStrictR=kTRUE ;
1290  if(fToUseCF){
1291  if(isOnFly)
1292  fConvCFCont->Fill(a,11) ;
1293  else
1294  fConvCFCont->Fill(a,12) ;
1295  }
1296 
1297 
1298  if((TMath::Abs(v0z)*zrSlope12)-zOffset > r ){ // cuts out regions where we do not reconstruct
1299  continue;
1300  }
1301  if(fToUseCF){
1302  if(isOnFly)
1303  fConvCFCont->Fill(a,13) ;
1304  else
1305  fConvCFCont->Fill(a,14) ;
1306  }
1307 
1308  if(TMath::Abs(v0z) > fmaxZ ){ // cuts out regions where we do not reconstruct
1309  continue;
1310  }
1311  Bool_t isStrictZ=kFALSE ;
1312  if((TMath::Abs(v0z)*zrSlope09)-zOffset < r )
1313  isStrictZ=kTRUE ;
1314 
1315  if(fToUseCF){
1316  if(isOnFly)
1317  fConvCFCont->Fill(a,15) ;
1318  else
1319  fConvCFCont->Fill(a,16) ;
1320  }
1321 
1322  if(photKF.GetNDF()<=0){
1323  continue;
1324  }
1325  if(fToUseCF){
1326  if(isOnFly)
1327  fConvCFCont->Fill(a,17) ;
1328  else
1329  fConvCFCont->Fill(a,18) ;
1330  }
1331 
1332  Double_t chi2V0 = photKF.GetChi2()/photKF.GetNDF();
1333  FillHistogram("All_chi2_eta_pt",chi2V0,photLV.Eta(),photLV.Pt()) ;
1334 
1335  if(chi2V0 > fchi2CutConversion || chi2V0 <=0){
1336  continue;
1337  }
1338  Bool_t isStrictChi=kFALSE ;
1339  if(chi2V0 < 0.7*fchi2CutConversion && chi2V0 >0){
1340  isStrictChi=kTRUE;
1341  }
1342 
1343  if(fToUseCF){
1344  if(isOnFly)
1345  fConvCFCont->Fill(a,19) ;
1346  else
1347  fConvCFCont->Fill(a,20) ;
1348  }
1349 
1350  const Double_t wideEtaCut=1.2 ;
1351  if(TMath::Abs(photLV.Eta())> wideEtaCut){
1352  continue;
1353  }
1354  if(TMath::Abs(paramPos->Eta())> wideEtaCut ||
1355  TMath::Abs(paramNeg->Eta())> wideEtaCut ){
1356  continue ;
1357  }
1358 
1359  Bool_t isWideEta=kTRUE ;
1360  if(TMath::Abs(photLV.Eta())< fetaCut &&
1361  TMath::Abs(paramPos->Eta())<fetaCut &&
1362  TMath::Abs(paramNeg->Eta()) < fetaCut){
1363  isWideEta=kFALSE;
1364  }
1365 
1366 
1367  if(photLV.Pt()<fptCut){
1368  continue;
1369  }
1370  if(fToUseCF){
1371  if(isOnFly)
1372  fConvCFCont->Fill(a,21) ;
1373  else
1374  fConvCFCont->Fill(a,22) ;
1375  }
1376 
1377 
1378  //Just QA plot
1379  if(photLV.Pt()>0.5){
1380  Double_t phi=photLV.Phi() ;
1381  while(phi<0.)phi+=TMath::TwoPi() ;
1382  while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
1383  FillHistogram("hQA_ConvPhiEta",phi,photLV.Eta()) ;
1384  }
1385 
1386  Double_t w=PlanarityAngle(paramPos,paramNeg) ;
1387  Bool_t isPlanarityCut = (0.08-0.22*w > m || 0.15*(w-2.4)>m) ;
1388  FillHistogram("All_w_vs_m",w,m) ;
1389 
1390  const Double_t armenterosAlphaCut=0.05 ;
1391  Double_t armenterosQtAlfa[2]={0.,0.} ;
1392  GetArmenterosQtAlfa(&posKF, &negKF, &photKF, armenterosQtAlfa ) ;
1393  Bool_t isArmQt=(armenterosQtAlfa[1]<armenterosAlphaCut) ;
1394 
1395  photLV.SetBit(kConvOnFly,isOnFly) ;
1396  photLV.SetBit(kConvArmQt,isArmQt) ;
1397  photLV.SetBit(kConvdEdx,isdEdx) ;
1398  photLV.SetBit(kConvProb,isProb) ;
1399  photLV.SetBit(kConvR,isStrictR) ;
1400  photLV.SetBit(kConvZR,isStrictZ) ;
1401  photLV.SetBit(kConvNDF,isStrictChi) ;
1402  photLV.SetBit(kConvEta,isWideEta) ;
1403  photLV.SetBit(kConvPlan,isPlanarityCut) ;
1404 
1405  new((*fConvEvent)[inConv]) TLorentzVector(photLV) ;
1406  fGammaV0s[inConv] = iv0 ;
1407  inConv++ ;
1408 
1409  //Single photon spectrum
1410  Double_t pt=photLV.Pt() ;
1411  if(isOnFly){
1412  //Default
1413  if(!isWideEta)
1414  FillHistogram("Single_conv_OnFly",pt) ;
1415  if(isdEdx && !isWideEta)
1416  FillHistogram("Single_conv_On_dEdx",pt) ;
1417  if(!isWideEta && isStrictR)
1418  FillHistogram("Single_conv_On_R120",pt) ;
1419  if( !isWideEta && isStrictZ)
1420  FillHistogram("Single_conv_On_Z",pt) ;
1421  if(!isWideEta && isStrictChi)
1422  FillHistogram("Single_conv_On_chi",pt) ;
1423  if(1)
1424  FillHistogram("Single_conv_On_Eta",pt) ;
1425  if(!isWideEta && isPlanarityCut)
1426  FillHistogram("Single_conv_On_Wcut",pt) ;
1427  if(!isWideEta && isArmQt)
1428  FillHistogram("Single_conv_On_ArmQt",pt) ;
1429  }
1430  else{
1431  if(!isWideEta)
1432  FillHistogram("Single_conv_Offline",pt) ;
1433  if(isdEdx && !isWideEta)
1434  FillHistogram("Single_conv_Off_dEdx",pt) ;
1435  if(!isWideEta && isStrictR)
1436  FillHistogram("Single_conv_Off_R120",pt) ;
1437  if(!isWideEta && isStrictZ)
1438  FillHistogram("Single_conv_Off_Z",pt) ;
1439  if(!isWideEta && isStrictChi)
1440  FillHistogram("Single_conv_Off_chi",pt) ;
1441  if(1)
1442  FillHistogram("Single_conv_Off_Eta",pt) ;
1443  if(!isWideEta && isPlanarityCut)
1444  FillHistogram("Single_conv_Off_Wcut",pt) ;
1445  if(!isWideEta && isArmQt)
1446  FillHistogram("Single_conv_Off_ArmQt",pt) ;
1447  }
1448 
1449  //Fill MC information
1450  if(fMCEvent){
1451  TParticle * negativeMC = fMCEvent->Particle(TMath::Abs(neg->GetLabel()));
1452  TParticle * positiveMC = fMCEvent->Particle(TMath::Abs(pos->GetLabel()));
1453 
1454  if(!negativeMC || !positiveMC)
1455  continue ;
1456 
1457  if(negativeMC->GetMother(0) != positiveMC->GetMother(0))
1458  continue ;
1459 
1460  if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){
1461  continue;
1462  }
1463  if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){
1464  continue;
1465  }
1466 
1467  TParticle * v0Gamma = fMCEvent->Particle(negativeMC->GetMother(0));
1468 
1469  if(negativeMC->GetUniqueID() != 5 || positiveMC->GetUniqueID() !=5){ // id 5 is conversion
1470  continue;
1471  }
1472  if(v0Gamma->GetPdgCode() == 22){
1473  FillHistogram("MC_V0_pt_eta_phi",v0Gamma->Pt(),v0Gamma->Eta(),v0Gamma->Phi()) ;
1474  FillHistogram("MC_V0_m_eta_pt",m,v0Gamma->Eta(),v0Gamma->Pt()) ;
1475  FillHistogram("MC_V0_chi2_eta_pt",chi2V0,v0Gamma->Eta(),v0Gamma->Pt()) ;
1476  FillHistogram("MC_V0_w_vs_m",w,m) ;
1477  }
1478  }
1479  }
1480 }
1481 //______________________________________________________________________
1482 void AliAnalysisTaskCaloConv::FillRealMixed(){
1483  // Fills real (same event) and Mixed (different events) inv.mass dsitributions
1484  // Moves current event to the list of stored events at the end
1485 
1486  Double_t vtx[3];
1487  vtx[0] = fESDEvent->GetPrimaryVertex()->GetX();
1488  vtx[1] = fESDEvent->GetPrimaryVertex()->GetY();
1489  vtx[2] = fESDEvent->GetPrimaryVertex()->GetZ();
1490 
1491  //Vtx class z-bin
1492  Int_t zvtx = (Int_t)((vtx[2]+10.)/2.) ;
1493  if(zvtx<0)zvtx=0 ;
1494  if(zvtx>9)zvtx=9 ;
1495 
1496  Double_t run = fESDEvent->GetRunNumber()+0.5;
1497  TString trigClasses = fESDEvent->GetFiredTriggerClasses();
1498  if(trigClasses.Contains("CINT1B-ABCE-NOPF-ALL"))
1499  FillHistogram("hRunTrigger",run,0.5) ;
1500  else
1501  FillHistogram("hRunTrigger",run,1.5) ;
1502 
1503  FillHistogram("hEvents",0.5) ;
1504  FillHistogram("hVtxBin",zvtx-0.5) ;
1505  FillHistogram("hRunEvents",run) ;
1506 
1507  //check if containers for mixed events exist
1508  //create new if necessary
1509  if(!fPHOSEvents[zvtx]) fPHOSEvents[zvtx]=new TList() ;
1510  if(!fEMCALEvents[zvtx]) fEMCALEvents[zvtx]=new TList() ;
1511  if(!fConvEvents[zvtx]) fConvEvents[zvtx]=new TList() ;
1512 
1513  Int_t nPHOS=fPHOSEvent->GetEntriesFast() ;
1514  Int_t nEMCAL=fEMCALEvent->GetEntriesFast() ;
1515  Int_t nConv = fConvEvent->GetEntriesFast() ;
1516  //Some QA histograms
1517  //Calculate number of good converion photons
1518  Int_t nConvGood=0 ;
1519  for(Int_t iConv = 0; iConv<nConv; iConv++){
1520  TLorentzVector * cnv = static_cast<TLorentzVector*>(fConvEvent->At(iConv)) ;
1521  if(cnv->TestBit(kConvOnFly) && !cnv->TestBit(kConvEta)){
1522  nConvGood++ ;
1523  }
1524  }
1525  FillHistogram("hRunConvs",run,double(nConvGood)) ;
1526  FillHistogram("hRunPHOS", run,double(nPHOS)) ;
1527  FillHistogram("hRunEMCAL",run,double(nEMCAL)) ;
1528 
1529  //Fill Real distributions
1530  for(Int_t iPHOS=0; iPHOS<nPHOS;iPHOS++){
1531  TLorentzVector * cal = static_cast<TLorentzVector*>(fPHOSEvent->At(iPHOS)) ;
1532  for(Int_t iConv = 0; iConv<nConv; iConv++){
1533  TLorentzVector * cnv = static_cast<TLorentzVector*>(fConvEvent->At(iConv)) ;
1534  TLorentzVector pi=*cal + *cnv ;
1535  Double_t alpha=TMath::Abs(cal->Energy()-cnv->Energy())/(cal->Energy()+cnv->Energy()) ;
1536  if(cnv->TestBit(kConvOnFly) && !cnv->TestBit(kConvEta)){
1537  FillHistogram("PHOS_Re_mvsPt_all",pi.M(),pi.Pt()) ;
1538  char keym[55] ;
1539  for(Int_t iw=0;iw<10;iw++){ //resolution
1540  for(Int_t in=0;in<10;in++){
1541  snprintf(keym,55,"hMC_nonlinearity_w%d_n%d",iw,in) ;
1542  Double_t mMod=0.,ptMod=0. ;
1543  Recalibrate(mMod, ptMod, cal, cnv, iw, in) ;
1544  FillHistogram(keym,mMod,ptMod) ;
1545  snprintf(keym,55,"hMC_nonlinearity_ConvPHOS_w%d_n%d",iw,in) ;
1546  RecalibrateConvPHOS(mMod, ptMod, cal, cnv, iw, in) ;
1547  FillHistogram(keym,mMod,ptMod) ;
1548 
1549  }
1550  }
1551  if(cal->TestBit(kCaloDistBad))
1552  FillHistogram("PHOS_Re_mvsPt_all_dist",pi.M(),pi.Pt()) ;
1553  FillHistogram("PHOS_Re_mvsPt_E",pi.M(),pi.Pt(),cal->Energy()) ;
1554  FillHistogram("PHOS_Re_mvsPt_alpha",pi.M(),pi.Pt(),alpha) ;
1555  if(cal->TestBit(kCaloPIDdisp))
1556  FillHistogram("PHOS_Re_mvsPt_Disp",pi.M(),pi.Pt()) ;
1557  if(cal->TestBit(kCaloPIDtof))
1558  FillHistogram("PHOS_Re_mvsPt_TOF",pi.M(),pi.Pt()) ;
1559  if(cal->TestBit(kCaloPIDneutral))
1560  FillHistogram("PHOS_Re_mvsPt_Neutral",pi.M(),pi.Pt()) ;
1561  if(cal->TestBit(kCaloPIDneutral) && cal->TestBit(kCaloPIDdisp))
1562  FillHistogram("PHOS_Re_mvsPt_DispNeutral",pi.M(),pi.Pt()) ;
1563  }
1564  //Vary Conversion cuts
1565  if(cnv->TestBit(kConvOnFly)){
1566  //Default
1567  if(!cnv->TestBit(kConvEta)){
1568  FillHistogram("PHOS_Re_mvsPt_OnFly",pi.M(),pi.Pt()) ;
1569  FillHistogram("PHOS_Re_mvsPt_OnFly_mult",pi.M(),pi.Pt(),fCentr) ;
1570  }
1571  if(cnv->TestBit(kConvdEdx) && !cnv->TestBit(kConvEta))
1572  FillHistogram("PHOS_Re_mvsPt_On_dEdx",pi.M(),pi.Pt()) ;
1573  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvR))
1574  FillHistogram("PHOS_Re_mvsPt_On_R120",pi.M(),pi.Pt()) ;
1575  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvZR))
1576  FillHistogram("PHOS_Re_mvsPt_On_Z",pi.M(),pi.Pt()) ;
1577  if( !cnv->TestBit(kConvEta) && cnv->TestBit(kConvNDF))
1578  FillHistogram("PHOS_Re_mvsPt_On_chi",pi.M(),pi.Pt()) ;
1579  if(1)
1580  FillHistogram("PHOS_Re_mvsPt_On_Eta",pi.M(),pi.Pt()) ;
1581  if( !cnv->TestBit(kConvEta) && cnv->TestBit(kConvPlan)){
1582  FillHistogram("PHOS_Re_mvsPt_On_Wcut",pi.M(),pi.Pt()) ;
1583  if(cal->TestBit(kCaloPIDneutral))
1584  FillHistogram("PHOS_Re_mvsPt_On_Wcut_Neu",pi.M(),pi.Pt()) ;
1585  }
1586  if( !cnv->TestBit(kConvEta) && cnv->TestBit(kConvArmQt))
1587  FillHistogram("PHOS_Re_mvsPt_On_ArmQt",pi.M(),pi.Pt()) ;
1588  }
1589  else{
1590  //Default
1591  if(!cnv->TestBit(kConvEta)){
1592  FillHistogram("PHOS_Re_mvsPt_Offline",pi.M(),pi.Pt()) ;
1593  FillHistogram("PHOS_Re_mvsPt_Offline_mult",pi.M(),pi.Pt(),fCentr) ;
1594  }
1595  if(cnv->TestBit(kConvdEdx) && !cnv->TestBit(kConvEta))
1596  FillHistogram("PHOS_Re_mvsPt_Off_dEdx",pi.M(),pi.Pt()) ;
1597  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvR))
1598  FillHistogram("PHOS_Re_mvsPt_Off_R120",pi.M(),pi.Pt()) ;
1599  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvZR))
1600  FillHistogram("PHOS_Re_mvsPt_Off_Z",pi.M(),pi.Pt()) ;
1601  if( !cnv->TestBit(kConvEta) && cnv->TestBit(kConvNDF))
1602  FillHistogram("PHOS_Re_mvsPt_Off_chi",pi.M(),pi.Pt()) ;
1603  if(1)
1604  FillHistogram("PHOS_Re_mvsPt_Off_Eta",pi.M(),pi.Pt()) ;
1605  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvPlan)){
1606  FillHistogram("PHOS_Re_mvsPt_Off_Wcut",pi.M(),pi.Pt()) ;
1607  if(cal->TestBit(kCaloPIDneutral))
1608  FillHistogram("PHOS_Re_mvsPt_Off_Wcut_Neu",pi.M(),pi.Pt()) ;
1609  }
1610  if( !cnv->TestBit(kConvEta) && cnv->TestBit(kConvArmQt))
1611  FillHistogram("PHOS_Re_mvsPt_Off_ArmQt",pi.M(),pi.Pt()) ;
1612  }
1613  }
1614  }
1615  //PHOS module-dependent histograms
1616  for(Int_t iPHOS=0; iPHOS<nPHOS;iPHOS++){
1617  TLorentzVector * cal = static_cast<TLorentzVector*>(fPHOSEvent->At(iPHOS)) ;
1618  Int_t mod=1;
1619  while(!cal->TestBit(BIT(17+mod)) && mod<5)mod++ ;
1620  TString base("PHOS_Re_mvsPt_mod") ; base+=mod ;
1621  TString full ;
1622  for(Int_t iConv = 0; iConv<nConv; iConv++){
1623  TLorentzVector * cnv = static_cast<TLorentzVector*>(fConvEvent->At(iConv)) ;
1624  TLorentzVector pi=*cal + *cnv ;
1625  full=base ; full+="_single" ;
1626  if(cnv->TestBit(kConvOnFly) && !cnv->TestBit(kConvEta)){
1627  FillHistogram(full,pi.M(),cal->Pt()) ;
1628  full=base ; full+="_all" ;
1629  FillHistogram(full,pi.M(),pi.Pt()) ;
1630  }
1631  }
1632  }
1633  for(Int_t iEMCAL=0; iEMCAL<nEMCAL;iEMCAL++){
1634  TLorentzVector * cal = static_cast<TLorentzVector*>(fEMCALEvent->At(iEMCAL)) ;
1635  Int_t mod=0;
1636  while(!cal->TestBit(BIT(17+mod)) && mod<6)mod++ ;
1637  TString base("EMCAL_Re_mvsPt_mod") ; base+=mod ;
1638  TString full ;
1639  for(Int_t iConv = 0; iConv<nConv; iConv++){
1640  TLorentzVector * cnv = static_cast<TLorentzVector*>(fConvEvent->At(iConv)) ;
1641  TLorentzVector pi=*cal + *cnv ;
1642 // Double_t alpha=TMath::Abs(cal->Energy()-cnv->Energy())/(cal->Energy()+cnv->Energy()) ;
1643  full=base+"_single" ;
1644  if(cnv->TestBit(kConvOnFly) && !cnv->TestBit(kConvEta)){
1645  FillHistogram(full,pi.M(),cal->Pt()) ;
1646  full=base+"_all" ;
1647  FillHistogram(full,pi.M(),pi.Pt()) ;
1648  char keym[55] ;
1649  for(Int_t iw=0;iw<10;iw++){ //resolution
1650  for(Int_t in=0;in<10;in++){
1651  snprintf(keym,55,"hMC_nonlinearity_EMCAL_w%d_n%d",iw,in) ;
1652  Double_t mMod=0.,ptMod=0. ;
1653  RecalibrateEMCAL(mMod, ptMod, cal, cnv, iw, in) ;
1654  FillHistogram(keym,mMod,ptMod) ;
1655  }
1656  }
1657  }
1658  //Vary Conversion cuts
1659  if(cnv->TestBit(kConvOnFly)){
1660  //Default
1661  if(!cnv->TestBit(kConvEta)){
1662  FillHistogram("EMCAL_Re_mvsPt_OnFly",pi.M(),pi.Pt()) ;
1663  FillHistogram("EMCAL_Re_mvsPt_OnFly_mult",pi.M(),pi.Pt(),fCentr) ;
1664  }
1665  if(cnv->TestBit(kConvdEdx) && !cnv->TestBit(kConvEta))
1666  FillHistogram("EMCAL_Re_mvsPt_On_dEdx",pi.M(),pi.Pt()) ;
1667  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvR))
1668  FillHistogram("EMCAL_Re_mvsPt_On_R120",pi.M(),pi.Pt()) ;
1669  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvZR))
1670  FillHistogram("EMCAL_Re_mvsPt_On_Z",pi.M(),pi.Pt()) ;
1671  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvNDF))
1672  FillHistogram("EMCAL_Re_mvsPt_On_chi",pi.M(),pi.Pt()) ;
1673  if(1)
1674  FillHistogram("EMCAL_Re_mvsPt_On_Eta",pi.M(),pi.Pt()) ;
1675  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvPlan))
1676  FillHistogram("EMCAL_Re_mvsPt_On_Wcut",pi.M(),pi.Pt()) ;
1677  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvArmQt))
1678  FillHistogram("EMCAL_Re_mvsPt_On_ArmQt",pi.M(),pi.Pt()) ;
1679  }
1680  else{
1681  //Default
1682  if(!cnv->TestBit(kConvEta)){
1683  FillHistogram("EMCAL_Re_mvsPt_Offline",pi.M(),pi.Pt()) ;
1684  FillHistogram("EMCAL_Re_mvsPt_Offline_mult",pi.M(),pi.Pt(),fCentr) ;
1685  }
1686  if(cnv->TestBit(kConvdEdx) && !cnv->TestBit(kConvEta))
1687  FillHistogram("EMCAL_Re_mvsPt_Off_dEdx",pi.M(),pi.Pt()) ;
1688  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvR))
1689  FillHistogram("EMCAL_Re_mvsPt_Off_R120",pi.M(),pi.Pt()) ;
1690  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvZR))
1691  FillHistogram("EMCAL_Re_mvsPt_Off_Z",pi.M(),pi.Pt()) ;
1692  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvNDF))
1693  FillHistogram("EMCAL_Re_mvsPt_Off_chi",pi.M(),pi.Pt()) ;
1694  if(1)
1695  FillHistogram("EMCAL_Re_mvsPt_Off_Eta",pi.M(),pi.Pt()) ;
1696  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvPlan))
1697  FillHistogram("EMCAL_Re_mvsPt_Off_Wcut",pi.M(),pi.Pt()) ;
1698  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvArmQt))
1699  FillHistogram("EMCAL_Re_mvsPt_Off_ArmQt",pi.M(),pi.Pt()) ;
1700  }
1701  }
1702  }
1703  //Now fill mixed
1704  TList * prevPHOS = fPHOSEvents[zvtx] ;
1705  TList * prevEMCAL = fEMCALEvents[zvtx] ;
1706  TList * prevConv = fConvEvents[zvtx] ;
1707 
1708  //PHOS
1709  for(Int_t iPHOS=0; iPHOS<nPHOS;iPHOS++){
1710  TLorentzVector * cal = static_cast<TLorentzVector*>(fPHOSEvent->At(iPHOS)) ;
1711  for(Int_t ev=0; ev<prevConv->GetSize();ev++){
1712  TClonesArray * mixConv = static_cast<TClonesArray*>(prevConv->At(ev)) ;
1713  for(Int_t iConv = 0; iConv<mixConv->GetEntriesFast(); iConv++){
1714  TLorentzVector * cnv = static_cast<TLorentzVector*>(mixConv->At(iConv)) ;
1715  TLorentzVector pi=*cal + *cnv ;
1716  if(cnv->TestBit(kConvOnFly) && !cnv->TestBit(kConvEta)){
1717  FillHistogram("PHOS_Mi_mvsPt_all",pi.M(),pi.Pt()) ;
1718  if(cal->TestBit(kCaloPIDdisp))
1719  FillHistogram("PHOS_Mi_mvsPt_Disp",pi.M(),pi.Pt()) ;
1720  if(cal->TestBit(kCaloPIDtof))
1721  FillHistogram("PHOS_Mi_mvsPt_TOF",pi.M(),pi.Pt()) ;
1722  if(cal->TestBit(kCaloPIDneutral))
1723  FillHistogram("PHOS_Mi_mvsPt_Neutral",pi.M(),pi.Pt()) ;
1724  if(cal->TestBit(kCaloPIDneutral) && cal->TestBit(kCaloPIDdisp))
1725  FillHistogram("PHOS_Mi_mvsPt_DispNeutral",pi.M(),pi.Pt()) ;
1726  }
1727  //Vary Conversion cuts
1728  if(cnv->TestBit(kConvOnFly)){
1729  //Default
1730  if(!cnv->TestBit(kConvEta))
1731  FillHistogram("PHOS_Mi_mvsPt_OnFly",pi.M(),pi.Pt()) ;
1732  if( cnv->TestBit(kConvdEdx) && !cnv->TestBit(kConvEta))
1733  FillHistogram("PHOS_Mi_mvsPt_On_dEdx",pi.M(),pi.Pt()) ;
1734  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvR))
1735  FillHistogram("PHOS_Mi_mvsPt_On_R120",pi.M(),pi.Pt()) ;
1736  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvZR))
1737  FillHistogram("PHOS_Mi_mvsPt_On_Z",pi.M(),pi.Pt()) ;
1738  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvNDF))
1739  FillHistogram("PHOS_Mi_mvsPt_On_chi",pi.M(),pi.Pt()) ;
1740  if(1)
1741  FillHistogram("PHOS_Mi_mvsPt_On_Eta",pi.M(),pi.Pt()) ;
1742  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvPlan))
1743  FillHistogram("PHOS_Mi_mvsPt_On_Wcut",pi.M(),pi.Pt()) ;
1744  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvArmQt))
1745  FillHistogram("PHOS_Mi_mvsPt_On_ArmQt",pi.M(),pi.Pt()) ;
1746  }
1747  else{
1748  //Default
1749  if(!cnv->TestBit(kConvEta))
1750  FillHistogram("PHOS_Mi_mvsPt_Offline",pi.M(),pi.Pt()) ;
1751  if(cnv->TestBit(kConvdEdx) && !cnv->TestBit(kConvEta))
1752  FillHistogram("PHOS_Mi_mvsPt_Off_dEdx",pi.M(),pi.Pt()) ;
1753  if( !cnv->TestBit(kConvEta) && cnv->TestBit(kConvR))
1754  FillHistogram("PHOS_Mi_mvsPt_Off_R120",pi.M(),pi.Pt()) ;
1755  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvZR))
1756  FillHistogram("PHOS_Mi_mvsPt_Off_Z",pi.M(),pi.Pt()) ;
1757  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvNDF))
1758  FillHistogram("PHOS_Mi_mvsPt_Off_chi",pi.M(),pi.Pt()) ;
1759  if(1)
1760  FillHistogram("PHOS_Mi_mvsPt_Off_Eta",pi.M(),pi.Pt()) ;
1761  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvPlan))
1762  FillHistogram("PHOS_Mi_mvsPt_Off_Wcut",pi.M(),pi.Pt()) ;
1763  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvArmQt))
1764  FillHistogram("PHOS_Mi_mvsPt_Off_ArmQt",pi.M(),pi.Pt()) ;
1765  }
1766  }
1767  }
1768  }
1769  for(Int_t iConv = 0; iConv<nConv; iConv++){
1770  TLorentzVector * cnv = static_cast<TLorentzVector*>(fConvEvent->At(iConv)) ;
1771  for(Int_t ev=0; ev<prevPHOS->GetSize();ev++){
1772  TClonesArray * mixPHOS = static_cast<TClonesArray*>(prevPHOS->At(ev)) ;
1773  for(Int_t iPHOS=0; iPHOS<mixPHOS->GetEntriesFast();iPHOS++){
1774  TLorentzVector * cal = static_cast<TLorentzVector*>(mixPHOS->At(iPHOS)) ;
1775  TLorentzVector pi=*cal + *cnv ;
1776  if(cnv->TestBit(kConvOnFly) && !cnv->TestBit(kConvEta)){
1777  FillHistogram("PHOS_Mi_mvsPt_all",pi.M(),pi.Pt()) ;
1778  if(cal->TestBit(kCaloPIDdisp))
1779  FillHistogram("PHOS_Mi_mvsPt_Disp",pi.M(),pi.Pt()) ;
1780  if(cal->TestBit(kCaloPIDtof))
1781  FillHistogram("PHOS_Mi_mvsPt_TOF",pi.M(),pi.Pt()) ;
1782  if(cal->TestBit(kCaloPIDneutral))
1783  FillHistogram("PHOS_Mi_mvsPt_Neutral",pi.M(),pi.Pt()) ;
1784  if(cal->TestBit(kCaloPIDneutral) && cal->TestBit(kCaloPIDdisp))
1785  FillHistogram("PHOS_Mi_mvsPt_DispNeutral",pi.M(),pi.Pt()) ;
1786  }
1787  //Vary Conversion cuts
1788  if(cnv->TestBit(kConvOnFly)){
1789  //Default
1790  if(!cnv->TestBit(kConvEta))
1791  FillHistogram("PHOS_Mi_mvsPt_OnFly",pi.M(),pi.Pt()) ;
1792  if(cnv->TestBit(kConvdEdx) && !cnv->TestBit(kConvEta))
1793  FillHistogram("PHOS_Mi_mvsPt_On_dEdx",pi.M(),pi.Pt()) ;
1794  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvR))
1795  FillHistogram("PHOS_Mi_mvsPt_On_R120",pi.M(),pi.Pt()) ;
1796  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvZR))
1797  FillHistogram("PHOS_Mi_mvsPt_On_Z",pi.M(),pi.Pt()) ;
1798  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvNDF))
1799  FillHistogram("PHOS_Mi_mvsPt_On_chi",pi.M(),pi.Pt()) ;
1800  if(1)
1801  FillHistogram("PHOS_Mi_mvsPt_On_Eta",pi.M(),pi.Pt()) ;
1802  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvPlan))
1803  FillHistogram("PHOS_Mi_mvsPt_On_Wcut",pi.M(),pi.Pt()) ;
1804  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvArmQt))
1805  FillHistogram("PHOS_Mi_mvsPt_On_ArmQt",pi.M(),pi.Pt()) ;
1806  }
1807  else{
1808  //Default
1809  if(!cnv->TestBit(kConvEta))
1810  FillHistogram("PHOS_Mi_mvsPt_Offline",pi.M(),pi.Pt()) ;
1811  if(cnv->TestBit(kConvdEdx) && !cnv->TestBit(kConvEta))
1812  FillHistogram("PHOS_Mi_mvsPt_Off_dEdx",pi.M(),pi.Pt()) ;
1813  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvR))
1814  FillHistogram("PHOS_Mi_mvsPt_Off_R120",pi.M(),pi.Pt()) ;
1815  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvZR))
1816  FillHistogram("PHOS_Mi_mvsPt_Off_Z",pi.M(),pi.Pt()) ;
1817  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvNDF))
1818  FillHistogram("PHOS_Mi_mvsPt_Off_chi",pi.M(),pi.Pt()) ;
1819  if(1)
1820  FillHistogram("PHOS_Mi_mvsPt_Off_Eta",pi.M(),pi.Pt()) ;
1821  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvPlan))
1822  FillHistogram("PHOS_Mi_mvsPt_Off_Wcut",pi.M(),pi.Pt()) ;
1823  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvArmQt))
1824  FillHistogram("PHOS_Mi_mvsPt_Off_ArmQt",pi.M(),pi.Pt()) ;
1825  }
1826  }
1827  }
1828  }
1829 
1830 /*
1831  //PHOS module dependent
1832  for(Int_t iPHOS=0; iPHOS<nPHOS;iPHOS++){
1833  TLorentzVector * cal = static_cast<TLorentzVector*>(fPHOSEvent->At(iPHOS)) ;
1834  Int_t mod=1;
1835  while(!cal->TestBit(BIT(17+mod)) && mod<5)mod++ ;
1836  TString base("PHOS_Mi_mvsPt_mod") ; base+=mod ;
1837  TString full ;
1838  for(Int_t ev=0; ev<prevConv->GetSize();ev++){
1839  TClonesArray * mixConv = static_cast<TClonesArray*>(prevConv->At(ev)) ;
1840  for(Int_t iConv = 0; iConv<mixConv->GetEntriesFast(); iConv++){
1841  TLorentzVector * cnv = static_cast<TLorentzVector*>(mixConv->At(iConv)) ;
1842  TLorentzVector pi=*cal + *cnv ;
1843  full=base+"_all" ;
1844  if(cnv->TestBit(kConvOnFly) && !cnv->TestBit(kConvEta)){
1845  FillHistogram(full,pi.M(),pi.Pt()) ;
1846  }
1847  }
1848  }
1849  }
1850  for(Int_t iConv = 0; iConv<nConv; iConv++){
1851  TLorentzVector * cnv = static_cast<TLorentzVector*>(fConvEvent->At(iConv)) ;
1852  for(Int_t ev=0; ev<prevPHOS->GetSize();ev++){
1853  TClonesArray * mixPHOS = static_cast<TClonesArray*>(prevPHOS->At(ev)) ;
1854  for(Int_t iPHOS=0; iPHOS<mixPHOS->GetEntriesFast();iPHOS++){
1855  TLorentzVector * cal = static_cast<TLorentzVector*>(mixPHOS->At(iPHOS)) ;
1856  Int_t mod=1;
1857  while(!cal->TestBit(BIT(17+mod)) && mod<5)mod++ ;
1858  TString base("PHOS_Mi_mvsPt_mod") ; base+=mod ;
1859  TString full ;
1860  TLorentzVector pi=*cal + *cnv ;
1861  full=base+"_all" ;
1862  if(cnv->TestBit(kConvOnFly) && !cnv->TestBit(kConvEta)){
1863  FillHistogram(full,pi.M(),pi.Pt()) ;
1864  }
1865  }
1866  }
1867  }
1868 */
1869 
1870  //EMCAL
1871  for(Int_t iEMCAL=0; iEMCAL<nEMCAL;iEMCAL++){
1872  TLorentzVector * cal = static_cast<TLorentzVector*>(fEMCALEvent->At(iEMCAL)) ;
1873  Int_t mod=0;
1874  while(!cal->TestBit(BIT(17+mod)) && mod<6)mod++ ;
1875  TString base("EMCAL_Mi_mvsPt_mod") ; base+=mod ;
1876  TString full ;
1877  for(Int_t ev=0; ev<prevConv->GetSize();ev++){
1878  TClonesArray * mixConv = static_cast<TClonesArray*>(prevConv->At(ev)) ;
1879  for(Int_t iConv = 0; iConv<mixConv->GetEntriesFast(); iConv++){
1880  TLorentzVector * cnv = static_cast<TLorentzVector*>(mixConv->At(iConv)) ;
1881  TLorentzVector pi=*cal + *cnv ;
1882  full=base+"_all" ;
1883  if(cnv->TestBit(kConvOnFly) && !cnv->TestBit(kConvEta)){
1884  FillHistogram(full,pi.M(),pi.Pt()) ;
1885  }
1886  if(cnv->TestBit(kConvOnFly)){
1887  //Default
1888  if(!cnv->TestBit(kConvEta))
1889  FillHistogram("EMCAL_Mi_mvsPt_OnFly",pi.M(),pi.Pt()) ;
1890  if(cnv->TestBit(kConvdEdx) && !cnv->TestBit(kConvEta))
1891  FillHistogram("EMCAL_Mi_mvsPt_On_dEdx",pi.M(),pi.Pt()) ;
1892  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvR))
1893  FillHistogram("EMCAL_Mi_mvsPt_On_R120",pi.M(),pi.Pt()) ;
1894  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvZR))
1895  FillHistogram("EMCAL_Mi_mvsPt_On_Z",pi.M(),pi.Pt()) ;
1896  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvNDF))
1897  FillHistogram("EMCAL_Mi_mvsPt_On_chi",pi.M(),pi.Pt()) ;
1898  if(1)
1899  FillHistogram("EMCAL_Mi_mvsPt_On_Eta",pi.M(),pi.Pt()) ;
1900  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvPlan))
1901  FillHistogram("EMCAL_Mi_mvsPt_On_Wcut",pi.M(),pi.Pt()) ;
1902  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvArmQt))
1903  FillHistogram("EMCAL_Mi_mvsPt_On_ArmQt",pi.M(),pi.Pt()) ;
1904  }
1905  else{
1906  //Default
1907  if(!cnv->TestBit(kConvEta))
1908  FillHistogram("EMCAL_Mi_mvsPt_Offline",pi.M(),pi.Pt()) ;
1909  if( cnv->TestBit(kConvdEdx) && !cnv->TestBit(kConvEta))
1910  FillHistogram("EMCAL_Mi_mvsPt_Off_dEdx",pi.M(),pi.Pt()) ;
1911  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvR))
1912  FillHistogram("EMCAL_Mi_mvsPt_Off_R120",pi.M(),pi.Pt()) ;
1913  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvZR))
1914  FillHistogram("EMCAL_Mi_mvsPt_Off_Z",pi.M(),pi.Pt()) ;
1915  if( !cnv->TestBit(kConvEta) && cnv->TestBit(kConvNDF))
1916  FillHistogram("EMCAL_Mi_mvsPt_Off_chi",pi.M(),pi.Pt()) ;
1917  if(1)
1918  FillHistogram("EMCAL_Mi_mvsPt_Off_Eta",pi.M(),pi.Pt()) ;
1919  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvPlan))
1920  FillHistogram("EMCAL_Mi_mvsPt_Off_Wcut",pi.M(),pi.Pt()) ;
1921  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvArmQt))
1922  FillHistogram("EMCAL_Mi_mvsPt_Off_ArmQt",pi.M(),pi.Pt()) ;
1923  }
1924  }
1925  }
1926  }
1927 /*
1928  for(Int_t iConv = 0; iConv<nConv; iConv++){
1929  TLorentzVector * cnv = static_cast<TLorentzVector*>(fConvEvent->At(iConv)) ;
1930  for(Int_t ev=0; ev<prevEMCAL->GetSize();ev++){
1931  TClonesArray * mixEMCAL = static_cast<TClonesArray*>(prevEMCAL->At(ev)) ;
1932  for(Int_t iEMCAL=0; iEMCAL<mixEMCAL->GetEntriesFast();iEMCAL++){
1933  TLorentzVector * cal = static_cast<TLorentzVector*>(mixEMCAL->At(iEMCAL)) ;
1934  Int_t mod=0;
1935  while(!cal->TestBit(BIT(17+mod)) && mod<6)mod++ ;
1936  TString base("EMCAL_Mi_mvsPt_mod") ; base+=mod ;
1937  TString full ;
1938  TLorentzVector pi=*cal + *cnv ;
1939  full=base+"_all" ;
1940  if(cnv->TestBit(kConvOnFly) && !cnv->TestBit(kConvEta)){
1941  FillHistogram(full,pi.M(),pi.Pt()) ;
1942  if(cal->TestBit(kCaloPIDdisp)){
1943  full=base+"_Disp" ;
1944  FillHistogram(full,pi.M(),pi.Pt()) ;
1945  }
1946  if(cal->TestBit(kCaloPIDtof)){
1947  full=base+"_TOF" ;
1948  FillHistogram(full,pi.M(),pi.Pt()) ;
1949  }
1950  if(cal->TestBit(kCaloPIDneutral)){
1951  full=base+"_Neutral" ;
1952  FillHistogram(full,pi.M(),pi.Pt()) ;
1953  }
1954  if(cal->TestBit(kCaloPIDneutral) && cal->TestBit(kCaloPIDdisp)){
1955  full=base+"_DispNeutral" ;
1956  FillHistogram(full,pi.M(),pi.Pt()) ;
1957  }
1958  }
1959  if(cnv->TestBit(kConvOnFly)){
1960  //Default
1961  if(!cnv->TestBit(kConvEta))
1962  FillHistogram("EMCAL_Mi_mvsPt_OnFly",pi.M(),pi.Pt()) ;
1963  if(cnv->TestBit(kConvdEdx) && !cnv->TestBit(kConvEta))
1964  FillHistogram("EMCAL_Mi_mvsPt_On_dEdx",pi.M(),pi.Pt()) ;
1965  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvR))
1966  FillHistogram("EMCAL_Mi_mvsPt_On_R120",pi.M(),pi.Pt()) ;
1967  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvZR))
1968  FillHistogram("EMCAL_Mi_mvsPt_On_Z",pi.M(),pi.Pt()) ;
1969  if( !cnv->TestBit(kConvEta) && cnv->TestBit(kConvNDF))
1970  FillHistogram("EMCAL_Mi_mvsPt_On_chi",pi.M(),pi.Pt()) ;
1971  if(1)
1972  FillHistogram("EMCAL_Mi_mvsPt_On_Eta",pi.M(),pi.Pt()) ;
1973  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvPlan))
1974  FillHistogram("EMCAL_Mi_mvsPt_On_Wcut",pi.M(),pi.Pt()) ;
1975  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvArmQt))
1976  FillHistogram("EMCAL_Mi_mvsPt_On_ArmQt",pi.M(),pi.Pt()) ;
1977  }
1978  else{
1979  //Default
1980  if(!cnv->TestBit(kConvEta))
1981  FillHistogram("EMCAL_Mi_mvsPt_Offline",pi.M(),pi.Pt()) ;
1982  if(cnv->TestBit(kConvdEdx) && !cnv->TestBit(kConvEta))
1983  FillHistogram("EMCAL_Mi_mvsPt_Off_dEdx",pi.M(),pi.Pt()) ;
1984  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvR))
1985  FillHistogram("EMCAL_Mi_mvsPt_Off_R120",pi.M(),pi.Pt()) ;
1986  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvZR))
1987  FillHistogram("EMCAL_Mi_mvsPt_Off_Z",pi.M(),pi.Pt()) ;
1988  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvNDF))
1989  FillHistogram("EMCAL_Mi_mvsPt_Off_chi",pi.M(),pi.Pt()) ;
1990  if(1)
1991  FillHistogram("EMCAL_Mi_mvsPt_Off_Eta",pi.M(),pi.Pt()) ;
1992  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvPlan))
1993  FillHistogram("EMCAL_Mi_mvsPt_Off_Wcut",pi.M(),pi.Pt()) ;
1994  if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvArmQt))
1995  FillHistogram("EMCAL_Mi_mvsPt_Off_ArmQt",pi.M(),pi.Pt()) ;
1996  }
1997  }
1998  }
1999  }
2000 */
2001 
2002 
2003 
2004  //Now we either add current events to stack or remove
2005  //Here we have some difficulty: conversion and calorimeter photons have different spectra.
2006  //So to correctly reproduce combinatorial background we have to preserve average number of
2007  //photons of both kinds per event. Therefore we should not reject empty PHOS/EMCAL events
2008  //though it will cost some memory. Reject only those events where no photons anywhere
2009 
2010  if((fPHOSEvent->GetEntriesFast()>0 || fEMCALEvent->GetEntriesFast()>0) && fConvEvent->GetEntriesFast()>0){
2011  prevPHOS->AddFirst(fPHOSEvent) ;
2012  fPHOSEvent=0;
2013  prevEMCAL->AddFirst(fEMCALEvent) ;
2014  fEMCALEvent=0 ;
2015  prevConv->AddFirst(fConvEvent) ;
2016  fConvEvent=0 ;
2017  if(prevPHOS->GetSize()>100){//Remove redundant events
2018  TClonesArray * tmp = static_cast<TClonesArray*>(prevPHOS->Last()) ;
2019  prevPHOS->RemoveLast() ;
2020  delete tmp ;
2021  tmp = static_cast<TClonesArray*>(prevEMCAL->Last()) ;
2022  prevEMCAL->RemoveLast() ;
2023  delete tmp ;
2024  tmp = static_cast<TClonesArray*>(prevConv->Last()) ;
2025  prevConv->RemoveLast() ;
2026  delete tmp ;
2027  }
2028  }
2029 
2030 }
2031 //___________________________________________________________________________
2032 void AliAnalysisTaskCaloConv::ProcessMC(){
2033  //ProcessMC
2034  //fill histograms for efficiensy etc. calculation
2035  if(!fMCEvent) return ;
2036 
2037  const Double_t rcut = 1. ; //cut for primary particles
2038  Double_t vtx[3];
2039  vtx[0] = fESDEvent->GetPrimaryVertex()->GetX();
2040  vtx[1] = fESDEvent->GetPrimaryVertex()->GetY();
2041  vtx[2] = fESDEvent->GetPrimaryVertex()->GetZ();
2042 
2043  Int_t nPHOS=fPHOSEvent->GetEntriesFast() ;
2044  Int_t nEMCAL=fEMCALEvent->GetEntriesFast() ;
2045  Int_t nConv = fConvEvent->GetEntriesFast() ;
2046 
2047  //---------First pi0/eta-----------------------------
2048  char partName[10] ;
2049  char hkey[55] ;
2050  for (Int_t iTracks = 0; iTracks < fMCEvent->GetNumberOfTracks(); iTracks++) {
2051  TParticle* particle = (TParticle *)fMCEvent->Particle(iTracks);
2052  if(particle->GetPdgCode() == 111)
2053  snprintf(partName,10,"pi0") ;
2054  else
2055  if(particle->GetPdgCode() == 221)
2056  snprintf(partName,10,"eta") ;
2057  else
2058  continue ;
2059 
2060  //Primary particle
2061  if(particle->R() >rcut)
2062  continue ;
2063 
2064  Double_t pt = particle->Pt() ;
2065  //Total number of pi0 with creation radius <1 cm
2066  snprintf(hkey,55,"hMC_CaloConv_all%s",partName) ;
2067  FillHistogram(hkey,pt) ;
2068  if(TMath::Abs(particle->Y())<1.){
2069  snprintf(hkey,55,"hMC_CaloConv_%s_unitEta",partName) ;
2070  FillHistogram(hkey,pt,fCentr) ;
2071  }
2072 
2073  //Check if one of photons converted
2074  if(particle->GetNDaughters()!=2)
2075  continue ; //Do not account Dalitz decays
2076 
2077  TParticle * gamma1 = fMCEvent->Particle(particle->GetFirstDaughter());
2078  TParticle * gamma2 = fMCEvent->Particle(particle->GetLastDaughter());
2079  //Number of pi0s decayed into acceptance
2080  Bool_t inAcc1 = (TMath::Abs(gamma1->Eta())<0.9) ;
2081  Bool_t inAcc2 = (TMath::Abs(gamma2->Eta())<0.9) ;
2082  Int_t mod1,mod2 ;
2083  Double_t x=0.,z=0. ;
2084  Bool_t hitPHOS1 = fPHOSgeom->ImpactOnEmc(gamma1, mod1, z,x) ;
2085  Bool_t hitPHOS2 = fPHOSgeom->ImpactOnEmc(gamma2, mod2, z,x) ;
2086  Bool_t hitEMCAL1= fEMCALgeom->Impact(gamma1) ;
2087  Bool_t hitEMCAL2= fEMCALgeom->Impact(gamma2) ;
2088 
2089  Bool_t goodPair=kFALSE ;
2090  if((inAcc1 && hitPHOS2) || (inAcc2 && hitPHOS1)){
2091  snprintf(hkey,55,"hMC_CaloConv_%sPHOSacc",partName) ;
2092  FillHistogram(hkey,pt) ;
2093  goodPair=kTRUE ;
2094  }
2095  if((inAcc1 && hitEMCAL2) || (inAcc2 && hitEMCAL1)){
2096  snprintf(hkey,55,"hMC_CaloConv_%sEMCALacc",partName) ;
2097  FillHistogram(hkey,pt) ;
2098  goodPair=kTRUE ;
2099  }
2100  if(!goodPair){
2101  continue ;
2102  }
2103 
2104  Bool_t converted1 = kFALSE ;
2105  if(gamma1->GetNDaughters()==2){
2106  TParticle * e1=fMCEvent->Particle(gamma1->GetFirstDaughter()) ;
2107  TParticle * e2=fMCEvent->Particle(gamma1->GetLastDaughter()) ;
2108  if(TMath::Abs(e1->GetPdgCode())==11 && TMath::Abs(e2->GetPdgCode())==11){ //conversion
2109  if(e1->R()<180.)
2110  converted1 = kTRUE ;
2111  }
2112  }
2113  Bool_t converted2 = kFALSE ;
2114  if(gamma2->GetNDaughters()==2){
2115  TParticle * e1=fMCEvent->Particle(gamma2->GetFirstDaughter()) ;
2116  TParticle * e2=fMCEvent->Particle(gamma2->GetLastDaughter()) ;
2117  if(TMath::Abs(e1->GetPdgCode())==11 && TMath::Abs(e2->GetPdgCode())==11){ //conversion
2118  if(e1->R()<180.)
2119  converted2 = kTRUE ;
2120  }
2121  }
2122 
2123  //Number of pi0s with one photon converted
2124  if((converted1 && !converted2 && hitPHOS2) || (!converted1 && hitPHOS1 && converted2)) {
2125  snprintf(hkey,55,"hMC_CaloConv_%s_PHOS_conv",partName) ;
2126  FillHistogram(hkey,pt) ;
2127  }
2128 
2129  if((converted1 && !converted2 && hitEMCAL2) || (!converted1 && hitEMCAL1 && converted2)) {
2130  snprintf(hkey,55,"hMC_CaloConv_%s_EMCAL_conv",partName) ;
2131  FillHistogram(hkey,pt) ;
2132  }
2133 
2134  //Both converted
2135  if(converted1 && converted2) {
2136  snprintf(hkey,55,"hMC_CaloConv_%s_bothphot_conv",partName) ;
2137  FillHistogram(hkey,pt) ;
2138  continue ;
2139  }
2140 
2141  //photon pointing calorimeter converted
2142  if((converted1 && hitPHOS1 && !hitEMCAL2) || (converted2 && hitPHOS2 && !hitEMCAL1) ||
2143  (converted1 && hitEMCAL1 && !hitPHOS2) || (converted2 && hitEMCAL2 && !hitPHOS1)){
2144  snprintf(hkey,55,"hMC_CaloConv_%s_convPhotInCalo",partName) ;
2145  FillHistogram(hkey,pt) ;
2146  continue ;
2147  }
2148 
2149  //Converted pi0 with v0 and photon PHOS or EMCAL
2150  Bool_t foundV01onfly=kFALSE, foundV01offline=kFALSE, foundV02onfly=kFALSE, foundV02offline=kFALSE ;
2151  Bool_t foundV01onflyPID=kFALSE, foundV01offlinePID=kFALSE, foundV02onflyPID=kFALSE, foundV02offlinePID=kFALSE ;
2152  TLorentzVector pConvOn,pConvOff ;
2153  for(Int_t iv0=0; iv0<fESDEvent->GetNumberOfV0s();iv0++){
2154  AliESDv0 * v0 = fESDEvent->GetV0(iv0) ;
2155 
2156  TParticle * negativeMC = fMCEvent->Particle(TMath::Abs(fESDEvent->GetTrack(v0->GetNindex())->GetLabel()));
2157  TParticle * positiveMC = fMCEvent->Particle(TMath::Abs(fESDEvent->GetTrack(v0->GetPindex())->GetLabel()));
2158 
2159  if(negativeMC->GetMother(0) != positiveMC->GetMother(0))
2160  continue ;
2161 
2162  if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){
2163  continue;
2164  }
2165  if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){
2166  continue;
2167  }
2168 
2169  TParticle * v0Gamma = fMCEvent->Particle(negativeMC->GetMother(0));
2170  Bool_t same = (v0Gamma == gamma1) ;
2171  TParticle * tmp = v0Gamma ;
2172  while(!same && tmp->GetFirstMother()>=0){
2173  tmp = fMCEvent->Particle(tmp->GetFirstMother());
2174  same = (tmp == gamma1) ;
2175  }
2176  if(same){
2177  if(v0->GetOnFlyStatus())
2178  foundV01onfly = kTRUE ;
2179  else
2180  foundV01offline= kTRUE ;
2181  for(Int_t iconv=0; iconv<nConv;iconv++){
2182  if(fGammaV0s[iconv] == iv0){
2183  TLorentzVector * cnv = static_cast<TLorentzVector*>(fConvEvent->At(iconv)) ;
2184  //default cuts
2185  if(!cnv->TestBit(kConvEta)){
2186  if(v0->GetOnFlyStatus()){
2187  pConvOn= *cnv ;
2188  foundV01onflyPID = kTRUE ;
2189  }
2190  else{
2191  pConvOff= *cnv ;
2192  foundV01offlinePID = kTRUE ;
2193  }
2194  }
2195  break ;
2196  }
2197  }
2198  continue ;
2199  }
2200  same = (v0Gamma == gamma2) ;
2201  tmp = v0Gamma ;
2202  while(!same && tmp->GetFirstMother()>=0){
2203  tmp = fMCEvent->Particle(tmp->GetFirstMother());
2204  same = (tmp == gamma2) ;
2205  }
2206  if(same){
2207  if(v0->GetOnFlyStatus())
2208  foundV02onfly = kTRUE ;
2209  else
2210  foundV02offline = kTRUE ;
2211  for(Int_t iconv=0; iconv<nConv;iconv++){
2212  if(fGammaV0s[iconv] == iv0){
2213  TLorentzVector * cnv = static_cast<TLorentzVector*>(fConvEvent->At(iconv)) ;
2214  //default cuts
2215  if(!cnv->TestBit(kConvEta)){
2216  if(v0->GetOnFlyStatus()){
2217  pConvOn= *cnv ;
2218  foundV02onflyPID = kTRUE ;
2219  }
2220  else{
2221  pConvOff= *cnv ;
2222  foundV02offlinePID = kTRUE ;
2223  }
2224  }
2225  break ;
2226  }
2227  }
2228  }
2229  }
2230 
2231  goodPair=kFALSE ;
2232  if((foundV01onfly && hitPHOS2) || (foundV02onfly && hitPHOS1)){
2233  snprintf(hkey,55,"hMC_CaloConv_%s_v0onfly_PHOSacc",partName) ;
2234  FillHistogram(hkey,pt) ;
2235  goodPair=kTRUE;
2236  }
2237  if((foundV01offline && hitPHOS2) || (foundV02offline && hitPHOS1)){
2238  snprintf(hkey,55,"hMC_CaloConv_%s_v0offline_PHOSacc",partName) ;
2239  FillHistogram(hkey,pt) ;
2240  goodPair=kTRUE;
2241  }
2242  if((foundV01onfly && hitEMCAL2) || (foundV02onfly && hitEMCAL1)){
2243  snprintf(hkey,55,"hMC_CaloConv_%s_v0onfly_EMCALacc",partName) ;
2244  FillHistogram(hkey,pt) ;
2245  goodPair=kTRUE;
2246  }
2247  if((foundV01offline && hitEMCAL2) || (foundV02offline && hitEMCAL1)){
2248  snprintf(hkey,55,"hMC_CaloConv_%s_v0offline_EMCALacc",partName) ;
2249  FillHistogram(hkey,pt) ;
2250  goodPair=kTRUE;
2251  }
2252  if(!goodPair){
2253  continue ;
2254  }
2255 
2256  //Converted pi0 with v0 and cluster in PHOS/EMCAL
2257  Bool_t cluInPHOS = kFALSE,cluInEMCAL=kFALSE ;
2258  Bool_t cluInPHOSpid = kFALSE,cluInEMCALpid=kFALSE ;
2259  Bool_t closeToBad= kFALSE ;
2260  TLorentzVector pCalo ;
2261  for (Int_t i=0; i<fESDEvent->GetNumberOfCaloClusters(); i++) {
2262  AliESDCaloCluster * clu = fESDEvent->GetCaloCluster(i);
2263  Int_t iprim = clu->GetLabel() ; //# of particle hit PHOS/EMCAL
2264  Bool_t matched = kFALSE ;
2265  while(iprim>=0 ) {
2266  if(iprim==particle->GetFirstDaughter() || iprim==particle->GetLastDaughter()){
2267  matched=kTRUE ;
2268  break ;
2269  }
2270  else{
2271  iprim=fMCEvent->Particle(iprim)->GetFirstMother() ;
2272  }
2273  }
2274  if(!matched)
2275  continue ;
2276  if(clu->IsPHOS() && (hitPHOS1 || hitPHOS2)){
2277  cluInPHOS=kTRUE ;
2278  //Check if cluster passed PID
2279  for(Int_t inPHOS=0; inPHOS<nPHOS;inPHOS++){
2280  if(fGammaPHOS[inPHOS] == i){
2281  cluInPHOSpid=kTRUE ;
2282  break ;
2283  }
2284  }
2285  clu->GetMomentum(pCalo ,vtx);
2286  if(clu->GetDistanceToBadChannel()<fBadDistCutPHOS)
2287  closeToBad=kTRUE ;
2288  break ;
2289  }
2290  if(!clu->IsPHOS() && (hitEMCAL1 || hitEMCAL2)){
2291  cluInEMCAL=kTRUE ;
2292  //Check if cluster passed PID
2293  for(Int_t inEMCAL=0; inEMCAL<nEMCAL;inEMCAL++){
2294  if(fGammaEMCAL[inEMCAL] == i){
2295  cluInPHOSpid=kTRUE ;
2296  break ;
2297  }
2298  }
2299  clu->GetMomentum(pCalo ,vtx);
2300  if(clu->GetDistanceToBadChannel()<fBadDistCutEMCAL)
2301  closeToBad=kTRUE ;
2302  break ;
2303  }
2304  }
2305 
2306  if(cluInPHOS){
2307  //OnFly
2308  if(foundV01onfly ||foundV02onfly){
2309  snprintf(hkey,55,"hMC_CaloConv_%s_v0onfly_PHOSclu",partName) ;
2310  FillHistogram(hkey,pt) ;
2311 
2312  if((foundV01onflyPID ||foundV02onflyPID) && cluInPHOSpid){
2313  snprintf(hkey,55,"hMC_CaloConv_%s_v0onfly_PHOSclu_pid",partName) ;
2314  FillHistogram(hkey,pt) ;
2315  for(Int_t iw=0;iw<10;iw++){ //resolution
2316  for(Int_t in=0;in<10;in++){
2317  char keym[55] ;
2318  snprintf(keym,55,"hMC_CaloConv_%s_v0onfly_PHOSclu_w%d_n%d",partName,iw,in) ;
2319  Double_t mMod=0.,ptMod=0. ;
2320  Recalibrate(mMod, ptMod, &pCalo, &pConvOn, iw, in) ;
2321  FillHistogram(keym,ptMod) ;
2322  snprintf(keym,55,"hMC_CaloConv_%s_v0onfly_ConvPHOSclu_w%d_n%d",partName,iw,in) ;
2323  RecalibrateConvPHOS(mMod, ptMod, &pCalo, &pConvOn, iw, in) ;
2324  FillHistogram(keym,ptMod) ;
2325  }
2326  }
2327  if(!closeToBad){
2328  snprintf(hkey,55,"hMC_CaloConv_%s_v0onfly_PHOSclu_good",partName) ;
2329  FillHistogram(hkey,pt) ;
2330  }
2331  Double_t m=(pCalo+pConvOn).M() ;
2332  Double_t ptm=(pCalo+pConvOn).Pt() ;
2333  snprintf(hkey,55,"hMC_CaloConv_%s_v0on_PHOSclu_ptRec",partName) ;
2334  FillHistogram(hkey,ptm) ;
2335  snprintf(hkey,55,"hMC_CaloConv_%s_v0on_PHOSclu_mvsPt",partName) ;
2336  FillHistogram(hkey,m,ptm) ;
2337  }
2338  }
2339 
2340  //Offline
2341  if(foundV01offline ||foundV02offline){
2342  snprintf(hkey,55,"hMC_CaloConv_%s_v0offline_PHOSclu",partName) ;
2343  FillHistogram(hkey,pt) ;
2344  if((foundV01offlinePID ||foundV02offlinePID) && cluInPHOSpid){
2345  snprintf(hkey,55,"hMC_CaloConv_%s_v0offline_PHOSclu_pid",partName) ;
2346  FillHistogram(hkey,pt) ;
2347  Double_t m=(pCalo+pConvOff).M() ;
2348  Double_t ptm=(pCalo+pConvOff).Pt() ;
2349  snprintf(hkey,55,"hMC_CaloConv_%s_v0off_PHOSclu_ptRec",partName) ;
2350  FillHistogram(hkey,ptm) ;
2351  snprintf(hkey,55,"hMC_CaloConv_%s_v0off_PHOSclu_mvsPt",partName) ;
2352  FillHistogram(hkey,m,ptm) ;
2353  if(!closeToBad){
2354  snprintf(hkey,55,"hMC_CaloConv_%s_v0offline_PHOSclu_good",partName) ;
2355  FillHistogram(hkey,pt) ;
2356  }
2357  }
2358  }
2359 
2360  if((foundV01onflyPID ||foundV02onflyPID) && cluInPHOSpid){
2361  TString base("hMC_CaloConv_") ; base+=partName; base+="_v0onfly_PHOSclu_mod" ;
2362  if(hitPHOS1)
2363  base+=mod1 ;
2364  else
2365  base+=mod2 ;
2366  FillHistogram(base.Data(),pt) ;
2367  }
2368 
2369  if((foundV01offlinePID ||foundV02offlinePID) && cluInPHOSpid){
2370  TString base("hMC_CaloConv_") ; base+=partName; base+="_v0offline_PHOSclu_mod" ;
2371  if(hitPHOS1)
2372  base+=mod1 ;
2373  else
2374  base+=mod2 ;
2375  FillHistogram(base.Data(),pt) ;
2376  }
2377  }
2378  if(cluInEMCAL && strcmp(partName,"pi0")==0){
2379  //OnFly
2380  if(foundV01onfly ||foundV02onfly){
2381  FillHistogram("hMC_CaloConv_pi0_v0onfly_EMCALclu",pt) ;
2382 
2383  if((foundV01onflyPID ||foundV02onflyPID) && cluInEMCALpid){
2384  FillHistogram("hMC_CaloConv_pi0_v0onfly_EMCALclu_pid",pt) ;
2385  for(Int_t iw=0;iw<10;iw++){ //resolution
2386  for(Int_t in=0;in<10;in++){
2387  char keym[55] ;
2388  snprintf(keym,55,"hMC_CaloConv_pi0_v0onfly_EMCALclu_w%d_n%d",iw,in) ;
2389  Double_t mMod=0.,ptMod=0. ;
2390  RecalibrateEMCAL(mMod, ptMod, &pCalo, &pConvOn, iw, in) ;
2391  FillHistogram(keym,ptMod) ;
2392  }
2393  }
2394  if(!closeToBad)
2395  FillHistogram("hMC_CaloConv_pi0_v0onfly_EMCALclu_good",pt) ;
2396  Double_t m=(pCalo+pConvOn).M() ;
2397  Double_t ptm=(pCalo+pConvOn).Pt() ;
2398  FillHistogram("hMC_CaloConv_pi0_v0on_EMCALclu_ptRec",ptm) ;
2399  FillHistogram("hMC_CaloConv_pi0_v0on_EMCALclu_mvsPt",m,ptm) ;
2400  }
2401  }
2402 
2403  //Offline
2404  if(foundV01offline ||foundV02offline){
2405  FillHistogram("hMC_CaloConv_pi0_v0offline_EMCALclu",pt) ;
2406  if((foundV01offlinePID ||foundV02offlinePID) && cluInEMCALpid){
2407  FillHistogram("hMC_CaloConv_pi0_v0offline_EMCALclu_pid",pt) ;
2408  Double_t m=(pCalo+pConvOff).M() ;
2409  Double_t ptm=(pCalo+pConvOff).Pt() ;
2410  FillHistogram("hMC_CaloConv_pi0_v0off_EMCALclu_ptRec",ptm) ;
2411  FillHistogram("hMC_CaloConv_pi0_v0off_EMCALclu_mvsPt",m,ptm) ;
2412  if(!closeToBad)
2413  FillHistogram("hMC_CaloConv_pi0_v0offline_EMCALclu_good",pt) ;
2414  }
2415  }
2416  }
2417  }
2418 
2419  //Construct Inv mass distributions for residual correlations
2420  if(fPHOSEvent && fConvEvent){
2421  for(Int_t iPHOS=0; iPHOS<fPHOSEvent->GetEntriesFast();iPHOS++){
2422  TLorentzVector * cal = static_cast<TLorentzVector*>(fPHOSEvent->At(iPHOS)) ;
2423  Int_t iclu=fGammaPHOS[iPHOS] ;
2424  AliESDCaloCluster * clu = fESDEvent->GetCaloCluster(iclu);
2425  Int_t iprimPHOS = clu->GetLabel() ; //# of particle hit PHOS/EMCAL
2426  for(Int_t iConv = 0; iConv<fConvEvent->GetEntriesFast(); iConv++){
2427  TLorentzVector * cnv = static_cast<TLorentzVector*>(fConvEvent->At(iConv)) ;
2428  if(!cnv->TestBit(kConvOnFly) || cnv->TestBit(kConvEta))
2429  continue;
2430 
2431  Int_t iv0=fGammaV0s[iConv] ;
2432  AliESDv0 * v0 = fESDEvent->GetV0(iv0) ;
2433  Int_t iprimNeg = TMath::Abs(fESDEvent->GetTrack(v0->GetNindex())->GetLabel()) ;
2434  Int_t iprimPos = TMath::Abs(fESDEvent->GetTrack(v0->GetPindex())->GetLabel()) ;
2435 
2436  //Check if there was a common ancistor
2437  Bool_t found = kFALSE ;
2438  Int_t curPHOS=iprimPHOS ;
2439  Int_t commonA=-1 ;
2440  while(!found && curPHOS>-1){
2441  Int_t curNeg=iprimNeg ;
2442  while(!found && curNeg>-1){
2443  if(curNeg==curPHOS){
2444  found=kTRUE ;
2445  commonA=curPHOS ;
2446  }
2447  else{
2448  curNeg=fMCEvent->Particle(curNeg)->GetFirstMother() ;
2449  }
2450  }
2451  curPHOS=fMCEvent->Particle(curPHOS)->GetFirstMother() ;
2452  }
2453  found = kFALSE ;
2454  curPHOS=iprimPHOS ;
2455  Int_t commonB=-1 ;
2456  while(!found && curPHOS>-1){
2457  Int_t curPos=iprimPos ;
2458  while(!found && curPos>-1){
2459  if(curPos==curPHOS){
2460  found=kTRUE ;
2461  commonB=curPHOS ;
2462  }
2463  else{
2464  curPos=fMCEvent->Particle(curPos)->GetFirstMother() ;
2465  }
2466  }
2467  curPHOS=fMCEvent->Particle(curPHOS)->GetFirstMother() ;
2468  }
2469  if(commonA != commonB){
2470  //Strange
2471  AliInfo(Form("CommonA=%d, commonB=%d",commonA,commonB)) ;
2472  }
2473  if(commonA>-1){//There was common particles
2474  Int_t pdg = fMCEvent->Particle(commonA)->GetPdgCode() ;
2475  TLorentzVector pi=*cal + *cnv ;
2476  Double_t m=pi.M() ;
2477  Double_t pt=pi.Pt() ;
2478  Double_t alpha=TMath::Abs(cal->Energy()-cnv->Energy())/(cal->Energy()+cnv->Energy()) ;
2479  switch(pdg){
2480  case 11:
2481  case -11:
2482  case 22: //conversion
2483  FillHistogram("hMC_Resid_PHOS_Phot_mvsPt",m,pt,alpha) ;
2484  break ;
2485  case 111: //pi0
2486  FillHistogram("hMC_Resid_PHOS_Pi0_mvsPt",m,pt,alpha) ;
2487  break ;
2488  case 221: //eta
2489  FillHistogram("hMC_Resid_PHOS_eta_mvsPt",m,pt,alpha) ;
2490  break ;
2491  case 321: //K+
2492  case -321: //K-
2493  case 310: //K0s
2494  case 130: //K0L
2495  FillHistogram("hMC_Resid_PHOS_K_mvsPt",m,pt,alpha) ;
2496  break ;
2497  case 211:
2498  case -211:
2499  FillHistogram("hMC_Resid_PHOS_pi_mvsPt",m,pt,alpha) ;
2500  break ;
2501  case -2212: //pbar
2502  case -2112: //nbar
2503  FillHistogram("hMC_Resid_PHOS_pbar_mvsPt",m,pt,alpha) ;
2504  break ;
2505  default: //else
2506  FillHistogram("hMC_Resid_PHOS_other_mvsPt",m,pt,alpha) ;
2507  break ;
2508  }
2509  }
2510  }
2511  }
2512  }
2513 
2514 
2515  if(fEMCALEvent && fConvEvent){
2516  for(Int_t iEMCAL=0; iEMCAL<fEMCALEvent->GetEntriesFast();iEMCAL++){
2517  TLorentzVector * cal = static_cast<TLorentzVector*>(fEMCALEvent->At(iEMCAL)) ;
2518  Int_t iclu=fGammaEMCAL[iEMCAL] ;
2519  AliESDCaloCluster * clu = fESDEvent->GetCaloCluster(iclu);
2520  Int_t iprimEMCAL = clu->GetLabel() ; //# of particle hit EMCAL
2521  for(Int_t iConv = 0; iConv<fConvEvent->GetEntriesFast(); iConv++){
2522  TLorentzVector * cnv = static_cast<TLorentzVector*>(fConvEvent->At(iConv)) ;
2523  if(!cnv->TestBit(kConvOnFly) || cnv->TestBit(kConvEta))
2524  continue;
2525  Int_t iv0=fGammaV0s[iConv] ;
2526  AliESDv0 * v0 = fESDEvent->GetV0(iv0) ;
2527  Int_t iprimNeg = TMath::Abs(fESDEvent->GetTrack(v0->GetNindex())->GetLabel()) ;
2528  Int_t iprimPos = TMath::Abs(fESDEvent->GetTrack(v0->GetPindex())->GetLabel()) ;
2529 
2530  //Check if there was a common ancistor
2531  Bool_t found = kFALSE ;
2532  Int_t curEMCAL=iprimEMCAL ;
2533  Int_t commonA=-1 ;
2534  while(!found && curEMCAL>-1){
2535  Int_t curNeg=iprimNeg ;
2536  while(!found && curNeg>-1){
2537  if(curNeg==curEMCAL){
2538  found=kTRUE ;
2539  commonA=curEMCAL ;
2540  }
2541  else{
2542  curNeg=fMCEvent->Particle(curNeg)->GetFirstMother() ;
2543  }
2544  }
2545  curEMCAL=fMCEvent->Particle(curEMCAL)->GetFirstMother() ;
2546  }
2547  found = kFALSE ;
2548  curEMCAL=iprimEMCAL ;
2549  Int_t commonB=-1 ;
2550  while(!found && curEMCAL>-1){
2551  Int_t curPos=iprimPos ;
2552  while(!found && curPos>-1){
2553  if(curPos==curEMCAL){
2554  found=kTRUE ;
2555  commonB=curEMCAL ;
2556  }
2557  else{
2558  curPos=fMCEvent->Particle(curPos)->GetFirstMother() ;
2559  }
2560  }
2561  curEMCAL=fMCEvent->Particle(curEMCAL)->GetFirstMother() ;
2562  }
2563 
2564  if(commonA != commonB){
2565  //Strange
2566  AliInfo(Form("CommonA=%d, commonB=%d",commonA,commonB)) ;
2567  }
2568  if(commonA>-1){//There was common particles
2569  Int_t pdg = fMCEvent->Particle(commonA)->GetPdgCode() ;
2570  TLorentzVector pi=*cal + *cnv ;
2571  Double_t m=pi.M() ;
2572  Double_t pt=pi.Pt() ;
2573  Double_t alpha=TMath::Abs(cal->Energy()-cnv->Energy())/(cal->Energy()+cnv->Energy()) ;
2574  switch(pdg){
2575  case 11:
2576  case -11:
2577  case 22: //conversion
2578  FillHistogram("hMC_Resid_EMCAL_Phot_mvsPt",m,pt,alpha) ;
2579  break ;
2580  case 111: //pi0
2581  FillHistogram("hMC_Resid_EMCAL_Pi0_mvsPt",m,pt,alpha) ;
2582  break ;
2583  case 221: //eta
2584  FillHistogram("hMC_Resid_EMCAL_eta_mvsPt",m,pt,alpha) ;
2585  break ;
2586  case 321: //K+
2587  case -321: //K-
2588  case 310: //K0s
2589  case 130: //K0L
2590  FillHistogram("hMC_Resid_EMCAL_K_mvsPt",m,pt,alpha) ;
2591  break ;
2592  case 211:
2593  case -211:
2594  FillHistogram("hMC_Resid_EMCAL_pi_mvsPt",m,pt,alpha) ;
2595  break ;
2596  case -2212: //pbar
2597  case -2112: //nbar
2598  FillHistogram("hMC_Resid_EMCAL_pbar_mvsPt",m,pt,alpha) ;
2599  break ;
2600  default: //else
2601  FillHistogram("hMC_Resid_EMCAL_other_mvsPt",m,pt,alpha) ;
2602  break ;
2603  }
2604  }
2605  }
2606  }
2607  }
2608 
2609 
2610  //------------- now photons ----------------
2611  for (Int_t iTracks = 0; iTracks < fMCEvent->GetNumberOfTracks(); iTracks++) {
2612  TParticle* particle = (TParticle *)fMCEvent->Particle(iTracks);
2613  if(particle->GetPdgCode() != 22)
2614  continue ;
2615 
2616  if(particle->R() >rcut)
2617  continue ;
2618 
2619  if(TMath::Abs(particle->Eta())>0.9)
2620  continue ;
2621 
2622  Double_t pt = particle->Pt() ;
2623  //Total number of pi0 with creation radius <1 cm
2624  FillHistogram("hMC_CaloConv_phot",pt) ;
2625 
2626  Int_t mod ;
2627  Double_t x=0.,z=0. ;
2628  Bool_t hitPHOS = fPHOSgeom->ImpactOnEmc(particle, mod, z,x) ;
2629  Bool_t hitEMCAL= fEMCALgeom->Impact(particle) ;
2630 
2631  //Photons in PHOS/EMCAL acceptance
2632  if(hitPHOS)
2633  FillHistogram("hMC_CaloConv_gammaPHOSacc",pt) ;
2634  if(hitEMCAL)
2635  FillHistogram("hMC_CaloConv_gammaEMCALacc",pt) ;
2636 
2637  //number of photons converted
2638  Bool_t converted = kFALSE ;
2639  if(particle->GetNDaughters()==2){
2640  TParticle * e1=fMCEvent->Particle(particle->GetFirstDaughter()) ;
2641  TParticle * e2=fMCEvent->Particle(particle->GetLastDaughter()) ;
2642  if(TMath::Abs(e1->GetPdgCode())==11 && TMath::Abs(e2->GetPdgCode())==11){ //conversion
2643  if(e1->R()<180.)
2644  converted = kTRUE ;
2645  }
2646  }
2647  if(converted)
2648  FillHistogram("hMC_CaloConv_gamma_conv",pt) ;
2649 
2650  //Converted photons with V0
2651  TLorentzVector pConv ;
2652  Bool_t foundV0=kFALSE ;
2653  for(Int_t iv0=0; iv0<fESDEvent->GetNumberOfV0s();iv0++){
2654  AliESDv0 * v0 = fESDEvent->GetV0(iv0) ;
2655 
2656  TParticle * negativeMC = fMCEvent->Particle(TMath::Abs(fESDEvent->GetTrack(v0->GetNindex())->GetLabel()));
2657  TParticle * positiveMC = fMCEvent->Particle(TMath::Abs(fESDEvent->GetTrack(v0->GetPindex())->GetLabel()));
2658 
2659  if(negativeMC->GetMother(0) != positiveMC->GetMother(0))
2660  continue ;
2661 
2662  if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){
2663  continue;
2664  }
2665  if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){
2666  continue;
2667  }
2668 
2669  TParticle * v0Gamma = fMCEvent->Particle(negativeMC->GetMother(0));
2670  Bool_t same = (v0Gamma == particle) ;
2671  TParticle * tmp = v0Gamma ;
2672  while(!same && tmp->GetFirstMother()>=0){
2673  tmp = fMCEvent->Particle(tmp->GetFirstMother());
2674  same = (tmp == particle) ;
2675  }
2676  if(same){
2677  foundV0 = kTRUE ;
2678  const AliExternalTrackParam * paramPos = v0->GetParamP() ;
2679  const AliExternalTrackParam * paramNeg = v0->GetParamN() ;
2680  AliKFParticle negKF(*paramNeg,11);
2681  AliKFParticle posKF(*paramPos,-11);
2682  pConv.SetXYZM(negKF.Px()+posKF.Px(),negKF.Py()+posKF.Py(),negKF.Pz()+negKF.Pz(),0.) ;
2683  break ;
2684  }
2685  }
2686  if(foundV0){
2687  FillHistogram("hMC_CaloConv_gamma_v0",pt) ;
2688  FillHistogram("hMC_CaloConv_gamma_v0_devsE",(particle->Energy()-pConv.E())/particle->Energy(),particle->Energy()) ;
2689  }
2690 
2691  //Registered in PHOS/EMCAL
2692  Bool_t cluInPHOS = kFALSE,cluInEMCAL=kFALSE ;
2693  TLorentzVector pCalo ;
2694  Bool_t dist1=kFALSE, dist2=kFALSE ;
2695  for (Int_t i=0; i<fESDEvent->GetNumberOfCaloClusters(); i++) {
2696  AliESDCaloCluster * clu = fESDEvent->GetCaloCluster(i);
2697  Int_t iprim = clu->GetLabel() ; //# of particle hit PHOS/EMCAL
2698  Bool_t matched = kFALSE ;
2699  while(iprim>=0 ) {
2700  if(iprim==iTracks){
2701  matched=kTRUE ;
2702  break ;
2703  }
2704  else{
2705  iprim=fMCEvent->Particle(iprim)->GetFirstMother() ;
2706  }
2707  }
2708  if(!matched)
2709  continue ;
2710  if(clu->IsPHOS() && hitPHOS){
2711  cluInPHOS=kTRUE ;
2712  clu->GetMomentum(pCalo ,vtx);
2713  if(clu->GetDistanceToBadChannel()<fBadDistCutPHOS)
2714  dist1=kTRUE ;
2715  if(clu->GetDistanceToBadChannel()<2.*fBadDistCutPHOS)
2716  dist2=kTRUE ;
2717  break ;
2718  }
2719  if(!clu->IsPHOS() && hitEMCAL){
2720  cluInEMCAL=kTRUE ;
2721  clu->GetMomentum(pCalo ,vtx);
2722  if(clu->GetDistanceToBadChannel()<fBadDistCutEMCAL)
2723  dist1=kTRUE ;
2724  if(clu->GetDistanceToBadChannel()<2.*fBadDistCutEMCAL)
2725  dist2=kTRUE ;
2726  break ;
2727  }
2728  }
2729 
2730  if(cluInPHOS){
2731  FillHistogram("hMC_CaloConv_gamma_PHOSclu",pt) ;
2732  if(!dist1)
2733  FillHistogram("hMC_CaloConv_gamma_PHOSclu_dist1",pt) ;
2734  if(!dist2)
2735  FillHistogram("hMC_CaloConv_gamma_PHOSclu_dist2",pt) ;
2736  FillHistogram("hMC_CaloConv_gamma_PHOSclu_recE",pCalo.E()) ;
2737  FillHistogram("hMC_CaloConv_gamma_PHOSclu_devsE",(particle->Energy()-pCalo.E())/particle->Energy(),particle->Energy()) ;
2738  }
2739  if(cluInEMCAL){
2740  FillHistogram("hMC_CaloConv_gamma_EMCALclu",pt) ;
2741  if(!dist1)
2742  FillHistogram("hMC_CaloConv_gamma_EMCALclu_dist1",pt) ;
2743  if(!dist2)
2744  FillHistogram("hMC_CaloConv_gamma_EMCALclu_dist2",pt) ;
2745  FillHistogram("hMC_CaloConv_gamma_EMCALclu_recE",pCalo.E()) ;
2746  FillHistogram("hMC_CaloConv_gamma_EMCALclu_devsE",(particle->Energy()-pCalo.E())/particle->Energy(),particle->Energy()) ;
2747  }
2748  }
2749 }
2750 //_____________________________________________________________________________
2751 void AliAnalysisTaskCaloConv::FillHistogram(const char * key,Double_t x)const{
2752  //FillHistogram
2753  TH1F * tmp = dynamic_cast<TH1F*>(fOutputContainer->FindObject(key)) ;
2754  if(!tmp){
2755  AliInfo(Form("can not find histogram <%s> ",key)) ;
2756  return ;
2757  }
2758  tmp->Fill(x) ;
2759 }
2760 //_____________________________________________________________________________
2761 void AliAnalysisTaskCaloConv::FillHistogram(const char * key,Double_t x,Double_t y)const{
2762  //FillHistogram
2763  TObject * tmp = fOutputContainer->FindObject(key) ;
2764  if(!tmp){
2765  AliInfo(Form("can not find histogram <%s> ",key)) ;
2766  return ;
2767  }
2768  if(tmp->IsA() == TClass::GetClass("TH1F")){
2769  ((TH1F*)tmp)->Fill(x,y) ;
2770  return ;
2771  }
2772  if(tmp->IsA() == TClass::GetClass("TH2F")){
2773  ((TH2F*)tmp)->Fill(x,y) ;
2774  return ;
2775  }
2776  AliError(Form("Calling FillHistogram with 2 parameters for histo <%s> of type %s",key,tmp->IsA()->GetName())) ;
2777 }
2778 
2779 //_____________________________________________________________________________
2780 void AliAnalysisTaskCaloConv::FillHistogram(const char * key,Double_t x,Double_t y, Double_t z) const{
2781  //Fills 1D histograms with key
2782  TObject * tmp = fOutputContainer->FindObject(key) ;
2783  if(!tmp){
2784  AliInfo(Form("can not find histogram <%s> ",key)) ;
2785  return ;
2786  }
2787  if(tmp->IsA() == TClass::GetClass("TH2F")){
2788  ((TH2F*)tmp)->Fill(x,y,z) ;
2789  return ;
2790  }
2791  if(tmp->IsA() == TClass::GetClass("TH3F")){
2792  ((TH3F*)tmp)->Fill(x,y,z) ;
2793  return ;
2794  }
2795 }
2796 //______________________________________________________________________________
2797 Double_t AliAnalysisTaskCaloConv::PlanarityAngle(const AliExternalTrackParam * paramPos,const AliExternalTrackParam * paramNeg)const {
2798  //calculate angle between e+e- plain and perpendicular to MF
2799  //We need sign of MagField to calculate orienation
2800 
2801  TVector3 u(paramPos->Px()+paramNeg->Px(),paramPos->Py()+paramNeg->Py(),paramPos->Pz()+paramNeg->Pz()) ;
2802  u.Unit() ;
2803  TVector3 vPos(paramPos->Px(),paramPos->Py(),paramPos->Pz()) ;
2804  TVector3 vNeg(paramNeg->Px(),paramNeg->Py(),paramNeg->Pz()) ;
2805  TVector3 v=vPos.Cross(vNeg) ;
2806  TVector3 w = u.Cross(v);
2807  TVector3 z(0,0,1.);
2808  TVector3 ua=u.Cross(z);
2809  Double_t wa = w.Angle(ua);
2810  Double_t mfield=fESDEvent->GetMagneticField() ;
2811  if(mfield>0.)
2812  return wa; //normal field
2813  else
2814  return TMath::Pi()-wa ; //reverse field
2815 
2816 }
2817 //______________________________________________________________________________
2818 Bool_t AliAnalysisTaskCaloConv::IsGoodChannel(const char * det, Int_t mod, Int_t ix, Int_t iz){
2819 //Check if this channel belogs to the good ones
2820 
2821  if(strcmp(det,"PHOS")==0){
2822  if(mod>5 || mod<1){
2823  AliError(Form("No bad map for PHOS module %d ",mod)) ;
2824  return kTRUE ;
2825  }
2826  if(!fPHOSBadMap[mod]){
2827  AliError(Form("No Bad map for PHOS module %d",mod)) ;
2828  return kTRUE ;
2829  }
2830  if(fPHOSBadMap[mod]->GetBinContent(ix,iz)>0)
2831  return kFALSE ;
2832  else
2833  return kTRUE ;
2834  }
2835  else{
2836  if(strcmp(det,"EMCAL")==0){
2837  if(mod>9 || mod<0){
2838  AliError(Form("No bad map for EMCAL module %d ",mod)) ;
2839  return kTRUE ;
2840  }
2841  if(!fEMCALBadMap[mod]){
2842  AliError(Form("No bad map for EMCAL module %d ",mod)) ;
2843  return kTRUE ;
2844  }
2845  if(fEMCALBadMap[mod]->GetBinContent(ix,iz)>0)
2846  return kFALSE ;
2847  else
2848  return kTRUE ;
2849  }
2850  else{
2851  AliError(Form("Can not find bad channels for detector %s ",det)) ;
2852  }
2853  }
2854 
2855  return kTRUE ;
2856 }
2857 //______________________________________________________________________________
2858 void AliAnalysisTaskCaloConv::Recalibrate(Double_t &m, Double_t &pt, const TLorentzVector *calo, const TLorentzVector * conv, Int_t iw, Int_t in) {
2859  //Apply decalibration and non-linearity
2860  TLorentzVector calo2(*calo) ;
2861  Double_t en=calo2.E() ;
2862 
2863  Double_t sigma=0.06+0.005*iw ; //additional smearing
2864  //Nonlinearity
2865  Double_t a=0.02*(in%6-2.5) ;
2866  Double_t b=0.5+1.*((Int_t)in/6) ;
2867  Double_t enNew=1.-a*TMath::Exp(-en/b) ;
2868  Double_t corr=gRandom->Gaus(enNew,sigma) ;
2869  calo2*=corr ;
2870 
2871  m=(calo2+ *conv).M() ;
2872  pt=(calo2+ *conv).Pt() ;
2873 
2874 }
2875 //______________________________________________________________________________
2876 void AliAnalysisTaskCaloConv::RecalibrateEMCAL(Double_t &m, Double_t &pt, const TLorentzVector *calo, const TLorentzVector * conv, Int_t iw, Int_t in) {
2877  //Apply decalibration and non-linearity
2878  TLorentzVector calo2(*calo) ;
2879  Double_t en=calo2.E() ;
2880 
2881  Double_t sigma=0.04+0.005*iw ; //additional smearing
2882  //Nonlinearity
2883  Double_t a=0.02*(in%6-2.5) ;
2884  Double_t b=0.25+0.5*((Int_t)in/6) ;
2885  Double_t enNew=1.-a*TMath::Exp(-en/b) ;
2886  Double_t corr=gRandom->Gaus(enNew,sigma) ;
2887  calo2*=corr ;
2888 
2889  m=(calo2+ *conv).M() ;
2890  pt=(calo2+ *conv).Pt() ;
2891 
2892 }
2893 //______________________________________________________________________________
2894 void AliAnalysisTaskCaloConv::RecalibrateConvPHOS(Double_t &m, Double_t &pt, const TLorentzVector *calo, const TLorentzVector * conv, Int_t iw, Int_t in) {
2895  //Apply decalibration and non-linearity
2896 
2897  //First default PHOS smearing
2898  TLorentzVector calo2(*calo) ;
2899  Double_t en=calo2.E() ;
2900 
2901  Double_t sigma=0.065 ; //additional smearing
2902  //Nonlinearity
2903  Double_t a=0.15 ;
2904  Double_t b=0.45 ;
2905  Double_t enNew=1.+a*TMath::Exp(-en/b) ;
2906  Double_t corr=gRandom->Gaus(enNew,sigma) ;
2907  calo2*=corr ;
2908 
2909  //Now conversion photon
2910  TLorentzVector conv2(*conv) ;
2911  //linear offset in z:
2912  Double_t eta=conv2.Eta() ;
2913  Double_t c=1.e-3*iw ;
2914  eta+= c *TMath::Sign(0.9-TMath::Abs(eta),eta) ;
2915 
2916  //Smear energy and add nonlinearity
2917  //Nonlinearity
2918  Double_t enConv=conv2.E() ;
2919  Double_t ac=0.02*(in%5) ;
2920  Double_t bc=0.25+0.5*((Int_t)in/5) ;
2921  Double_t enNewc=1.+ac*TMath::Exp(-enConv/bc) ;
2922  corr=gRandom->Gaus(enNewc,0.01) ;
2923  Double_t ptc=conv2.Pt()*corr ;
2924  conv2.SetPtEtaPhiM(ptc,eta,conv2.Phi(),0.) ;
2925 
2926  m =(calo2 + conv2).M() ;
2927  pt=(calo2 + conv2).Pt() ;
2928 
2929 }
2930 //______________________________________________________________________________
2931 void AliAnalysisTaskCaloConv::GetArmenterosQtAlfa(AliKFParticle* positiveKFParticle, AliKFParticle * negativeKFParticle, AliKFParticle * gammaKFCandidate, Double_t armenterosQtAlfa[2] ){
2932  //see header file for documentation
2933 
2934  TVector3 momentumVectorPositiveKF(positiveKFParticle->GetPx(),positiveKFParticle->GetPy(),positiveKFParticle->GetPz());
2935  TVector3 momentumVectorNegativeKF(negativeKFParticle->GetPx(),negativeKFParticle->GetPy(),negativeKFParticle->GetPz());
2936  TVector3 vecV0(gammaKFCandidate->GetPx(),gammaKFCandidate->GetPy(),gammaKFCandidate->GetPz());
2937 
2938  Float_t thetaV0pos=TMath::ACos(( momentumVectorPositiveKF* vecV0)/(momentumVectorPositiveKF.Mag() * vecV0.Mag()));
2939  Float_t thetaV0neg=TMath::ACos(( momentumVectorNegativeKF* vecV0)/(momentumVectorNegativeKF.Mag() * vecV0.Mag()));
2940 
2941  Float_t alfa =((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)-(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg))/
2942  ((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)+(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)) ;
2943 
2944 
2945  Float_t qt = momentumVectorPositiveKF.Mag()*TMath::Sin(thetaV0pos);
2946 
2947  armenterosQtAlfa[0]=qt;
2948  armenterosQtAlfa[1]=alfa;
2949 
2950 }
2951 
2952 
2953 
2954 
2955 
2956 
2957 
2958 
pdg
Int_t pdg
Definition: charmCutsOptimization.C:34
Double_t
double Double_t
Definition: External.C:58
TH3F
Definition: External.C:260
TH2F
Definition: External.C:236
AliAnalysisTaskCaloConv::UserExec
virtual void UserExec(Option_t *option)
Definition: AliAnalysisTaskCaloConv.cxx:237
AliAnalysisTaskCaloConv::FillHistogram
void FillHistogram(const char *key, Double_t x) const
Definition: AliAnalysisTaskCaloConv.cxx:2751
AliAnalysisTaskCaloConv::IsGoodChannel
Bool_t IsGoodChannel(const char *det="PHOS", Int_t mod=1, Int_t ix=1, Int_t iz=1)
Definition: AliAnalysisTaskCaloConv.cxx:2818
c
TCanvas * c
Definition: TestFitELoss.C:172
AliAnalysisTaskCaloConv::fPi0Thresh1
Double_t fPi0Thresh1
EMCAL geometry.
Definition: AliAnalysisTaskCaloConv.h:131
AliAnalysisTaskCaloConv::RecalibrateConvPHOS
void RecalibrateConvPHOS(Double_t &m, Double_t &pt, const TLorentzVector *calo, const TLorentzVector *conv, Int_t iw, Int_t in)
Definition: AliAnalysisTaskCaloConv.cxx:2894
gRandom
TRandom * gRandom
kProton
Definition: FastCentEstimators.C:25
sigma
Double_t * sigma
Definition: CalculateAveragePt.C:58
AliAnalysisTaskCaloConv::GetArmenterosQtAlfa
void GetArmenterosQtAlfa(AliKFParticle *positiveKFParticle, AliKFParticle *negativeKFParticle, AliKFParticle *gammaKFCandidate, Double_t armenterosQtAlfa[2])
Definition: AliAnalysisTaskCaloConv.cxx:2931
Int_t
int Int_t
Definition: External.C:63
Float_t
float Float_t
Definition: External.C:68
AliAnalysisTaskCaloConv::fEMCALgeom
AliEMCALGeometry * fEMCALgeom
PHOS geometry.
Definition: AliAnalysisTaskCaloConv.h:130
ptmax
const Double_t ptmax
Definition: AddTaskCFDStar.C:39
AliAnalysisTaskCaloConv::PlanarityAngle
Double_t PlanarityAngle(const AliExternalTrackParam *pos, const AliExternalTrackParam *neg) const
Definition: AliAnalysisTaskCaloConv.cxx:2797
AliAnalysisTaskCaloConv::RecalibrateEMCAL
void RecalibrateEMCAL(Double_t &m, Double_t &pt, const TLorentzVector *calo, const TLorentzVector *conv, Int_t iw, Int_t in)
Definition: AliAnalysisTaskCaloConv.cxx:2876
AliAnalysisTaskCaloConv::fESDpid
AliESDpid * fESDpid
pointer to the ESDEvent
Definition: AliAnalysisTaskCaloConv.h:112
AliAnalysisTaskCaloConv::ConnectInputData
virtual void ConnectInputData(Option_t *option)
Definition: AliAnalysisTaskCaloConv.cxx:378
Option_t
const char Option_t
Definition: External.C:48
AliAnalysisTaskCaloConv::Recalibrate
void Recalibrate(Double_t &m, Double_t &pt, const TLorentzVector *calo, const TLorentzVector *conv, Int_t iw, Int_t in)
Definition: AliAnalysisTaskCaloConv.cxx:2858
pi
const Double_t pi
Definition: AliAnalysisTaskDijetHadron.cxx:72
Bool_t
bool Bool_t
Definition: External.C:53
AliAnalysisTaskCaloConv.h