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