Fixing bug preventing pi0s to be added to tclonesarray of reconstructed pi0 if aod...
[u/mrichter/AliRoot.git] / PWG4 / GammaConv / AliAnalysisTaskGammaConversion.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: Ana Marin, Kathrin Koch, Kenneth Aamodt                        *
5  * Version 1.1                                                            *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 ////////////////////////////////////////////////
17 //--------------------------------------------- 
18 // Class used to do analysis on conversion pairs
19 //---------------------------------------------
20 ////////////////////////////////////////////////
21
22 // root
23 #include <TChain.h>
24
25 // analysis
26 #include "AliAnalysisTaskGammaConversion.h"
27 #include "AliStack.h"
28 #include "AliLog.h"
29 #include "AliESDtrackCuts.h"
30 #include "TNtuple.h"
31 //#include "AliCFManager.h"  // for CF
32 //#include "AliCFContainer.h"   // for CF
33 #include "AliESDInputHandler.h"
34 #include "AliAnalysisManager.h"
35 #include "AliAODPWG4Particle.h"
36 #include "AliAODPWG4ParticleCorrelation.h"
37 #include "AliGammaConversionAODObject.h"
38 #include "AliAODConversionParticle.h"
39 #include "AliGammaConversionBGHandler.h"
40 #include "AliESDCaloCluster.h" // for combining PHOS and GammaConv
41 #include "AliKFVertex.h"
42 #include "AliGenPythiaEventHeader.h"
43 #include "AliGenDPMjetEventHeader.h"
44 #include "AliGenEventHeader.h"
45 #include <AliMCEventHandler.h>
46 #include "TRandom3.h"
47 #include "AliTriggerAnalysis.h"
48 #include "AliCentrality.h"
49
50 class AliESDTrackCuts;
51 class AliCFContainer;
52 class AliCFManager;
53 class AliKFVertex;
54 class AliAODHandler;
55 class AliAODEvent;
56 class ALiESDEvent;
57 class AliMCEvent;
58 class AliMCEventHandler;
59 class AliESDInputHandler;
60 class AliAnalysisManager;
61 class Riostream;
62 class TFile;
63 class TInterpreter;
64 class TSystem;
65 class TROOT;
66
67 ClassImp(AliAnalysisTaskGammaConversion)
68
69
70 AliAnalysisTaskGammaConversion::AliAnalysisTaskGammaConversion():
71 AliAnalysisTaskSE(),
72   fV0Reader(NULL),
73   fStack(NULL),
74   fMCTruth(NULL),    // for CF
75   fGCMCEvent(NULL),    // for CF
76   fESDEvent(NULL),      
77   fOutputContainer(NULL),
78   fCFManager(0x0),   // for CF
79   fHistograms(NULL),
80   fTriggerCINT1B(kFALSE),
81   fDoMCTruth(kFALSE),
82   fDoNeutralMeson(kFALSE),
83   fDoOmegaMeson(kFALSE),
84   fDoJet(kFALSE),
85   fDoChic(kFALSE),
86   fRecalculateV0ForGamma(kFALSE),
87   fKFReconstructedGammasTClone(NULL),
88   fKFReconstructedPi0sTClone(NULL),
89   fKFRecalculatedGammasTClone(NULL),
90   fCurrentEventPosElectronTClone(NULL),
91   fCurrentEventNegElectronTClone(NULL),
92   fKFReconstructedGammasCutTClone(NULL),
93   fPreviousEventTLVNegElectronTClone(NULL),
94   fPreviousEventTLVPosElectronTClone(NULL),     
95   fElectronv1(),
96   fElectronv2(),
97   fGammav1(),
98   fGammav2(),
99   fElectronRecalculatedv1(),
100   fElectronRecalculatedv2(),
101   fElectronMass(-1),
102   fGammaMass(-1),
103   fPi0Mass(-1),
104   fEtaMass(-1),
105   fGammaWidth(-1),
106   fPi0Width(-1),
107   fEtaWidth(-1),
108   fMinOpeningAngleGhostCut(0.),
109   fEsdTrackCuts(NULL),
110   fCalculateBackground(kFALSE),
111   fWriteNtuple(kFALSE),
112   fGammaNtuple(NULL),
113   fNeutralMesonNtuple(NULL),
114   fTotalNumberOfAddedNtupleEntries(0),
115   fChargedParticles(NULL),
116   fChargedParticlesId(),
117   fGammaPtHighest(0.),
118   fMinPtForGammaJet(1.),
119   fMinIsoConeSize(0.2),
120   fMinPtIsoCone(0.7),
121   fMinPtGamChargedCorr(0.5),
122   fMinPtJetCone(0.5),
123   fLeadingChargedIndex(-1),
124   fLowPtMapping(1.),
125   fHighPtMapping(3.),
126   fDoCF(kFALSE),
127   fAODGamma(NULL),
128   fAODPi0(NULL),
129   fAODOmega(NULL),
130   fAODBranchName("GammaConv"),
131   fOutputAODClassName("AliAODConversionParticle"),
132   fKFCreateAOD(kTRUE),
133   fKFForceAOD(kFALSE),
134   fKFDeltaAODFileName(""),
135   fDoNeutralMesonV0MCCheck(kFALSE),
136   fUseTrackMultiplicityForBG(kTRUE),
137   fMoveParticleAccordingToVertex(kFALSE),
138   fApplyChi2Cut(kFALSE),
139   fNRandomEventsForBG(15),
140   fNDegreesPMBackground(15),
141   fDoRotation(kTRUE),
142   fCheckBGProbability(kTRUE),
143   fKFReconstructedGammasV0Index(),
144   fRemovePileUp(kFALSE),
145   fSelectV0AND(kFALSE),
146   fTriggerAnalysis(NULL),
147   fMultiplicity(0),
148   fUseMultiplicity(0), 
149   fUseMultiplicityBin(0),
150   fUseCentrality(0), 
151   fUseCentralityBin(0)
152 {
153   // Default constructor
154
155   /*   Kenneth: the default constructor should not have any define input/output or the call to SetESDtrackCuts
156   // Common I/O in slot 0
157   DefineInput (0, TChain::Class());
158   DefineOutput(0, TTree::Class());
159         
160   // Your private output
161   DefineOutput(1, TList::Class());
162         
163   // Define standard ESD track cuts for Gamma-hadron correlation 
164   SetESDtrackCuts();
165   */
166 }
167
168 AliAnalysisTaskGammaConversion::AliAnalysisTaskGammaConversion(const char* name):
169   AliAnalysisTaskSE(name),
170   fV0Reader(NULL),
171   fStack(NULL),
172   fMCTruth(NULL),    // for CF
173   fGCMCEvent(NULL),    // for CF
174   fESDEvent(NULL),      
175   fOutputContainer(0x0),
176   fCFManager(0x0),   // for CF
177   fHistograms(NULL),
178   fTriggerCINT1B(kFALSE),
179   fDoMCTruth(kFALSE),
180   fDoNeutralMeson(kFALSE),
181   fDoOmegaMeson(kFALSE),
182   fDoJet(kFALSE),
183   fDoChic(kFALSE),
184   fRecalculateV0ForGamma(kFALSE),
185   fKFReconstructedGammasTClone(NULL),
186   fKFReconstructedPi0sTClone(NULL),
187   fKFRecalculatedGammasTClone(NULL),
188   fCurrentEventPosElectronTClone(NULL),
189   fCurrentEventNegElectronTClone(NULL),
190   fKFReconstructedGammasCutTClone(NULL),
191   fPreviousEventTLVNegElectronTClone(NULL),
192   fPreviousEventTLVPosElectronTClone(NULL),     
193   fElectronv1(),
194   fElectronv2(),
195   fGammav1(),
196   fGammav2(),
197   fElectronRecalculatedv1(),
198   fElectronRecalculatedv2(),
199   fElectronMass(-1),
200   fGammaMass(-1),
201   fPi0Mass(-1),
202   fEtaMass(-1),
203   fGammaWidth(-1),
204   fPi0Width(-1),
205   fEtaWidth(-1),
206   fMinOpeningAngleGhostCut(0.),
207   fEsdTrackCuts(NULL),
208   fCalculateBackground(kFALSE),
209   fWriteNtuple(kFALSE),
210   fGammaNtuple(NULL),
211   fNeutralMesonNtuple(NULL),
212   fTotalNumberOfAddedNtupleEntries(0),
213   fChargedParticles(NULL),
214   fChargedParticlesId(),
215   fGammaPtHighest(0.),
216   fMinPtForGammaJet(1.),
217   fMinIsoConeSize(0.2),
218   fMinPtIsoCone(0.7),
219   fMinPtGamChargedCorr(0.5),
220   fMinPtJetCone(0.5),
221   fLeadingChargedIndex(-1),
222   fLowPtMapping(1.),
223   fHighPtMapping(3.),
224   fDoCF(kFALSE),
225   fAODGamma(NULL),
226   fAODPi0(NULL),
227   fAODOmega(NULL),
228   fAODBranchName("GammaConv"),
229   fOutputAODClassName("AliAODConversionParticle"),
230   fKFCreateAOD(kTRUE),
231   fKFForceAOD(kFALSE),
232   fKFDeltaAODFileName(""),
233   fDoNeutralMesonV0MCCheck(kFALSE),
234   fUseTrackMultiplicityForBG(kTRUE),
235   fMoveParticleAccordingToVertex(kFALSE),
236   fApplyChi2Cut(kFALSE),
237   fNRandomEventsForBG(15),
238   fNDegreesPMBackground(15),
239   fDoRotation(kTRUE),
240   fCheckBGProbability(kTRUE),
241   fKFReconstructedGammasV0Index(),
242   fRemovePileUp(kFALSE),
243   fSelectV0AND(kFALSE),
244   fTriggerAnalysis(NULL),
245   fMultiplicity(0),
246   fUseMultiplicity(0), 
247   fUseMultiplicityBin(0),
248   fUseCentrality(0), 
249   fUseCentralityBin(0)
250 {
251   // Common I/O in slot 0
252   DefineInput (0, TChain::Class());
253   DefineOutput(0, TTree::Class());
254         
255   // Your private output
256   DefineOutput(1, TList::Class());
257   DefineOutput(2, AliCFContainer::Class());  // for CF
258         
259         
260   // Define standard ESD track cuts for Gamma-hadron correlation 
261   SetESDtrackCuts();
262
263 }
264
265 AliAnalysisTaskGammaConversion::~AliAnalysisTaskGammaConversion() 
266 {
267   // Remove all pointers
268         
269   if(fOutputContainer){
270     fOutputContainer->Clear() ; 
271     delete fOutputContainer ;
272   }
273   if(fHistograms){
274     delete fHistograms;
275   }
276   if(fV0Reader){
277     delete fV0Reader;
278   }
279         
280   // for CF
281   if(fCFManager){
282     delete fCFManager;
283   }
284
285   if(fEsdTrackCuts){
286     delete fEsdTrackCuts;
287   }
288
289   //Delete AODs
290   if (fAODGamma) {
291     fAODGamma->Clear();
292     delete fAODGamma;
293   }
294   fAODGamma = NULL;
295
296   if (fAODPi0) {
297     fAODPi0->Clear();
298     delete fAODPi0;
299   }
300   fAODPi0 = NULL;
301
302   if (fAODOmega) {
303     fAODOmega->Clear();
304     delete fAODOmega;
305   }
306   fAODOmega = NULL;
307
308   if(fTriggerAnalysis) {
309     delete fTriggerAnalysis;
310   }
311
312
313 }
314
315
316 void AliAnalysisTaskGammaConversion::Init()
317 {
318   // Initialization
319   // AliLog::SetGlobalLogLevel(AliLog::kError);
320 }
321 void AliAnalysisTaskGammaConversion::SetESDtrackCuts()
322 {
323   // SetESDtrackCuts
324   if (fEsdTrackCuts!=NULL){
325     delete fEsdTrackCuts;
326   }
327   fEsdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts");
328   //standard cuts from:
329   //http://aliceinfo.cern.ch/alicvs/viewvc/PWG0/dNdEta/CreateCuts.C?revision=1.4&view=markup
330
331   // Cuts used up to 3rd of March
332
333  //  fEsdTrackCuts->SetMinNClustersTPC(50);
334 //   fEsdTrackCuts->SetMaxChi2PerClusterTPC(3.5);
335 //   fEsdTrackCuts->SetRequireTPCRefit(kTRUE);
336 //   fEsdTrackCuts->SetRequireITSRefit(kTRUE);
337 //   fEsdTrackCuts->SetMaxNsigmaToVertex(3);
338 //   fEsdTrackCuts->SetRequireSigmaToVertex(kTRUE);
339
340   //------- To be tested-----------
341   // Cuts used up to 26th of Agost
342 //    Int_t minNClustersTPC = 70;
343 //    Double_t maxChi2PerClusterTPC = 4.0;
344 //    Double_t maxDCAtoVertexXY = 2.4; // cm
345 //    Double_t maxDCAtoVertexZ  = 3.2; // cm
346 //    fEsdTrackCuts->SetRequireSigmaToVertex(kFALSE);
347 //    fEsdTrackCuts->SetRequireTPCRefit(kTRUE);
348 //    fEsdTrackCuts->SetRequireITSRefit(kTRUE);
349 //    //   fEsdTrackCuts->SetRequireTPCStandAlone(kTRUE);
350 //    fEsdTrackCuts->SetAcceptKinkDaughters(kFALSE);
351 //    fEsdTrackCuts->SetMinNClustersTPC(minNClustersTPC);
352 //    fEsdTrackCuts->SetMaxChi2PerClusterTPC(maxChi2PerClusterTPC);
353 //    fEsdTrackCuts->SetMaxDCAToVertexXY(maxDCAtoVertexXY);
354 //    fEsdTrackCuts->SetMaxDCAToVertexZ(maxDCAtoVertexZ);
355 //    fEsdTrackCuts->SetDCAToVertex2D(kTRUE);
356 //    fEsdTrackCuts->SetEtaRange(-0.8, 0.8);
357 //    fEsdTrackCuts->SetPtRange(0.15);
358
359 //    fEsdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kFirst);
360
361
362 // Using standard function  for setting Cuts
363   Bool_t selectPrimaries=kTRUE;
364   fEsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(selectPrimaries);
365   fEsdTrackCuts->SetMaxDCAToVertexZ(2);
366   fEsdTrackCuts->SetEtaRange(-0.8, 0.8);
367   fEsdTrackCuts->SetPtRange(0.15);
368   
369   //-----  From Jacek 10.03.03 ------------------/
370 //     minNClustersTPC = 70;
371 //     maxChi2PerClusterTPC = 4.0;
372 //     maxDCAtoVertexXY = 2.4; // cm
373 //     maxDCAtoVertexZ  = 3.2; // cm
374
375 //     esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
376 //     esdTrackCuts->SetRequireTPCRefit(kFALSE);
377 //     esdTrackCuts->SetRequireTPCStandAlone(kTRUE);
378 //     esdTrackCuts->SetAcceptKinkDaughters(kFALSE);
379 //     esdTrackCuts->SetMinNClustersTPC(minNClustersTPC);
380 //     esdTrackCuts->SetMaxChi2PerClusterTPC(maxChi2PerClusterTPC);
381 //     esdTrackCuts->SetMaxDCAToVertexXY(maxDCAtoVertexXY);
382 //     esdTrackCuts->SetMaxDCAToVertexZ(maxDCAtoVertexZ);
383 //     esdTrackCuts->SetDCAToVertex2D(kTRUE);
384   
385
386
387   //  fEsdTrackCuts->SetAcceptKinkDaughters(kFALSE);
388   //  fV0Reader->SetESDtrackCuts(fEsdTrackCuts);
389 }
390
391 void AliAnalysisTaskGammaConversion::UserExec(Option_t */*option*/)
392 {
393   // Execute analysis for current event
394
395   //  Load the esdpid from the esdhandler if exists (tender was applied) otherwise set the Bethe Bloch parameters
396   Int_t eventQuality=-1;
397
398   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
399   AliESDInputHandler *esdHandler=0x0;
400   if ( (esdHandler=dynamic_cast<AliESDInputHandler*>(man->GetInputEventHandler())) && esdHandler->GetESDpid() ){
401     AliV0Reader::SetESDpid(esdHandler->GetESDpid());
402   } else {
403     //load esd pid bethe bloch parameters depending on the existance of the MC handler
404     // yes: MC parameters
405     // no:  data parameters
406     if (!AliV0Reader::GetESDpid()){
407       if (fMCEvent ) {
408         AliV0Reader::InitESDpid();
409       } else {
410         AliV0Reader::InitESDpid(1);
411       }
412     }
413   } 
414
415   //Must set fForceAOD to true for the AOD to get filled. Should only be done when running independent chain / train. 
416   if(fKFForceAOD) {
417     if (!AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()) AliFatal("Cannot run ESD filter without an output event handler");
418     AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kTRUE);
419   }
420
421   if(fV0Reader == NULL){
422     // Write warning here cuts and so on are default if this ever happens
423   }
424
425   if (fMCEvent ) {
426     // To avoid crashes due to unzip errors. Sometimes the trees are not there.
427
428     AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
429     if (!mcHandler){ 
430       AliError("Could not retrive MC event handler!"); 
431       return; 
432
433       eventQuality=0;
434       fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
435     }
436     if (!mcHandler->InitOk() ){
437       eventQuality=0;
438       fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
439       return;
440     }
441     if (!mcHandler->TreeK() ){
442       eventQuality=0;
443       fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
444       return;
445     }
446     if (!mcHandler->TreeTR() ) {
447       eventQuality=0;
448       fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
449       return;
450     }
451   }
452
453   fV0Reader->SetInputAndMCEvent(InputEvent(), MCEvent());
454
455   fV0Reader->Initialize();
456   fDoMCTruth = fV0Reader->GetDoMCTruth();
457
458   if(fAODGamma) fAODGamma->Delete();
459   if(fAODPi0) fAODPi0->Delete();
460   if(fAODOmega) fAODOmega->Delete();
461         
462   if(fKFReconstructedGammasTClone == NULL){
463     fKFReconstructedGammasTClone = new TClonesArray("AliKFParticle",0);
464   }
465   if(fCurrentEventPosElectronTClone == NULL){
466     fCurrentEventPosElectronTClone = new TClonesArray("AliESDtrack",0);
467   }
468   if(fCurrentEventNegElectronTClone == NULL){
469     fCurrentEventNegElectronTClone = new TClonesArray("AliESDtrack",0);
470   }
471   if(fKFReconstructedGammasCutTClone == NULL){
472     fKFReconstructedGammasCutTClone = new TClonesArray("AliKFParticle",0);
473   }
474   if(fPreviousEventTLVNegElectronTClone == NULL){
475     fPreviousEventTLVNegElectronTClone = new TClonesArray("TLorentzVector",0);
476   }
477   if(fPreviousEventTLVPosElectronTClone == NULL){
478     fPreviousEventTLVPosElectronTClone  = new TClonesArray("TLorentzVector",0);
479   }
480   if(fChargedParticles == NULL){
481     fChargedParticles = new TClonesArray("AliESDtrack",0);
482   }
483
484   if(fKFReconstructedPi0sTClone == NULL){
485     fKFReconstructedPi0sTClone = new TClonesArray("AliKFParticle",0);
486   }
487  
488   if(fKFRecalculatedGammasTClone == NULL){
489     fKFRecalculatedGammasTClone = new TClonesArray("AliKFParticle",0);
490   }
491
492   if(fTriggerAnalysis== NULL){
493     fTriggerAnalysis = new AliTriggerAnalysis;
494   }
495
496   //clear TClones
497   fKFReconstructedGammasTClone->Delete();
498   fCurrentEventPosElectronTClone->Delete();
499   fCurrentEventNegElectronTClone->Delete();
500   fKFReconstructedGammasCutTClone->Delete();
501   fPreviousEventTLVNegElectronTClone->Delete();
502   fPreviousEventTLVPosElectronTClone->Delete();
503   fKFReconstructedPi0sTClone->Delete();
504   fKFRecalculatedGammasTClone->Delete();
505
506   //clear vectors
507   fElectronv1.clear();
508   fElectronv2.clear();
509
510   fGammav1.clear();
511   fGammav2.clear();
512
513   fElectronRecalculatedv1.clear();
514   fElectronRecalculatedv2.clear();
515
516         
517   fChargedParticles->Delete();  
518
519   fChargedParticlesId.clear();  
520
521   fKFReconstructedGammasV0Index.clear();
522         
523   //Clear the data in the v0Reader
524   //  fV0Reader->UpdateEventByEventData();
525
526   //Take Only events with proper trigger
527   /*
528   if(fTriggerCINT1B){
529     if(!fV0Reader->GetESDEvent()->IsTriggerClassFired("CINT1B-ABCE-NOPF-ALL")) return;
530   }
531   */
532
533   if(fV0Reader->CheckForPrimaryVertex() == kFALSE){
534     //    cout<< "Event not taken"<< endl;
535     eventQuality=1;
536     fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
537     if(fDoMCTruth){
538       CheckMesonProcessTypeEventQuality(eventQuality);
539     }
540     return; // aborts if the primary vertex does not have contributors.
541   }
542
543
544   if(!fV0Reader->CheckForPrimaryVertexZ() ){
545     eventQuality=2;
546     fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
547     if(fDoMCTruth){
548       CheckMesonProcessTypeEventQuality(eventQuality);
549     }
550     return;
551   }
552
553   if(fRemovePileUp && fV0Reader->GetESDEvent()->IsPileupFromSPD()) {
554     eventQuality=4;
555     fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
556     return;
557   }
558
559   Bool_t v0A       = fTriggerAnalysis->IsOfflineTriggerFired(fV0Reader->GetESDEvent(), AliTriggerAnalysis::kV0A);
560   Bool_t v0C       = fTriggerAnalysis->IsOfflineTriggerFired(fV0Reader->GetESDEvent(), AliTriggerAnalysis::kV0C);
561   Bool_t v0AND = v0A && v0C;
562
563   if(fSelectV0AND && !v0AND){
564     eventQuality=5;
565     fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
566     return;
567   }
568   fMultiplicity = fEsdTrackCuts->CountAcceptedTracks(fV0Reader->GetESDEvent());
569
570   if( CalculateMultiplicityBin() != fUseMultiplicityBin){
571     eventQuality=6;
572     fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
573     return;
574   }
575
576   if(fV0Reader->GetIsHeavyIon()){
577     if(fUseCentrality>0){
578       AliCentrality *esdCentrality = fV0Reader->GetESDEvent()->GetCentrality();
579       Int_t centralityC = -1;
580
581       if(fUseCentrality==1){
582         centralityC = esdCentrality->GetCentralityClass10("V0M");
583         if( centralityC != fUseCentralityBin ){
584           eventQuality=7;
585           fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
586           return;
587         }
588       }
589
590       if(fUseCentrality==2){
591         centralityC = esdCentrality->GetCentralityClass10("CL1");
592         if( centralityC != fUseCentralityBin ){
593           eventQuality=7;
594           fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
595           return;
596         }
597       }
598     }
599   }
600   eventQuality=3;
601   fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
602
603
604
605   fHistograms->FillHistogram("ESD_NumberOfGoodESDTracks",fMultiplicity);
606   if (fV0Reader->GetNumberOfContributorsVtx()>=1){
607     fHistograms->FillHistogram("ESD_NumberOfGoodESDTracksVtx",fMultiplicity);
608   } 
609
610
611
612   // Process the MC information
613   if(fDoMCTruth){
614     ProcessMCData();
615   }
616         
617   //Process the v0 information with no cuts
618   ProcessV0sNoCut();
619
620   // Process the v0 information
621   ProcessV0s();
622   
623
624   //Fill Gamma AOD
625   if(fKFCreateAOD) {
626     FillAODWithConversionGammas() ; 
627   }
628
629
630   // Process reconstructed gammas
631   if(fDoNeutralMeson == kTRUE){
632     ProcessGammasForNeutralMesonAnalysis();
633
634   }
635   
636   if(fDoMCTruth == kTRUE){
637     CheckV0Efficiency();
638   }
639   //Process reconstructed gammas electrons for Chi_c Analysis
640   if(fDoChic == kTRUE){
641     ProcessGammaElectronsForChicAnalysis();
642   }
643   // Process reconstructed gammas for gamma Jet/hadron correlations
644   if(fDoJet == kTRUE){
645     ProcessGammasForGammaJetAnalysis();
646   }
647         
648   //calculate background if flag is set
649   if(fCalculateBackground){
650     CalculateBackground();
651   }
652
653    if(fDoNeutralMeson == kTRUE){
654 //     ProcessConvPHOSGammasForNeutralMesonAnalysis();
655      if(fDoOmegaMeson == kTRUE){
656        ProcessGammasForOmegaMesonAnalysis();
657      }
658    }
659
660   //Clear the data in the v0Reader
661   fV0Reader->UpdateEventByEventData();
662   if(fRecalculateV0ForGamma==kTRUE){
663     RecalculateV0ForGamma();
664   }
665   PostData(1, fOutputContainer);
666   PostData(2, fCFManager->GetParticleContainer());  // for CF
667         
668 }
669
670 // void AliAnalysisTaskGammaConversion::ConnectInputData(Option_t *option){
671 //   // see header file for documentation
672 //   //  printf("   ConnectInputData %s\n", GetName());
673
674 //   AliAnalysisTaskSE::ConnectInputData(option);
675
676 //   if(fV0Reader == NULL){
677 //     // Write warning here cuts and so on are default if this ever happens
678 //   }
679 //   fV0Reader->Initialize();
680 //   fDoMCTruth = fV0Reader->GetDoMCTruth();
681 // }
682
683 void AliAnalysisTaskGammaConversion::CheckMesonProcessTypeEventQuality(Int_t evtQ){
684   // Check meson process type event quality
685   fStack= MCEvent()->Stack();
686   fGCMCEvent=MCEvent();
687
688   for (Int_t iTracks = 0; iTracks < fStack->GetNprimary(); iTracks++) {
689     TParticle* particle = (TParticle *)fStack->Particle(iTracks);
690     if (!particle) {
691       //print warning here
692       continue;
693     }
694     if(particle->GetPdgCode()!=111){     //Pi0
695       continue;
696     }
697     if(TMath::Abs(particle->Eta())> fV0Reader->GetEtaCut() )    continue;
698     if(evtQ==1){
699       switch(GetProcessType(fGCMCEvent)){
700       case  kProcSD:
701         fHistograms->FillHistogram("MC_SD_EvtQ1_Pi0_Pt", particle->Pt());
702         break;
703       case  kProcDD:
704         fHistograms->FillHistogram("MC_DD_EvtQ1_Pi0_Pt", particle->Pt());
705         break;
706       case  kProcND:
707         fHistograms->FillHistogram("MC_ND_EvtQ1_Pi0_Pt", particle->Pt());
708         break;
709       default:
710         AliError("Unknown Process");
711       }
712     }
713     if(evtQ==2){
714       switch(GetProcessType(fGCMCEvent)){
715       case  kProcSD:
716         fHistograms->FillHistogram("MC_SD_EvtQ2_Pi0_Pt", particle->Pt());
717         break;
718       case  kProcDD:
719         fHistograms->FillHistogram("MC_DD_EvtQ2_Pi0_Pt", particle->Pt());
720         break;
721       case  kProcND:
722         fHistograms->FillHistogram("MC_ND_EvtQ2_Pi0_Pt", particle->Pt());
723         break;
724       default:
725         AliError("Unknown Process");
726       }
727     }
728
729
730   }
731
732 }
733
734 void AliAnalysisTaskGammaConversion::ProcessMCData(){
735   // see header file for documentation
736   //InputEvent(), MCEvent());
737   /* TestAnaMarin
738   fStack = fV0Reader->GetMCStack();
739   fMCTruth = fV0Reader->GetMCTruth();  // for CF
740   fGCMCEvent = fV0Reader->GetMCEvent();  // for CF
741   */
742   fStack= MCEvent()->Stack();
743   fGCMCEvent=MCEvent();
744         
745   // for CF
746   Double_t containerInput[3];
747   if(fDoCF){
748     if(!fGCMCEvent) cout << "NO MC INFO FOUND" << endl;
749     fCFManager->SetEventInfo(fGCMCEvent);
750   } 
751   // end for CF
752         
753         
754   if(fV0Reader->CheckForPrimaryVertex() == kFALSE){
755     return; // aborts if the primary vertex does not have contributors.
756   }
757   for (Int_t iTracks = 0; iTracks < fStack->GetNprimary(); iTracks++) {
758     //  for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++) {
759     TParticle* particle = (TParticle *)fStack->Particle(iTracks);
760
761
762
763     if (!particle) {
764       //print warning here
765       continue;
766     }
767                 
768     ///////////////////////Begin Chic Analysis/////////////////////////////
769     if(fDoChic) {
770       if(particle->GetPdgCode() == 443){//Is JPsi       
771         if(particle->GetNDaughters()==2){
772           if(TMath::Abs(fStack->Particle(particle->GetFirstDaughter())->GetPdgCode()) == 11 &&
773              TMath::Abs(fStack->Particle(particle->GetLastDaughter())->GetPdgCode()) == 11){
774
775             TParticle* daug0 = fStack->Particle(particle->GetFirstDaughter());
776             TParticle* daug1 = fStack->Particle(particle->GetLastDaughter());
777             if(TMath::Abs(daug0->Eta()) < 0.9 && TMath::Abs(daug1->Eta()) < 0.9)
778               fHistograms->FillTable("Table_Electrons",3);//e+ e-  from J/Psi inside acceptance
779                                         
780             if( TMath::Abs(daug0->Eta()) < 0.9){
781               if(daug0->GetPdgCode() == -11)
782                 fHistograms->FillTable("Table_Electrons",1);//e+  from J/Psi inside acceptance
783               else
784                 fHistograms->FillTable("Table_Electrons",2);//e-   from J/Psi inside acceptance
785                                                 
786             }
787             if(TMath::Abs(daug1->Eta()) < 0.9){
788               if(daug1->GetPdgCode() == -11)
789                 fHistograms->FillTable("Table_Electrons",1);//e+  from J/Psi inside acceptance
790               else
791                 fHistograms->FillTable("Table_Electrons",2);//e-   from J/Psi inside acceptance
792             }
793           }
794         }
795       }
796       //              const int CHI_C0   = 10441;
797       //              const int CHI_C1   = 20443;
798       //              const int CHI_C2   = 445
799       if(particle->GetPdgCode() == 22){//gamma from JPsi
800         if(particle->GetMother(0) > -1){
801           if(fStack->Particle(particle->GetMother(0))->GetPdgCode() == 10441 ||
802              fStack->Particle(particle->GetMother(0))->GetPdgCode() == 20443 ||
803              fStack->Particle(particle->GetMother(0))->GetPdgCode() == 445){
804             if(TMath::Abs(particle->Eta()) < 1.2)
805               fHistograms->FillTable("Table_Electrons",17);// gamma from chic inside accptance
806           }
807         }
808       }
809       if(particle->GetPdgCode() == 10441 || particle->GetPdgCode() == 20443 || particle->GetPdgCode() == 445){
810         if( particle->GetNDaughters() == 2){
811           TParticle* daug0 = fStack->Particle(particle->GetFirstDaughter());
812           TParticle* daug1 = fStack->Particle(particle->GetLastDaughter());
813                                 
814           if( (daug0->GetPdgCode() == 443 || daug0->GetPdgCode() == 22) && (daug1->GetPdgCode() == 443 || daug1->GetPdgCode() == 22) ){
815             if( daug0->GetPdgCode() == 443){
816               TParticle* daugE0 = fStack->Particle(daug0->GetFirstDaughter());
817               TParticle* daugE1 = fStack->Particle(daug0->GetLastDaughter());
818               if( TMath::Abs(daug1->Eta()) < 1.2 && TMath::Abs(daugE0->Eta()) < 0.9 && TMath::Abs(daugE1->Eta()) < 0.9 )
819                 fHistograms->FillTable("Table_Electrons",18);
820                                                 
821             }//if
822             else if (daug1->GetPdgCode() == 443){
823               TParticle* daugE0 = fStack->Particle(daug1->GetFirstDaughter());
824               TParticle* daugE1 = fStack->Particle(daug1->GetLastDaughter());
825               if( TMath::Abs(daug0->Eta()) < 1.2 && TMath::Abs(daugE0->Eta()) < 0.9 && TMath::Abs(daugE1->Eta()) < 0.9 )
826                 fHistograms->FillTable("Table_Electrons",18);
827             }//else if
828           }//gamma o Jpsi
829         }//GetNDaughters
830       }
831     }
832                 
833     /////////////////////End Chic Analysis////////////////////////////
834                 
835                 
836     //    if(TMath::Abs(particle->Eta())> fV0Reader->GetEtaCut() )      continue;
837                 
838     if(particle->R()>fV0Reader->GetMaxRCut())   continue; // cuts on distance from collision point
839                 
840     Double_t tmpPhi=particle->Phi();
841                 
842     if(particle->Phi()> TMath::Pi()){
843       tmpPhi = particle->Phi()-(2*TMath::Pi());
844     }
845                 
846     Double_t rapidity;
847     if(particle->Energy() - particle->Pz() == 0 || particle->Energy() + particle->Pz() == 0){
848       rapidity=0;
849     }
850     else{
851       rapidity = 0.5*(TMath::Log((particle->Energy()+particle->Pz()) / (particle->Energy()-particle->Pz())));
852     }   
853                 
854
855
856     if(iTracks<=fStack->GetNprimary() ){
857       if ( particle->GetPdgCode()== -211 ||  particle->GetPdgCode()== 211 ||
858            particle->GetPdgCode()== 2212 ||  particle->GetPdgCode()==-2212 ||
859            particle->GetPdgCode()== 321  ||  particle->GetPdgCode()==-321 ){
860         if(TMath::Abs(particle->Eta())> fV0Reader->GetEtaCut() )        continue;
861         fHistograms->FillHistogram("MC_PhysicalPrimaryCharged_Pt", particle->Pt());
862       }
863     }
864
865  
866
867     //process the gammas
868     if (particle->GetPdgCode() == 22){
869       if(TMath::Abs(particle->Eta())> fV0Reader->GetEtaCut() )  continue;      
870
871       if(particle->GetMother(0) >-1 && fStack->Particle(particle->GetMother(0))->GetPdgCode() == 22){
872         continue; // no photon as mothers!
873       }
874       
875       if(particle->GetMother(0) >= fStack->GetNprimary()){
876         continue; // the gamma has a mother, and it is not a primary particle
877       }
878                 
879       if(particle->GetMother(0) >-1){
880          fHistograms->FillHistogram("MC_DecayAllGamma_Pt", particle->Pt()); // All
881          switch(fStack->Particle(particle->GetMother(0))->GetPdgCode()){
882          case 111: // Pi0
883             fHistograms->FillHistogram("MC_DecayPi0Gamma_Pt", particle->Pt());
884             break;
885          case 113: // Rho0
886             fHistograms->FillHistogram("MC_DecayRho0Gamma_Pt", particle->Pt());
887             break;
888          case 221: // Eta
889             fHistograms->FillHistogram("MC_DecayEtaGamma_Pt", particle->Pt());
890             break;
891          case 223: // Omega
892             fHistograms->FillHistogram("MC_DecayOmegaGamma_Pt", particle->Pt());
893             break;
894          case 310: // K_s0
895             fHistograms->FillHistogram("MC_DecayK0sGamma_Pt", particle->Pt());
896             break;
897          case 331: // Eta'
898             fHistograms->FillHistogram("MC_DecayEtapGamma_Pt", particle->Pt());
899             break;
900          }
901       }
902         
903       fHistograms->FillHistogram("MC_allGamma_Energy", particle->Energy());
904       fHistograms->FillHistogram("MC_allGamma_Pt", particle->Pt());
905       fHistograms->FillHistogram("MC_allGamma_Eta", particle->Eta());
906       fHistograms->FillHistogram("MC_allGamma_Phi", tmpPhi);
907       fHistograms->FillHistogram("MC_allGamma_Rapid", rapidity);
908                         
909       // for CF
910       if(fDoCF){
911         containerInput[0] = particle->Pt();
912         containerInput[1] = particle->Eta();
913         if(particle->GetMother(0) >=0){
914           containerInput[2] = fStack->Particle(particle->GetMother(0))->GetMass();
915         }
916         else{
917           containerInput[2]=-1;
918         }
919       
920         fCFManager->GetParticleContainer()->Fill(containerInput,kStepGenerated);                                        // generated gamma
921       }         
922       if(particle->GetMother(0) < 0){   // direct gamma
923         fHistograms->FillHistogram("MC_allDirectGamma_Energy",particle->Energy());
924         fHistograms->FillHistogram("MC_allDirectGamma_Pt", particle->Pt());
925         fHistograms->FillHistogram("MC_allDirectGamma_Eta", particle->Eta());
926         fHistograms->FillHistogram("MC_allDirectGamma_Phi", tmpPhi);
927         fHistograms->FillHistogram("MC_allDirectGamma_Rapid", rapidity);                                
928       }
929                         
930       // looking for conversion (electron + positron from pairbuilding (= 5) )
931       TParticle* ePos = NULL;
932       TParticle* eNeg = NULL;
933                         
934       if(particle->GetNDaughters() >= 2){
935         for(Int_t daughterIndex=particle->GetFirstDaughter();daughterIndex<=particle->GetLastDaughter();daughterIndex++){
936           TParticle *tmpDaughter = fStack->Particle(daughterIndex);
937           if(tmpDaughter->GetUniqueID() == 5){
938             if(tmpDaughter->GetPdgCode() == 11){
939               eNeg = tmpDaughter;
940             }
941             else if(tmpDaughter->GetPdgCode() == -11){
942               ePos = tmpDaughter;
943             }
944           }
945         }
946       }
947                         
948                         
949       if(ePos == NULL || eNeg == NULL){ // means we do not have two daughters from pair production
950         continue;
951       }
952                         
953                         
954       Double_t ePosPhi = ePos->Phi();
955       if(ePos->Phi()> TMath::Pi()) ePosPhi = ePos->Phi()-(2*TMath::Pi());
956                         
957       Double_t eNegPhi = eNeg->Phi();
958       if(eNeg->Phi()> TMath::Pi()) eNegPhi = eNeg->Phi()-(2*TMath::Pi());
959                         
960                         
961       if(ePos->Pt()<fV0Reader->GetPtCut() || eNeg->Pt()<fV0Reader->GetPtCut()){
962         continue; // no reconstruction below the Pt cut
963       }
964                         
965       if(TMath::Abs(ePos->Eta())> fV0Reader->GetEtaCut() || TMath::Abs(eNeg->Eta())> fV0Reader->GetEtaCut()){
966         continue;
967       } 
968                         
969       if(ePos->R()>fV0Reader->GetMaxRCut()){
970         continue; // cuts on distance from collision point
971       }
972
973       if(TMath::Abs(ePos->Vz()) > fV0Reader->GetMaxZCut()){
974         continue;   // outside material
975       }
976                         
977                         
978       if((TMath::Abs(ePos->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()  > ePos->R()){
979         continue;               // line cut to exclude regions where we do not reconstruct
980       }         
981                 
982                         
983       // for CF
984       if(fDoCF){
985         fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructable);  // reconstructable gamma        
986       }
987       fHistograms->FillHistogram("MC_ConvGamma_Energy", particle->Energy());
988       fHistograms->FillHistogram("MC_ConvGamma_Pt", particle->Pt());
989       fHistograms->FillHistogram("MC_ConvGamma_Eta", particle->Eta());
990       fHistograms->FillHistogram("MC_ConvGamma_Phi", tmpPhi);
991       fHistograms->FillHistogram("MC_ConvGamma_Rapid", rapidity);
992       fHistograms->FillHistogram("MC_ConvGamma_Pt_Eta", particle->Pt(),particle->Eta());
993                         
994       fHistograms->FillHistogram("MC_E_Energy", eNeg->Energy());
995       fHistograms->FillHistogram("MC_E_Pt", eNeg->Pt());
996       fHistograms->FillHistogram("MC_E_Eta", eNeg->Eta());
997       fHistograms->FillHistogram("MC_E_Phi", eNegPhi);
998                         
999       fHistograms->FillHistogram("MC_P_Energy", ePos->Energy());
1000       fHistograms->FillHistogram("MC_P_Pt", ePos->Pt());
1001       fHistograms->FillHistogram("MC_P_Eta", ePos->Eta());
1002       fHistograms->FillHistogram("MC_P_Phi", ePosPhi);
1003                         
1004                         
1005       // begin Mapping 
1006       Int_t rBin    = fHistograms->GetRBin(ePos->R());
1007       Int_t zBin    = fHistograms->GetZBin(ePos->Vz());
1008       Int_t phiBin  = fHistograms->GetPhiBin(particle->Phi());
1009       Double_t rFMD=30;
1010       Double_t rITSTPCMin=50;
1011       Double_t rITSTPCMax=80;
1012
1013       TVector3 vtxPos(ePos->Vx(),ePos->Vy(),ePos->Vz());        
1014       
1015       TString nameMCMappingPhiR="";
1016       nameMCMappingPhiR.Form("MC_Conversion_Mapping_Phi%02d_R%02d",phiBin,rBin);
1017       // fHistograms->FillHistogram(nameMCMappingPhiR, ePos->Vz(), particle->Eta());
1018                         
1019       TString nameMCMappingPhi="";
1020       nameMCMappingPhi.Form("MC_Conversion_Mapping_Phi%02d",phiBin);
1021       //      fHistograms->FillHistogram(nameMCMappingPhi, particle->Eta());
1022       //fHistograms->FillHistogram(nameMCMappingPhi, ePos->Vz(), particle->Eta());
1023                         
1024       TString nameMCMappingR="";
1025       nameMCMappingR.Form("MC_Conversion_Mapping_R%02d",rBin);
1026       //      fHistograms->FillHistogram(nameMCMappingR, particle->Eta());
1027       //fHistograms->FillHistogram(nameMCMappingR,ePos->Vz(), particle->Eta());
1028                         
1029       TString nameMCMappingPhiInR="";
1030       nameMCMappingPhiInR.Form("MC_Conversion_Mapping_Phi_in_R_%02d",rBin);
1031       //      fHistograms->FillHistogram(nameMCMappingPhiInR, tmpPhi);
1032       fHistograms->FillHistogram(nameMCMappingPhiInR, vtxPos.Phi());
1033
1034       TString nameMCMappingZInR="";
1035       nameMCMappingZInR.Form("MC_Conversion_Mapping_Z_in_R_%02d",rBin);
1036       fHistograms->FillHistogram(nameMCMappingZInR,ePos->Vz() );
1037
1038
1039       TString nameMCMappingPhiInZ="";
1040       nameMCMappingPhiInZ.Form("MC_Conversion_Mapping_Phi_in_Z_%02d",zBin);
1041       //      fHistograms->FillHistogram(nameMCMappingPhiInR, tmpPhi);
1042       fHistograms->FillHistogram(nameMCMappingPhiInZ, vtxPos.Phi());
1043
1044
1045       if(ePos->R()<rFMD){
1046         TString nameMCMappingFMDPhiInZ="";
1047         nameMCMappingFMDPhiInZ.Form("MC_Conversion_Mapping_FMD_Phi_in_Z_%02d",zBin);
1048         fHistograms->FillHistogram(nameMCMappingFMDPhiInZ, vtxPos.Phi());
1049       }
1050
1051       if(ePos->R()>rITSTPCMin  && ePos->R()<rITSTPCMax){
1052         TString nameMCMappingITSTPCPhiInZ="";
1053         nameMCMappingITSTPCPhiInZ.Form("MC_Conversion_Mapping_ITSTPC_Phi_in_Z_%02d",zBin);
1054         fHistograms->FillHistogram(nameMCMappingITSTPCPhiInZ, vtxPos.Phi());
1055       }
1056
1057       TString nameMCMappingRInZ="";
1058       nameMCMappingRInZ.Form("MC_Conversion_Mapping_R_in_Z_%02d",zBin);
1059       fHistograms->FillHistogram(nameMCMappingRInZ,ePos->R() );
1060
1061       if(particle->Pt() > fLowPtMapping && particle->Pt()< fHighPtMapping){
1062         TString nameMCMappingMidPtPhiInR="";
1063         nameMCMappingMidPtPhiInR.Form("MC_Conversion_Mapping_MidPt_Phi_in_R_%02d",rBin);
1064         fHistograms->FillHistogram(nameMCMappingMidPtPhiInR, vtxPos.Phi());
1065         
1066         TString nameMCMappingMidPtZInR="";
1067         nameMCMappingMidPtZInR.Form("MC_Conversion_Mapping_MidPt_Z_in_R_%02d",rBin);
1068         fHistograms->FillHistogram(nameMCMappingMidPtZInR,ePos->Vz() );
1069         
1070         
1071         TString nameMCMappingMidPtPhiInZ="";
1072         nameMCMappingMidPtPhiInZ.Form("MC_Conversion_Mapping_MidPt_Phi_in_Z_%02d",zBin);
1073         fHistograms->FillHistogram(nameMCMappingMidPtPhiInZ, vtxPos.Phi());
1074
1075
1076         if(ePos->R()<rFMD){
1077           TString nameMCMappingMidPtFMDPhiInZ="";
1078           nameMCMappingMidPtFMDPhiInZ.Form("MC_Conversion_Mapping_MidPt_FMD_Phi_in_Z_%02d",zBin);
1079           fHistograms->FillHistogram(nameMCMappingMidPtFMDPhiInZ, vtxPos.Phi());
1080         }
1081
1082         
1083         TString nameMCMappingMidPtRInZ="";
1084         nameMCMappingMidPtRInZ.Form("MC_Conversion_Mapping_MidPt_R_in_Z_%02d",zBin);
1085         fHistograms->FillHistogram(nameMCMappingMidPtRInZ,ePos->R() );
1086
1087       }
1088
1089       //end mapping
1090                         
1091       fHistograms->FillHistogram("MC_Conversion_R",ePos->R());
1092       fHistograms->FillHistogram("MC_Conversion_ZR",ePos->Vz(),ePos->R());
1093       fHistograms->FillHistogram("MC_Conversion_XY",ePos->Vx(),ePos->Vy());
1094       fHistograms->FillHistogram("MC_Conversion_OpeningAngle",GetMCOpeningAngle(ePos, eNeg));
1095       fHistograms->FillHistogram("MC_ConvGamma_E_AsymmetryP",particle->P(),eNeg->P()/particle->P());
1096       fHistograms->FillHistogram("MC_ConvGamma_P_AsymmetryP",particle->P(),ePos->P()/particle->P());
1097
1098
1099       if(particle->GetMother(0) < 0){ // no mother = direct gamma, still inside converted
1100         fHistograms->FillHistogram("MC_ConvDirectGamma_Energy",particle->Energy());
1101         fHistograms->FillHistogram("MC_ConvDirectGamma_Pt", particle->Pt());
1102         fHistograms->FillHistogram("MC_ConvDirectGamma_Eta", particle->Eta());
1103         fHistograms->FillHistogram("MC_ConvDirectGamma_Phi", tmpPhi);
1104         fHistograms->FillHistogram("MC_ConvDirectGamma_Rapid", rapidity);
1105                                 
1106       } // end direct gamma
1107       else{   // mother exits 
1108         /*      if( fStack->Particle(particle->GetMother(0))->GetPdgCode()==10441 ||//chic0 
1109                 fStack->Particle(particle->GetMother(0))->GetPdgCode()==20443 ||//psi2S
1110                 fStack->Particle(particle->GetMother(0))->GetPdgCode()==445  //chic2
1111                 ){ 
1112                 fMCGammaChic.push_back(particle);
1113                 }
1114         */
1115       }  // end if mother exits
1116     } // end if particle is a photon
1117                 
1118                 
1119                 
1120     // process motherparticles (2 gammas as daughters)
1121     // the motherparticle had already to pass the R and the eta cut, but no line cut.
1122     // the line cut is just valid for the conversions!
1123                 
1124     if(particle->GetNDaughters() == 2){
1125                         
1126       TParticle* daughter0 = (TParticle*)fStack->Particle(particle->GetFirstDaughter());
1127       TParticle* daughter1 = (TParticle*)fStack->Particle(particle->GetLastDaughter());
1128                         
1129       if(daughter0->GetPdgCode() != 22 || daughter1->GetPdgCode() != 22) continue; //check for gamma gamma daughters
1130
1131       if(TMath::Abs(rapidity) > fV0Reader->GetRapidityMesonCut() ) continue; 
1132
1133       // Check the acceptance for both gammas
1134       Bool_t gammaEtaCut = kTRUE;
1135       if(TMath::Abs(daughter0->Eta()) > fV0Reader->GetEtaCut() || TMath::Abs(daughter1->Eta()) > fV0Reader->GetEtaCut()  ) gammaEtaCut = kFALSE;
1136       Bool_t gammaRCut = kTRUE;
1137       if(daughter0->R() > fV0Reader->GetMaxRCut() || daughter1->R() > fV0Reader->GetMaxRCut()  ) gammaRCut = kFALSE;
1138                         
1139                         
1140                         
1141       // check for conversions now -> have to pass eta, R and line cut!
1142       Bool_t daughter0Electron = kFALSE;
1143       Bool_t daughter0Positron = kFALSE;
1144       Bool_t daughter1Electron = kFALSE;
1145       Bool_t daughter1Positron = kFALSE;
1146                         
1147       if(daughter0->GetNDaughters() >= 2){  // first gamma
1148         for(Int_t TrackIndex=daughter0->GetFirstDaughter();TrackIndex<=daughter0->GetLastDaughter();TrackIndex++){
1149           TParticle *tmpDaughter = fStack->Particle(TrackIndex);
1150           if(tmpDaughter->GetUniqueID() == 5){
1151             if(tmpDaughter->GetPdgCode() == 11){
1152               if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){
1153                 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()  < tmpDaughter->R() ){
1154                   if(tmpDaughter->R()< fV0Reader->GetMaxRCut()){
1155                     daughter0Electron = kTRUE;
1156                   }
1157                 }
1158                                                                 
1159               }
1160             }
1161             else if(tmpDaughter->GetPdgCode() == -11){
1162               if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){
1163                 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()  < tmpDaughter->R() ){
1164                   if(tmpDaughter->R()< fV0Reader->GetMaxRCut()){
1165                     daughter0Positron = kTRUE;
1166                   }
1167                 }                                               
1168               }                                 
1169             }
1170           }
1171         }
1172       }
1173                         
1174                         
1175       if(daughter1->GetNDaughters() >= 2){   // second gamma
1176         for(Int_t TrackIndex=daughter1->GetFirstDaughter();TrackIndex<=daughter1->GetLastDaughter();TrackIndex++){
1177           TParticle *tmpDaughter = fStack->Particle(TrackIndex);
1178           if(tmpDaughter->GetUniqueID() == 5){
1179             if(tmpDaughter->GetPdgCode() == 11){
1180               if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){
1181                 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()  < tmpDaughter->R() ){
1182                   if(tmpDaughter->R()< fV0Reader->GetMaxRCut()){
1183                     daughter1Electron = kTRUE;
1184                   }
1185                 }
1186                                                                 
1187               }
1188             }
1189             else if(tmpDaughter->GetPdgCode() == -11){
1190               if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){
1191                 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()  < tmpDaughter->R() ){
1192                   if(tmpDaughter->R()< fV0Reader->GetMaxRCut()){
1193                     daughter1Positron = kTRUE;
1194                   }
1195                 }
1196                                                                 
1197               }
1198                                                         
1199             }
1200           }
1201         }
1202       }
1203                         
1204                         
1205       if(particle->GetPdgCode()==111){     //Pi0
1206         if( iTracks >= fStack->GetNprimary()){
1207           fHistograms->FillHistogram("MC_Pi0_Secondaries_Eta", particle->Eta());
1208           fHistograms->FillHistogram("MC_Pi0_Secondaries_Rapid", rapidity);
1209           fHistograms->FillHistogram("MC_Pi0_Secondaries_Phi", tmpPhi);
1210           fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt", particle->Pt());
1211           fHistograms->FillHistogram("MC_Pi0_Secondaries_Energy", particle->Energy());
1212           fHistograms->FillHistogram("MC_Pi0_Secondaries_R", particle->R());
1213           fHistograms->FillHistogram("MC_Pi0_Secondaries_ZR", particle->Vz(),particle->R());
1214           fHistograms->FillHistogram("MC_Pi0_Secondaries_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));
1215           fHistograms->FillHistogram("MC_Pi0_Secondaries_XY", particle->Vx(),particle->Vy());//only fill from one daughter to avoid multiple filling
1216                                         
1217           if(gammaEtaCut && gammaRCut){  
1218             //if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){
1219             fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());
1220             fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);
1221             if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){
1222               fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());
1223               fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);
1224             }
1225           }
1226         }
1227         else{
1228           fHistograms->FillHistogram("MC_Pi0_Eta", particle->Eta());    
1229           fHistograms->FillHistogram("MC_Pi0_Rapid", rapidity);
1230           fHistograms->FillHistogram("MC_Pi0_Phi", tmpPhi);
1231           fHistograms->FillHistogram("MC_Pi0_Pt", particle->Pt());
1232           fHistograms->FillHistogram("MC_Pi0_Pt_vs_Rapid", particle->Pt(),rapidity);
1233           fHistograms->FillHistogram("MC_Pi0_Energy", particle->Energy());
1234           fHistograms->FillHistogram("MC_Pi0_R", particle->R());
1235           fHistograms->FillHistogram("MC_Pi0_ZR", particle->Vz(),particle->R());
1236           fHistograms->FillHistogram("MC_Pi0_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));
1237           fHistograms->FillHistogram("MC_Pi0_XY", particle->Vx(), particle->Vy());//only fill from one daughter to avoid multiple filling
1238           if(TMath::Abs(particle->Eta())<0.9)fHistograms->FillHistogram("MC_Pi0_Pt_Fiducial", particle->Pt());
1239
1240           switch(GetProcessType(fGCMCEvent)){
1241           case  kProcSD:
1242             fHistograms->FillHistogram("MC_SD_EvtQ3_Pi0_Pt", particle->Pt());
1243             break;
1244           case  kProcDD:
1245             fHistograms->FillHistogram("MC_DD_EvtQ3_Pi0_Pt", particle->Pt());
1246             break;
1247           case  kProcND:
1248             fHistograms->FillHistogram("MC_ND_EvtQ3_Pi0_Pt", particle->Pt());
1249             break;
1250           default:
1251             AliError("Unknown Process");
1252           }
1253                         
1254           if(gammaEtaCut && gammaRCut){
1255             //    if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){
1256             fHistograms->FillHistogram("MC_Pi0_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());
1257             fHistograms->FillHistogram("MC_Pi0_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);
1258             if(TMath::Abs(particle->Eta())<0.9)fHistograms->FillHistogram("MC_Pi0_Pt_withinAcceptance_Fiducial", particle->Pt());
1259
1260             if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){
1261               fHistograms->FillHistogram("MC_Pi0_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());
1262               fHistograms->FillHistogram("MC_Pi0_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);
1263               fHistograms->FillHistogram("MC_Pi0_ZR_ConvGamma_withinAcceptance", particle->Vz(),particle->R());
1264               fHistograms->FillHistogram("MC_Pi0_ConvGamma_OpeningAngle_Pt", particle->Pt(),GetMCOpeningAngle(daughter0,daughter1));
1265               fHistograms->FillHistogram("MC_Pi0_ConvGamma_PtGamma_Pt", particle->Pt(),daughter0->Pt());
1266               fHistograms->FillHistogram("MC_Pi0_ConvGamma_PtGamma_Pt", particle->Pt(),daughter1->Pt());
1267
1268               Double_t alfa=0.;
1269               if((daughter0->Energy()+daughter1->Energy()) > 0.){
1270                 alfa= TMath::Abs((daughter0->Energy()-daughter1->Energy())/(daughter0->Energy()+daughter1->Energy()));
1271               }
1272               fHistograms->FillHistogram("MC_Pi0_alpha",alfa);
1273               if(TMath::Abs(particle->Eta())<0.9)fHistograms->FillHistogram("MC_Pi0_Pt_ConvGamma_withinAcceptance_Fiducial", particle->Pt());
1274
1275             }
1276           }
1277         }
1278       }
1279                         
1280       if(particle->GetPdgCode()==221){   //Eta
1281         fHistograms->FillHistogram("MC_Eta_Eta", particle->Eta());
1282         fHistograms->FillHistogram("MC_Eta_Rapid", rapidity);
1283         fHistograms->FillHistogram("MC_Eta_Phi",tmpPhi);
1284         fHistograms->FillHistogram("MC_Eta_Pt", particle->Pt());
1285         fHistograms->FillHistogram("MC_Eta_Pt_vs_Rapid", particle->Pt(),rapidity);
1286         fHistograms->FillHistogram("MC_Eta_Energy", particle->Energy());
1287         fHistograms->FillHistogram("MC_Eta_R", particle->R());
1288         fHistograms->FillHistogram("MC_Eta_ZR", particle->Vz(),particle->R());
1289         fHistograms->FillHistogram("MC_Eta_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));
1290         fHistograms->FillHistogram("MC_Eta_XY", particle->Vx(), particle->Vy());//only fill from one daughter to avoid multiple filling
1291                                 
1292         if(gammaEtaCut && gammaRCut){  
1293           //    if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){
1294           fHistograms->FillHistogram("MC_Eta_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());
1295           fHistograms->FillHistogram("MC_Eta_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);
1296           if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){
1297             fHistograms->FillHistogram("MC_Eta_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());
1298             fHistograms->FillHistogram("MC_Eta_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);
1299             fHistograms->FillHistogram("MC_Eta_ZR_ConvGamma_withinAcceptance", particle->Vz(),particle->R());
1300             fHistograms->FillHistogram("MC_Eta_ConvGamma_OpeningAngle_Pt", particle->Pt(),GetMCOpeningAngle(daughter0,daughter1));
1301             fHistograms->FillHistogram("MC_Eta_ConvGamma_PtGamma_Pt", particle->Pt(),daughter0->Pt());
1302             fHistograms->FillHistogram("MC_Eta_ConvGamma_PtGamma_Pt", particle->Pt(),daughter1->Pt());
1303
1304             Double_t alfa=0.;
1305             if((daughter0->Energy()+daughter1->Energy()) > 0.){
1306               alfa= TMath::Abs((daughter0->Energy()-daughter1->Energy())/(daughter0->Energy()+daughter1->Energy()));
1307               }
1308             fHistograms->FillHistogram("MC_Eta_alpha",alfa);
1309
1310           }
1311                                         
1312         }
1313                                 
1314       }
1315                         
1316       // all motherparticles with 2 gammas as daughters
1317       fHistograms->FillHistogram("MC_Mother_R", particle->R());
1318       fHistograms->FillHistogram("MC_Mother_ZR", particle->Vz(),particle->R());
1319       fHistograms->FillHistogram("MC_Mother_XY", particle->Vx(),particle->Vy());
1320       fHistograms->FillHistogram("MC_Mother_Mass", particle->GetCalcMass());
1321       fHistograms->FillHistogram("MC_Mother_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));
1322       fHistograms->FillHistogram("MC_Mother_Energy", particle->Energy());
1323       fHistograms->FillHistogram("MC_Mother_Pt", particle->Pt());
1324       fHistograms->FillHistogram("MC_Mother_Eta", particle->Eta());
1325       fHistograms->FillHistogram("MC_Mother_Rapid", rapidity);
1326       fHistograms->FillHistogram("MC_Mother_Phi",tmpPhi);
1327       fHistograms->FillHistogram("MC_Mother_InvMass_vs_Pt",particle->GetMass(),particle->Pt());                 
1328                         
1329       if(gammaEtaCut && gammaRCut){  
1330         //      if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){
1331         fHistograms->FillHistogram("MC_Mother_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());
1332         fHistograms->FillHistogram("MC_Mother_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);
1333         fHistograms->FillHistogram("MC_Mother_InvMass_vs_Pt_withinAcceptance",particle->GetMass(),particle->Pt());                      
1334         if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){
1335           fHistograms->FillHistogram("MC_Mother_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());
1336           fHistograms->FillHistogram("MC_Mother_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);
1337           fHistograms->FillHistogram("MC_Mother_InvMass_vs_Pt_ConvGamma_withinAcceptance",particle->GetMass(),particle->Pt());                  
1338                                         
1339         }
1340                                 
1341                                 
1342       }  // end passed R and eta cut
1343                         
1344     } // end if(particle->GetNDaughters() == 2)
1345                 
1346   }// end for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++)
1347
1348 } // end ProcessMCData
1349
1350
1351
1352 void AliAnalysisTaskGammaConversion::FillNtuple(){
1353   //Fills the ntuple with the different values
1354         
1355   if(fGammaNtuple == NULL){
1356     return;
1357   }
1358   Int_t numberOfV0s = fV0Reader->GetNumberOfV0s();
1359   for(Int_t i=0;i<numberOfV0s;i++){
1360     Float_t values[27] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
1361     AliESDv0 * cV0 = fV0Reader->GetV0(i);
1362     Double_t negPID=0;
1363     Double_t posPID=0;
1364     fV0Reader->GetPIDProbability(negPID,posPID);
1365     values[0]=cV0->GetOnFlyStatus();
1366     values[1]=fV0Reader->CheckForPrimaryVertex();
1367     values[2]=negPID;
1368     values[3]=posPID;
1369     values[4]=fV0Reader->GetX();
1370     values[5]=fV0Reader->GetY();
1371     values[6]=fV0Reader->GetZ();
1372     values[7]=fV0Reader->GetXYRadius();
1373     values[8]=fV0Reader->GetMotherCandidateNDF();
1374     values[9]=fV0Reader->GetMotherCandidateChi2();
1375     values[10]=fV0Reader->GetMotherCandidateEnergy();
1376     values[11]=fV0Reader->GetMotherCandidateEta();
1377     values[12]=fV0Reader->GetMotherCandidatePt();
1378     values[13]=fV0Reader->GetMotherCandidateMass();
1379     values[14]=fV0Reader->GetMotherCandidateWidth();
1380     //    values[15]=fV0Reader->GetMotherMCParticle()->Pt();   MOVED TO THE END, HAS TO BE CALLED AFTER HasSameMother NB: still has the same entry in the array
1381     values[16]=fV0Reader->GetOpeningAngle();
1382     values[17]=fV0Reader->GetNegativeTrackEnergy();
1383     values[18]=fV0Reader->GetNegativeTrackPt();
1384     values[19]=fV0Reader->GetNegativeTrackEta();
1385     values[20]=fV0Reader->GetNegativeTrackPhi();
1386     values[21]=fV0Reader->GetPositiveTrackEnergy();
1387     values[22]=fV0Reader->GetPositiveTrackPt();
1388     values[23]=fV0Reader->GetPositiveTrackEta();
1389     values[24]=fV0Reader->GetPositiveTrackPhi();
1390     values[25]=fV0Reader->HasSameMCMother();
1391     if(values[25] != 0){
1392       values[26]=fV0Reader->GetMotherMCParticlePDGCode();
1393       values[15]=fV0Reader->GetMotherMCParticle()->Pt();
1394     }
1395     fTotalNumberOfAddedNtupleEntries++;
1396     fGammaNtuple->Fill(values);
1397   }
1398   fV0Reader->ResetV0IndexNumber();
1399         
1400 }
1401
1402 void AliAnalysisTaskGammaConversion::ProcessV0sNoCut(){
1403   // Process all the V0's without applying any cuts to it
1404         
1405   Int_t numberOfV0s = fV0Reader->GetNumberOfV0s();
1406   for(Int_t i=0;i<numberOfV0s;i++){
1407     /*AliESDv0 * cV0 = */fV0Reader->GetV0(i);
1408                 
1409     if(fV0Reader->CheckForPrimaryVertex() == kFALSE){
1410       continue;
1411     }
1412
1413     //    if( !fV0Reader->GetV0(i)->GetOnFlyStatus()){
1414     if( !fV0Reader->CheckV0FinderStatus(i)){
1415       continue;
1416     }
1417
1418
1419     if( !((fV0Reader->GetNegativeESDTrack())->GetStatus() & AliESDtrack::kTPCrefit) || 
1420         !((fV0Reader->GetPositiveESDTrack())->GetStatus() & AliESDtrack::kTPCrefit) ){
1421       continue;
1422     }
1423
1424     if( fV0Reader->GetNegativeESDTrack()->GetSign()== fV0Reader->GetPositiveESDTrack()->GetSign()){
1425       continue;
1426     }
1427
1428     if( fV0Reader->GetNegativeESDTrack()->GetKinkIndex(0) > 0 || 
1429         fV0Reader->GetPositiveESDTrack()->GetKinkIndex(0) > 0) {                        
1430       continue;
1431     }
1432     if(TMath::Abs(fV0Reader->GetMotherCandidateEta())> fV0Reader->GetEtaCut()){
1433       continue;
1434     }
1435     if(TMath::Abs(fV0Reader->GetPositiveTrackEta())> fV0Reader->GetEtaCut()){
1436       continue;
1437     }
1438     if(TMath::Abs(fV0Reader->GetNegativeTrackEta())> fV0Reader->GetEtaCut()){
1439       continue;
1440     }
1441     if((TMath::Abs(fV0Reader->GetZ())*fV0Reader->GetLineCutZRSlope())-fV0Reader->GetLineCutZValue() > fV0Reader->GetXYRadius() ){ // cuts out regions where we do not reconstruct
1442       continue; 
1443     }
1444     if(fDoMCTruth){
1445                         
1446       if(fV0Reader->HasSameMCMother() == kFALSE){
1447         continue;
1448       }
1449                         
1450       TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
1451       TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
1452                         
1453       if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){
1454         continue;
1455       }
1456       if(negativeMC->GetPdgCode() == positiveMC->GetPdgCode()){
1457         continue;
1458       }
1459
1460       if(negativeMC->GetUniqueID() != 5 || positiveMC->GetUniqueID() !=5){ // id 5 is conversion
1461         continue;
1462       }
1463                         
1464       if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
1465                                 
1466         fHistograms->FillHistogram("ESD_NoCutConvGamma_Pt", fV0Reader->GetMotherCandidatePt());
1467         fHistograms->FillHistogram("ESD_NoCutConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());
1468         fHistograms->FillHistogram("ESD_NoCutConvGamma_Eta", fV0Reader->GetMotherCandidateEta());                               
1469         fHistograms->FillHistogram("ESD_NoCutConvGamma_Phi", fV0Reader->GetMotherCandidatePhi());
1470         fHistograms->FillHistogram("ESD_NoCutConvGamma_Mass", fV0Reader->GetMotherCandidateMass());
1471         fHistograms->FillHistogram("ESD_NoCutConvGamma_Width", fV0Reader->GetMotherCandidateWidth());
1472         fHistograms->FillHistogram("ESD_NoCutConvGamma_Chi2", fV0Reader->GetMotherCandidateChi2());
1473         fHistograms->FillHistogram("ESD_NoCutConvGamma_NDF", fV0Reader->GetMotherCandidateNDF());
1474         fHistograms->FillHistogram("ESD_NoCutConvGamma_Rapid", fV0Reader->GetMotherCandidateRapidity());
1475         fHistograms->FillHistogram("ESD_NoCutConvGamma_Pt_Eta", fV0Reader->GetMotherCandidatePt(),fV0Reader->GetMotherCandidateEta());
1476                                 
1477         fHistograms->FillHistogram("ESD_NoCutConvGamma_Pt_Chi2", fV0Reader->GetMotherCandidatePt(), fV0Reader->GetMotherCandidateChi2());
1478         fHistograms->FillHistogram("ESD_NoCutConvGamma_Eta_Chi2", fV0Reader->GetMotherCandidateEta(), fV0Reader->GetMotherCandidateChi2());
1479                                 
1480         fHistograms->FillHistogram("ESD_NoCutConversion_XY", fV0Reader->GetX(),fV0Reader->GetY());
1481         fHistograms->FillHistogram("ESD_NoCutConversion_R", fV0Reader->GetXYRadius());
1482         fHistograms->FillHistogram("ESD_NoCutConversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());
1483         fHistograms->FillHistogram("ESD_NoCutConversion_OpeningAngle", fV0Reader->GetOpeningAngle());
1484         fHistograms->FillHistogram("ESD_NoCutConvGamma_CosPointingAngle", fV0Reader->GetCosPointingAngle());
1485         fHistograms->FillHistogram("ESD_NoCutConvGamma_DcaDaughters", fV0Reader->GetDcaDaughters());
1486         fHistograms->FillHistogram("ESD_NoCutConvGamma_NormDcaDistDaughters", fV0Reader->GetNormDcaDistDaughters());
1487         fHistograms->FillHistogram("ESD_NoCutConvGamma_LikelihoodAP", fV0Reader->GetLikelihoodAP());
1488
1489         fHistograms->FillHistogram("ESD_NoCutConvGamma_E_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetNegativeTrackP()/fV0Reader->GetMotherCandidateP());
1490         fHistograms->FillHistogram("ESD_NoCutConvGamma_P_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetPositiveTrackP()/fV0Reader->GetMotherCandidateP());
1491         fHistograms->FillHistogram("ESD_NoCutConvGamma_E_dEdxP",fV0Reader->GetNegativeTrackP(),fV0Reader->GetNegativeTrackTPCdEdx());
1492         fHistograms->FillHistogram("ESD_NoCutConvGamma_P_dEdxP",fV0Reader->GetPositiveTrackP(),fV0Reader->GetPositiveTrackTPCdEdx());
1493
1494         //store MCTruth properties
1495         fHistograms->FillHistogram("ESD_NoCutConvGamma_MC_Pt_Eta", fV0Reader->GetMotherMCParticle()->Pt(),fV0Reader->GetMotherMCParticle()->Eta());
1496         fHistograms->FillHistogram("ESD_NoCutConversion_MC_ZR", negativeMC->Vz(),negativeMC->R());
1497         fHistograms->FillHistogram("ESD_NoCutConversion_MC_XY", negativeMC->Vx(),negativeMC->Vy());
1498       }
1499     }
1500   }
1501   fV0Reader->ResetV0IndexNumber();
1502 }
1503
1504 void AliAnalysisTaskGammaConversion::ProcessV0s(){
1505   // see header file for documentation
1506
1507
1508   if(fWriteNtuple == kTRUE){
1509     FillNtuple();
1510   }
1511         
1512   Int_t nSurvivingV0s=0;
1513   fV0Reader->ResetNGoodV0s();
1514   while(fV0Reader->NextV0()){
1515     nSurvivingV0s++;
1516                 
1517
1518     TVector3 vtxConv(fV0Reader->GetX(),fV0Reader->GetY(), fV0Reader->GetZ());
1519         
1520     //-------------------------- filling v0 information -------------------------------------
1521     fHistograms->FillHistogram("ESD_Conversion_R", fV0Reader->GetXYRadius());
1522     fHistograms->FillHistogram("ESD_Conversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());
1523     fHistograms->FillHistogram("ESD_Conversion_XY", fV0Reader->GetX(),fV0Reader->GetY());
1524     fHistograms->FillHistogram("ESD_Conversion_OpeningAngle", fV0Reader->GetOpeningAngle());    
1525
1526     // Specific histograms for beam pipe studies
1527     if( TMath::Abs(fV0Reader->GetZ()) < fV0Reader->GetLineCutZValue() ){
1528       fHistograms->FillHistogram("ESD_Conversion_XY_BeamPipe", fV0Reader->GetX(),fV0Reader->GetY());
1529       fHistograms->FillHistogram("ESD_Conversion_RPhi_BeamPipe", vtxConv.Phi(),fV0Reader->GetXYRadius());
1530     }
1531
1532         
1533     fHistograms->FillHistogram("ESD_E_Energy", fV0Reader->GetNegativeTrackEnergy());
1534     fHistograms->FillHistogram("ESD_E_Pt", fV0Reader->GetNegativeTrackPt());
1535     fHistograms->FillHistogram("ESD_E_Eta", fV0Reader->GetNegativeTrackEta());
1536     fHistograms->FillHistogram("ESD_E_Phi", fV0Reader->GetNegativeTrackPhi());
1537     fHistograms->FillHistogram("ESD_E_nTPCClusters", fV0Reader->GetNegativeTracknTPCClusters());
1538     fHistograms->FillHistogram("ESD_E_nITSClusters", fV0Reader->GetNegativeTracknITSClusters());
1539     if(fV0Reader->GetNegativeTracknTPCFClusters()!=0 && fV0Reader->GetNegativeTracknTPCClusters()!=0 ){
1540       Double_t eClsToF= (Double_t)fV0Reader->GetNegativeTracknTPCClusters()/(Double_t)fV0Reader->GetNegativeTracknTPCFClusters();
1541       fHistograms->FillHistogram("ESD_E_nTPCClustersToFP", fV0Reader->GetNegativeTrackP(),eClsToF );
1542       fHistograms->FillHistogram("ESD_E_TPCchi2", fV0Reader->GetNegativeTrackTPCchi2()/(Double_t)fV0Reader->GetNegativeTracknTPCClusters());
1543     }
1544
1545
1546
1547     fHistograms->FillHistogram("ESD_P_Energy", fV0Reader->GetPositiveTrackEnergy());
1548     fHistograms->FillHistogram("ESD_P_Pt", fV0Reader->GetPositiveTrackPt());
1549     fHistograms->FillHistogram("ESD_P_Eta", fV0Reader->GetPositiveTrackEta());
1550     fHistograms->FillHistogram("ESD_P_Phi", fV0Reader->GetPositiveTrackPhi());
1551     fHistograms->FillHistogram("ESD_P_nTPCClusters", fV0Reader->GetPositiveTracknTPCClusters());
1552     fHistograms->FillHistogram("ESD_P_nITSClusters", fV0Reader->GetPositiveTracknITSClusters());
1553     if(fV0Reader->GetPositiveTracknTPCFClusters()!=0 && (Double_t)fV0Reader->GetPositiveTracknTPCClusters()!=0 ){
1554       Double_t pClsToF= (Double_t)fV0Reader->GetPositiveTracknTPCClusters()/(Double_t)fV0Reader->GetPositiveTracknTPCFClusters();
1555       fHistograms->FillHistogram("ESD_P_nTPCClustersToFP",fV0Reader->GetPositiveTrackP(), pClsToF);
1556       fHistograms->FillHistogram("ESD_P_TPCchi2", fV0Reader->GetPositiveTrackTPCchi2()/(Double_t)fV0Reader->GetPositiveTracknTPCClusters());
1557     }
1558
1559
1560     fHistograms->FillHistogram("ESD_ConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());
1561     fHistograms->FillHistogram("ESD_ConvGamma_Pt", fV0Reader->GetMotherCandidatePt());
1562     fHistograms->FillHistogram("ESD_ConvGamma_Eta", fV0Reader->GetMotherCandidateEta());
1563     fHistograms->FillHistogram("ESD_ConvGamma_Phi", fV0Reader->GetMotherCandidatePhi());
1564     fHistograms->FillHistogram("ESD_ConvGamma_Mass", fV0Reader->GetMotherCandidateMass());
1565     fHistograms->FillHistogram("ESD_ConvGamma_Width", fV0Reader->GetMotherCandidateWidth());
1566     fHistograms->FillHistogram("ESD_ConvGamma_Chi2", fV0Reader->GetMotherCandidateChi2());
1567     fHistograms->FillHistogram("ESD_ConvGamma_NDF", fV0Reader->GetMotherCandidateNDF());
1568     fHistograms->FillHistogram("ESD_ConvGamma_Rapid", fV0Reader->GetMotherCandidateRapidity());
1569     fHistograms->FillHistogram("ESD_ConvGamma_Pt_Eta", fV0Reader->GetMotherCandidatePt(),fV0Reader->GetMotherCandidateEta());
1570                 
1571     fHistograms->FillHistogram("ESD_ConvGamma_Pt_Chi2", fV0Reader->GetMotherCandidatePt(), fV0Reader->GetMotherCandidateChi2());
1572     fHistograms->FillHistogram("ESD_ConvGamma_Eta_Chi2", fV0Reader->GetMotherCandidateEta(), fV0Reader->GetMotherCandidateChi2());
1573     
1574     fHistograms->FillHistogram("ESD_ConvGamma_CosPointingAngle", fV0Reader->GetCosPointingAngle());
1575     fHistograms->FillHistogram("ESD_ConvGamma_DcaDaughters", fV0Reader->GetDcaDaughters());
1576     fHistograms->FillHistogram("ESD_ConvGamma_NormDcaDistDaughters", fV0Reader->GetNormDcaDistDaughters());
1577     fHistograms->FillHistogram("ESD_ConvGamma_LikelihoodAP", fV0Reader->GetLikelihoodAP());
1578     
1579     fHistograms->FillHistogram("ESD_ConvGamma_E_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetNegativeTrackP()/fV0Reader->GetMotherCandidateP());
1580     fHistograms->FillHistogram("ESD_ConvGamma_P_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetPositiveTrackP()/fV0Reader->GetMotherCandidateP());
1581     fHistograms->FillHistogram("ESD_ConvGamma_E_dEdxP",fV0Reader->GetNegativeTrackP(),fV0Reader->GetNegativeTrackTPCdEdx());
1582     fHistograms->FillHistogram("ESD_ConvGamma_P_dEdxP",fV0Reader->GetPositiveTrackP(),fV0Reader->GetPositiveTrackTPCdEdx());
1583     
1584     Double_t negPID=0;
1585     Double_t posPID=0;
1586     fV0Reader->GetPIDProbability(negPID,posPID);
1587     fHistograms->FillHistogram("ESD_ConvGamma_E_EProbP",fV0Reader->GetNegativeTrackP(),negPID);
1588     fHistograms->FillHistogram("ESD_ConvGamma_P_EProbP",fV0Reader->GetPositiveTrackP(),posPID);
1589
1590     Double_t negPIDmupi=0;
1591     Double_t posPIDmupi=0;
1592     fV0Reader->GetPIDProbabilityMuonPion(negPIDmupi,posPIDmupi);
1593     fHistograms->FillHistogram("ESD_ConvGamma_E_mupiProbP",fV0Reader->GetNegativeTrackP(),negPIDmupi);
1594     fHistograms->FillHistogram("ESD_ConvGamma_P_mupiProbP",fV0Reader->GetPositiveTrackP(),posPIDmupi);
1595
1596     Double_t armenterosQtAlfa[2];
1597     fV0Reader->GetArmenterosQtAlfa(fV0Reader-> GetNegativeKFParticle(), 
1598                                    fV0Reader-> GetPositiveKFParticle(), 
1599                                    fV0Reader->GetMotherCandidateKFCombination(),
1600                                    armenterosQtAlfa);
1601    
1602     fHistograms->FillHistogram("ESD_ConvGamma_alfa_qt",armenterosQtAlfa[1],armenterosQtAlfa[0]);
1603  
1604
1605     // begin mapping
1606     Int_t rBin    = fHistograms->GetRBin(fV0Reader->GetXYRadius());
1607     Int_t zBin    = fHistograms->GetZBin(fV0Reader->GetZ());
1608     Int_t phiBin  = fHistograms->GetPhiBin(fV0Reader->GetNegativeTrackPhi());
1609     Double_t rFMD=30;
1610     Double_t rITSTPCMin=50;
1611     Double_t rITSTPCMax=80;
1612
1613
1614     //    Double_t motherCandidateEta= fV0Reader->GetMotherCandidateEta();
1615                 
1616     TString nameESDMappingPhiR="";
1617     nameESDMappingPhiR.Form("ESD_Conversion_Mapping_Phi%02d_R%02d",phiBin,rBin);
1618     //fHistograms->FillHistogram(nameESDMappingPhiR, fV0Reader->GetZ(), motherCandidateEta);
1619                 
1620     TString nameESDMappingPhi="";
1621     nameESDMappingPhi.Form("ESD_Conversion_Mapping_Phi%02d",phiBin);
1622     //fHistograms->FillHistogram(nameESDMappingPhi, fV0Reader->GetZ(), motherCandidateEta);
1623                 
1624     TString nameESDMappingR="";
1625     nameESDMappingR.Form("ESD_Conversion_Mapping_R%02d",rBin);
1626     //fHistograms->FillHistogram(nameESDMappingR, fV0Reader->GetZ(), motherCandidateEta);  
1627                 
1628     TString nameESDMappingPhiInR="";
1629     nameESDMappingPhiInR.Form("ESD_Conversion_Mapping_Phi_in_R_%02d",rBin);
1630     //    fHistograms->FillHistogram(nameESDMappingPhiInR, fV0Reader->GetMotherCandidatePhi());
1631     fHistograms->FillHistogram(nameESDMappingPhiInR, vtxConv.Phi());
1632
1633     TString nameESDMappingZInR="";
1634     nameESDMappingZInR.Form("ESD_Conversion_Mapping_Z_in_R_%02d",rBin);
1635     fHistograms->FillHistogram(nameESDMappingZInR, fV0Reader->GetZ());
1636
1637     TString nameESDMappingPhiInZ="";
1638     nameESDMappingPhiInZ.Form("ESD_Conversion_Mapping_Phi_in_Z_%02d",zBin);
1639     //    fHistograms->FillHistogram(nameESDMappingPhiInR, fV0Reader->GetMotherCandidatePhi());
1640     fHistograms->FillHistogram(nameESDMappingPhiInZ, vtxConv.Phi());
1641
1642     if(fV0Reader->GetXYRadius()<rFMD){
1643       TString nameESDMappingFMDPhiInZ="";
1644       nameESDMappingFMDPhiInZ.Form("ESD_Conversion_Mapping_FMD_Phi_in_Z_%02d",zBin);
1645       fHistograms->FillHistogram(nameESDMappingFMDPhiInZ, vtxConv.Phi());
1646     }
1647
1648     if(fV0Reader->GetXYRadius()>rITSTPCMin && fV0Reader->GetXYRadius()<rITSTPCMax){
1649       TString nameESDMappingITSTPCPhiInZ="";
1650       nameESDMappingITSTPCPhiInZ.Form("ESD_Conversion_Mapping_ITSTPC_Phi_in_Z_%02d",zBin);
1651       fHistograms->FillHistogram(nameESDMappingITSTPCPhiInZ, vtxConv.Phi());
1652     }
1653
1654     TString nameESDMappingRInZ="";
1655     nameESDMappingRInZ.Form("ESD_Conversion_Mapping_R_in_Z_%02d",zBin);
1656     fHistograms->FillHistogram(nameESDMappingRInZ, fV0Reader->GetXYRadius());
1657
1658     if(fV0Reader->GetMotherCandidatePt() > fLowPtMapping && fV0Reader->GetMotherCandidatePt()< fHighPtMapping){
1659       TString nameESDMappingMidPtPhiInR="";
1660       nameESDMappingMidPtPhiInR.Form("ESD_Conversion_Mapping_MidPt_Phi_in_R_%02d",rBin);
1661       fHistograms->FillHistogram(nameESDMappingMidPtPhiInR, vtxConv.Phi());
1662       
1663       TString nameESDMappingMidPtZInR="";
1664       nameESDMappingMidPtZInR.Form("ESD_Conversion_Mapping_MidPt_Z_in_R_%02d",rBin);
1665       fHistograms->FillHistogram(nameESDMappingMidPtZInR, fV0Reader->GetZ());
1666       
1667       TString nameESDMappingMidPtPhiInZ="";
1668       nameESDMappingMidPtPhiInZ.Form("ESD_Conversion_Mapping_MidPt_Phi_in_Z_%02d",zBin);
1669       fHistograms->FillHistogram(nameESDMappingMidPtPhiInZ, vtxConv.Phi());
1670       if(fV0Reader->GetXYRadius()<rFMD){
1671         TString nameESDMappingMidPtFMDPhiInZ="";
1672         nameESDMappingMidPtFMDPhiInZ.Form("ESD_Conversion_Mapping_MidPt_FMD_Phi_in_Z_%02d",zBin);
1673         fHistograms->FillHistogram(nameESDMappingMidPtFMDPhiInZ, vtxConv.Phi());
1674       }
1675
1676       
1677       TString nameESDMappingMidPtRInZ="";
1678       nameESDMappingMidPtRInZ.Form("ESD_Conversion_Mapping_MidPt_R_in_Z_%02d",zBin);
1679       fHistograms->FillHistogram(nameESDMappingMidPtRInZ, fV0Reader->GetXYRadius());
1680
1681     }
1682
1683
1684     // end mapping
1685                 
1686     new((*fKFReconstructedGammasTClone)[fKFReconstructedGammasTClone->GetEntriesFast()])  AliKFParticle(*fV0Reader->GetMotherCandidateKFCombination());
1687     fKFReconstructedGammasV0Index.push_back(fV0Reader->GetCurrentV0IndexNumber()-1);
1688     //    fKFReconstructedGammas.push_back(*fV0Reader->GetMotherCandidateKFCombination());
1689     fElectronv1.push_back(fV0Reader->GetCurrentV0()->GetPindex());
1690     fElectronv2.push_back(fV0Reader->GetCurrentV0()->GetNindex());
1691                 
1692     //----------------------------------- checking for "real" conversions (MC match) --------------------------------------
1693     if(fDoMCTruth){
1694       TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
1695       TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
1696                         
1697       if(fV0Reader->HasSameMCMother() == kFALSE){
1698         fHistograms->FillHistogram("ESD_TrueConvCombinatorial_R", fV0Reader->GetXYRadius());
1699         fHistograms->FillHistogram("ESD_TrueConvCombinatorial_Pt", fV0Reader->GetMotherCandidatePt());
1700         if(TMath::Abs(negativeMC->GetPdgCode())==11 && TMath::Abs(positiveMC->GetPdgCode())==11){
1701           fHistograms->FillHistogram("ESD_TrueConvCombinatorialElec_R", fV0Reader->GetXYRadius());
1702           fHistograms->FillHistogram("ESD_TrueConvCombinatorialElec_Pt", fV0Reader->GetMotherCandidatePt());
1703         }
1704         continue;
1705       }
1706       // Moved up to check true electron background
1707       //      TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
1708       //      TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
1709
1710       if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){
1711         fHistograms->FillHistogram("ESD_TrueConvHadronicBck_R", fV0Reader->GetXYRadius());
1712         fHistograms->FillHistogram("ESD_TrueConvHadronicBck_Pt", fV0Reader->GetMotherCandidatePt());
1713         continue;
1714       }
1715
1716       if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){
1717         continue;
1718       }
1719
1720       UInt_t statusPos = fV0Reader->GetPositiveESDTrack()->GetStatus(); 
1721       UInt_t statusNeg = fV0Reader->GetNegativeESDTrack()->GetStatus(); 
1722       UChar_t itsPixelPos = fV0Reader->GetPositiveESDTrack()->GetITSClusterMap();
1723       UChar_t itsPixelNeg = fV0Reader->GetNegativeESDTrack()->GetITSClusterMap();
1724
1725       // Using the UniqueID Phojet does not get the Dalitz right
1726       //      if( (negativeMC->GetUniqueID() == 4 && positiveMC->GetUniqueID() ==4) ||
1727       //  (negativeMC->GetUniqueID() == 0 && positiveMC->GetUniqueID() ==0) ){// fill r distribution for Dalitz decays 
1728       if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 111){ //pi0
1729         fHistograms->FillHistogram("ESD_TrueDalitzContamination_R", fV0Reader->GetXYRadius());
1730         fHistograms->FillHistogram("ESD_TrueConvDalitzPi0_R", fV0Reader->GetXYRadius());
1731         //--------Histos for HFE 
1732
1733         if(statusPos & AliESDtrack::kTOFpid){
1734           fHistograms->FillHistogram("ESD_TrueConvDalitzPi0_SinglePos_R", fV0Reader->GetXYRadius());
1735           if( TESTBIT(itsPixelPos, 0) ){ 
1736             fHistograms->FillHistogram("ESD_TrueConvDalitzPi0_SinglePos_kFirst_R", fV0Reader->GetXYRadius());
1737           }
1738         }
1739         if(statusNeg & AliESDtrack::kTOFpid){
1740           fHistograms->FillHistogram("ESD_TrueConvDalitzPi0_SingleNeg_R", fV0Reader->GetXYRadius());
1741           if( TESTBIT(itsPixelNeg, 0) ){ 
1742             fHistograms->FillHistogram("ESD_TrueConvDalitzPi0_SingleNeg_kFirst_R", fV0Reader->GetXYRadius());
1743           }
1744         }
1745         //--------------------------------------------------------
1746
1747       }
1748       if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 221){ //eta
1749         fHistograms->FillHistogram("ESD_TrueConvDalitzEta_R", fV0Reader->GetXYRadius());
1750       }
1751
1752       //}
1753
1754       if(negativeMC->GetUniqueID() != 5 || positiveMC->GetUniqueID() !=5){// check if the daughters come from a conversion
1755         continue;
1756       }
1757                         
1758       if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
1759
1760         if(fDoCF){
1761           Double_t containerInput[3];
1762           containerInput[0] = fV0Reader->GetMotherCandidatePt();
1763           containerInput[1] = fV0Reader->GetMotherCandidateEta();
1764           containerInput[2] = fV0Reader->GetMotherCandidateMass();
1765           fCFManager->GetParticleContainer()->Fill(containerInput,kStepTrueGamma); // for CF 
1766         }
1767
1768         fHistograms->FillHistogram("ESD_TrueConvGamma_Pt", fV0Reader->GetMotherCandidatePt());
1769         fHistograms->FillHistogram("ESD_TrueConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());
1770         fHistograms->FillHistogram("ESD_TrueConvGamma_Eta", fV0Reader->GetMotherCandidateEta());                                
1771         fHistograms->FillHistogram("ESD_TrueConvGamma_Phi", fV0Reader->GetMotherCandidatePhi());
1772         fHistograms->FillHistogram("ESD_TrueConvGamma_Mass", fV0Reader->GetMotherCandidateMass());
1773         fHistograms->FillHistogram("ESD_TrueConvGamma_Width", fV0Reader->GetMotherCandidateWidth());
1774         fHistograms->FillHistogram("ESD_TrueConvGamma_Chi2", fV0Reader->GetMotherCandidateChi2());
1775         fHistograms->FillHistogram("ESD_TrueConvGamma_NDF", fV0Reader->GetMotherCandidateNDF());
1776         fHistograms->FillHistogram("ESD_TrueConvGamma_Pt_Eta", fV0Reader->GetMotherCandidatePt(),fV0Reader->GetMotherCandidateEta());
1777         fHistograms->FillHistogram("ESD_TrueConvGamma_Rapid", fV0Reader->GetMotherCandidateRapidity());
1778         fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLength", /*fV0Reader->GetNegativeTrackLength()*/fV0Reader->GetNegativeNTPCClusters());
1779         fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLength", /*fV0Reader->GetPositiveTrackLength()*/fV0Reader->GetPositiveNTPCClusters());
1780         fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLengthVSInvMass",/*fV0Reader->GetNegativeTrackLength()*/fV0Reader->GetNegativeNTPCClusters(),fV0Reader->GetMotherCandidateMass());
1781         fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLengthVSInvMass",/*fV0Reader->GetPositiveTrackLength()*/fV0Reader->GetPositiveNTPCClusters(),fV0Reader->GetMotherCandidateMass());
1782                                 
1783         fHistograms->FillHistogram("ESD_TrueConvGamma_Pt_Chi2", fV0Reader->GetMotherCandidatePt(), fV0Reader->GetMotherCandidateChi2());
1784         fHistograms->FillHistogram("ESD_TrueConvGamma_Eta_Chi2", fV0Reader->GetMotherCandidateEta(), fV0Reader->GetMotherCandidateChi2());
1785                                 
1786                                 
1787         fHistograms->FillHistogram("ESD_TrueConversion_XY", fV0Reader->GetX(),fV0Reader->GetY());
1788         fHistograms->FillHistogram("ESD_TrueConversion_R", fV0Reader->GetXYRadius());
1789         fHistograms->FillHistogram("ESD_TrueConversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());
1790         fHistograms->FillHistogram("ESD_TrueConversion_OpeningAngle", fV0Reader->GetOpeningAngle());
1791
1792         //----Histos for HFE-------------------------------------- 
1793
1794         if(statusPos & AliESDtrack::kTOFpid){
1795           fHistograms->FillHistogram("ESD_TrueConversion_SinglePos_R", positiveMC->R(),fV0Reader->GetPositiveMCParticle()->Pt());
1796           if( TESTBIT(itsPixelPos, 0) ){ 
1797             fHistograms->FillHistogram("ESD_TrueConversion_SinglePos_kFirst_R", positiveMC->R(),fV0Reader->GetPositiveMCParticle()->Pt());
1798           }
1799         }
1800         if(statusNeg & AliESDtrack::kTOFpid){
1801           fHistograms->FillHistogram("ESD_TrueConversion_SingleNeg_R", negativeMC->R(),fV0Reader->GetNegativeMCParticle()->Pt());
1802           if( TESTBIT(itsPixelNeg, 0) ){ 
1803              fHistograms->FillHistogram("ESD_TrueConversion_SingleNeg_kFirst_R", negativeMC->R(),fV0Reader->GetNegativeMCParticle()->Pt());
1804            }
1805         }
1806         //--------------------------------------------------------
1807
1808         fHistograms->FillHistogram("ESD_TrueConvGamma_CosPointingAngle", fV0Reader->GetCosPointingAngle());
1809         fHistograms->FillHistogram("ESD_TrueConvGamma_DcaDaughters", fV0Reader->GetDcaDaughters());
1810         fHistograms->FillHistogram("ESD_TrueConvGamma_NormDcaDistDaughters", fV0Reader->GetNormDcaDistDaughters());
1811         fHistograms->FillHistogram("ESD_TrueConvGamma_LikelihoodAP", fV0Reader->GetLikelihoodAP());
1812         if (fV0Reader->GetMotherCandidateP() != 0) {
1813           fHistograms->FillHistogram("ESD_TrueConvGamma_E_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetNegativeTrackP()/fV0Reader->GetMotherCandidateP());
1814           fHistograms->FillHistogram("ESD_TrueConvGamma_P_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetPositiveTrackP()/fV0Reader->GetMotherCandidateP());
1815         } else { cout << "Error::fV0Reader->GetNegativeTrackP() == 0 !!!" << endl; }
1816         
1817         fHistograms->FillHistogram("ESD_TrueConvGamma_E_dEdxP",fV0Reader->GetNegativeTrackP(),fV0Reader->GetNegativeTrackTPCdEdx());
1818         fHistograms->FillHistogram("ESD_TrueConvGamma_P_dEdxP",fV0Reader->GetPositiveTrackP(),fV0Reader->GetPositiveTrackTPCdEdx());
1819         fHistograms->FillHistogram("ESD_TrueConvGamma_alfa_qt",armenterosQtAlfa[1],armenterosQtAlfa[0]);
1820  
1821
1822                                 
1823         //store MCTruth properties
1824         fHistograms->FillHistogram("ESD_TrueConvGamma_MC_Pt_Eta", fV0Reader->GetMotherMCParticle()->Pt(),fV0Reader->GetMotherMCParticle()->Eta());
1825         fHistograms->FillHistogram("ESD_TrueConversion_MC_ZR", negativeMC->Vz(),negativeMC->R());
1826         fHistograms->FillHistogram("ESD_TrueConversion_MC_XY", negativeMC->Vx(),negativeMC->Vy());
1827         
1828         //resolution
1829         Double_t mcpt   = fV0Reader->GetMotherMCParticle()->Pt();
1830         Double_t esdpt  = fV0Reader->GetMotherCandidatePt();
1831         Double_t resdPt = 0.;
1832         if(mcpt > 0){
1833           resdPt = ((esdpt - mcpt)/mcpt)*100.;
1834         }
1835         else if(mcpt < 0){
1836           cout<<"Pt of MC particle is negative, this will cause wrong calculation of resPt"<<endl; 
1837         }
1838                                 
1839         fHistograms->FillHistogram("Resolution_Gamma_dPt_Pt", mcpt, resdPt);
1840         fHistograms->FillHistogram("Resolution_MC_Pt", mcpt);
1841         fHistograms->FillHistogram("Resolution_ESD_Pt", esdpt);
1842         fHistograms->FillHistogram("Resolution_Gamma_dPt_Phi", fV0Reader->GetMotherCandidatePhi(), resdPt);
1843                                 
1844         Double_t resdZ = 0.;
1845         if(fV0Reader->GetNegativeMCParticle()->Vz() != 0){
1846           resdZ = ((fV0Reader->GetZ() -fV0Reader->GetNegativeMCParticle()->Vz())/fV0Reader->GetNegativeMCParticle()->Vz())*100.;
1847         }
1848         Double_t resdZAbs = 0.;
1849         resdZAbs = (fV0Reader->GetZ() -fV0Reader->GetNegativeMCParticle()->Vz());
1850
1851         fHistograms->FillHistogram("Resolution_dZAbs_VS_R", fV0Reader->GetNegativeMCParticle()->R(), resdZAbs);
1852         fHistograms->FillHistogram("Resolution_dZ", fV0Reader->GetNegativeMCParticle()->Vz(), resdZ);
1853         fHistograms->FillHistogram("Resolution_MC_Z", fV0Reader->GetNegativeMCParticle()->Vz());
1854         fHistograms->FillHistogram("Resolution_ESD_Z", fV0Reader->GetZ());
1855                 
1856         // new for dPt_Pt-histograms for Electron and Positron
1857         Double_t mcEpt = fV0Reader->GetNegativeMCParticle()->Pt();
1858         Double_t resEdPt = 0.;
1859         if (mcEpt > 0){ 
1860                 resEdPt = ((fV0Reader->GetNegativeTrackPt()-mcEpt)/mcEpt)*100.;
1861         }
1862         UInt_t statusN = fV0Reader->GetNegativeESDTrack()->GetStatus(); 
1863  //    AliESDtrack * negTrk = fV0Reader->GetNegativeESDTrack();
1864         UInt_t kTRDoutN =  (statusN & AliESDtrack::kTRDout);
1865
1866         Int_t nITSclsE= fV0Reader->GetNegativeTracknITSClusters();
1867         // filling Resolution_Pt_dPt with respect to the Number of ITS clusters for Positrons
1868          switch(nITSclsE){
1869           case 0: // 0 ITS clusters
1870                 fHistograms->FillHistogram("Resolution_E_dPt_Pt_ITS0", mcEpt, resEdPt);
1871             break;
1872           case 1:  // 1 ITS cluster
1873                 fHistograms->FillHistogram("Resolution_E_dPt_Pt_ITS1", mcEpt, resEdPt);
1874             break;
1875           case 2:  // 2 ITS clusters
1876                 fHistograms->FillHistogram("Resolution_E_dPt_Pt_ITS2", mcEpt, resEdPt);
1877             break;
1878           case 3: // 3 ITS clusters
1879                 fHistograms->FillHistogram("Resolution_E_dPt_Pt_ITS3", mcEpt, resEdPt);
1880             break;
1881           case 4: // 4 ITS clusters
1882                 fHistograms->FillHistogram("Resolution_E_dPt_Pt_ITS4", mcEpt, resEdPt);
1883             break;
1884           case 5: // 5 ITS clusters
1885                 fHistograms->FillHistogram("Resolution_E_dPt_Pt_ITS5", mcEpt, resEdPt);
1886             break;
1887           case 6: // 6 ITS clusters
1888                 fHistograms->FillHistogram("Resolution_E_dPt_Pt_ITS6", mcEpt, resEdPt);
1889             break;
1890           }
1891         //Filling histograms with respect to Electron resolution
1892         fHistograms->FillHistogram("Resolution_E_dPt_Pt", mcEpt, resEdPt);
1893         fHistograms->FillHistogram("Resolution_E_dPt_Phi", fV0Reader->GetNegativeTrackPhi(), resEdPt);
1894         if(kTRDoutN){
1895         fHistograms->FillHistogram("Resolution_E_nTRDtracklets_ESDPt", fV0Reader->GetNegativeTrackPt(), fV0Reader->GetNegativeESDTrack()->GetTRDntracklets());
1896         fHistograms->FillHistogram("Resolution_E_nTRDtracklets_MCPt", mcEpt, fV0Reader->GetNegativeESDTrack()->GetTRDntracklets());     
1897         fHistograms->FillHistogram("Resolution_E_nTRDclusters_ESDPt",fV0Reader->GetNegativeTrackPt(), fV0Reader->GetNegativeESDTrack()->GetTRDncls());
1898         fHistograms->FillHistogram("Resolution_E_nTRDclusters_MCPt",mcEpt, fV0Reader->GetNegativeESDTrack()->GetTRDncls());
1899         fHistograms->FillHistogram("Resolution_E_TRDsignal_ESDPt", fV0Reader->GetNegativeTrackPt(), fV0Reader->GetNegativeESDTrack()->GetTRDsignal());
1900         }
1901
1902         Double_t mcPpt = fV0Reader->GetPositiveMCParticle()->Pt();
1903         Double_t resPdPt = 0;
1904         if (mcPpt > 0){ 
1905                 resPdPt = ((fV0Reader->GetPositiveTrackPt()-mcPpt)/mcPpt)*100.;
1906         }
1907
1908         UInt_t statusP = fV0Reader->GetPositiveESDTrack()->GetStatus(); 
1909 //     AliESDtrack * posTr= fV0Reader->GetPositiveESDTrack();
1910         UInt_t kTRDoutP =  (statusP & AliESDtrack::kTRDout);
1911         
1912         Int_t nITSclsP = fV0Reader->GetPositiveTracknITSClusters();
1913         // filling Resolution_Pt_dPt with respect to the Number of ITS clusters for Positrons
1914          switch(nITSclsP){
1915           case 0: // 0 ITS clusters
1916                 fHistograms->FillHistogram("Resolution_P_dPt_Pt_ITS0", mcPpt, resPdPt);
1917             break;
1918           case 1:  // 1 ITS cluster
1919                 fHistograms->FillHistogram("Resolution_P_dPt_Pt_ITS1", mcPpt, resPdPt);
1920             break;
1921           case 2:  // 2 ITS clusters
1922                 fHistograms->FillHistogram("Resolution_P_dPt_Pt_ITS2", mcPpt, resPdPt);
1923             break;
1924           case 3: // 3 ITS clusters
1925                 fHistograms->FillHistogram("Resolution_P_dPt_Pt_ITS3", mcPpt, resPdPt);
1926             break;
1927           case 4: // 4 ITS clusters
1928                 fHistograms->FillHistogram("Resolution_P_dPt_Pt_ITS4", mcPpt, resPdPt);
1929             break;
1930           case 5: // 5 ITS clusters
1931                 fHistograms->FillHistogram("Resolution_P_dPt_Pt_ITS5", mcPpt, resPdPt);
1932             break;
1933           case 6: // 6 ITS clusters
1934                 fHistograms->FillHistogram("Resolution_P_dPt_Pt_ITS6", mcPpt, resPdPt);
1935             break;
1936           }
1937         //Filling histograms with respect to Positron resolution
1938         fHistograms->FillHistogram("Resolution_P_dPt_Pt", mcPpt, resPdPt);
1939         fHistograms->FillHistogram("Resolution_P_dPt_Phi", fV0Reader->GetPositiveTrackPhi(), resPdPt);
1940         if(kTRDoutP){
1941         fHistograms->FillHistogram("Resolution_P_nTRDtracklets_ESDPt", fV0Reader->GetPositiveTrackPt(), fV0Reader->GetPositiveESDTrack()->GetTRDntracklets());
1942         fHistograms->FillHistogram("Resolution_P_nTRDtracklets_MCPt", mcPpt, fV0Reader->GetPositiveESDTrack()->GetTRDntracklets());
1943         fHistograms->FillHistogram("Resolution_P_nTRDclusters_ESDPt",fV0Reader->GetPositiveTrackPt(), fV0Reader->GetPositiveESDTrack()->GetTRDncls());
1944         fHistograms->FillHistogram("Resolution_P_nTRDclusters_MCPt",mcPpt, fV0Reader->GetPositiveESDTrack()->GetTRDncls());
1945         fHistograms->FillHistogram("Resolution_P_TRDsignal_ESDPt", fV0Reader->GetPositiveTrackPt(), fV0Reader->GetPositiveESDTrack()->GetTRDsignal());
1946         }
1947
1948                 
1949         Double_t resdR = 0.;
1950         if(fV0Reader->GetNegativeMCParticle()->R() != 0){
1951           resdR = ((fV0Reader->GetXYRadius() - fV0Reader->GetNegativeMCParticle()->R())/fV0Reader->GetNegativeMCParticle()->R())*100.;
1952         }
1953         Double_t resdRAbs = 0.;
1954         resdRAbs = (fV0Reader->GetXYRadius() - fV0Reader->GetNegativeMCParticle()->R());
1955
1956         fHistograms->FillHistogram("Resolution_dRAbs_VS_R", fV0Reader->GetNegativeMCParticle()->R(), resdRAbs);
1957         fHistograms->FillHistogram("Resolution_dR", fV0Reader->GetNegativeMCParticle()->R(), resdR);
1958         fHistograms->FillHistogram("Resolution_MC_R", fV0Reader->GetNegativeMCParticle()->R());
1959         fHistograms->FillHistogram("Resolution_ESD_R", fV0Reader->GetXYRadius());
1960         fHistograms->FillHistogram("Resolution_R_dPt", fV0Reader->GetNegativeMCParticle()->R(), resdPt);
1961  
1962         Double_t resdPhiAbs=0.;
1963         resdPhiAbs=0.;
1964         resdPhiAbs= (fV0Reader->GetMotherCandidatePhi()-fV0Reader->GetNegativeMCParticle()->Phi());
1965         fHistograms->FillHistogram("Resolution_dPhiAbs_VS_R", fV0Reader->GetNegativeMCParticle()->R(), resdPhiAbs);
1966
1967       }//if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22)
1968     }//if(fDoMCTruth)
1969   }//while(fV0Reader->NextV0)
1970   fHistograms->FillHistogram("ESD_NumberOfSurvivingV0s", nSurvivingV0s);
1971   fHistograms->FillHistogram("ESD_NumberOfV0s", fV0Reader->GetNumberOfV0s());
1972   fHistograms->FillHistogram("ESD_NumberOfContributorsVtx", fV0Reader->GetNumberOfContributorsVtx());
1973
1974   fV0Reader->ResetV0IndexNumber();
1975 }
1976
1977 void AliAnalysisTaskGammaConversion::AddToAODBranch(TClonesArray * branch, AliAODPWG4Particle & particle) {
1978   //See header file for documentation
1979
1980   Int_t i = branch->GetEntriesFast();
1981   if(! (fOutputAODClassName.Contains("Correlation")) ) {
1982     new((*branch)[i])  AliAODPWG4Particle(particle);
1983   } else {
1984     new((*branch)[i])  AliAODPWG4ParticleCorrelation(particle);
1985   }
1986 }
1987
1988 void AliAnalysisTaskGammaConversion::AddToAODBranch(TClonesArray * branch, AliGammaConversionAODObject & particle) {
1989   //See header file for documentation
1990
1991   Int_t i = branch->GetEntriesFast();
1992   new((*branch)[i])  AliGammaConversionAODObject(particle);
1993 }
1994
1995 void AliAnalysisTaskGammaConversion::AddToAODBranch(TClonesArray * branch, AliAODConversionParticle & particle) {
1996   //See header file for documentation
1997
1998   Int_t i = branch->GetEntriesFast();
1999   new((*branch)[i])  AliAODConversionParticle(particle);
2000 }
2001
2002
2003 void AliAnalysisTaskGammaConversion::FillAODWithConversionGammas(){
2004   // Fill AOD with reconstructed Gamma
2005         
2006   for(Int_t gammaIndex=0;gammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();gammaIndex++){
2007     AliKFParticle * gammakf = dynamic_cast<AliKFParticle*>(fKFReconstructedGammasTClone->At(gammaIndex));
2008     
2009     if(fOutputAODClassName.Contains("AliAODPWG4Particle")) {
2010       AliAODPWG4Particle gamma = AliAODPWG4Particle(gammakf->Px(),gammakf->Py(),gammakf->Pz(), gammakf->E());
2011       //gamma.SetLabel(-1);//How to get the MC label of the reconstructed gamma?
2012       gamma.SetTrackLabel( fElectronv1[gammaIndex], fElectronv2[gammaIndex] ); //How to get the MC label of the 2 electrons that form the gamma?
2013       gamma.SetDetector("CTS"); //tag the gamma as reconstructed in the central barrel
2014       gamma.SetPdg(AliPID::kEleCon); //photon id
2015       gamma.SetTag(-1); //Here I usually put a flag saying that montecarlo says it is prompt, decay fragmentation photon, or hadrons or whatever
2016       //PH      gamma.SetChi2(gammakf->Chi2());
2017       
2018       AddToAODBranch(fAODGamma, gamma);
2019       
2020     } else if(fOutputAODClassName.Contains("ConversionParticle")) {
2021       TLorentzVector momentum(gammakf->Px(),gammakf->Py(),gammakf->Pz(), gammakf->E());
2022       AliAODConversionParticle gamma = AliAODConversionParticle(momentum);
2023       //gamma.SetLabel(-1);//How to get the MC label of the reconstructed gamma?
2024       gamma.SetTrackLabels( fElectronv1[gammaIndex], fElectronv2[gammaIndex] ); //How to get the MC label of the 2 electrons that form the gamma?
2025       //gamma.SetPdg(AliPID::kEleCon); //photon id
2026       //gamma.SetTag(-1); //Here I usually put a flag saying that montecarlo says it is prompt, decay fragmentation photon, or hadrons or whatever
2027       //PH      gamma.SetChi2(gammakf->Chi2());
2028       gamma.SetTrackLabels( fElectronv1[gammaIndex], fElectronv2[gammaIndex] );
2029       gamma.SetESDEvent(dynamic_cast<AliESDEvent*>(InputEvent()));    
2030       AddToAODBranch(fAODGamma, gamma);
2031       
2032       
2033
2034     } else {
2035       AliGammaConversionAODObject gamma;
2036       gamma.SetPx(gammakf->GetPx());
2037       gamma.SetPy(gammakf->GetPy());
2038       gamma.SetPz(gammakf->GetPz());
2039       gamma.SetLabel1(fElectronv1[gammaIndex]);
2040       gamma.SetLabel2(fElectronv2[gammaIndex]);
2041       gamma.SetChi2(gammakf->Chi2());
2042       gamma.SetE(gammakf->E());
2043       gamma.SetESDEvent(dynamic_cast<AliESDEvent*>(InputEvent()));
2044       AddToAODBranch(fAODGamma, gamma);
2045     }
2046         
2047   }
2048 }
2049 void AliAnalysisTaskGammaConversion::ProcessGammasForOmegaMesonAnalysis(){
2050   // omega meson analysis pi0+gamma decay
2051   for(Int_t firstPi0Index=0;firstPi0Index<fKFReconstructedPi0sTClone->GetEntriesFast();firstPi0Index++){
2052     AliKFParticle * omegaCandidatePi0Daughter = (AliKFParticle *)fKFReconstructedPi0sTClone->At(firstPi0Index);
2053     for(Int_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();firstGammaIndex++){
2054
2055       AliKFParticle * omegaCandidateGammaDaughter = (AliKFParticle *)fKFReconstructedGammasTClone->At(firstGammaIndex);
2056       if(fGammav1[firstPi0Index]==firstGammaIndex || fGammav2[firstPi0Index]==firstGammaIndex){
2057         continue;
2058       }
2059
2060       AliKFParticle omegaCandidate(*omegaCandidatePi0Daughter,*omegaCandidateGammaDaughter);
2061       Double_t massOmegaCandidate = 0.;
2062       Double_t widthOmegaCandidate = 0.;
2063
2064       omegaCandidate.GetMass(massOmegaCandidate,widthOmegaCandidate);
2065
2066       if ( massOmegaCandidate > 733 && massOmegaCandidate < 833 ) {
2067         //AddOmegaToAOD(&omegaCandidate, massOmegaCandidate, firstPi0Index, firstGammaIndex);
2068       }
2069       
2070       fHistograms->FillHistogram("ESD_Omega_InvMass_vs_Pt",massOmegaCandidate ,omegaCandidate.GetPt());
2071       fHistograms->FillHistogram("ESD_Omega_InvMass",massOmegaCandidate);
2072  
2073       //delete omegaCandidate;
2074
2075     }// end of omega reconstruction in pi0+gamma channel
2076
2077     if(fDoJet == kTRUE){
2078       AliKFParticle* negPiKF=NULL;
2079       AliKFParticle* posPiKF=NULL;
2080       
2081       // look at the pi+pi+pi0 channel 
2082       for(Int_t iCh=0;iCh<fChargedParticles->GetEntriesFast();iCh++){
2083         AliESDtrack* posTrack = (AliESDtrack*)(fChargedParticles->At(iCh));
2084         if (posTrack->GetSign()<0) continue;
2085         if(TMath::Abs(fV0Reader->GetESDpid()->NumberOfSigmasTPC(posTrack,AliPID::kPion))>2.) continue;
2086         if (posPiKF) delete posPiKF; posPiKF=NULL;
2087         posPiKF = new AliKFParticle( *(posTrack) ,211);
2088         
2089         for(Int_t jCh=0;jCh<fChargedParticles->GetEntriesFast();jCh++){
2090           AliESDtrack* negTrack = (AliESDtrack*)(fChargedParticles->At(jCh));
2091           if( negTrack->GetSign()>0) continue;
2092           if(TMath::Abs(fV0Reader->GetESDpid()->NumberOfSigmasTPC(negTrack,AliPID::kPion))>2.) continue;
2093           if (negPiKF) delete negPiKF; negPiKF=NULL;
2094           negPiKF = new AliKFParticle( *(negTrack) ,-211);
2095           AliKFParticle omegaCandidatePipPinPi0(*omegaCandidatePi0Daughter,*posPiKF,*negPiKF);
2096           Double_t massOmegaCandidatePipPinPi0 = 0.;
2097           Double_t widthOmegaCandidatePipPinPi0 = 0.;
2098           
2099           omegaCandidatePipPinPi0.GetMass(massOmegaCandidatePipPinPi0,widthOmegaCandidatePipPinPi0);
2100
2101           if ( massOmegaCandidatePipPinPi0 > 733 && massOmegaCandidatePipPinPi0 < 833 ) {
2102             // AddOmegaToAOD(&omegaCandidatePipPinPi0, massOmegaCandidatePipPinPi0, -1, -1);
2103           }
2104           
2105           fHistograms->FillHistogram("ESD_OmegaPipPinPi0_InvMass_vs_Pt",massOmegaCandidatePipPinPi0 ,omegaCandidatePipPinPi0.GetPt());
2106           fHistograms->FillHistogram("ESD_OmegaPipPinPi0_InvMass",massOmegaCandidatePipPinPi0);
2107
2108           //  delete omegaCandidatePipPinPi0;
2109         }
2110       }
2111
2112       if (posPiKF) delete posPiKF; posPiKF=NULL;     if (negPiKF) delete negPiKF; negPiKF=NULL;
2113
2114     } // checking ig gammajet because in that case the chargedparticle list is created
2115
2116   }
2117
2118   if(fCalculateBackground){
2119
2120     AliGammaConversionBGHandler * bgHandler = fV0Reader->GetBGHandler();
2121     
2122     Int_t zbin= bgHandler->GetZBinIndex(fV0Reader->GetVertexZ());
2123     Int_t mbin = 0;
2124     if(fUseTrackMultiplicityForBG == kTRUE){
2125       mbin = bgHandler->GetMultiplicityBinIndex(fV0Reader->CountESDTracks());
2126     }
2127     else{
2128       mbin = bgHandler->GetMultiplicityBinIndex(fV0Reader->GetNGoodV0s());
2129     }
2130     
2131     AliGammaConversionBGHandler::GammaConversionVertex *bgEventVertex = NULL;
2132
2133     // Background calculation for the omega
2134     for(Int_t nEventsInBG=0;nEventsInBG <fV0Reader->GetNBGEvents();nEventsInBG++){
2135       AliGammaConversionKFVector * previousEventV0s = bgHandler->GetBGGoodV0s(zbin,mbin,nEventsInBG);
2136       
2137       if(fMoveParticleAccordingToVertex == kTRUE){
2138         bgEventVertex = bgHandler->GetBGEventVertex(zbin,mbin,nEventsInBG);
2139       }
2140       for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
2141         AliKFParticle previousGoodV0 = (AliKFParticle)(*(previousEventV0s->at(iPrevious)));
2142
2143         if(fMoveParticleAccordingToVertex == kTRUE){
2144           MoveParticleAccordingToVertex(&previousGoodV0,bgEventVertex);
2145         }
2146
2147         for(Int_t firstPi0Index=0;firstPi0Index<fKFReconstructedPi0sTClone->GetEntriesFast();firstPi0Index++){
2148           AliKFParticle * omegaCandidatePi0Daughter = (AliKFParticle *)fKFReconstructedPi0sTClone->At(firstPi0Index);
2149           AliKFParticle * omegaBckCandidate = new AliKFParticle(*omegaCandidatePi0Daughter,previousGoodV0);
2150           Double_t massOmegaBckCandidate = 0.;
2151           Double_t widthOmegaBckCandidate = 0.;
2152           
2153           omegaBckCandidate->GetMass(massOmegaBckCandidate,widthOmegaBckCandidate);
2154
2155
2156           fHistograms->FillHistogram("ESD_Omega_Bck_InvMass_vs_Pt",massOmegaBckCandidate ,omegaBckCandidate->GetPt());
2157           fHistograms->FillHistogram("ESD_Omega_Bck_InvMass",massOmegaBckCandidate);
2158
2159           delete omegaBckCandidate;
2160         }
2161       }
2162     }
2163   } // end of checking if background calculation is available
2164 }
2165
2166
2167 void AliAnalysisTaskGammaConversion::AddOmegaToAOD(const AliKFParticle * const omegakf, Double_t mass, Int_t omegaDaughter, Int_t gammaDaughter) {
2168   //See header file for documentation
2169   AliGammaConversionAODObject omega;
2170   omega.SetPx(omegakf->GetPx());
2171   omega.SetPy(omegakf->GetPy());
2172   omega.SetPz(omegakf->GetPz());
2173   omega.SetChi2(omegakf->GetChi2());
2174   omega.SetE(omegakf->GetE());
2175   omega.SetIMass(mass);
2176   omega.SetLabel1(omegaDaughter);
2177   // //dynamic_cast<AliAODPWG4Particle*>(fAODBranch->At(daughter1))->SetTagged(kTRUE);
2178   omega.SetLabel2(gammaDaughter);
2179   AddToAODBranch(fAODOmega, omega);
2180 }
2181
2182
2183
2184 void AliAnalysisTaskGammaConversion::ProcessGammasForNeutralMesonAnalysis(){
2185   // see header file for documentation
2186         
2187   //  for(UInt_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammas.size();firstGammaIndex++){
2188   //    for(UInt_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fKFReconstructedGammas.size();secondGammaIndex++){
2189
2190   fESDEvent = fV0Reader->GetESDEvent();
2191
2192   if(fKFReconstructedGammasTClone->GetEntriesFast()>fV0Reader->GetNumberOfV0s()){
2193     cout<<"Warning, number of entries in the tclone is bigger than number of v0s"<<endl;
2194   }
2195
2196   for(Int_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();firstGammaIndex++){
2197     for(Int_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();secondGammaIndex++){
2198                         
2199       //      AliKFParticle * twoGammaDecayCandidateDaughter0 = &fKFReconstructedGammas[firstGammaIndex];
2200       //      AliKFParticle * twoGammaDecayCandidateDaughter1 = &fKFReconstructedGammas[secondGammaIndex];
2201                         
2202       AliKFParticle * twoGammaDecayCandidateDaughter0 = (AliKFParticle *)fKFReconstructedGammasTClone->At(firstGammaIndex);
2203       AliKFParticle * twoGammaDecayCandidateDaughter1 = (AliKFParticle *)fKFReconstructedGammasTClone->At(secondGammaIndex);
2204                         
2205       if(fElectronv1[firstGammaIndex]==fElectronv1[secondGammaIndex] || fElectronv1[firstGammaIndex]==fElectronv2[secondGammaIndex]){
2206         continue;
2207       }
2208       if(fElectronv2[firstGammaIndex]==fElectronv1[secondGammaIndex] || fElectronv2[firstGammaIndex]==fElectronv2[secondGammaIndex]){
2209         continue;
2210       }
2211                         
2212       AliKFParticle *twoGammaCandidate = new AliKFParticle(*twoGammaDecayCandidateDaughter0,*twoGammaDecayCandidateDaughter1);
2213                         
2214       Double_t massTwoGammaCandidate = 0.;
2215       Double_t widthTwoGammaCandidate = 0.;
2216       Double_t chi2TwoGammaCandidate =10000.;   
2217       twoGammaCandidate->GetMass(massTwoGammaCandidate,widthTwoGammaCandidate);
2218       //      if(twoGammaCandidate->GetNDF()>0){
2219       //        chi2TwoGammaCandidate = twoGammaCandidate->GetChi2()/twoGammaCandidate->GetNDF();
2220       chi2TwoGammaCandidate = twoGammaCandidate->GetChi2();
2221                                 
2222         fHistograms->FillHistogram("ESD_Mother_Chi2",chi2TwoGammaCandidate);
2223         if((chi2TwoGammaCandidate>0 && chi2TwoGammaCandidate<fV0Reader->GetChi2CutMeson()) || fApplyChi2Cut == kFALSE){
2224                                         
2225           TVector3 momentumVectorTwoGammaCandidate(twoGammaCandidate->GetPx(),twoGammaCandidate->GetPy(),twoGammaCandidate->GetPz());
2226           TVector3 spaceVectorTwoGammaCandidate(twoGammaCandidate->GetX(),twoGammaCandidate->GetY(),twoGammaCandidate->GetZ());
2227                                         
2228           Double_t openingAngleTwoGammaCandidate = twoGammaDecayCandidateDaughter0->GetAngle(*twoGammaDecayCandidateDaughter1);                                 
2229           Double_t rapidity;
2230           if(twoGammaCandidate->GetE() - twoGammaCandidate->GetPz() <= 0 || twoGammaCandidate->GetE() + twoGammaCandidate->GetPz() <= 0){
2231             cout << "Error: |Pz| > E !!!! " << endl;
2232             rapidity=0;
2233           }
2234           else{
2235             rapidity = 0.5*(TMath::Log((twoGammaCandidate->GetE() +twoGammaCandidate->GetPz()) / (twoGammaCandidate->GetE()-twoGammaCandidate->GetPz())));
2236           }
2237
2238           if(TMath::Abs(rapidity) > fV0Reader->GetRapidityMesonCut()){
2239             delete twoGammaCandidate;
2240             continue;   // rapidity cut
2241           }
2242
2243
2244           Double_t alfa=0.0;
2245           if( (twoGammaDecayCandidateDaughter0->GetE()+twoGammaDecayCandidateDaughter1->GetE()) != 0){
2246             alfa=TMath::Abs((twoGammaDecayCandidateDaughter0->GetE()-twoGammaDecayCandidateDaughter1->GetE())
2247               /(twoGammaDecayCandidateDaughter0->GetE()+twoGammaDecayCandidateDaughter1->GetE()));
2248           }
2249                                         
2250           if(openingAngleTwoGammaCandidate < fMinOpeningAngleGhostCut){
2251             delete twoGammaCandidate;
2252             continue;   // minimum opening angle to avoid using ghosttracks
2253           }
2254                         
2255           if(alfa>fV0Reader->GetAlphaMinCutMeson() && alfa<fV0Reader->GetAlphaCutMeson()){
2256             fHistograms->FillHistogram("ESD_Mother_GammaDaughter_OpeningAngle", openingAngleTwoGammaCandidate);
2257             fHistograms->FillHistogram("ESD_Mother_Energy", twoGammaCandidate->GetE());
2258             fHistograms->FillHistogram("ESD_Mother_Pt", momentumVectorTwoGammaCandidate.Pt());
2259             fHistograms->FillHistogram("ESD_Mother_Eta", momentumVectorTwoGammaCandidate.Eta());
2260             fHistograms->FillHistogram("ESD_Mother_Rapid", rapidity);                                   
2261             fHistograms->FillHistogram("ESD_Mother_Phi", spaceVectorTwoGammaCandidate.Phi());
2262             fHistograms->FillHistogram("ESD_Mother_Mass", massTwoGammaCandidate);
2263             fHistograms->FillHistogram("ESD_Mother_alfa", alfa);
2264             if(massTwoGammaCandidate>0.1 && massTwoGammaCandidate<0.15){
2265               fHistograms->FillHistogram("ESD_Mother_alfa_Pi0", alfa);
2266             }
2267             if(massTwoGammaCandidate>0.5 && massTwoGammaCandidate<0.57){
2268               fHistograms->FillHistogram("ESD_Mother_alfa_Eta", alfa);
2269             }
2270
2271             fHistograms->FillHistogram("ESD_Mother_R", spaceVectorTwoGammaCandidate.Pt());    // Pt in Space == R!!!
2272             fHistograms->FillHistogram("ESD_Mother_ZR", twoGammaCandidate->GetZ(), spaceVectorTwoGammaCandidate.Pt());
2273             fHistograms->FillHistogram("ESD_Mother_XY", twoGammaCandidate->GetX(), twoGammaCandidate->GetY());
2274             fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2275             fHistograms->FillHistogram("ESD_Mother_InvMass",massTwoGammaCandidate);         
2276             fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_alpha",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2277           }
2278           if(alfa<0.1){
2279             fHistograms->FillHistogram("ESD_Mother_InvMass_vs_E_alpha",massTwoGammaCandidate ,twoGammaCandidate->GetE());
2280           }
2281
2282
2283           if(fCalculateBackground){
2284             /* Kenneth, just for testing*/
2285             AliGammaConversionBGHandler * bgHandlerTest = fV0Reader->GetBGHandler();
2286             
2287             Int_t zbin= bgHandlerTest->GetZBinIndex(fV0Reader->GetVertexZ());
2288             Int_t mbin=0;
2289             Int_t multKAA=0;
2290             if(fUseTrackMultiplicityForBG == kTRUE){
2291               multKAA=fV0Reader->CountESDTracks();
2292               mbin = bgHandlerTest->GetMultiplicityBinIndex(fV0Reader->CountESDTracks());
2293             }
2294             else{// means we use #v0s for multiplicity
2295               multKAA=fV0Reader->GetNGoodV0s();
2296               mbin = bgHandlerTest->GetMultiplicityBinIndex(fV0Reader->GetNGoodV0s());
2297             }
2298             //      cout<<"Filling bin number "<<zbin<<" and "<<mbin<<endl;
2299             //      cout<<"Corresponding to z = "<<fV0Reader->GetVertexZ()<<" and m = "<<multKAA<<endl;
2300             if(alfa>fV0Reader->GetAlphaMinCutMeson() && alfa<fV0Reader->GetAlphaCutMeson()){
2301               fHistograms->FillHistogram(Form("%d%dESD_Mother_InvMass",zbin,mbin),massTwoGammaCandidate);
2302               fHistograms->FillHistogram(Form("%d%dESD_Mother_InvMass_vs_Pt",zbin,mbin),massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2303               /* end Kenneth, just for testing*/
2304               fHistograms->FillHistogram(Form("%dESD_Mother_InvMass_vs_Pt",mbin),massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2305             }
2306           }
2307           /*      if(fCalculateBackground){
2308             AliGammaConversionBGHandler * bgHandler = fV0Reader->GetBGHandler();
2309             Int_t mbin= bgHandler->GetMultiplicityBinIndex(fV0Reader->CountESDTracks());
2310             fHistograms->FillHistogram(Form("%dESD_Mother_InvMass_vs_Pt",mbin),massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2311             }*/
2312           //      if(fDoNeutralMesonV0MCCheck){
2313           if(fDoMCTruth){
2314             //Kenneth: Checking the eta of the gamma to check the difference between 0.9 and 1.2
2315             Int_t indexKF1 = fKFReconstructedGammasV0Index.at(firstGammaIndex);
2316             if(indexKF1<fV0Reader->GetNumberOfV0s()){
2317               fV0Reader->GetV0(indexKF1);//updates to the correct v0
2318               Double_t eta1 = fV0Reader->GetMotherCandidateEta();
2319               Bool_t isRealPi0=kFALSE;
2320               Bool_t isRealEta=kFALSE;
2321               Int_t gamma1MotherLabel=-1;
2322               if(fV0Reader->HasSameMCMother() == kTRUE){
2323                 //cout<<"This v0 is a real v0!!!!"<<endl;
2324                 TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
2325                 TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
2326                 if(TMath::Abs(negativeMC->GetPdgCode())==11 && TMath::Abs(positiveMC->GetPdgCode())==11){
2327                   if(negativeMC->GetUniqueID() == 5 && positiveMC->GetUniqueID() ==5){
2328                     if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
2329                       gamma1MotherLabel=fV0Reader->GetMotherMCParticle()->GetFirstMother();
2330                     }
2331                   }
2332                 }
2333               }
2334               Int_t indexKF2 = fKFReconstructedGammasV0Index.at(secondGammaIndex);
2335               if(indexKF1 == indexKF2){
2336                 cout<<"index of the two KF particles are the same.... should not happen"<<endl;
2337               }
2338               if(indexKF2<fV0Reader->GetNumberOfV0s()){
2339                 fV0Reader->GetV0(indexKF2);
2340                 Double_t eta2 = fV0Reader->GetMotherCandidateEta();
2341                 Int_t gamma2MotherLabel=-1;
2342                 if(fV0Reader->HasSameMCMother() == kTRUE){
2343                   TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
2344                   TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
2345                   if(TMath::Abs(negativeMC->GetPdgCode())==11 && TMath::Abs(positiveMC->GetPdgCode())==11){
2346                     if(negativeMC->GetUniqueID() == 5 && positiveMC->GetUniqueID() ==5){
2347                       if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
2348                         gamma2MotherLabel=fV0Reader->GetMotherMCParticle()->GetFirstMother();
2349                       }
2350                     }
2351                   }
2352                 }
2353                 if(gamma1MotherLabel>=0 && gamma1MotherLabel==gamma2MotherLabel){
2354                   if(fV0Reader->CheckIfPi0IsMother(gamma1MotherLabel)){
2355                     isRealPi0=kTRUE;
2356                   }
2357                   if(fV0Reader->CheckIfEtaIsMother(gamma1MotherLabel)){
2358                     isRealEta=kTRUE;
2359                   }
2360
2361                 }
2362                 if(alfa>fV0Reader->GetAlphaMinCutMeson() && alfa<fV0Reader->GetAlphaCutMeson()){
2363                   if(TMath::Abs(eta1)>0.9 && TMath::Abs(eta2)>0.9){
2364                     //            fHistograms->FillHistogram("ESD_Mother_InvMass_1212",massTwoGammaCandidate);
2365                     //            fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt1212",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
2366                     if(isRealPi0 || isRealEta){
2367                       fHistograms->FillHistogram("ESD_TruePi0_InvMass_1212",massTwoGammaCandidate);
2368                       fHistograms->FillHistogram("ESD_TruePi0_OpeningAngle_1212",openingAngleTwoGammaCandidate);
2369                       fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt1212",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
2370                       fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2371                       fHistograms->FillHistogram("ESD_TruePi0_InvMass",massTwoGammaCandidate);
2372                       fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt_alpha",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2373                     }
2374
2375                     if(!isRealPi0 && !isRealEta){
2376                       if(gamma1MotherLabel>-1 && gamma2MotherLabel>-1){
2377                         fHistograms->FillHistogram("ESD_TrueBckGG_InvMass",massTwoGammaCandidate);
2378                       }else{
2379                         fHistograms->FillHistogram("ESD_TrueBckCont_InvMass",massTwoGammaCandidate);
2380                       }
2381                     }
2382
2383                   }
2384                   else if(TMath::Abs(eta1)>0.9 || TMath::Abs(eta2)>0.9){
2385                   //              fHistograms->FillHistogram("ESD_Mother_InvMass_0912",massTwoGammaCandidate);
2386                   //              fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt0912",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
2387                   
2388                     if(isRealPi0 || isRealEta){
2389                       fHistograms->FillHistogram("ESD_TruePi0_InvMass_0912",massTwoGammaCandidate);
2390                       fHistograms->FillHistogram("ESD_TruePi0_OpeningAngle_0912",openingAngleTwoGammaCandidate);
2391                       fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt0912",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
2392                       fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2393                       fHistograms->FillHistogram("ESD_TruePi0_InvMass",massTwoGammaCandidate);
2394                       fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt_alpha",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2395                     }
2396                     if(!isRealPi0 && !isRealEta){
2397                       if(gamma1MotherLabel>-1 && gamma2MotherLabel>-1){
2398                         fHistograms->FillHistogram("ESD_TrueBckGG_InvMass",massTwoGammaCandidate);
2399                       }else{
2400                         fHistograms->FillHistogram("ESD_TrueBckCont_InvMass",massTwoGammaCandidate);
2401                       }
2402                     }
2403                   }
2404                   else{
2405                   //              fHistograms->FillHistogram("ESD_Mother_InvMass_0909",massTwoGammaCandidate);
2406                   //              fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt0909",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
2407                     if(isRealPi0 || isRealEta){
2408                       fHistograms->FillHistogram("ESD_TruePi0_InvMass_0909",massTwoGammaCandidate);
2409                       fHistograms->FillHistogram("ESD_TruePi0_OpeningAngle_0909",openingAngleTwoGammaCandidate);
2410                       fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt0909",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
2411                       fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2412                       fHistograms->FillHistogram("ESD_TruePi0_InvMass",massTwoGammaCandidate);
2413                       fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt_alpha",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2414                       if(gamma1MotherLabel > fV0Reader->GetMCStack()->GetNprimary()){
2415                         fHistograms->FillHistogram("ESD_TruePi0Sec_InvMass",massTwoGammaCandidate);
2416                       }
2417                     }
2418                     if(!isRealPi0 && !isRealEta){
2419                       if(gamma1MotherLabel>-1 && gamma2MotherLabel>-1){
2420                         fHistograms->FillHistogram("ESD_TrueBckGG_InvMass",massTwoGammaCandidate);
2421                       }else{
2422                         fHistograms->FillHistogram("ESD_TrueBckCont_InvMass",massTwoGammaCandidate);
2423                       }
2424                     }
2425                   }
2426                 }
2427               }
2428             }
2429           }
2430           if(alfa>fV0Reader->GetAlphaMinCutMeson() && alfa<fV0Reader->GetAlphaCutMeson()){
2431             if ( TMath::Abs(twoGammaDecayCandidateDaughter0->GetEta())<0.9 &&  TMath::Abs(twoGammaDecayCandidateDaughter1->GetEta())<0.9 ){
2432               fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_Fiducial",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2433               fHistograms->FillHistogram("ESD_Mother_InvMass_Fiducial",massTwoGammaCandidate);
2434             }
2435             
2436             if(TMath::Abs(twoGammaDecayCandidateDaughter0->GetEta())>0.9 && TMath::Abs(twoGammaDecayCandidateDaughter1->GetEta())>0.9){
2437               fHistograms->FillHistogram("ESD_Mother_InvMass_1212",massTwoGammaCandidate);
2438               fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt1212",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
2439             }
2440             else if(TMath::Abs(twoGammaDecayCandidateDaughter0->GetEta())>0.9 || TMath::Abs(twoGammaDecayCandidateDaughter1->GetEta())>0.9){
2441               fHistograms->FillHistogram("ESD_Mother_InvMass_0912",massTwoGammaCandidate);
2442               fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt0912",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
2443             }
2444             else{
2445               fHistograms->FillHistogram("ESD_Mother_InvMass_0909",massTwoGammaCandidate);
2446               fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt0909",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
2447             }
2448
2449             Double_t lowMassPi0=0.1;
2450             Double_t highMassPi0=0.15;
2451             if ( ( massTwoGammaCandidate > lowMassPi0) && (massTwoGammaCandidate < highMassPi0) ){
2452               new((*fKFReconstructedPi0sTClone)[fKFReconstructedPi0sTClone->GetEntriesFast()])  AliKFParticle(*twoGammaCandidate);
2453               fGammav1.push_back(firstGammaIndex);
2454               fGammav2.push_back(secondGammaIndex);
2455               if( fKFCreateAOD ) {
2456                 AddPionToAOD(twoGammaCandidate, massTwoGammaCandidate, firstGammaIndex, secondGammaIndex);
2457               }
2458             }
2459           }
2460
2461         }
2462         delete twoGammaCandidate;
2463     }
2464   }
2465 }
2466
2467 void AliAnalysisTaskGammaConversion::AddPionToAOD(AliKFParticle * pionkf, Double_t mass, Int_t daughter1, Int_t daughter2) {
2468   //See header file for documentation
2469   if(fOutputAODClassName.Contains("AODObject")) {
2470     AliGammaConversionAODObject pion;
2471     pion.SetPx(pionkf->GetPx());
2472     pion.SetPy(pionkf->GetPy());
2473     pion.SetPz(pionkf->GetPz());
2474     pion.SetChi2(pionkf->GetChi2());
2475     pion.SetE(pionkf->GetE());
2476     pion.SetIMass(mass);
2477     pion.SetLabel1(daughter1);
2478     pion.SetLabel2(daughter2);
2479     AddToAODBranch(fAODPi0, pion);
2480   } else {
2481     TLorentzVector momentum(pionkf->Px(),pionkf->Py(),pionkf->Pz(), pionkf->E());
2482     AliAODConversionParticle pion = AliAODConversionParticle(momentum);
2483     pion.SetTrackLabels( daughter1, daughter2 );
2484     pion.SetChi2(pionkf->GetChi2());
2485     AddToAODBranch(fAODPi0, pion);
2486     
2487   }
2488 }
2489
2490
2491 /*
2492   void AliAnalysisTaskGammaConversion::ProcessConvPHOSGammasForNeutralMesonAnalysis(){
2493
2494   // see header file for documentation
2495   // Analyse Pi0 with one photon from Phos and 1 photon from conversions
2496         
2497
2498
2499   Double_t vtx[3];
2500   vtx[0] = fV0Reader->GetPrimaryVertex()->GetX();
2501   vtx[1] = fV0Reader->GetPrimaryVertex()->GetY();
2502   vtx[2] = fV0Reader->GetPrimaryVertex()->GetZ();
2503
2504
2505   // Loop over all CaloClusters and consider only the PHOS ones:
2506   AliESDCaloCluster *clu;
2507   TLorentzVector pPHOS;
2508   TLorentzVector gammaPHOS;
2509   TLorentzVector gammaGammaConv;
2510   TLorentzVector pi0GammaConvPHOS;
2511   TLorentzVector gammaGammaConvBck;
2512   TLorentzVector pi0GammaConvPHOSBck;
2513
2514
2515   for (Int_t i=0; i<fV0Reader->GetESDEvent()->GetNumberOfCaloClusters(); i++) {
2516   clu = fV0Reader->GetESDEvent()->GetCaloCluster(i);
2517   if ( !clu->IsPHOS() || clu->E()<0.1 ) continue;
2518   clu ->GetMomentum(pPHOS ,vtx);
2519   for(Int_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();firstGammaIndex++){
2520   AliKFParticle * twoGammaDecayCandidateDaughter0 = (AliKFParticle *)fKFReconstructedGammasTClone->At(firstGammaIndex);
2521   gammaGammaConv.SetXYZM(twoGammaDecayCandidateDaughter0->Px(),twoGammaDecayCandidateDaughter0->Py(),twoGammaDecayCandidateDaughter0->Pz(),0.);
2522   gammaPHOS.SetXYZM(pPHOS.Px(),pPHOS.Py(),pPHOS.Pz(),0.);
2523   pi0GammaConvPHOS=gammaGammaConv+gammaPHOS;
2524   fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvPHOS",pi0GammaConvPHOS.M());
2525   fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_GammaConvPHOS",pi0GammaConvPHOS.M(),pi0GammaConvPHOS.Pt());
2526
2527   TVector3 v3D0(twoGammaDecayCandidateDaughter0->Px(),twoGammaDecayCandidateDaughter0->Py(),twoGammaDecayCandidateDaughter0->Pz());
2528   TVector3 v3D1(gammaPHOS.Px(),gammaPHOS.Py(),gammaPHOS.Pz());
2529   Double_t opanConvPHOS= v3D0.Angle(v3D1);
2530   if ( opanConvPHOS < 0.35){
2531   fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvPHOS_OpanLow",pi0GammaConvPHOS.M());
2532   }else{
2533   fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvPHOS_OpanHigh",pi0GammaConvPHOS.M());
2534   }
2535
2536   }
2537
2538   //  Now the LorentVector pPHOS is obtained and can be paired with the converted proton
2539   }
2540   //==== End of the PHOS cluster selection ============
2541   TLorentzVector pEMCAL;
2542   TLorentzVector gammaEMCAL;
2543   TLorentzVector pi0GammaConvEMCAL;
2544   TLorentzVector pi0GammaConvEMCALBck;
2545
2546   for (Int_t i=0; i<fV0Reader->GetESDEvent()->GetNumberOfCaloClusters(); i++) {
2547   clu = fV0Reader->GetESDEvent()->GetCaloCluster(i);
2548   if ( !clu->IsEMCAL()  || clu->E()<0.1 ) continue;
2549   if (clu->GetNCells() <= 1) continue;
2550   if ( clu->GetTOF()*1e9 < 550  || clu->GetTOF()*1e9 > 750) continue;
2551
2552   clu ->GetMomentum(pEMCAL ,vtx);
2553   for(Int_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();firstGammaIndex++){
2554   AliKFParticle * twoGammaDecayCandidateDaughter0 = (AliKFParticle *)fKFReconstructedGammasTClone->At(firstGammaIndex);
2555   gammaGammaConv.SetXYZM(twoGammaDecayCandidateDaughter0->Px(),
2556   twoGammaDecayCandidateDaughter0->Py(),
2557   twoGammaDecayCandidateDaughter0->Pz(),0.);
2558   gammaEMCAL.SetXYZM(pEMCAL.Px(),pEMCAL.Py(),pEMCAL.Pz(),0.);
2559   pi0GammaConvEMCAL=gammaGammaConv+gammaEMCAL;
2560   fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvEMCAL",pi0GammaConvEMCAL.M());
2561   fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_GammaConvEMCAL",pi0GammaConvEMCAL.M(),pi0GammaConvEMCAL.Pt());
2562   TVector3 v3D0(twoGammaDecayCandidateDaughter0->Px(),
2563   twoGammaDecayCandidateDaughter0->Py(),
2564   twoGammaDecayCandidateDaughter0->Pz());
2565   TVector3 v3D1(gammaEMCAL.Px(),gammaEMCAL.Py(),gammaEMCAL.Pz());
2566
2567
2568   Double_t opanConvEMCAL= v3D0.Angle(v3D1);
2569   if ( opanConvEMCAL < 0.35){
2570   fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvEMCAL_OpanLow",pi0GammaConvEMCAL.M());
2571   }else{
2572   fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvEMCAL_OpanHigh",pi0GammaConvEMCAL.M());
2573   }
2574
2575   }
2576   if(fCalculateBackground){
2577   for(Int_t nEventsInBG=0;nEventsInBG <fV0Reader->GetNBGEvents();nEventsInBG++){
2578   AliGammaConversionKFVector * previousEventV0s = fV0Reader->GetBGGoodV0s(nEventsInBG);
2579   for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
2580   AliKFParticle previousGoodV0 = (AliKFParticle)(*(previousEventV0s->at(iPrevious)));
2581   gammaGammaConvBck.SetXYZM(previousGoodV0.Px(),
2582   previousGoodV0.Py(),
2583   previousGoodV0.Pz(),0.);
2584   pi0GammaConvEMCALBck=gammaGammaConvBck+gammaEMCAL;
2585   fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvEMCAL_Bck",pi0GammaConvEMCALBck.M());
2586   fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_GammaConvEMCAL_Bck",pi0GammaConvEMCALBck.M(),
2587   pi0GammaConvEMCALBck.Pt());
2588   }
2589   }
2590       
2591   //  Now the LorentVector pEMCAL is obtained and can be paired with the converted proton
2592   } // end of checking if background photons are available
2593   }
2594   //==== End of the PHOS cluster selection ============
2595
2596   }
2597 */
2598
2599 void AliAnalysisTaskGammaConversion::MoveParticleAccordingToVertex(AliKFParticle * particle,const AliGammaConversionBGHandler::GammaConversionVertex *vertex){
2600   //see header file for documentation
2601
2602    Double_t dx = vertex->fX - fESDEvent->GetPrimaryVertex()->GetX();
2603    Double_t dy = vertex->fY - fESDEvent->GetPrimaryVertex()->GetY();
2604    Double_t dz = vertex->fZ - fESDEvent->GetPrimaryVertex()->GetZ();
2605   
2606   //  cout<<"dx, dy, dz: ["<<dx<<","<<dy<<","<<dz<<"]"<<endl;
2607    particle->X() = particle->GetX() - dx;
2608    particle->Y() = particle->GetY() - dy;
2609    particle->Z() = particle->GetZ() - dz;
2610 }
2611
2612 void AliAnalysisTaskGammaConversion::RotateKFParticle(AliKFParticle * kfParticle,Double_t angle){
2613   // Rotate the kf particle
2614   Double_t c = cos(angle);
2615   Double_t s = sin(angle);
2616   
2617   Double_t mA[7][ 7];
2618   for( Int_t i=0; i<7; i++ ){
2619     for( Int_t j=0; j<7; j++){
2620       mA[i][j] = 0;
2621     }
2622   }
2623   for( int i=0; i<7; i++ ){
2624     mA[i][i] = 1;
2625   }
2626   mA[0][0] =  c;  mA[0][1] = s;
2627   mA[1][0] = -s;  mA[1][1] = c;
2628   mA[3][3] =  c;  mA[3][4] = s;
2629   mA[4][3] = -s;  mA[4][4] = c;
2630   
2631   Double_t mAC[7][7];
2632   Double_t mAp[7];
2633   
2634   for( Int_t i=0; i<7; i++ ){
2635     mAp[i] = 0;
2636     for( Int_t k=0; k<7; k++){
2637       mAp[i]+=mA[i][k] * kfParticle->GetParameter(k);
2638     }
2639   }
2640   
2641   for( Int_t i=0; i<7; i++){
2642     kfParticle->Parameter(i) = mAp[i];
2643   }
2644
2645   for( Int_t i=0; i<7; i++ ){
2646     for( Int_t j=0; j<7; j++ ){
2647       mAC[i][j] = 0;
2648       for( Int_t k=0; k<7; k++ ){
2649         mAC[i][j]+= mA[i][k] * kfParticle->GetCovariance(k,j);
2650       }
2651     }
2652   }
2653
2654   for( Int_t i=0; i<7; i++ ){
2655     for( Int_t j=0; j<=i; j++ ){
2656       Double_t xx = 0;
2657       for( Int_t k=0; k<7; k++){
2658         xx+= mAC[i][k]*mA[j][k];
2659       }
2660       kfParticle->Covariance(i,j) = xx;
2661     }
2662   }
2663 }
2664
2665
2666 void AliAnalysisTaskGammaConversion::CalculateBackground(){
2667   // see header file for documentation
2668
2669
2670   TClonesArray * currentEventV0s = fV0Reader->GetCurrentEventGoodV0s();
2671
2672   AliGammaConversionBGHandler * bgHandler = fV0Reader->GetBGHandler();
2673   
2674   Int_t zbin= bgHandler->GetZBinIndex(fV0Reader->GetVertexZ());
2675   Int_t mbin = 0;
2676   if(fUseTrackMultiplicityForBG == kTRUE){
2677     mbin = bgHandler->GetMultiplicityBinIndex(fV0Reader->CountESDTracks());
2678   }
2679   else{
2680     mbin = bgHandler->GetMultiplicityBinIndex(fV0Reader->GetNGoodV0s());
2681   }
2682
2683   if(fDoRotation == kTRUE){
2684     TRandom3 *random = new TRandom3();
2685     for(Int_t iCurrent=0;iCurrent<currentEventV0s->GetEntriesFast();iCurrent++){
2686       AliKFParticle currentEventGoodV0 = *(AliKFParticle *)(currentEventV0s->At(iCurrent)); 
2687       for(Int_t iCurrent2=iCurrent+1;iCurrent2<currentEventV0s->GetEntriesFast();iCurrent2++){
2688         for(Int_t nRandom=0;nRandom<fNRandomEventsForBG;nRandom++){
2689         
2690           AliKFParticle currentEventGoodV02 = *(AliKFParticle *)(currentEventV0s->At(iCurrent2));
2691
2692           if(fCheckBGProbability == kTRUE){
2693             Double_t massBGprob =0.;
2694             Double_t widthBGprob = 0.;
2695             AliKFParticle *backgroundCandidateProb = new AliKFParticle(currentEventGoodV0,currentEventGoodV02);
2696             backgroundCandidateProb->GetMass(massBGprob,widthBGprob);
2697             if(massBGprob>0.1 && massBGprob<0.14){
2698               if(random->Rndm()>bgHandler->GetBGProb(zbin,mbin)){
2699                 delete backgroundCandidateProb;
2700                 continue;
2701               }
2702             }
2703             delete backgroundCandidateProb;
2704           }
2705         
2706           Double_t nRadiansPM = fNDegreesPMBackground*TMath::Pi()/180;
2707           
2708           Double_t rotationValue = random->Rndm()*2*nRadiansPM + TMath::Pi()-nRadiansPM;
2709
2710           RotateKFParticle(&currentEventGoodV02,rotationValue);
2711
2712           AliKFParticle *backgroundCandidate = new AliKFParticle(currentEventGoodV0,currentEventGoodV02);
2713
2714           Double_t massBG =0.;
2715           Double_t widthBG = 0.;
2716           Double_t chi2BG =10000.;      
2717           backgroundCandidate->GetMass(massBG,widthBG);
2718
2719           //      if(backgroundCandidate->GetNDF()>0){
2720           chi2BG = backgroundCandidate->GetChi2();
2721           if((chi2BG>0 && chi2BG<fV0Reader->GetChi2CutMeson())  || fApplyChi2Cut == kFALSE){
2722           
2723             TVector3 momentumVectorbackgroundCandidate(backgroundCandidate->GetPx(),backgroundCandidate->GetPy(),backgroundCandidate->GetPz());
2724             TVector3 spaceVectorbackgroundCandidate(backgroundCandidate->GetX(),backgroundCandidate->GetY(),backgroundCandidate->GetZ());
2725           
2726             Double_t openingAngleBG = currentEventGoodV0.GetAngle(currentEventGoodV02);
2727           
2728             Double_t rapidity;
2729             if(backgroundCandidate->GetE() - backgroundCandidate->GetPz() == 0 || backgroundCandidate->GetE() + backgroundCandidate->GetPz() == 0) rapidity=0;
2730             else rapidity = 0.5*(TMath::Log((backgroundCandidate->GetE() +backgroundCandidate->GetPz()) / (backgroundCandidate->GetE()-backgroundCandidate->GetPz())));
2731           
2732             if(TMath::Abs(rapidity) > fV0Reader->GetRapidityMesonCut() ){
2733               delete backgroundCandidate;   
2734               continue;   // rapidity cut
2735             }                   
2736                                         
2737           
2738             Double_t alfa=0.0;
2739             if( (currentEventGoodV0.GetE()+currentEventGoodV02.GetE()) != 0){
2740               alfa=TMath::Abs((currentEventGoodV0.GetE()-currentEventGoodV02.GetE())
2741                               /(currentEventGoodV0.GetE()+currentEventGoodV02.GetE()));
2742             }
2743           
2744           
2745             if(openingAngleBG < fMinOpeningAngleGhostCut ){
2746               delete backgroundCandidate;   
2747               continue;   // minimum opening angle to avoid using ghosttracks
2748             }                   
2749           
2750             // original
2751             if(alfa>fV0Reader->GetAlphaMinCutMeson() &&   alfa<fV0Reader->GetAlphaCutMeson()){
2752               fHistograms->FillHistogram("ESD_Background_GammaDaughter_OpeningAngle", openingAngleBG);
2753               fHistograms->FillHistogram("ESD_Background_Energy", backgroundCandidate->GetE());
2754               fHistograms->FillHistogram("ESD_Background_Pt",  momentumVectorbackgroundCandidate.Pt());
2755               fHistograms->FillHistogram("ESD_Background_Eta", momentumVectorbackgroundCandidate.Eta());
2756               fHistograms->FillHistogram("ESD_Background_Rapidity", rapidity);
2757               fHistograms->FillHistogram("ESD_Background_Phi", spaceVectorbackgroundCandidate.Phi());
2758               fHistograms->FillHistogram("ESD_Background_Mass", massBG);
2759               fHistograms->FillHistogram("ESD_Background_R", spaceVectorbackgroundCandidate.Pt());  // Pt in Space == R!!!!
2760               fHistograms->FillHistogram("ESD_Background_ZR", backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
2761               fHistograms->FillHistogram("ESD_Background_XY", backgroundCandidate->GetX(), backgroundCandidate->GetY());
2762               fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt",massBG,momentumVectorbackgroundCandidate.Pt());
2763               fHistograms->FillHistogram("ESD_Background_InvMass",massBG);
2764               fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_alpha",massBG,momentumVectorbackgroundCandidate.Pt());
2765
2766               if(massBG>0.1 && massBG<0.15){
2767                 fHistograms->FillHistogram("ESD_Background_alfa_Pi0", alfa);
2768               }
2769               if(massBG>0.5 && massBG<0.57){
2770                 fHistograms->FillHistogram("ESD_Background_alfa_Eta", alfa);
2771               }
2772
2773               if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 &&  TMath::Abs(currentEventGoodV02.GetEta())<0.9 ){
2774                 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_Fiducial",massBG,momentumVectorbackgroundCandidate.Pt());
2775                 fHistograms->FillHistogram("ESD_Background_InvMass_Fiducial",massBG);
2776               }
2777               
2778               fHistograms->FillHistogram(Form("%d%dESD_Background_GammaDaughter_OpeningAngle",zbin,mbin), openingAngleBG);
2779               fHistograms->FillHistogram(Form("%d%dESD_Background_Energy",zbin,mbin), backgroundCandidate->GetE());
2780               fHistograms->FillHistogram(Form("%d%dESD_Background_Pt",zbin,mbin),  momentumVectorbackgroundCandidate.Pt());
2781               fHistograms->FillHistogram(Form("%d%dESD_Background_Eta",zbin,mbin), momentumVectorbackgroundCandidate.Eta());
2782               fHistograms->FillHistogram(Form("%d%dESD_Background_Rapidity",zbin,mbin), rapidity);
2783               fHistograms->FillHistogram(Form("%d%dESD_Background_Phi",zbin,mbin), spaceVectorbackgroundCandidate.Phi());
2784               fHistograms->FillHistogram(Form("%d%dESD_Background_Mass",zbin,mbin), massBG);
2785               fHistograms->FillHistogram(Form("%d%dESD_Background_R",zbin,mbin), spaceVectorbackgroundCandidate.Pt());  // Pt in Space == R!!!!
2786               fHistograms->FillHistogram(Form("%d%dESD_Background_ZR",zbin,mbin), backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
2787               fHistograms->FillHistogram(Form("%d%dESD_Background_XY",zbin,mbin), backgroundCandidate->GetX(), backgroundCandidate->GetY());
2788               fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_vs_Pt",zbin,mbin),massBG,momentumVectorbackgroundCandidate.Pt());
2789               fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass",zbin,mbin),massBG);
2790               
2791               if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 &&  TMath::Abs(currentEventGoodV02.GetEta())<0.9 ){
2792                 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_vs_Pt_Fiducial",zbin,mbin),massBG,momentumVectorbackgroundCandidate.Pt());
2793                 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_Fiducial",zbin,mbin),massBG);
2794               }
2795             }
2796             if(alfa<0.1){
2797               fHistograms->FillHistogram("ESD_Background_InvMass_vs_E_alpha",massBG ,backgroundCandidate->GetE());
2798             }
2799
2800           }
2801           //}
2802           delete backgroundCandidate;      
2803         }
2804       }
2805     }
2806     delete random;
2807   }
2808   else{ // means no rotation
2809     AliGammaConversionBGHandler::GammaConversionVertex *bgEventVertex = NULL;
2810             
2811     if(fUseTrackMultiplicityForBG){
2812       //    cout<<"Using charged track multiplicity for background calculation"<<endl;
2813       for(Int_t nEventsInBG=0;nEventsInBG <fV0Reader->GetNBGEvents();nEventsInBG++){
2814
2815         AliGammaConversionKFVector * previousEventV0s = bgHandler->GetBGGoodV0s(zbin,mbin,nEventsInBG);//fV0Reader->GetBGGoodV0s(nEventsInBG);
2816       
2817         if(fMoveParticleAccordingToVertex == kTRUE){
2818           bgEventVertex = bgHandler->GetBGEventVertex(zbin,mbin,nEventsInBG);
2819         }
2820
2821         for(Int_t iCurrent=0;iCurrent<currentEventV0s->GetEntriesFast();iCurrent++){
2822           AliKFParticle currentEventGoodV0 = *(AliKFParticle *)(currentEventV0s->At(iCurrent)); 
2823           for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
2824             AliKFParticle previousGoodV0 = (AliKFParticle)(*(previousEventV0s->at(iPrevious)));
2825             AliKFParticle previousGoodV0test = (AliKFParticle)(*(previousEventV0s->at(iPrevious)));
2826
2827             //cout<<"Primary Vertex event: ["<<fESDEvent->GetPrimaryVertex()->GetX()<<","<<fESDEvent->GetPrimaryVertex()->GetY()<<","<<fESDEvent->GetPrimaryVertex()->GetZ()<<"]"<<endl;
2828             //cout<<"BG prim Vertex event: ["<<bgEventVertex->fX<<","<<bgEventVertex->fY<<","<<bgEventVertex->fZ<<"]"<<endl;
2829           
2830             //cout<<"XYZ of particle before transport: ["<<previousGoodV0.X()<<","<<previousGoodV0.Y()<<","<<previousGoodV0.Z()<<"]"<<endl;
2831             if(fMoveParticleAccordingToVertex == kTRUE){
2832               MoveParticleAccordingToVertex(&previousGoodV0,bgEventVertex);
2833             }
2834             //cout<<"XYZ of particle after transport: ["<<previousGoodV0.X()<<","<<previousGoodV0.Y()<<","<<previousGoodV0.Z()<<"]"<<endl;
2835
2836             AliKFParticle *backgroundCandidate = new AliKFParticle(currentEventGoodV0,previousGoodV0);
2837         
2838             Double_t massBG =0.;
2839             Double_t widthBG = 0.;
2840             Double_t chi2BG =10000.;    
2841             backgroundCandidate->GetMass(massBG,widthBG);
2842             //  if(backgroundCandidate->GetNDF()>0){
2843             //    chi2BG = backgroundCandidate->GetChi2()/backgroundCandidate->GetNDF();
2844             chi2BG = backgroundCandidate->GetChi2();
2845             if((chi2BG>0 && chi2BG<fV0Reader->GetChi2CutMeson()) || fApplyChi2Cut == kFALSE){
2846                                         
2847               TVector3 momentumVectorbackgroundCandidate(backgroundCandidate->GetPx(),backgroundCandidate->GetPy(),backgroundCandidate->GetPz());
2848               TVector3 spaceVectorbackgroundCandidate(backgroundCandidate->GetX(),backgroundCandidate->GetY(),backgroundCandidate->GetZ());
2849                                         
2850               Double_t openingAngleBG = currentEventGoodV0.GetAngle(previousGoodV0);
2851                                         
2852               Double_t rapidity;
2853             
2854               if(backgroundCandidate->GetE() - backgroundCandidate->GetPz() <= 0 || backgroundCandidate->GetE() + backgroundCandidate->GetPz() <= 0){
2855                 cout << "Error: |Pz| > E !!!! " << endl;
2856                 rapidity=0;
2857               } else {
2858                 rapidity = 0.5*(TMath::Log((backgroundCandidate->GetE() +backgroundCandidate->GetPz()) / (backgroundCandidate->GetE()-backgroundCandidate->GetPz())));
2859               }                         
2860               if(TMath::Abs(rapidity) > fV0Reader->GetRapidityMesonCut() ){
2861                 delete backgroundCandidate;   
2862                 continue;   // rapidity cut
2863               }                 
2864                                                         
2865         
2866               Double_t alfa=0.0;
2867               if( (currentEventGoodV0.GetE()+previousGoodV0.GetE()) != 0){
2868                 alfa=TMath::Abs((currentEventGoodV0.GetE()-previousGoodV0.GetE())
2869                                 /(currentEventGoodV0.GetE()+previousGoodV0.GetE()));
2870               }
2871                         
2872                                         
2873               if(openingAngleBG < fMinOpeningAngleGhostCut ){
2874                 delete backgroundCandidate;   
2875                 continue;   // minimum opening angle to avoid using ghosttracks
2876               }                 
2877
2878               // original
2879               if(alfa>fV0Reader->GetAlphaMinCutMeson() &&   alfa<fV0Reader->GetAlphaCutMeson()){
2880                 fHistograms->FillHistogram("ESD_Background_GammaDaughter_OpeningAngle", openingAngleBG);
2881                 fHistograms->FillHistogram("ESD_Background_Energy", backgroundCandidate->GetE());
2882                 fHistograms->FillHistogram("ESD_Background_Pt",  momentumVectorbackgroundCandidate.Pt());
2883                 fHistograms->FillHistogram("ESD_Background_Eta", momentumVectorbackgroundCandidate.Eta());
2884                 fHistograms->FillHistogram("ESD_Background_Rapidity", rapidity);
2885                 fHistograms->FillHistogram("ESD_Background_Phi", spaceVectorbackgroundCandidate.Phi());
2886                 fHistograms->FillHistogram("ESD_Background_Mass", massBG);
2887                 fHistograms->FillHistogram("ESD_Background_R", spaceVectorbackgroundCandidate.Pt());  // Pt in Space == R!!!!
2888                 fHistograms->FillHistogram("ESD_Background_ZR", backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
2889                 fHistograms->FillHistogram("ESD_Background_XY", backgroundCandidate->GetX(), backgroundCandidate->GetY());
2890                 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt",massBG,momentumVectorbackgroundCandidate.Pt());
2891                 fHistograms->FillHistogram("ESD_Background_InvMass",massBG);
2892                 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_alpha",massBG,momentumVectorbackgroundCandidate.Pt());
2893
2894                 if(massBG>0.1 && massBG<0.15){
2895                   fHistograms->FillHistogram("ESD_Background_alfa_Pi0", alfa);
2896                 }
2897                 if(massBG>0.5 && massBG<0.57){
2898                   fHistograms->FillHistogram("ESD_Background_alfa_Eta", alfa);
2899                 }
2900
2901                 if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 &&  TMath::Abs(previousGoodV0.GetEta())<0.9 ){
2902                   fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_Fiducial",massBG,momentumVectorbackgroundCandidate.Pt());
2903                   fHistograms->FillHistogram("ESD_Background_InvMass_Fiducial",massBG);
2904                 }
2905
2906                 // test
2907                 fHistograms->FillHistogram(Form("%d%dESD_Background_GammaDaughter_OpeningAngle",zbin,mbin), openingAngleBG);
2908                 fHistograms->FillHistogram(Form("%d%dESD_Background_Energy",zbin,mbin), backgroundCandidate->GetE());
2909                 fHistograms->FillHistogram(Form("%d%dESD_Background_Pt",zbin,mbin),  momentumVectorbackgroundCandidate.Pt());
2910                 fHistograms->FillHistogram(Form("%d%dESD_Background_Eta",zbin,mbin), momentumVectorbackgroundCandidate.Eta());
2911                 fHistograms->FillHistogram(Form("%d%dESD_Background_Rapidity",zbin,mbin), rapidity);
2912                 fHistograms->FillHistogram(Form("%d%dESD_Background_Phi",zbin,mbin), spaceVectorbackgroundCandidate.Phi());
2913                 fHistograms->FillHistogram(Form("%d%dESD_Background_Mass",zbin,mbin), massBG);
2914                 fHistograms->FillHistogram(Form("%d%dESD_Background_R",zbin,mbin), spaceVectorbackgroundCandidate.Pt());  // Pt in Space == R!!!!
2915                 fHistograms->FillHistogram(Form("%d%dESD_Background_ZR",zbin,mbin), backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
2916                 fHistograms->FillHistogram(Form("%d%dESD_Background_XY",zbin,mbin), backgroundCandidate->GetX(), backgroundCandidate->GetY());
2917                 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_vs_Pt",zbin,mbin),massBG,momentumVectorbackgroundCandidate.Pt());
2918                 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass",zbin,mbin),massBG);
2919                 
2920                 if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 &&  TMath::Abs(previousGoodV0.GetEta())<0.9 ){
2921                   fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_vs_Pt_Fiducial",zbin,mbin),massBG,momentumVectorbackgroundCandidate.Pt());
2922                   fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_Fiducial",zbin,mbin),massBG);
2923                 }
2924                 //        }
2925               }
2926               if(alfa<0.1){
2927                 fHistograms->FillHistogram("ESD_Background_InvMass_vs_E_alpha",massBG ,backgroundCandidate->GetE());
2928               }
2929
2930             }
2931             delete backgroundCandidate;      
2932           }
2933         }
2934       }
2935     }
2936     else{ // means using #V0s for multiplicity
2937
2938       //    cout<<"Using the v0 multiplicity to calculate background"<<endl;
2939     
2940       fHistograms->FillHistogram("ESD_Background_z_m",zbin,mbin);
2941       fHistograms->FillHistogram("ESD_Mother_multpilicityVSv0s",fV0Reader->CountESDTracks(),fV0Reader->GetNumberOfV0s());
2942
2943       for(Int_t nEventsInBG=0;nEventsInBG <fV0Reader->GetNBGEvents();nEventsInBG++){
2944         AliGammaConversionKFVector * previousEventV0s = bgHandler->GetBGGoodV0s(zbin,mbin,nEventsInBG);// fV0Reader->GetBGGoodV0s(nEventsInBG);
2945         if(previousEventV0s){
2946         
2947           if(fMoveParticleAccordingToVertex == kTRUE){
2948             bgEventVertex = bgHandler->GetBGEventVertex(zbin,mbin,nEventsInBG);
2949           }
2950
2951           for(Int_t iCurrent=0;iCurrent<currentEventV0s->GetEntriesFast();iCurrent++){
2952             AliKFParticle currentEventGoodV0 = *(AliKFParticle *)(currentEventV0s->At(iCurrent)); 
2953             for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
2954               AliKFParticle previousGoodV0 = (AliKFParticle)(*(previousEventV0s->at(iPrevious)));
2955
2956               if(fMoveParticleAccordingToVertex == kTRUE){
2957                 MoveParticleAccordingToVertex(&previousGoodV0,bgEventVertex);
2958               }
2959
2960               AliKFParticle *backgroundCandidate = new AliKFParticle(currentEventGoodV0,previousGoodV0);
2961               Double_t massBG =0.;
2962               Double_t widthBG = 0.;
2963               Double_t chi2BG =10000.;  
2964               backgroundCandidate->GetMass(massBG,widthBG);
2965               /*            if(backgroundCandidate->GetNDF()>0){
2966                             chi2BG = backgroundCandidate->GetChi2()/backgroundCandidate->GetNDF();
2967                             {//remember to remove
2968                             TVector3 momentumVectorbackgroundCandidate(backgroundCandidate->GetPx(),backgroundCandidate->GetPy(),backgroundCandidate->GetPz());
2969                             TVector3 spaceVectorbackgroundCandidate(backgroundCandidate->GetX(),backgroundCandidate->GetY(),backgroundCandidate->GetZ());
2970               
2971                             Double_t openingAngleBG = currentEventGoodV0.GetAngle(previousGoodV0);
2972                             fHistograms->FillHistogram("ESD_Background_GammaDaughter_OpeningAngle_nochi2", openingAngleBG);
2973                             }
2974               */
2975               chi2BG = backgroundCandidate->GetChi2();
2976               if((chi2BG>0 && chi2BG<fV0Reader->GetChi2CutMeson()) || fApplyChi2Cut == kFALSE){
2977                 TVector3 momentumVectorbackgroundCandidate(backgroundCandidate->GetPx(),backgroundCandidate->GetPy(),backgroundCandidate->GetPz());
2978                 TVector3 spaceVectorbackgroundCandidate(backgroundCandidate->GetX(),backgroundCandidate->GetY(),backgroundCandidate->GetZ());
2979                                         
2980                 Double_t openingAngleBG = currentEventGoodV0.GetAngle(previousGoodV0);
2981                                         
2982                 Double_t rapidity;
2983                 if(backgroundCandidate->GetE() - backgroundCandidate->GetPz() == 0 || backgroundCandidate->GetE() + backgroundCandidate->GetPz() == 0) rapidity=0;
2984                 else rapidity = 0.5*(TMath::Log((backgroundCandidate->GetE() +backgroundCandidate->GetPz()) / (backgroundCandidate->GetE()-backgroundCandidate->GetPz())));
2985                                         
2986                 if(TMath::Abs(rapidity) > fV0Reader->GetRapidityMesonCut() ){
2987                   delete backgroundCandidate;   
2988                   continue;   // rapidity cut
2989                 }                       
2990                                                                 
2991
2992                 Double_t alfa=0.0;
2993                 if( (currentEventGoodV0.GetE()+previousGoodV0.GetE()) != 0){