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