eb2f7b38467c8b6a9253e992db00a7c6af860c8f
[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 ( fKFCreateAOD && (massTwoGammaCandidate > lowMassPi0) && (massTwoGammaCandidate < highMassPi0) ){
2452               new((*fKFReconstructedPi0sTClone)[fKFReconstructedPi0sTClone->GetEntriesFast()])  AliKFParticle(*twoGammaCandidate);
2453               fGammav1.push_back(firstGammaIndex);
2454               fGammav2.push_back(secondGammaIndex);
2455               AddPionToAOD(twoGammaCandidate, massTwoGammaCandidate, firstGammaIndex, secondGammaIndex);
2456             }
2457           }
2458
2459         }
2460           //}
2461         delete twoGammaCandidate;
2462     }
2463   }
2464 }
2465
2466 void AliAnalysisTaskGammaConversion::AddPionToAOD(AliKFParticle * pionkf, Double_t mass, Int_t daughter1, Int_t daughter2) {
2467   //See header file for documentation
2468   if(fOutputAODClassName.Contains("AODObject")) {
2469     AliGammaConversionAODObject pion;
2470     pion.SetPx(pionkf->GetPx());
2471     pion.SetPy(pionkf->GetPy());
2472     pion.SetPz(pionkf->GetPz());
2473     pion.SetChi2(pionkf->GetChi2());
2474     pion.SetE(pionkf->GetE());
2475     pion.SetIMass(mass);
2476     pion.SetLabel1(daughter1);
2477     pion.SetLabel2(daughter2);
2478     AddToAODBranch(fAODPi0, pion);
2479   } else {
2480     TLorentzVector momentum(pionkf->Px(),pionkf->Py(),pionkf->Pz(), pionkf->E());
2481     AliAODConversionParticle pion = AliAODConversionParticle(momentum);
2482     pion.SetTrackLabels( daughter1, daughter2 );
2483     pion.SetChi2(pionkf->GetChi2());
2484     AddToAODBranch(fAODPi0, pion);
2485     
2486   }
2487 }
2488
2489
2490 /*
2491   void AliAnalysisTaskGammaConversion::ProcessConvPHOSGammasForNeutralMesonAnalysis(){
2492
2493   // see header file for documentation
2494   // Analyse Pi0 with one photon from Phos and 1 photon from conversions
2495         
2496
2497
2498   Double_t vtx[3];
2499   vtx[0] = fV0Reader->GetPrimaryVertex()->GetX();
2500   vtx[1] = fV0Reader->GetPrimaryVertex()->GetY();
2501   vtx[2] = fV0Reader->GetPrimaryVertex()->GetZ();
2502
2503
2504   // Loop over all CaloClusters and consider only the PHOS ones:
2505   AliESDCaloCluster *clu;
2506   TLorentzVector pPHOS;
2507   TLorentzVector gammaPHOS;
2508   TLorentzVector gammaGammaConv;
2509   TLorentzVector pi0GammaConvPHOS;
2510   TLorentzVector gammaGammaConvBck;
2511   TLorentzVector pi0GammaConvPHOSBck;
2512
2513
2514   for (Int_t i=0; i<fV0Reader->GetESDEvent()->GetNumberOfCaloClusters(); i++) {
2515   clu = fV0Reader->GetESDEvent()->GetCaloCluster(i);
2516   if ( !clu->IsPHOS() || clu->E()<0.1 ) continue;
2517   clu ->GetMomentum(pPHOS ,vtx);
2518   for(Int_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();firstGammaIndex++){
2519   AliKFParticle * twoGammaDecayCandidateDaughter0 = (AliKFParticle *)fKFReconstructedGammasTClone->At(firstGammaIndex);
2520   gammaGammaConv.SetXYZM(twoGammaDecayCandidateDaughter0->Px(),twoGammaDecayCandidateDaughter0->Py(),twoGammaDecayCandidateDaughter0->Pz(),0.);
2521   gammaPHOS.SetXYZM(pPHOS.Px(),pPHOS.Py(),pPHOS.Pz(),0.);
2522   pi0GammaConvPHOS=gammaGammaConv+gammaPHOS;
2523   fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvPHOS",pi0GammaConvPHOS.M());
2524   fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_GammaConvPHOS",pi0GammaConvPHOS.M(),pi0GammaConvPHOS.Pt());
2525
2526   TVector3 v3D0(twoGammaDecayCandidateDaughter0->Px(),twoGammaDecayCandidateDaughter0->Py(),twoGammaDecayCandidateDaughter0->Pz());
2527   TVector3 v3D1(gammaPHOS.Px(),gammaPHOS.Py(),gammaPHOS.Pz());
2528   Double_t opanConvPHOS= v3D0.Angle(v3D1);
2529   if ( opanConvPHOS < 0.35){
2530   fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvPHOS_OpanLow",pi0GammaConvPHOS.M());
2531   }else{
2532   fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvPHOS_OpanHigh",pi0GammaConvPHOS.M());
2533   }
2534
2535   }
2536
2537   //  Now the LorentVector pPHOS is obtained and can be paired with the converted proton
2538   }
2539   //==== End of the PHOS cluster selection ============
2540   TLorentzVector pEMCAL;
2541   TLorentzVector gammaEMCAL;
2542   TLorentzVector pi0GammaConvEMCAL;
2543   TLorentzVector pi0GammaConvEMCALBck;
2544
2545   for (Int_t i=0; i<fV0Reader->GetESDEvent()->GetNumberOfCaloClusters(); i++) {
2546   clu = fV0Reader->GetESDEvent()->GetCaloCluster(i);
2547   if ( !clu->IsEMCAL()  || clu->E()<0.1 ) continue;
2548   if (clu->GetNCells() <= 1) continue;
2549   if ( clu->GetTOF()*1e9 < 550  || clu->GetTOF()*1e9 > 750) continue;
2550
2551   clu ->GetMomentum(pEMCAL ,vtx);
2552   for(Int_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();firstGammaIndex++){
2553   AliKFParticle * twoGammaDecayCandidateDaughter0 = (AliKFParticle *)fKFReconstructedGammasTClone->At(firstGammaIndex);
2554   gammaGammaConv.SetXYZM(twoGammaDecayCandidateDaughter0->Px(),
2555   twoGammaDecayCandidateDaughter0->Py(),
2556   twoGammaDecayCandidateDaughter0->Pz(),0.);
2557   gammaEMCAL.SetXYZM(pEMCAL.Px(),pEMCAL.Py(),pEMCAL.Pz(),0.);
2558   pi0GammaConvEMCAL=gammaGammaConv+gammaEMCAL;
2559   fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvEMCAL",pi0GammaConvEMCAL.M());
2560   fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_GammaConvEMCAL",pi0GammaConvEMCAL.M(),pi0GammaConvEMCAL.Pt());
2561   TVector3 v3D0(twoGammaDecayCandidateDaughter0->Px(),
2562   twoGammaDecayCandidateDaughter0->Py(),
2563   twoGammaDecayCandidateDaughter0->Pz());
2564   TVector3 v3D1(gammaEMCAL.Px(),gammaEMCAL.Py(),gammaEMCAL.Pz());
2565
2566
2567   Double_t opanConvEMCAL= v3D0.Angle(v3D1);
2568   if ( opanConvEMCAL < 0.35){
2569   fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvEMCAL_OpanLow",pi0GammaConvEMCAL.M());
2570   }else{
2571   fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvEMCAL_OpanHigh",pi0GammaConvEMCAL.M());
2572   }
2573
2574   }
2575   if(fCalculateBackground){
2576   for(Int_t nEventsInBG=0;nEventsInBG <fV0Reader->GetNBGEvents();nEventsInBG++){
2577   AliGammaConversionKFVector * previousEventV0s = fV0Reader->GetBGGoodV0s(nEventsInBG);
2578   for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
2579   AliKFParticle previousGoodV0 = (AliKFParticle)(*(previousEventV0s->at(iPrevious)));
2580   gammaGammaConvBck.SetXYZM(previousGoodV0.Px(),
2581   previousGoodV0.Py(),
2582   previousGoodV0.Pz(),0.);
2583   pi0GammaConvEMCALBck=gammaGammaConvBck+gammaEMCAL;
2584   fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvEMCAL_Bck",pi0GammaConvEMCALBck.M());
2585   fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_GammaConvEMCAL_Bck",pi0GammaConvEMCALBck.M(),
2586   pi0GammaConvEMCALBck.Pt());
2587   }
2588   }
2589       
2590   //  Now the LorentVector pEMCAL is obtained and can be paired with the converted proton
2591   } // end of checking if background photons are available
2592   }
2593   //==== End of the PHOS cluster selection ============
2594
2595   }
2596 */
2597
2598 void AliAnalysisTaskGammaConversion::MoveParticleAccordingToVertex(AliKFParticle * particle,const AliGammaConversionBGHandler::GammaConversionVertex *vertex){
2599   //see header file for documentation
2600
2601    Double_t dx = vertex->fX - fESDEvent->GetPrimaryVertex()->GetX();
2602    Double_t dy = vertex->fY - fESDEvent->GetPrimaryVertex()->GetY();
2603    Double_t dz = vertex->fZ - fESDEvent->GetPrimaryVertex()->GetZ();
2604   
2605   //  cout<<"dx, dy, dz: ["<<dx<<","<<dy<<","<<dz<<"]"<<endl;
2606    particle->X() = particle->GetX() - dx;
2607    particle->Y() = particle->GetY() - dy;
2608    particle->Z() = particle->GetZ() - dz;
2609 }
2610
2611 void AliAnalysisTaskGammaConversion::RotateKFParticle(AliKFParticle * kfParticle,Double_t angle){
2612   // Rotate the kf particle
2613   Double_t c = cos(angle);
2614   Double_t s = sin(angle);
2615   
2616   Double_t mA[7][ 7];
2617   for( Int_t i=0; i<7; i++ ){
2618     for( Int_t j=0; j<7; j++){
2619       mA[i][j] = 0;
2620     }
2621   }
2622   for( int i=0; i<7; i++ ){
2623     mA[i][i] = 1;
2624   }
2625   mA[0][0] =  c;  mA[0][1] = s;
2626   mA[1][0] = -s;  mA[1][1] = c;
2627   mA[3][3] =  c;  mA[3][4] = s;
2628   mA[4][3] = -s;  mA[4][4] = c;
2629   
2630   Double_t mAC[7][7];
2631   Double_t mAp[7];
2632   
2633   for( Int_t i=0; i<7; i++ ){
2634     mAp[i] = 0;
2635     for( Int_t k=0; k<7; k++){
2636       mAp[i]+=mA[i][k] * kfParticle->GetParameter(k);
2637     }
2638   }
2639   
2640   for( Int_t i=0; i<7; i++){
2641     kfParticle->Parameter(i) = mAp[i];
2642   }
2643
2644   for( Int_t i=0; i<7; i++ ){
2645     for( Int_t j=0; j<7; j++ ){
2646       mAC[i][j] = 0;
2647       for( Int_t k=0; k<7; k++ ){
2648         mAC[i][j]+= mA[i][k] * kfParticle->GetCovariance(k,j);
2649       }
2650     }
2651   }
2652
2653   for( Int_t i=0; i<7; i++ ){
2654     for( Int_t j=0; j<=i; j++ ){
2655       Double_t xx = 0;
2656       for( Int_t k=0; k<7; k++){
2657         xx+= mAC[i][k]*mA[j][k];
2658       }
2659       kfParticle->Covariance(i,j) = xx;
2660     }
2661   }
2662 }
2663
2664
2665 void AliAnalysisTaskGammaConversion::CalculateBackground(){
2666   // see header file for documentation
2667
2668
2669   TClonesArray * currentEventV0s = fV0Reader->GetCurrentEventGoodV0s();
2670
2671   AliGammaConversionBGHandler * bgHandler = fV0Reader->GetBGHandler();
2672   
2673   Int_t zbin= bgHandler->GetZBinIndex(fV0Reader->GetVertexZ());
2674   Int_t mbin = 0;
2675   if(fUseTrackMultiplicityForBG == kTRUE){
2676     mbin = bgHandler->GetMultiplicityBinIndex(fV0Reader->CountESDTracks());
2677   }
2678   else{
2679     mbin = bgHandler->GetMultiplicityBinIndex(fV0Reader->GetNGoodV0s());
2680   }
2681
2682   if(fDoRotation == kTRUE){
2683     TRandom3 *random = new TRandom3();
2684     for(Int_t iCurrent=0;iCurrent<currentEventV0s->GetEntriesFast();iCurrent++){
2685       AliKFParticle currentEventGoodV0 = *(AliKFParticle *)(currentEventV0s->At(iCurrent)); 
2686       for(Int_t iCurrent2=iCurrent+1;iCurrent2<currentEventV0s->GetEntriesFast();iCurrent2++){
2687         for(Int_t nRandom=0;nRandom<fNRandomEventsForBG;nRandom++){
2688         
2689           AliKFParticle currentEventGoodV02 = *(AliKFParticle *)(currentEventV0s->At(iCurrent2));
2690
2691           if(fCheckBGProbability == kTRUE){
2692             Double_t massBGprob =0.;
2693             Double_t widthBGprob = 0.;
2694             AliKFParticle *backgroundCandidateProb = new AliKFParticle(currentEventGoodV0,currentEventGoodV02);
2695             backgroundCandidateProb->GetMass(massBGprob,widthBGprob);
2696             if(massBGprob>0.1 && massBGprob<0.14){
2697               if(random->Rndm()>bgHandler->GetBGProb(zbin,mbin)){
2698                 delete backgroundCandidateProb;
2699                 continue;
2700               }
2701             }
2702             delete backgroundCandidateProb;
2703           }
2704         
2705           Double_t nRadiansPM = fNDegreesPMBackground*TMath::Pi()/180;
2706           
2707           Double_t rotationValue = random->Rndm()*2*nRadiansPM + TMath::Pi()-nRadiansPM;
2708
2709           RotateKFParticle(&currentEventGoodV02,rotationValue);
2710
2711           AliKFParticle *backgroundCandidate = new AliKFParticle(currentEventGoodV0,currentEventGoodV02);
2712
2713           Double_t massBG =0.;
2714           Double_t widthBG = 0.;
2715           Double_t chi2BG =10000.;      
2716           backgroundCandidate->GetMass(massBG,widthBG);
2717
2718           //      if(backgroundCandidate->GetNDF()>0){
2719           chi2BG = backgroundCandidate->GetChi2();
2720           if((chi2BG>0 && chi2BG<fV0Reader->GetChi2CutMeson())  || fApplyChi2Cut == kFALSE){
2721           
2722             TVector3 momentumVectorbackgroundCandidate(backgroundCandidate->GetPx(),backgroundCandidate->GetPy(),backgroundCandidate->GetPz());
2723             TVector3 spaceVectorbackgroundCandidate(backgroundCandidate->GetX(),backgroundCandidate->GetY(),backgroundCandidate->GetZ());
2724           
2725             Double_t openingAngleBG = currentEventGoodV0.GetAngle(currentEventGoodV02);
2726           
2727             Double_t rapidity;
2728             if(backgroundCandidate->GetE() - backgroundCandidate->GetPz() == 0 || backgroundCandidate->GetE() + backgroundCandidate->GetPz() == 0) rapidity=0;
2729             else rapidity = 0.5*(TMath::Log((backgroundCandidate->GetE() +backgroundCandidate->GetPz()) / (backgroundCandidate->GetE()-backgroundCandidate->GetPz())));
2730           
2731             if(TMath::Abs(rapidity) > fV0Reader->GetRapidityMesonCut() ){
2732               delete backgroundCandidate;   
2733               continue;   // rapidity cut
2734             }                   
2735                                         
2736           
2737             Double_t alfa=0.0;
2738             if( (currentEventGoodV0.GetE()+currentEventGoodV02.GetE()) != 0){
2739               alfa=TMath::Abs((currentEventGoodV0.GetE()-currentEventGoodV02.GetE())
2740                               /(currentEventGoodV0.GetE()+currentEventGoodV02.GetE()));
2741             }
2742           
2743           
2744             if(openingAngleBG < fMinOpeningAngleGhostCut ){
2745               delete backgroundCandidate;   
2746               continue;   // minimum opening angle to avoid using ghosttracks
2747             }                   
2748           
2749             // original
2750             if(alfa>fV0Reader->GetAlphaMinCutMeson() &&   alfa<fV0Reader->GetAlphaCutMeson()){
2751               fHistograms->FillHistogram("ESD_Background_GammaDaughter_OpeningAngle", openingAngleBG);
2752               fHistograms->FillHistogram("ESD_Background_Energy", backgroundCandidate->GetE());
2753               fHistograms->FillHistogram("ESD_Background_Pt",  momentumVectorbackgroundCandidate.Pt());
2754               fHistograms->FillHistogram("ESD_Background_Eta", momentumVectorbackgroundCandidate.Eta());
2755               fHistograms->FillHistogram("ESD_Background_Rapidity", rapidity);
2756               fHistograms->FillHistogram("ESD_Background_Phi", spaceVectorbackgroundCandidate.Phi());
2757               fHistograms->FillHistogram("ESD_Background_Mass", massBG);
2758               fHistograms->FillHistogram("ESD_Background_R", spaceVectorbackgroundCandidate.Pt());  // Pt in Space == R!!!!
2759               fHistograms->FillHistogram("ESD_Background_ZR", backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
2760               fHistograms->FillHistogram("ESD_Background_XY", backgroundCandidate->GetX(), backgroundCandidate->GetY());
2761               fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt",massBG,momentumVectorbackgroundCandidate.Pt());
2762               fHistograms->FillHistogram("ESD_Background_InvMass",massBG);
2763               fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_alpha",massBG,momentumVectorbackgroundCandidate.Pt());
2764
2765               if(massBG>0.1 && massBG<0.15){
2766                 fHistograms->FillHistogram("ESD_Background_alfa_Pi0", alfa);
2767               }
2768               if(massBG>0.5 && massBG<0.57){
2769                 fHistograms->FillHistogram("ESD_Background_alfa_Eta", alfa);
2770               }
2771
2772               if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 &&  TMath::Abs(currentEventGoodV02.GetEta())<0.9 ){
2773                 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_Fiducial",massBG,momentumVectorbackgroundCandidate.Pt());
2774                 fHistograms->FillHistogram("ESD_Background_InvMass_Fiducial",massBG);
2775               }
2776               
2777               fHistograms->FillHistogram(Form("%d%dESD_Background_GammaDaughter_OpeningAngle",zbin,mbin), openingAngleBG);
2778               fHistograms->FillHistogram(Form("%d%dESD_Background_Energy",zbin,mbin), backgroundCandidate->GetE());
2779               fHistograms->FillHistogram(Form("%d%dESD_Background_Pt",zbin,mbin),  momentumVectorbackgroundCandidate.Pt());
2780               fHistograms->FillHistogram(Form("%d%dESD_Background_Eta",zbin,mbin), momentumVectorbackgroundCandidate.Eta());
2781               fHistograms->FillHistogram(Form("%d%dESD_Background_Rapidity",zbin,mbin), rapidity);
2782               fHistograms->FillHistogram(Form("%d%dESD_Background_Phi",zbin,mbin), spaceVectorbackgroundCandidate.Phi());
2783               fHistograms->FillHistogram(Form("%d%dESD_Background_Mass",zbin,mbin), massBG);
2784               fHistograms->FillHistogram(Form("%d%dESD_Background_R",zbin,mbin), spaceVectorbackgroundCandidate.Pt());  // Pt in Space == R!!!!
2785               fHistograms->FillHistogram(Form("%d%dESD_Background_ZR",zbin,mbin), backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
2786               fHistograms->FillHistogram(Form("%d%dESD_Background_XY",zbin,mbin), backgroundCandidate->GetX(), backgroundCandidate->GetY());
2787               fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_vs_Pt",zbin,mbin),massBG,momentumVectorbackgroundCandidate.Pt());
2788               fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass",zbin,mbin),massBG);
2789               
2790               if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 &&  TMath::Abs(currentEventGoodV02.GetEta())<0.9 ){
2791                 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_vs_Pt_Fiducial",zbin,mbin),massBG,momentumVectorbackgroundCandidate.Pt());
2792                 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_Fiducial",zbin,mbin),massBG);
2793               }
2794             }
2795             if(alfa<0.1){
2796               fHistograms->FillHistogram("ESD_Background_InvMass_vs_E_alpha",massBG ,backgroundCandidate->GetE());
2797             }
2798
2799           }
2800           //}
2801           delete backgroundCandidate;      
2802         }
2803       }
2804     }
2805     delete random;
2806   }
2807   else{ // means no rotation
2808     AliGammaConversionBGHandler::GammaConversionVertex *bgEventVertex = NULL;
2809             
2810     if(fUseTrackMultiplicityForBG){
2811       //    cout<<"Using charged track multiplicity for background calculation"<<endl;
2812       for(Int_t nEventsInBG=0;nEventsInBG <fV0Reader->GetNBGEvents();nEventsInBG++){
2813
2814         AliGammaConversionKFVector * previousEventV0s = bgHandler->GetBGGoodV0s(zbin,mbin,nEventsInBG);//fV0Reader->GetBGGoodV0s(nEventsInBG);
2815       
2816         if(fMoveParticleAccordingToVertex == kTRUE){
2817           bgEventVertex = bgHandler->GetBGEventVertex(zbin,mbin,nEventsInBG);
2818         }
2819
2820         for(Int_t iCurrent=0;iCurrent<currentEventV0s->GetEntriesFast();iCurrent++){
2821           AliKFParticle currentEventGoodV0 = *(AliKFParticle *)(currentEventV0s->At(iCurrent)); 
2822           for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
2823             AliKFParticle previousGoodV0 = (AliKFParticle)(*(previousEventV0s->at(iPrevious)));
2824             AliKFParticle previousGoodV0test = (AliKFParticle)(*(previousEventV0s->at(iPrevious)));
2825
2826             //cout<<"Primary Vertex event: ["<<fESDEvent->GetPrimaryVertex()->GetX()<<","<<fESDEvent->GetPrimaryVertex()->GetY()<<","<<fESDEvent->GetPrimaryVertex()->GetZ()<<"]"<<endl;
2827             //cout<<"BG prim Vertex event: ["<<bgEventVertex->fX<<","<<bgEventVertex->fY<<","<<bgEventVertex->fZ<<"]"<<endl;
2828           
2829             //cout<<"XYZ of particle before transport: ["<<previousGoodV0.X()<<","<<previousGoodV0.Y()<<","<<previousGoodV0.Z()<<"]"<<endl;
2830             if(fMoveParticleAccordingToVertex == kTRUE){
2831               MoveParticleAccordingToVertex(&previousGoodV0,bgEventVertex);
2832             }
2833             //cout<<"XYZ of particle after transport: ["<<previousGoodV0.X()<<","<<previousGoodV0.Y()<<","<<previousGoodV0.Z()<<"]"<<endl;
2834
2835             AliKFParticle *backgroundCandidate = new AliKFParticle(currentEventGoodV0,previousGoodV0);
2836         
2837             Double_t massBG =0.;
2838             Double_t widthBG = 0.;
2839             Double_t chi2BG =10000.;    
2840             backgroundCandidate->GetMass(massBG,widthBG);
2841             //  if(backgroundCandidate->GetNDF()>0){
2842             //    chi2BG = backgroundCandidate->GetChi2()/backgroundCandidate->GetNDF();
2843             chi2BG = backgroundCandidate->GetChi2();
2844             if((chi2BG>0 && chi2BG<fV0Reader->GetChi2CutMeson()) || fApplyChi2Cut == kFALSE){
2845                                         
2846               TVector3 momentumVectorbackgroundCandidate(backgroundCandidate->GetPx(),backgroundCandidate->GetPy(),backgroundCandidate->GetPz());
2847               TVector3 spaceVectorbackgroundCandidate(backgroundCandidate->GetX(),backgroundCandidate->GetY(),backgroundCandidate->GetZ());
2848                                         
2849               Double_t openingAngleBG = currentEventGoodV0.GetAngle(previousGoodV0);
2850                                         
2851               Double_t rapidity;
2852             
2853               if(backgroundCandidate->GetE() - backgroundCandidate->GetPz() <= 0 || backgroundCandidate->GetE() + backgroundCandidate->GetPz() <= 0){
2854                 cout << "Error: |Pz| > E !!!! " << endl;
2855                 rapidity=0;
2856               } else {
2857                 rapidity = 0.5*(TMath::Log((backgroundCandidate->GetE() +backgroundCandidate->GetPz()) / (backgroundCandidate->GetE()-backgroundCandidate->GetPz())));
2858               }                         
2859               if(TMath::Abs(rapidity) > fV0Reader->GetRapidityMesonCut() ){
2860                 delete backgroundCandidate;   
2861                 continue;   // rapidity cut
2862               }                 
2863                                                         
2864         
2865               Double_t alfa=0.0;
2866               if( (currentEventGoodV0.GetE()+previousGoodV0.GetE()) != 0){
2867                 alfa=TMath::Abs((currentEventGoodV0.GetE()-previousGoodV0.GetE())
2868                                 /(currentEventGoodV0.GetE()+previousGoodV0.GetE()));
2869               }
2870                         
2871                                         
2872               if(openingAngleBG < fMinOpeningAngleGhostCut ){
2873                 delete backgroundCandidate;   
2874                 continue;   // minimum opening angle to avoid using ghosttracks
2875               }                 
2876
2877               // original
2878               if(alfa>fV0Reader->GetAlphaMinCutMeson() &&   alfa<fV0Reader->GetAlphaCutMeson()){
2879                 fHistograms->FillHistogram("ESD_Background_GammaDaughter_OpeningAngle", openingAngleBG);
2880                 fHistograms->FillHistogram("ESD_Background_Energy", backgroundCandidate->GetE());
2881                 fHistograms->FillHistogram("ESD_Background_Pt",  momentumVectorbackgroundCandidate.Pt());
2882                 fHistograms->FillHistogram("ESD_Background_Eta", momentumVectorbackgroundCandidate.Eta());
2883                 fHistograms->FillHistogram("ESD_Background_Rapidity", rapidity);
2884                 fHistograms->FillHistogram("ESD_Background_Phi", spaceVectorbackgroundCandidate.Phi());
2885                 fHistograms->FillHistogram("ESD_Background_Mass", massBG);
2886                 fHistograms->FillHistogram("ESD_Background_R", spaceVectorbackgroundCandidate.Pt());  // Pt in Space == R!!!!
2887                 fHistograms->FillHistogram("ESD_Background_ZR", backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
2888                 fHistograms->FillHistogram("ESD_Background_XY", backgroundCandidate->GetX(), backgroundCandidate->GetY());
2889                 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt",massBG,momentumVectorbackgroundCandidate.Pt());
2890                 fHistograms->FillHistogram("ESD_Background_InvMass",massBG);
2891                 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_alpha",massBG,momentumVectorbackgroundCandidate.Pt());
2892
2893                 if(massBG>0.1 && massBG<0.15){
2894                   fHistograms->FillHistogram("ESD_Background_alfa_Pi0", alfa);
2895                 }
2896                 if(massBG>0.5 && massBG<0.57){
2897                   fHistograms->FillHistogram("ESD_Background_alfa_Eta", alfa);
2898                 }
2899
2900                 if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 &&  TMath::Abs(previousGoodV0.GetEta())<0.9 ){
2901                   fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_Fiducial",massBG,momentumVectorbackgroundCandidate.Pt());
2902                   fHistograms->FillHistogram("ESD_Background_InvMass_Fiducial",massBG);
2903                 }
2904
2905                 // test
2906                 fHistograms->FillHistogram(Form("%d%dESD_Background_GammaDaughter_OpeningAngle",zbin,mbin), openingAngleBG);
2907                 fHistograms->FillHistogram(Form("%d%dESD_Background_Energy",zbin,mbin), backgroundCandidate->GetE());
2908                 fHistograms->FillHistogram(Form("%d%dESD_Background_Pt",zbin,mbin),  momentumVectorbackgroundCandidate.Pt());
2909                 fHistograms->FillHistogram(Form("%d%dESD_Background_Eta",zbin,mbin), momentumVectorbackgroundCandidate.Eta());
2910                 fHistograms->FillHistogram(Form("%d%dESD_Background_Rapidity",zbin,mbin), rapidity);
2911                 fHistograms->FillHistogram(Form("%d%dESD_Background_Phi",zbin,mbin), spaceVectorbackgroundCandidate.Phi());
2912                 fHistograms->FillHistogram(Form("%d%dESD_Background_Mass",zbin,mbin), massBG);
2913                 fHistograms->FillHistogram(Form("%d%dESD_Background_R",zbin,mbin), spaceVectorbackgroundCandidate.Pt());  // Pt in Space == R!!!!
2914                 fHistograms->FillHistogram(Form("%d%dESD_Background_ZR",zbin,mbin), backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
2915                 fHistograms->FillHistogram(Form("%d%dESD_Background_XY",zbin,mbin), backgroundCandidate->GetX(), backgroundCandidate->GetY());
2916                 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_vs_Pt",zbin,mbin),massBG,momentumVectorbackgroundCandidate.Pt());
2917                 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass",zbin,mbin),massBG);
2918                 
2919                 if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 &&  TMath::Abs(previousGoodV0.GetEta())<0.9 ){
2920                   fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_vs_Pt_Fiducial",zbin,mbin),massBG,momentumVectorbackgroundCandidate.Pt());
2921                   fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_Fiducial",zbin,mbin),massBG);
2922                 }
2923                 //        }
2924               }
2925               if(alfa<0.1){
2926                 fHistograms->FillHistogram("ESD_Background_InvMass_vs_E_alpha",massBG ,backgroundCandidate->GetE());
2927               }
2928
2929             }
2930             delete backgroundCandidate;      
2931           }
2932         }
2933       }
2934     }
2935     else{ // means using #V0s for multiplicity
2936
2937       //    cout<<"Using the v0 multiplicity to calculate background"<<endl;
2938     
2939       fHistograms->FillHistogram("ESD_Background_z_m",zbin,mbin);
2940       fHistograms->FillHistogram("ESD_Mother_multpilicityVSv0s",fV0Reader->CountESDTracks(),fV0Reader->GetNumberOfV0s());
2941
2942       for(Int_t nEventsInBG=0;nEventsInBG <fV0Reader->GetNBGEvents();nEventsInBG++){
2943         AliGammaConversionKFVector * previousEventV0s = bgHandler->GetBGGoodV0s(zbin,mbin,nEventsInBG);// fV0Reader->GetBGGoodV0s(nEventsInBG);
2944         if(previousEventV0s){
2945         
2946           if(fMoveParticleAccordingToVertex == kTRUE){
2947             bgEventVertex = bgHandler->GetBGEventVertex(zbin,mbin,nEventsInBG);
2948           }
2949
2950           for(Int_t iCurrent=0;iCurrent<currentEventV0s->GetEntriesFast();iCurrent++){
2951             AliKFParticle currentEventGoodV0 = *(AliKFParticle *)(currentEventV0s->At(iCurrent)); 
2952             for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
2953               AliKFParticle previousGoodV0 = (AliKFParticle)(*(previousEventV0s->at(iPrevious)));
2954
2955               if(fMoveParticleAccordingToVertex == kTRUE){
2956                 MoveParticleAccordingToVertex(&previousGoodV0,bgEventVertex);
2957               }
2958
2959               AliKFParticle *backgroundCandidate = new AliKFParticle(currentEventGoodV0,previousGoodV0);
2960               Double_t massBG =0.;
2961               Double_t widthBG = 0.;
2962               Double_t chi2BG =10000.;  
2963               backgroundCandidate->GetMass(massBG,widthBG);
2964               /*            if(backgroundCandidate->GetNDF()>0){
2965                             chi2BG = backgroundCandidate->GetChi2()/backgroundCandidate->GetNDF();
2966                             {//remember to remove
2967                             TVector3 momentumVectorbackgroundCandidate(backgroundCandidate->GetPx(),backgroundCandidate->GetPy(),backgroundCandidate->GetPz());
2968                             TVector3 spaceVectorbackgroundCandidate(backgroundCandidate->GetX(),backgroundCandidate->GetY(),backgroundCandidate->GetZ());
2969               
2970                             Double_t openingAngleBG = currentEventGoodV0.GetAngle(previousGoodV0);
2971                             fHistograms->FillHistogram("ESD_Background_GammaDaughter_OpeningAngle_nochi2", openingAngleBG);
2972                             }
2973               */
2974               chi2BG = backgroundCandidate->GetChi2();
2975               if((chi2BG>0 && chi2BG<fV0Reader->GetChi2CutMeson()) || fApplyChi2Cut == kFALSE){
2976                 TVector3 momentumVectorbackgroundCandidate(backgroundCandidate->GetPx(),backgroundCandidate->GetPy(),backgroundCandidate->GetPz());
2977                 TVector3 spaceVectorbackgroundCandidate(backgroundCandidate->GetX(),backgroundCandidate->GetY(),backgroundCandidate->GetZ());
2978                                         
2979                 Double_t openingAngleBG = currentEventGoodV0.GetAngle(previousGoodV0);
2980                                         
2981                 Double_t rapidity;
2982                 if(backgroundCandidate->GetE() - backgroundCandidate->GetPz() == 0 || backgroundCandidate->GetE() + backgroundCandidate->GetPz() == 0) rapidity=0;
2983                 else rapidity = 0.5*(TMath::Log((backgroundCandidate->GetE() +backgroundCandidate->GetPz()) / (backgroundCandidate->GetE()-backgroundCandidate->GetPz())));
2984                                         
2985                 if(TMath::Abs(rapidity) > fV0Reader->GetRapidityMesonCut() ){
2986                   delete backgroundCandidate;   
2987                   continue;   // rapidity cut
2988                 }                       
2989                                                                 
2990
2991                 Double_t alfa=0.0;
2992                 if( (currentEventGoodV0.GetE()+previousGoodV0.GetE()) != 0){
2993                   alfa=TMath::Abs((currentEventGoodV0.GetE()-previousGoodV0.GetE())
2994                                   /(currentEventGoodV0.GetE()+previousGoodV0.GetE()));
2995                 }
2996                         
2997                                         
2998                 if(openingAngleBG < fMinOpeningAngleGhostCut ){
2999                   delete backgroundCandidate;   
3000                   continue;   // minimum opening angle to avoid using ghosttracks
3001                 }                       
3002
3003                 if(alfa>fV0Reader->GetAlphaMinCutMeson() &&   alfa<fV0Reader->GetAlphaCutMeson()){
3004                   fHistograms->FillHistogram("ESD_Background_GammaDaughter_OpeningAngle", openingAngleBG);
3005                   fHistograms->FillHistogram("ESD_Background_Energy", backgroundCandidate->GetE());
3006                   fHistograms->FillHistogram("ESD_Background_Pt",  momentumVectorbackgroundCandidate.Pt());
3007                   fHistograms->FillHistogram("ESD_Background_Eta", momentumVectorbackgroundCandidate.Eta());
3008                   fHistograms->FillHistogram("ESD_Background_Rapidity", rapidity);
3009                   fHistograms->FillHistogram("ESD_Background_Phi", spaceVectorbackgroundCandidate.Phi());
3010                   fHistograms->FillHistogram("ESD_Background_Mass", massBG);
3011                   fHistograms->FillHistogram("ESD_Background_R", spaceVectorbackgroundCandidate.Pt());  // Pt in Space == R!!!!
3012                   fHistograms->FillHistogram("ESD_Background_ZR", backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
3013                   fHistograms->FillHistogram("ESD_Background_XY", backgroundCandidate->GetX(), backgroundCandidate->GetY());
3014                   fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt",massBG,momentumVectorbackgroundCandidate.Pt());
3015                   fHistograms->FillHistogram("ESD_Background_InvMass",massBG);
3016                   
3017
3018                   fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_alpha",massBG,momentumVectorbackgroundCandidate.Pt());
3019
3020                   if(massBG>0.1 && massBG<0.15){
3021                     fHistograms->FillHistogram("ESD_Background_alfa_Pi0", alfa);
3022                   }
3023                   if(massBG>0.5 && massBG<0.57){
3024                     fHistograms->FillHistogram("ESD_Background_alfa_Eta", alfa);
3025                   }
3026
3027                   if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 &&  TMath::Abs(previousGoodV0.GetEta())<0.9 ){
3028                     fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_Fiducial",massBG,momentumVectorbackgroundCandidate.Pt());
3029                     fHistograms->FillHistogram("ESD_Background_InvMass_Fiducial",massBG);
3030                   }
3031                   
3032                   if(massBG>0.5 && massBG<0.6){
3033                     fHistograms->FillHistogram("ESD_Background_alfa_pt0506",momentumVectorbackgroundCandidate.Pt(),alfa);
3034                   }
3035                   if(massBG>0.3 && massBG<0.4){
3036                     fHistograms->FillHistogram("ESD_Background_alfa_pt0304",momentumVectorbackgroundCandidate.Pt(),alfa);
3037                   }
3038                   
3039                   // test
3040                   fHistograms->FillHistogram(Form("%d%dESD_Background_GammaDaughter_OpeningAngle",zbin,mbin), openingAngleBG);
3041                   fHistograms->FillHistogram(Form("%d%dESD_Background_Energy",zbin,mbin), backgroundCandidate->GetE());
3042                   fHistograms->FillHistogram(Form("%d%dESD_Background_Pt",zbin,mbin),  momentumVectorbackgroundCandidate.Pt());
3043                   fHistograms->FillHistogram(Form("%d%dESD_Background_Eta",zbin,mbin), momentumVectorbackgroundCandidate.Eta());
3044                   fHistograms->FillHistogram(Form("%d%dESD_Background_Rapidity",zbin,mbin), rapidity);
3045                   fHistograms->FillHistogram(Form("%d%dESD_Background_Phi",zbin,mbin), spaceVectorbackgroundCandidate.Phi());
3046                   fHistograms->FillHistogram(Form("%d%dESD_Background_Mass",zbin,mbin), massBG);
3047                   fHistograms->FillHistogram(Form("%d%dESD_Background_R",zbin,mbin), spaceVectorbackgroundCandidate.Pt());  // Pt in Space == R!!!!
3048                   fHistograms->FillHistogram(Form("%d%dESD_Background_ZR",zbin,mbin), backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
3049                   fHistograms->FillHistogram(Form("%d%dESD_Background_XY",zbin,mbin), backgroundCandidate->GetX(), backgroundCandidate->GetY());
3050                   fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_vs_Pt",zbin,mbin),massBG,momentumVectorbackgroundCandidate.Pt());
3051                   fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass",zbin,mbin),massBG);
3052                   
3053                   if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 &&  TMath::Abs(previousGoodV0.GetEta())<0.9 ){
3054                     fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_vs_Pt_Fiducial",zbin,mbin),massBG,momentumVectorbackgroundCandidate.Pt());
3055                     fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_Fiducial",zbin,mbin),massBG);
3056                   }
3057                 }
3058
3059                 if(alfa<0.1){
3060                   fHistograms->FillHistogram("ESD_Background_InvMass_vs_E_alpha",massBG ,backgroundCandidate->GetE());
3061                 }
3062                 //  }
3063               }
3064               delete backgroundCandidate;      
3065             }
3066           }
3067         }
3068       }
3069     } // end else (means use #v0s as multiplicity)
3070   } // end no rotation
3071 }
3072
3073
3074 void AliAnalysisTaskGammaConversion::ProcessGammasForGammaJetAnalysis(){
3075   //ProcessGammasForGammaJetAnalysis
3076         
3077   Double_t distIsoMin;
3078         
3079   CreateListOfChargedParticles();
3080         
3081         
3082   //  for(UInt_t gammaIndex=0;gammaIndex<fKFReconstructedGammas.size();gammaIndex++){
3083   for(Int_t gammaIndex=0;gammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();gammaIndex++){
3084     AliKFParticle * currentGamma = (AliKFParticle*)fKFReconstructedGammasTClone->At(gammaIndex);
3085     TVector3 momentumVectorCurrentGamma(currentGamma->GetPx(),currentGamma->GetPy(),currentGamma->GetPz());
3086                 
3087     if( momentumVectorCurrentGamma.Pt()> fMinPtForGammaJet){
3088       distIsoMin=GetMinimumDistanceToCharge(gammaIndex);
3089                         
3090       if (distIsoMin > fMinIsoConeSize && fLeadingChargedIndex>=0){
3091         CalculateJetCone(gammaIndex);
3092       }
3093     }
3094   }
3095 }
3096
3097 //____________________________________________________________________
3098 Bool_t AliAnalysisTaskGammaConversion::IsGoodImpPar(const AliESDtrack *const track)
3099 {
3100 //
3101 // check whether particle has good DCAr(Pt) impact
3102 // parameter. Only for TPC+ITS tracks (7*sigma cut)
3103 // Origin: Andrea Dainese
3104 //
3105
3106 Float_t d0z0[2],covd0z0[3];
3107 track->GetImpactParameters(d0z0,covd0z0);
3108 Float_t sigma= 0.0050+0.0060/TMath::Power(track->Pt(),0.9);
3109 Float_t d0max = 7.*sigma;
3110 if(TMath::Abs(d0z0[0]) < d0max) return kTRUE;
3111
3112 return kFALSE;
3113 }
3114
3115
3116 void AliAnalysisTaskGammaConversion::CreateListOfChargedParticles(){
3117   // CreateListOfChargedParticles
3118         
3119   fESDEvent = fV0Reader->GetESDEvent();
3120   Int_t numberOfESDTracks=0;
3121   for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
3122     AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);
3123                 
3124     if(!curTrack){
3125       continue;
3126     }
3127     // Not needed if Standard function used.
3128 //     if(!IsGoodImpPar(curTrack)){
3129 //       continue;
3130 //     }
3131                 
3132     if(fEsdTrackCuts->AcceptTrack(curTrack) ){
3133       new((*fChargedParticles)[fChargedParticles->GetEntriesFast()])  AliESDtrack(*curTrack);
3134       //      fChargedParticles.push_back(curTrack);
3135       fChargedParticlesId.push_back(iTracks);
3136       numberOfESDTracks++;
3137     }
3138   }
3139 // Moved to UserExec using CountAcceptedTracks function. runjet is not needed by default
3140 //   fHistograms->FillHistogram("ESD_NumberOfGoodESDTracks",numberOfESDTracks);
3141 //   cout<<"esdtracks::"<< numberOfESDTracks<<endl;
3142 //   if (fV0Reader->GetNumberOfContributorsVtx()>=1){
3143 //     fHistograms->FillHistogram("ESD_NumberOfGoodESDTracksVtx",numberOfESDTracks);
3144 //   } 
3145 }
3146 void AliAnalysisTaskGammaConversion::RecalculateV0ForGamma(){
3147   //recalculates v0 for gamma
3148
3149  Double_t massE=0.00051099892;
3150  TLorentzVector curElecPos;
3151  TLorentzVector curElecNeg;
3152  TLorentzVector curGamma;
3153
3154  TLorentzVector curGammaAt;
3155  TLorentzVector curElecPosAt;
3156  TLorentzVector curElecNegAt;
3157  AliKFVertex primVtxGamma(*(fESDEvent->GetPrimaryVertex()));
3158  AliKFVertex primVtxImprovedGamma = primVtxGamma;
3159
3160  const AliESDVertex *vtxT3D=fESDEvent->GetPrimaryVertex();
3161
3162  Double_t xPrimaryVertex=vtxT3D->GetXv();
3163  Double_t yPrimaryVertex=vtxT3D->GetYv();
3164  Double_t zPrimaryVertex=vtxT3D->GetZv();
3165  // Float_t primvertex[3]={xPrimaryVertex,yPrimaryVertex,zPrimaryVertex};
3166
3167  Float_t nsigmaTPCtrackPos;
3168  Float_t nsigmaTPCtrackNeg;
3169  Float_t nsigmaTPCtrackPosToPion;
3170  Float_t nsigmaTPCtrackNegToPion;
3171  AliKFParticle* negKF=NULL;
3172  AliKFParticle* posKF=NULL;
3173
3174  for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
3175    AliESDtrack* posTrack = fESDEvent->GetTrack(iTracks);
3176    if(!posTrack){
3177      continue;
3178    }
3179    if (posKF) delete posKF; posKF=NULL;
3180    if(posTrack->GetSign()<0) continue;
3181    if(!(posTrack->GetStatus() & AliESDtrack::kTPCrefit))continue;
3182    if(posTrack->GetKinkIndex(0)>0 ) continue;
3183    if(posTrack->GetNcls(1)<50)continue;
3184    Double_t momPos[3];
3185    //    posTrack->GetConstrainedPxPyPz(momPos);
3186    posTrack->GetPxPyPz(momPos);
3187    AliESDtrack *ptrk=fESDEvent->GetTrack(iTracks);
3188    curElecPos.SetXYZM(momPos[0],momPos[1],momPos[2],massE);
3189    if(TMath::Abs(curElecPos.Eta())<0.9) continue;
3190    posKF = new AliKFParticle( *(posTrack),-11);
3191
3192    nsigmaTPCtrackPos = fV0Reader->GetESDpid()->NumberOfSigmasTPC(posTrack,AliPID::kElectron);
3193    nsigmaTPCtrackPosToPion = fV0Reader->GetESDpid()->NumberOfSigmasTPC(posTrack,AliPID::kPion);
3194
3195    if ( nsigmaTPCtrackPos>5.|| nsigmaTPCtrackPos<-2.){
3196      continue;
3197    }
3198   
3199    if(pow((momPos[0]*momPos[0]+momPos[1]*momPos[1]+momPos[2]*momPos[2]),0.5)>0.5 && nsigmaTPCtrackPosToPion<1){
3200      continue;
3201    }
3202
3203
3204
3205    for(Int_t jTracks = 0; jTracks < fESDEvent->GetNumberOfTracks(); jTracks++){
3206      AliESDtrack* negTrack = fESDEvent->GetTrack(jTracks);
3207      if(!negTrack){
3208        continue;
3209      }
3210      if (negKF) delete negKF; negKF=NULL;
3211      if(negTrack->GetSign()>0) continue;
3212      if(!(negTrack->GetStatus() & AliESDtrack::kTPCrefit))continue;
3213      if(negTrack->GetKinkIndex(0)>0 ) continue;
3214      if(negTrack->GetNcls(1)<50)continue;
3215      Double_t momNeg[3];
3216      //    negTrack->GetConstrainedPxPyPz(momNeg);
3217      negTrack->GetPxPyPz(momNeg);
3218
3219      nsigmaTPCtrackNeg = fV0Reader->GetESDpid()->NumberOfSigmasTPC(negTrack,AliPID::kElectron);     
3220      nsigmaTPCtrackNegToPion = fV0Reader->GetESDpid()->NumberOfSigmasTPC(negTrack,AliPID::kPion);
3221      if ( nsigmaTPCtrackNeg>5. || nsigmaTPCtrackNeg<-2.){
3222        continue;
3223      }
3224      if(pow((momNeg[0]*momNeg[0]+momNeg[1]*momNeg[1]+momNeg[2]*momNeg[2]),0.5)>0.5 && nsigmaTPCtrackNegToPion<1){
3225        continue;
3226      }
3227      AliESDtrack *ntrk=fESDEvent->GetTrack(jTracks);
3228      curElecNeg.SetXYZM(momNeg[0],momNeg[1],momNeg[2],massE);
3229      if(TMath::Abs(curElecNeg.Eta())<0.9) continue;
3230      negKF = new AliKFParticle( *(negTrack) ,11);
3231
3232      Double_t b=fESDEvent->GetMagneticField();
3233      Double_t xn, xp, dca=ntrk->GetDCA(ptrk,b,xn,xp);
3234      AliExternalTrackParam nt(*ntrk), pt(*ptrk);
3235      nt.PropagateTo(xn,b); pt.PropagateTo(xp,b);
3236
3237
3238      //--- Like in ITSV0Finder
3239      AliExternalTrackParam ntAt0(*ntrk), ptAt0(*ptrk);
3240      Double_t xxP,yyP,alphaP;
3241      Double_t rP[3];
3242
3243      //     if (!ptAt0.GetGlobalXYZat(ptAt0->GetX(),xxP,yyP,zzP)) continue;
3244      if (!ptAt0.GetXYZAt(ptAt0.GetX(),b,rP)) continue;
3245      xxP=rP[0];
3246      yyP=rP[1];
3247      alphaP = TMath::ATan2(yyP,xxP);
3248
3249
3250      ptAt0.Propagate(alphaP,0,b);
3251      Float_t ptfacP  = (1.+100.*TMath::Abs(ptAt0.GetC(b)));
3252
3253      //     Double_t distP      = ptAt0.GetY();
3254      Double_t normP      = ptfacP*TMath::Sqrt(ptAt0.GetSigmaY2());
3255      Double_t normdist0P = TMath::Abs(ptAt0.GetY()/normP);
3256      Double_t normdist1P = TMath::Abs((ptAt0.GetZ()-zPrimaryVertex)/(ptfacP*TMath::Sqrt(ptAt0.GetSigmaZ2())));
3257      Double_t normdistP  = TMath::Sqrt(normdist0P*normdist0P+normdist1P*normdist1P);
3258   
3259
3260      Double_t xxN,yyN,alphaN;
3261      Double_t rN[3];
3262      //     if (!ntAt0.GetGlobalXYZat(ntAt0->GetX(),xxN,yyN,zzN)) continue;
3263      if (!ntAt0.GetXYZAt(ntAt0.GetX(),b,rN)) continue;
3264      xxN=rN[0];
3265      yyN=rN[1];
3266  
3267      alphaN = TMath::ATan2(yyN,xxN);
3268
3269      ntAt0.Propagate(alphaN,0,b);
3270
3271      Float_t ptfacN  = (1.+100.*TMath::Abs(ntAt0.GetC(b)));
3272      //     Double_t distN      = ntAt0.GetY();
3273      Double_t normN      = ptfacN*TMath::Sqrt(ntAt0.GetSigmaY2());
3274      Double_t normdist0N = TMath::Abs(ntAt0.GetY()/normN);
3275      Double_t normdist1N = TMath::Abs((ntAt0.GetZ()-zPrimaryVertex)/(ptfacN*TMath::Sqrt(ntAt0.GetSigmaZ2())));
3276      Double_t normdistN  = TMath::Sqrt(normdist0N*normdist0N+normdist1N*normdist1N);
3277   
3278      //-----------------------------
3279
3280      Double_t momNegAt[3];
3281      nt.GetPxPyPz(momNegAt);
3282      curElecNegAt.SetXYZM(momNegAt[0],momNegAt[1],momNegAt[2],massE);
3283
3284      Double_t momPosAt[3];
3285      pt.GetPxPyPz(momPosAt);
3286      curElecPosAt.SetXYZM(momPosAt[0],momPosAt[1],momPosAt[2],massE);
3287      if(dca>1){
3288        continue;
3289      }
3290
3291      //     Double_t dneg= negTrack->GetD(xPrimaryVertex,yPrimaryVertex,b);
3292      //     Double_t dpos= posTrack->GetD(xPrimaryVertex,yPrimaryVertex,b);
3293      AliESDv0 vertex(nt,jTracks,pt,iTracks);
3294     
3295
3296      Float_t cpa=vertex.GetV0CosineOfPointingAngle(xPrimaryVertex,yPrimaryVertex,zPrimaryVertex);
3297
3298  
3299
3300      //  cout<< "v0Rr::"<< v0Rr<<endl;
3301      // if (pvertex.GetRr()<0.5){
3302      // continue;
3303      //}
3304      //     cout<<"vertex.GetChi2V0()"<<vertex.GetChi2V0()<<endl;
3305      if(cpa<0.9)continue;
3306      //     if (vertex.GetChi2V0() > 30) continue;
3307      //     cout<<"xp+xn::"<<xp<<" "<<xn<<endl;
3308      if ((xn+xp) < 0.4) continue;
3309      if (TMath::Abs(ntrk->GetD(xPrimaryVertex,yPrimaryVertex,b))<0.05)
3310        if (TMath::Abs(ptrk->GetD(xPrimaryVertex,yPrimaryVertex,b))<0.05) continue;
3311
3312      //cout<<"pass"<<endl;
3313
3314      AliKFParticle v0GammaC;
3315      v0GammaC+=(*negKF);
3316      v0GammaC+=(*posKF);
3317      v0GammaC.SetMassConstraint(0,0.001);
3318      primVtxImprovedGamma+=v0GammaC;
3319      v0GammaC.SetProductionVertex(primVtxImprovedGamma);
3320
3321
3322      curGamma=curElecNeg+curElecPos;
3323      curGammaAt=curElecNegAt+curElecPosAt;
3324      
3325      // invariant mass versus pt of K0short
3326      
3327      Double_t chi2V0GammaC=100000.;
3328      if( v0GammaC.GetNDF() != 0) {
3329        chi2V0GammaC = v0GammaC.GetChi2()/v0GammaC.GetNDF();
3330      }else{
3331        cout<< "ERROR::v0K0C.GetNDF()" << endl;
3332      }
3333
3334      if(chi2V0GammaC<200 &&chi2V0GammaC>0 ){
3335        if(fHistograms != NULL){
3336          fHistograms->FillHistogram("ESD_RecalculateV0_InvMass",v0GammaC.GetMass());
3337          fHistograms->FillHistogram("ESD_RecalculateV0_Pt",v0GammaC.GetPt());
3338          fHistograms->FillHistogram("ESD_RecalculateV0_E_dEdxP",curElecNegAt.P(),negTrack->GetTPCsignal());
3339          fHistograms->FillHistogram("ESD_RecalculateV0_P_dEdxP",curElecPosAt.P(),posTrack->GetTPCsignal());
3340          fHistograms->FillHistogram("ESD_RecalculateV0_cpa",cpa);
3341          fHistograms->FillHistogram("ESD_RecalculateV0_dca",dca);
3342          fHistograms->FillHistogram("ESD_RecalculateV0_normdistP",normdistP);
3343          fHistograms->FillHistogram("ESD_RecalculateV0_normdistN",normdistN);
3344
3345          new((*fKFRecalculatedGammasTClone)[fKFRecalculatedGammasTClone->GetEntriesFast()])  AliKFParticle(v0GammaC);
3346          fElectronRecalculatedv1.push_back(iTracks);
3347          fElectronRecalculatedv2.push_back(jTracks);
3348        }
3349      }
3350    }
3351  }
3352  
3353  for(Int_t firstGammaIndex=0;firstGammaIndex<fKFRecalculatedGammasTClone->GetEntriesFast();firstGammaIndex++){
3354    for(Int_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fKFRecalculatedGammasTClone->GetEntriesFast();secondGammaIndex++){
3355       AliKFParticle * twoGammaDecayCandidateDaughter0 = (AliKFParticle *)fKFRecalculatedGammasTClone->At(firstGammaIndex);
3356       AliKFParticle * twoGammaDecayCandidateDaughter1 = (AliKFParticle *)fKFRecalculatedGammasTClone->At(secondGammaIndex);
3357                         
3358       if(fElectronRecalculatedv1[firstGammaIndex]==fElectronRecalculatedv1[secondGammaIndex]){
3359         continue;
3360       }
3361       if( fElectronRecalculatedv2[firstGammaIndex]==fElectronRecalculatedv2[secondGammaIndex]){
3362         continue;
3363       }
3364                         
3365       AliKFParticle twoGammaCandidate(*twoGammaDecayCandidateDaughter0,*twoGammaDecayCandidateDaughter1);
3366       fHistograms->FillHistogram("ESD_RecalculateGG_InvMass",twoGammaCandidate.GetMass());              
3367       fHistograms->FillHistogram("ESD_RecalculateGG_InvMass_vs_Pt",twoGammaCandidate.GetMass(),twoGammaCandidate.GetPt());              
3368    }
3369  }
3370 }
3371 void AliAnalysisTaskGammaConversion::CalculateJetCone(Int_t gammaIndex){
3372   // CaculateJetCone
3373         
3374   Double_t cone;
3375   Double_t coneSize=0.3;
3376   Double_t ptJet=0;
3377         
3378   //  AliKFParticle * currentGamma = &fKFReconstructedGammas[gammaIndex];
3379   AliKFParticle * currentGamma = (AliKFParticle*)fKFReconstructedGammasTClone->At(gammaIndex);
3380
3381   TVector3 momentumVectorCurrentGamma(currentGamma->GetPx(),currentGamma->GetPy(),currentGamma->GetPz());
3382         
3383   AliESDtrack* leadingCharged = (AliESDtrack*)(fChargedParticles->At(fLeadingChargedIndex));
3384
3385   Double_t momLeadingCharged[3];
3386   leadingCharged->GetConstrainedPxPyPz(momLeadingCharged);
3387         
3388   TVector3 momentumVectorLeadingCharged(momLeadingCharged[0],momLeadingCharged[1],momLeadingCharged[2]);
3389         
3390   Double_t phi1=momentumVectorLeadingCharged.Phi();
3391   Double_t eta1=momentumVectorLeadingCharged.Eta();
3392   Double_t phi3=momentumVectorCurrentGamma.Phi();
3393         
3394   for(Int_t iCh=0;iCh<fChargedParticles->GetEntriesFast();iCh++){
3395     AliESDtrack* curTrack = (AliESDtrack*)(fChargedParticles->At(iCh));
3396     Int_t chId = fChargedParticlesId[iCh];
3397     if(fLeadingChargedIndex==chId || fLeadingChargedIndex==chId) continue;
3398     Double_t mom[3];
3399     curTrack->GetConstrainedPxPyPz(mom);
3400     TVector3 momentumVectorChargedParticle(mom[0],mom[1],mom[2]);
3401     Double_t phi2=momentumVectorChargedParticle.Phi();
3402     Double_t eta2=momentumVectorChargedParticle.Eta();
3403                 
3404                 
3405     cone=100.;
3406     if( TMath::Abs(phi2 - phi1) <= ( TMath::TwoPi()-coneSize) ){
3407       cone = TMath::Sqrt(  TMath::Power((eta2-eta1),2)+ TMath::Power((phi2-phi1),2) );
3408     }else{
3409       if( (phi2 - phi1)> TMath::TwoPi()-coneSize ){
3410         cone = TMath::Sqrt(  TMath::Power((eta2-eta1),2)+ TMath::Power((phi2-TMath::TwoPi()-phi1),2) );
3411       }
3412       if( (phi2 - phi1)< -(TMath::TwoPi()-coneSize) ){
3413         cone = TMath::Sqrt(  TMath::Power((eta2-eta1),2)+ TMath::Power((phi2+TMath::TwoPi()-phi1),2) );
3414       }
3415     }
3416                 
3417     if(cone <coneSize&& momentumVectorChargedParticle.Pt()>fMinPtJetCone ){
3418       ptJet+= momentumVectorChargedParticle.Pt();
3419       Double_t ffzHdrGam = momentumVectorChargedParticle.Pt()/momentumVectorCurrentGamma.Pt();
3420       Double_t imbalanceHdrGam=-momentumVectorChargedParticle.Dot(momentumVectorCurrentGamma)/momentumVectorCurrentGamma.Mag2();
3421       fHistograms->FillHistogram("ESD_FFzHdrGam",ffzHdrGam);
3422       fHistograms->FillHistogram("ESD_ImbalanceHdrGam",imbalanceHdrGam);
3423                         
3424     }
3425                 
3426     Double_t dphiHdrGam=phi3-phi2;
3427     if ( dphiHdrGam < (-TMath::PiOver2())){
3428       dphiHdrGam+=(TMath::TwoPi());
3429     }
3430                 
3431     if ( dphiHdrGam > (3.*TMath::PiOver2()) ){
3432       dphiHdrGam-=(TMath::TwoPi());
3433     }
3434                 
3435     if (momentumVectorChargedParticle.Pt()>fMinPtGamChargedCorr){
3436       fHistograms->FillHistogram("ESD_dphiHdrGamIsolated",dphiHdrGam);
3437     }
3438   }//track loop
3439         
3440         
3441 }
3442
3443 Double_t AliAnalysisTaskGammaConversion::GetMinimumDistanceToCharge(Int_t indexHighestPtGamma){
3444   // GetMinimumDistanceToCharge
3445         
3446   Double_t fIsoMin=100.;
3447   Double_t ptLeadingCharged=-1.;
3448
3449   fLeadingChargedIndex=-1;
3450         
3451   AliKFParticle * gammaHighest