Added new background scheeme, did some cleanup. added new bethe block parameters...
[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 "AliGammaConversionAODObject.h"
34 #include "AliGammaConversionBGHandler.h"
35
36 class AliCFContainer;
37 class AliCFManager;
38 class AliKFVertex;
39 class AliAODHandler;
40 class AliAODEvent;
41 class ALiESDEvent;
42 class AliMCEvent;
43 class AliMCEventHandler;
44 class AliESDInputHandler;
45 class AliAnalysisManager;
46 class Riostream;
47 class TFile;
48 class TInterpreter;
49 class TSystem;
50 class TROOT;
51
52 ClassImp(AliAnalysisTaskGammaConversion)
53
54
55 AliAnalysisTaskGammaConversion::AliAnalysisTaskGammaConversion():
56 AliAnalysisTaskSE(),
57   fV0Reader(NULL),
58   fStack(NULL),
59   fMCTruth(NULL),    // for CF
60   fGCMCEvent(NULL),    // for CF
61   fESDEvent(NULL),      
62   fOutputContainer(NULL),
63   fCFManager(0x0),   // for CF
64   fHistograms(NULL),
65   fTriggerCINT1B(kFALSE),
66   fDoMCTruth(kFALSE),
67   fDoNeutralMeson(kFALSE),
68   fDoJet(kFALSE),
69   fDoChic(kFALSE),
70   fKFReconstructedGammasTClone(NULL),
71   fCurrentEventPosElectronTClone(NULL),
72   fCurrentEventNegElectronTClone(NULL),
73   fKFReconstructedGammasCutTClone(NULL),
74   fPreviousEventTLVNegElectronTClone(NULL),
75   fPreviousEventTLVPosElectronTClone(NULL),     
76 //  fKFReconstructedGammas(),
77   fElectronv1(),
78   fElectronv2(),
79 //  fCurrentEventPosElectron(),
80 //  fCurrentEventNegElectron(),
81 //  fKFReconstructedGammasCut(),                        
82 //  fPreviousEventTLVNegElectron(),
83 //  fPreviousEventTLVPosElectron(),     
84   fElectronMass(-1),
85   fGammaMass(-1),
86   fPi0Mass(-1),
87   fEtaMass(-1),
88   fGammaWidth(-1),
89   fPi0Width(-1),
90   fEtaWidth(-1),
91   fMinOpeningAngleGhostCut(0.),
92   fEsdTrackCuts(NULL),
93   fCalculateBackground(kFALSE),
94   fWriteNtuple(kFALSE),
95   fGammaNtuple(NULL),
96   fNeutralMesonNtuple(NULL),
97   fTotalNumberOfAddedNtupleEntries(0),
98   fChargedParticles(NULL),
99   fChargedParticlesId(),
100   fGammaPtHighest(0.),
101   fMinPtForGammaJet(1.),
102   fMinIsoConeSize(0.2),
103   fMinPtIsoCone(0.7),
104   fMinPtGamChargedCorr(0.5),
105   fMinPtJetCone(0.5),
106   fLeadingChargedIndex(-1),
107   fLowPtMapping(1.),
108   fHighPtMapping(3.),
109   fDoCF(kFALSE),
110   fAODBranch(NULL),
111   fAODBranchName("GammaConv")//,
112   //  fAODObjects(NULL)
113 {
114   // Default constructor
115
116   /*   Kenneth: the default constructor should not have any define input/output or the call to SetESDtrackCuts
117   // Common I/O in slot 0
118   DefineInput (0, TChain::Class());
119   DefineOutput(0, TTree::Class());
120         
121   // Your private output
122   DefineOutput(1, TList::Class());
123         
124   // Define standard ESD track cuts for Gamma-hadron correlation 
125   SetESDtrackCuts();
126   */
127 }
128
129 AliAnalysisTaskGammaConversion::AliAnalysisTaskGammaConversion(const char* name):
130   AliAnalysisTaskSE(name),
131   fV0Reader(NULL),
132   fStack(NULL),
133   fMCTruth(NULL),    // for CF
134   fGCMCEvent(NULL),    // for CF
135   fESDEvent(NULL),      
136   fOutputContainer(0x0),
137   fCFManager(0x0),   // for CF
138   fHistograms(NULL),
139   fTriggerCINT1B(kFALSE),
140   fDoMCTruth(kFALSE),
141   fDoNeutralMeson(kFALSE),
142   fDoJet(kFALSE),
143   fDoChic(kFALSE),
144   fKFReconstructedGammasTClone(NULL),
145   fCurrentEventPosElectronTClone(NULL),
146   fCurrentEventNegElectronTClone(NULL),
147   fKFReconstructedGammasCutTClone(NULL),
148   fPreviousEventTLVNegElectronTClone(NULL),
149   fPreviousEventTLVPosElectronTClone(NULL),     
150   //  fKFReconstructedGammas(),
151   fElectronv1(),
152   fElectronv2(),
153   //  fCurrentEventPosElectron(),
154   //  fCurrentEventNegElectron(),
155   //  fKFReconstructedGammasCut(),      
156   //  fPreviousEventTLVNegElectron(),
157   //  fPreviousEventTLVPosElectron(),
158   fElectronMass(-1),
159   fGammaMass(-1),
160   fPi0Mass(-1),
161   fEtaMass(-1),
162   fGammaWidth(-1),
163   fPi0Width(-1),
164   fEtaWidth(-1),
165   fMinOpeningAngleGhostCut(0.),
166   fEsdTrackCuts(NULL),
167   fCalculateBackground(kFALSE),
168   fWriteNtuple(kFALSE),
169   fGammaNtuple(NULL),
170   fNeutralMesonNtuple(NULL),
171   fTotalNumberOfAddedNtupleEntries(0),
172   fChargedParticles(NULL),
173   fChargedParticlesId(),
174   fGammaPtHighest(0.),
175   fMinPtForGammaJet(1.),
176   fMinIsoConeSize(0.2),
177   fMinPtIsoCone(0.7),
178   fMinPtGamChargedCorr(0.5),
179   fMinPtJetCone(0.5),
180   fLeadingChargedIndex(-1),
181   fLowPtMapping(1.),
182   fHighPtMapping(3.),
183   fDoCF(kFALSE),
184   fAODBranch(NULL),
185   fAODBranchName("GammaConv")//,
186   // fAODObjects(NULL)
187 {
188   // Common I/O in slot 0
189   DefineInput (0, TChain::Class());
190   DefineOutput(0, TTree::Class());
191         
192   // Your private output
193   DefineOutput(1, TList::Class());
194   DefineOutput(2, AliCFContainer::Class());  // for CF
195         
196         
197   // Define standard ESD track cuts for Gamma-hadron correlation 
198   SetESDtrackCuts();
199 }
200
201 AliAnalysisTaskGammaConversion::~AliAnalysisTaskGammaConversion() 
202 {
203   // Remove all pointers
204         
205   if(fOutputContainer){
206     fOutputContainer->Clear() ; 
207     delete fOutputContainer ;
208   }
209   if(fHistograms){
210     delete fHistograms;
211   }
212   if(fV0Reader){
213     delete fV0Reader;
214   }
215         
216   // for CF
217   if(fCFManager){
218     delete fCFManager;
219   }
220
221   if(fEsdTrackCuts){
222     delete fEsdTrackCuts;
223   }
224
225   if (fAODBranch) {
226     fAODBranch->Clear();
227     delete fAODBranch ;
228   }
229 }
230
231
232 void AliAnalysisTaskGammaConversion::Init()
233 {
234   // Initialization
235   // AliLog::SetGlobalLogLevel(AliLog::kError);
236 }
237 void AliAnalysisTaskGammaConversion::SetESDtrackCuts()
238 {
239   // SetESDtrackCuts
240   if (fEsdTrackCuts!=NULL){
241     delete fEsdTrackCuts;
242   }
243   fEsdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts");
244   //standard cuts from:
245   //http://aliceinfo.cern.ch/alicvs/viewvc/PWG0/dNdEta/CreateCuts.C?revision=1.4&view=markup
246   //fEsdTrackCuts->SetMinNClustersTPC(50);
247   //fEsdTrackCuts->SetMaxChi2PerClusterTPC(3.5);
248   //fEsdTrackCuts->SetMaxCovDiagonalElements(2,2,0.5,0.5,2);
249   fEsdTrackCuts->SetRequireTPCRefit(kTRUE);
250   fEsdTrackCuts->SetRequireITSRefit(kTRUE);
251   fEsdTrackCuts->SetMaxNsigmaToVertex(3);
252   fEsdTrackCuts->SetRequireSigmaToVertex(kTRUE);
253   //  fEsdTrackCuts->SetAcceptKinkDaughters(kFALSE);
254         
255 }
256
257 void AliAnalysisTaskGammaConversion::Exec(Option_t */*option*/)
258 {
259   // Execute analysis for current event
260   ConnectInputData("");
261         
262   //Each event needs an empty branch
263   //  fAODBranch->Clear();
264   fAODBranch->Delete();
265         
266   if(fKFReconstructedGammasTClone == NULL){
267     fKFReconstructedGammasTClone = new TClonesArray("AliKFParticle",0);
268   }
269   if(fCurrentEventPosElectronTClone == NULL){
270     fCurrentEventPosElectronTClone = new TClonesArray("AliESDtrack",0);
271   }
272   if(fCurrentEventNegElectronTClone == NULL){
273     fCurrentEventNegElectronTClone = new TClonesArray("AliESDtrack",0);
274   }
275   if(fKFReconstructedGammasCutTClone == NULL){
276     fKFReconstructedGammasCutTClone = new TClonesArray("AliKFParticle",0);
277   }
278   if(fPreviousEventTLVNegElectronTClone == NULL){
279     fPreviousEventTLVNegElectronTClone = new TClonesArray("TLorentzVector",0);
280   }
281   if(fPreviousEventTLVPosElectronTClone == NULL){
282     fPreviousEventTLVPosElectronTClone  = new TClonesArray("TLorentzVector",0);
283   }
284   if(fChargedParticles == NULL){
285     fChargedParticles = new TClonesArray("AliESDtrack",0);
286   }
287         
288   //clear TClones
289   fKFReconstructedGammasTClone->Delete();
290   fCurrentEventPosElectronTClone->Delete();
291   fCurrentEventNegElectronTClone->Delete();
292   fKFReconstructedGammasCutTClone->Delete();
293   fPreviousEventTLVNegElectronTClone->Delete();
294   fPreviousEventTLVPosElectronTClone->Delete();
295
296   //clear vectors
297   fElectronv1.clear();
298   fElectronv2.clear();
299         
300   fChargedParticles->Delete();  
301
302   fChargedParticlesId.clear();  
303         
304   //Clear the data in the v0Reader
305   //  fV0Reader->UpdateEventByEventData();
306
307   //Take Only events with proper trigger
308   /*
309   if(fTriggerCINT1B){
310     if(!fV0Reader->GetESDEvent()->IsTriggerClassFired("CINT1B-ABCE-NOPF-ALL")) return;
311   }
312   */
313         
314   // Process the MC information
315   if(fDoMCTruth){
316     ProcessMCData();
317   }
318         
319   //Process the v0 information with no cuts
320   ProcessV0sNoCut();
321
322   // Process the v0 information
323   ProcessV0s();
324   
325
326   //Fill Gamma AOD
327   FillAODWithConversionGammas() ; 
328         
329   // Process reconstructed gammas
330   if(fDoNeutralMeson == kTRUE){
331     ProcessGammasForNeutralMesonAnalysis();
332   }
333
334   if(fDoMCTruth == kTRUE){
335     CheckV0Efficiency();
336   }
337   //Process reconstructed gammas electrons for Chi_c Analysis
338   if(fDoChic == kTRUE){
339     ProcessGammaElectronsForChicAnalysis();
340   }
341   // Process reconstructed gammas for gamma Jet/hadron correlations
342   if(fDoJet == kTRUE){
343     ProcessGammasForGammaJetAnalysis();
344   }
345         
346   //calculate background if flag is set
347   if(fCalculateBackground){
348     CalculateBackground();
349   }
350         
351   //Clear the data in the v0Reader
352   fV0Reader->UpdateEventByEventData();
353
354   PostData(1, fOutputContainer);
355   PostData(2, fCFManager->GetParticleContainer());  // for CF
356         
357 }
358
359 void AliAnalysisTaskGammaConversion::ConnectInputData(Option_t *option){
360   // see header file for documentation
361
362   AliAnalysisTaskSE::ConnectInputData(option);
363
364   if(fV0Reader == NULL){
365     // Write warning here cuts and so on are default if this ever happens
366   }
367   fV0Reader->Initialize();
368   fDoMCTruth = fV0Reader->GetDoMCTruth();
369 }
370
371
372
373 void AliAnalysisTaskGammaConversion::ProcessMCData(){
374   // see header file for documentation
375         
376   fStack = fV0Reader->GetMCStack();
377   fMCTruth = fV0Reader->GetMCTruth();  // for CF
378   fGCMCEvent = fV0Reader->GetMCEvent();  // for CF
379         
380         
381   // for CF
382   Double_t containerInput[3];
383   if(fDoCF){
384     if(!fGCMCEvent) cout << "NO MC INFO FOUND" << endl;
385     fCFManager->SetEventInfo(fGCMCEvent);
386   } 
387   // end for CF
388         
389         
390   if(fV0Reader->CheckForPrimaryVertex() == kFALSE){
391     return; // aborts if the primary vertex does not have contributors.
392   }
393         
394   for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++) {
395     TParticle* particle = (TParticle *)fStack->Particle(iTracks);
396                 
397     if (!particle) {
398       //print warning here
399       continue;
400     }
401                 
402     ///////////////////////Begin Chic Analysis/////////////////////////////
403                 
404     if(particle->GetPdgCode() == 443){//Is JPsi 
405       if(particle->GetNDaughters()==2){
406         if(TMath::Abs(fStack->Particle(particle->GetFirstDaughter())->GetPdgCode()) == 11 &&
407            TMath::Abs(fStack->Particle(particle->GetLastDaughter())->GetPdgCode()) == 11){
408           TParticle* daug0 = fStack->Particle(particle->GetFirstDaughter());
409           TParticle* daug1 = fStack->Particle(particle->GetLastDaughter());
410           if(TMath::Abs(daug0->Eta()) < 0.9 && TMath::Abs(daug1->Eta()) < 0.9)
411             fHistograms->FillTable("Table_Electrons",3);//e+ e-  from J/Psi inside acceptance
412                                         
413           if( TMath::Abs(daug0->Eta()) < 0.9){
414             if(daug0->GetPdgCode() == -11)
415               fHistograms->FillTable("Table_Electrons",1);//e+  from J/Psi inside acceptance
416             else
417               fHistograms->FillTable("Table_Electrons",2);//e-   from J/Psi inside acceptance
418                                                 
419           }
420           if(TMath::Abs(daug1->Eta()) < 0.9){
421             if(daug1->GetPdgCode() == -11)
422               fHistograms->FillTable("Table_Electrons",1);//e+  from J/Psi inside acceptance
423             else
424               fHistograms->FillTable("Table_Electrons",2);//e-   from J/Psi inside acceptance
425           }
426         }
427       }
428     }
429     //              const int CHI_C0   = 10441;
430     //              const int CHI_C1   = 20443;
431     //              const int CHI_C2   = 445
432     if(particle->GetPdgCode() == 22){//gamma from JPsi
433       if(particle->GetMother(0) > -1){
434         if(fStack->Particle(particle->GetMother(0))->GetPdgCode() == 10441 ||
435            fStack->Particle(particle->GetMother(0))->GetPdgCode() == 20443 ||
436            fStack->Particle(particle->GetMother(0))->GetPdgCode() == 445){
437           if(TMath::Abs(particle->Eta()) < 1.2)
438             fHistograms->FillTable("Table_Electrons",17);// gamma from chic inside accptance
439         }
440       }
441     }
442     if(particle->GetPdgCode() == 10441 || particle->GetPdgCode() == 20443 || particle->GetPdgCode() == 445){
443       if( particle->GetNDaughters() == 2){
444         TParticle* daug0 = fStack->Particle(particle->GetFirstDaughter());
445         TParticle* daug1 = fStack->Particle(particle->GetLastDaughter());
446                                 
447         if( (daug0->GetPdgCode() == 443 || daug0->GetPdgCode() == 22) && (daug1->GetPdgCode() == 443 || daug1->GetPdgCode() == 22) ){
448           if( daug0->GetPdgCode() == 443){
449             TParticle* daugE0 = fStack->Particle(daug0->GetFirstDaughter());
450             TParticle* daugE1 = fStack->Particle(daug0->GetLastDaughter());
451             if( TMath::Abs(daug1->Eta()) < 1.2 && TMath::Abs(daugE0->Eta()) < 0.9 && TMath::Abs(daugE1->Eta()) < 0.9 )
452               fHistograms->FillTable("Table_Electrons",18);
453                                                 
454           }//if
455           else if (daug1->GetPdgCode() == 443){
456             TParticle* daugE0 = fStack->Particle(daug1->GetFirstDaughter());
457             TParticle* daugE1 = fStack->Particle(daug1->GetLastDaughter());
458             if( TMath::Abs(daug0->Eta()) < 1.2 && TMath::Abs(daugE0->Eta()) < 0.9 && TMath::Abs(daugE1->Eta()) < 0.9 )
459               fHistograms->FillTable("Table_Electrons",18);
460           }//else if
461         }//gamma o Jpsi
462       }//GetNDaughters
463     }
464                 
465                 
466     /////////////////////End Chic Analysis////////////////////////////
467                 
468                 
469     if(TMath::Abs(particle->Eta())> fV0Reader->GetEtaCut() )    continue;
470                 
471     if(particle->R()>fV0Reader->GetMaxRCut())   continue; // cuts on distance from collision point
472                 
473     Double_t tmpPhi=particle->Phi();
474                 
475     if(particle->Phi()> TMath::Pi()){
476       tmpPhi = particle->Phi()-(2*TMath::Pi());
477     }
478                 
479     Double_t rapidity;
480     if(particle->Energy() - particle->Pz() == 0 || particle->Energy() + particle->Pz() == 0){
481       rapidity=0;
482     }
483     else{
484       rapidity = 0.5*(TMath::Log((particle->Energy()+particle->Pz()) / (particle->Energy()-particle->Pz())));
485     }   
486                 
487     //process the gammas
488     if (particle->GetPdgCode() == 22){
489       
490
491       if(particle->GetMother(0) >-1 && fStack->Particle(particle->GetMother(0))->GetPdgCode() == 22){
492         continue; // no photon as mothers!
493       }
494       
495       if(particle->GetMother(0) >= fStack->GetNprimary()){
496         continue; // the gamma has a mother, and it is not a primary particle
497       }
498                         
499       fHistograms->FillHistogram("MC_allGamma_Energy", particle->Energy());
500       fHistograms->FillHistogram("MC_allGamma_Pt", particle->Pt());
501       fHistograms->FillHistogram("MC_allGamma_Eta", particle->Eta());
502       fHistograms->FillHistogram("MC_allGamma_Phi", tmpPhi);
503       fHistograms->FillHistogram("MC_allGamma_Rapid", rapidity);
504                         
505       // for CF
506       if(fDoCF){
507         containerInput[0] = particle->Pt();
508         containerInput[1] = particle->Eta();
509         if(particle->GetMother(0) >=0){
510           containerInput[2] = fStack->Particle(particle->GetMother(0))->GetMass();
511         }
512         else{
513           containerInput[2]=-1;
514         }
515       
516         fCFManager->GetParticleContainer()->Fill(containerInput,kStepGenerated);                                        // generated gamma
517       }         
518       if(particle->GetMother(0) < 0){   // direct gamma
519         fHistograms->FillHistogram("MC_allDirectGamma_Energy",particle->Energy());
520         fHistograms->FillHistogram("MC_allDirectGamma_Pt", particle->Pt());
521         fHistograms->FillHistogram("MC_allDirectGamma_Eta", particle->Eta());
522         fHistograms->FillHistogram("MC_allDirectGamma_Phi", tmpPhi);
523         fHistograms->FillHistogram("MC_allDirectGamma_Rapid", rapidity);                                
524       }
525                         
526       // looking for conversion (electron + positron from pairbuilding (= 5) )
527       TParticle* ePos = NULL;
528       TParticle* eNeg = NULL;
529                         
530       if(particle->GetNDaughters() >= 2){
531         for(Int_t daughterIndex=particle->GetFirstDaughter();daughterIndex<=particle->GetLastDaughter();daughterIndex++){
532           TParticle *tmpDaughter = fStack->Particle(daughterIndex);
533           if(tmpDaughter->GetUniqueID() == 5){
534             if(tmpDaughter->GetPdgCode() == 11){
535               eNeg = tmpDaughter;
536             }
537             else if(tmpDaughter->GetPdgCode() == -11){
538               ePos = tmpDaughter;
539             }
540           }
541         }
542       }
543                         
544                         
545       if(ePos == NULL || eNeg == NULL){ // means we do not have two daughters from pair production
546         continue;
547       }
548                         
549                         
550       Double_t ePosPhi = ePos->Phi();
551       if(ePos->Phi()> TMath::Pi()) ePosPhi = ePos->Phi()-(2*TMath::Pi());
552                         
553       Double_t eNegPhi = eNeg->Phi();
554       if(eNeg->Phi()> TMath::Pi()) eNegPhi = eNeg->Phi()-(2*TMath::Pi());
555                         
556                         
557       if(ePos->Pt()<fV0Reader->GetPtCut() || eNeg->Pt()<fV0Reader->GetPtCut()){
558         continue; // no reconstruction below the Pt cut
559       }
560                         
561       if(TMath::Abs(ePos->Eta())> fV0Reader->GetEtaCut() || TMath::Abs(eNeg->Eta())> fV0Reader->GetEtaCut()){
562         continue;
563       } 
564                         
565       if(ePos->R()>fV0Reader->GetMaxRCut()){
566         continue; // cuts on distance from collision point
567       }
568
569       if(TMath::Abs(ePos->Vz()) > fV0Reader->GetMaxZCut()){
570         continue;   // outside material
571       }
572                         
573                         
574       if((TMath::Abs(ePos->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()  > ePos->R()){
575         continue;               // line cut to exclude regions where we do not reconstruct
576       }         
577                 
578                         
579       // for CF
580       if(fDoCF){
581         fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructable);  // reconstructable gamma        
582       }
583       fHistograms->FillHistogram("MC_ConvGamma_Energy", particle->Energy());
584       fHistograms->FillHistogram("MC_ConvGamma_Pt", particle->Pt());
585       fHistograms->FillHistogram("MC_ConvGamma_Eta", particle->Eta());
586       fHistograms->FillHistogram("MC_ConvGamma_Phi", tmpPhi);
587       fHistograms->FillHistogram("MC_ConvGamma_Rapid", rapidity);
588       fHistograms->FillHistogram("MC_ConvGamma_Pt_Eta", particle->Pt(),particle->Eta());
589                         
590       fHistograms->FillHistogram("MC_E_Energy", eNeg->Energy());
591       fHistograms->FillHistogram("MC_E_Pt", eNeg->Pt());
592       fHistograms->FillHistogram("MC_E_Eta", eNeg->Eta());
593       fHistograms->FillHistogram("MC_E_Phi", eNegPhi);
594                         
595       fHistograms->FillHistogram("MC_P_Energy", ePos->Energy());
596       fHistograms->FillHistogram("MC_P_Pt", ePos->Pt());
597       fHistograms->FillHistogram("MC_P_Eta", ePos->Eta());
598       fHistograms->FillHistogram("MC_P_Phi", ePosPhi);
599                         
600                         
601       // begin Mapping 
602       Int_t rBin    = fHistograms->GetRBin(ePos->R());
603       Int_t zBin    = fHistograms->GetZBin(ePos->Vz());
604       Int_t phiBin  = fHistograms->GetPhiBin(particle->Phi());
605
606       TVector3 vtxPos(ePos->Vx(),ePos->Vy(),ePos->Vz());        
607       
608       TString nameMCMappingPhiR="";
609       nameMCMappingPhiR.Form("MC_Conversion_Mapping_Phi%02d_R%02d",phiBin,rBin);
610       fHistograms->FillHistogram(nameMCMappingPhiR, ePos->Vz(), particle->Eta());
611                         
612       TString nameMCMappingPhi="";
613       nameMCMappingPhi.Form("MC_Conversion_Mapping_Phi%02d",phiBin);
614       //      fHistograms->FillHistogram(nameMCMappingPhi, particle->Eta());
615       fHistograms->FillHistogram(nameMCMappingPhi, ePos->Vz(), particle->Eta());
616                         
617       TString nameMCMappingR="";
618       nameMCMappingR.Form("MC_Conversion_Mapping_R%02d",rBin);
619       //      fHistograms->FillHistogram(nameMCMappingR, particle->Eta());
620       fHistograms->FillHistogram(nameMCMappingR,ePos->Vz(), particle->Eta());
621                         
622       TString nameMCMappingPhiInR="";
623       nameMCMappingPhiInR.Form("MC_Conversion_Mapping_Phi_in_R_%02d",rBin);
624       //      fHistograms->FillHistogram(nameMCMappingPhiInR, tmpPhi);
625       fHistograms->FillHistogram(nameMCMappingPhiInR, vtxPos.Phi());
626
627       TString nameMCMappingZInR="";
628       nameMCMappingZInR.Form("MC_Conversion_Mapping_Z_in_R_%02d",rBin);
629       fHistograms->FillHistogram(nameMCMappingZInR,ePos->Vz() );
630
631
632       TString nameMCMappingPhiInZ="";
633       nameMCMappingPhiInZ.Form("MC_Conversion_Mapping_Phi_in_Z_%02d",zBin);
634       //      fHistograms->FillHistogram(nameMCMappingPhiInR, tmpPhi);
635       fHistograms->FillHistogram(nameMCMappingPhiInZ, vtxPos.Phi());
636
637       TString nameMCMappingRInZ="";
638       nameMCMappingRInZ.Form("MC_Conversion_Mapping_R_in_Z_%02d",zBin);
639       fHistograms->FillHistogram(nameMCMappingRInZ,ePos->R() );
640
641       if(particle->Pt() > fLowPtMapping && particle->Pt()< fHighPtMapping){
642         TString nameMCMappingMidPtPhiInR="";
643         nameMCMappingMidPtPhiInR.Form("MC_Conversion_Mapping_MidPt_Phi_in_R_%02d",rBin);
644         fHistograms->FillHistogram(nameMCMappingMidPtPhiInR, vtxPos.Phi());
645         
646         TString nameMCMappingMidPtZInR="";
647         nameMCMappingMidPtZInR.Form("MC_Conversion_Mapping_MidPt_Z_in_R_%02d",rBin);
648         fHistograms->FillHistogram(nameMCMappingMidPtZInR,ePos->Vz() );
649         
650         
651         TString nameMCMappingMidPtPhiInZ="";
652         nameMCMappingMidPtPhiInZ.Form("MC_Conversion_Mapping_MidPt_Phi_in_Z_%02d",zBin);
653         fHistograms->FillHistogram(nameMCMappingMidPtPhiInZ, vtxPos.Phi());
654         
655         TString nameMCMappingMidPtRInZ="";
656         nameMCMappingMidPtRInZ.Form("MC_Conversion_Mapping_MidPt_R_in_Z_%02d",zBin);
657         fHistograms->FillHistogram(nameMCMappingMidPtRInZ,ePos->R() );
658
659       }
660
661       //end mapping
662                         
663       fHistograms->FillHistogram("MC_Conversion_R",ePos->R());
664       fHistograms->FillHistogram("MC_Conversion_ZR",ePos->Vz(),ePos->R());
665       fHistograms->FillHistogram("MC_Conversion_XY",ePos->Vx(),ePos->Vy());
666       fHistograms->FillHistogram("MC_Conversion_OpeningAngle",GetMCOpeningAngle(ePos, eNeg));
667       fHistograms->FillHistogram("MC_ConvGamma_E_AsymmetryP",particle->P(),eNeg->P()/particle->P());
668       fHistograms->FillHistogram("MC_ConvGamma_P_AsymmetryP",particle->P(),ePos->P()/particle->P());
669
670
671       if(particle->GetMother(0) < 0){ // no mother = direct gamma, still inside converted
672         fHistograms->FillHistogram("MC_ConvDirectGamma_Energy",particle->Energy());
673         fHistograms->FillHistogram("MC_ConvDirectGamma_Pt", particle->Pt());
674         fHistograms->FillHistogram("MC_ConvDirectGamma_Eta", particle->Eta());
675         fHistograms->FillHistogram("MC_ConvDirectGamma_Phi", tmpPhi);
676         fHistograms->FillHistogram("MC_ConvDirectGamma_Rapid", rapidity);
677                                 
678       } // end direct gamma
679       else{   // mother exits 
680         /*      if( fStack->Particle(particle->GetMother(0))->GetPdgCode()==10441 ||//chic0 
681                 fStack->Particle(particle->GetMother(0))->GetPdgCode()==20443 ||//psi2S
682                 fStack->Particle(particle->GetMother(0))->GetPdgCode()==445  //chic2
683                 ){ 
684                 fMCGammaChic.push_back(particle);
685                 }
686         */
687       }  // end if mother exits
688     } // end if particle is a photon
689                 
690                 
691                 
692     // process motherparticles (2 gammas as daughters)
693     // the motherparticle had already to pass the R and the eta cut, but no line cut.
694     // the line cut is just valid for the conversions!
695                 
696     if(particle->GetNDaughters() == 2){
697                         
698       TParticle* daughter0 = (TParticle*)fStack->Particle(particle->GetFirstDaughter());
699       TParticle* daughter1 = (TParticle*)fStack->Particle(particle->GetLastDaughter());
700                         
701       if(daughter0->GetPdgCode() != 22 || daughter1->GetPdgCode() != 22) continue; //check for gamma gamma daughters
702                         
703       // Check the acceptance for both gammas
704       Bool_t gammaEtaCut = kTRUE;
705       if(TMath::Abs(daughter0->Eta()) > fV0Reader->GetEtaCut() || TMath::Abs(daughter1->Eta()) > fV0Reader->GetEtaCut()  ) gammaEtaCut = kFALSE;
706       Bool_t gammaRCut = kTRUE;
707       if(daughter0->R() > fV0Reader->GetMaxRCut() || daughter1->R() > fV0Reader->GetMaxRCut()  ) gammaRCut = kFALSE;
708                         
709                         
710                         
711       // check for conversions now -> have to pass eta, R and line cut!
712       Bool_t daughter0Electron = kFALSE;
713       Bool_t daughter0Positron = kFALSE;
714       Bool_t daughter1Electron = kFALSE;
715       Bool_t daughter1Positron = kFALSE;
716                         
717       if(daughter0->GetNDaughters() >= 2){  // first gamma
718         for(Int_t TrackIndex=daughter0->GetFirstDaughter();TrackIndex<=daughter0->GetLastDaughter();TrackIndex++){
719           TParticle *tmpDaughter = fStack->Particle(TrackIndex);
720           if(tmpDaughter->GetUniqueID() == 5){
721             if(tmpDaughter->GetPdgCode() == 11){
722               if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){
723                 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()  < tmpDaughter->R() ){
724                   if(tmpDaughter->R()< fV0Reader->GetMaxRCut()){
725                     daughter0Electron = kTRUE;
726                   }
727                 }
728                                                                 
729               }
730             }
731             else if(tmpDaughter->GetPdgCode() == -11){
732               if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){
733                 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()  < tmpDaughter->R() ){
734                   if(tmpDaughter->R()< fV0Reader->GetMaxRCut()){
735                     daughter0Positron = kTRUE;
736                   }
737                 }                                               
738               }                                 
739             }
740           }
741         }
742       }
743                         
744                         
745       if(daughter1->GetNDaughters() >= 2){   // second gamma
746         for(Int_t TrackIndex=daughter1->GetFirstDaughter();TrackIndex<=daughter1->GetLastDaughter();TrackIndex++){
747           TParticle *tmpDaughter = fStack->Particle(TrackIndex);
748           if(tmpDaughter->GetUniqueID() == 5){
749             if(tmpDaughter->GetPdgCode() == 11){
750               if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){
751                 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()  < tmpDaughter->R() ){
752                   if(tmpDaughter->R()< fV0Reader->GetMaxRCut()){
753                     daughter1Electron = kTRUE;
754                   }
755                 }
756                                                                 
757               }
758             }
759             else if(tmpDaughter->GetPdgCode() == -11){
760               if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){
761                 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()  < tmpDaughter->R() ){
762                   if(tmpDaughter->R()< fV0Reader->GetMaxRCut()){
763                     daughter1Positron = kTRUE;
764                   }
765                 }
766                                                                 
767               }
768                                                         
769             }
770           }
771         }
772       }
773                         
774                         
775                         
776                         
777       if(particle->GetPdgCode()==111){     //Pi0
778         if( iTracks >= fStack->GetNprimary()){
779           fHistograms->FillHistogram("MC_Pi0_Secondaries_Eta", particle->Eta());
780           fHistograms->FillHistogram("MC_Pi0_Secondaries_Rapid", rapidity);
781           fHistograms->FillHistogram("MC_Pi0_Secondaries_Phi", tmpPhi);
782           fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt", particle->Pt());
783           fHistograms->FillHistogram("MC_Pi0_Secondaries_Energy", particle->Energy());
784           fHistograms->FillHistogram("MC_Pi0_Secondaries_R", particle->R());
785           fHistograms->FillHistogram("MC_Pi0_Secondaries_ZR", particle->Vz(),particle->R());
786           fHistograms->FillHistogram("MC_Pi0_Secondaries_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));
787           fHistograms->FillHistogram("MC_Pi0_Secondaries_XY", particle->Vx(),particle->Vy());//only fill from one daughter to avoid multiple filling
788                                         
789           if(gammaEtaCut && gammaRCut){  
790             //if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){
791             fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());
792             fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);
793             if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){
794               fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());
795               fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);
796             }
797           }
798         }
799         else{
800           fHistograms->FillHistogram("MC_Pi0_Eta", particle->Eta());    
801           fHistograms->FillHistogram("MC_Pi0_Rapid", rapidity);
802           fHistograms->FillHistogram("MC_Pi0_Phi", tmpPhi);
803           fHistograms->FillHistogram("MC_Pi0_Pt", particle->Pt());
804           fHistograms->FillHistogram("MC_Pi0_Energy", particle->Energy());
805           fHistograms->FillHistogram("MC_Pi0_R", particle->R());
806           fHistograms->FillHistogram("MC_Pi0_ZR", particle->Vz(),particle->R());
807           fHistograms->FillHistogram("MC_Pi0_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));
808           fHistograms->FillHistogram("MC_Pi0_XY", particle->Vx(), particle->Vy());//only fill from one daughter to avoid multiple filling
809                                         
810           if(gammaEtaCut && gammaRCut){
811             //    if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){
812             fHistograms->FillHistogram("MC_Pi0_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());
813             fHistograms->FillHistogram("MC_Pi0_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);
814             if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){
815               fHistograms->FillHistogram("MC_Pi0_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());
816               fHistograms->FillHistogram("MC_Pi0_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);
817               fHistograms->FillHistogram("MC_Pi0_ZR_ConvGamma_withinAcceptance", particle->Vz(),particle->R());
818             }
819           }
820         }
821       }
822                         
823       if(particle->GetPdgCode()==221){   //Eta
824         fHistograms->FillHistogram("MC_Eta_Eta", particle->Eta());
825         fHistograms->FillHistogram("MC_Eta_Rapid", rapidity);
826         fHistograms->FillHistogram("MC_Eta_Phi",tmpPhi);
827         fHistograms->FillHistogram("MC_Eta_Pt", particle->Pt());
828         fHistograms->FillHistogram("MC_Eta_Energy", particle->Energy());
829         fHistograms->FillHistogram("MC_Eta_R", particle->R());
830         fHistograms->FillHistogram("MC_Eta_ZR", particle->Vz(),particle->R());
831         fHistograms->FillHistogram("MC_Eta_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));
832         fHistograms->FillHistogram("MC_Eta_XY", particle->Vx(), particle->Vy());//only fill from one daughter to avoid multiple filling
833                                 
834         if(gammaEtaCut && gammaRCut){  
835           //    if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){
836           fHistograms->FillHistogram("MC_Eta_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());
837           fHistograms->FillHistogram("MC_Eta_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);
838           if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){
839             fHistograms->FillHistogram("MC_Eta_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());
840             fHistograms->FillHistogram("MC_Eta_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);
841             fHistograms->FillHistogram("MC_Eta_ZR_ConvGamma_withinAcceptance", particle->Vz(),particle->R());
842           }
843                                         
844         }
845                                 
846       }
847                         
848       // all motherparticles with 2 gammas as daughters
849       fHistograms->FillHistogram("MC_Mother_R", particle->R());
850       fHistograms->FillHistogram("MC_Mother_ZR", particle->Vz(),particle->R());
851       fHistograms->FillHistogram("MC_Mother_XY", particle->Vx(),particle->Vy());
852       fHistograms->FillHistogram("MC_Mother_Mass", particle->GetCalcMass());
853       fHistograms->FillHistogram("MC_Mother_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));
854       fHistograms->FillHistogram("MC_Mother_Energy", particle->Energy());
855       fHistograms->FillHistogram("MC_Mother_Pt", particle->Pt());
856       fHistograms->FillHistogram("MC_Mother_Eta", particle->Eta());
857       fHistograms->FillHistogram("MC_Mother_Rapid", rapidity);
858       fHistograms->FillHistogram("MC_Mother_Phi",tmpPhi);
859       fHistograms->FillHistogram("MC_Mother_InvMass_vs_Pt",particle->GetMass(),particle->Pt());                 
860                         
861       if(gammaEtaCut && gammaRCut){  
862         //      if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){
863         fHistograms->FillHistogram("MC_Mother_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());
864         fHistograms->FillHistogram("MC_Mother_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);
865         fHistograms->FillHistogram("MC_Mother_InvMass_vs_Pt_withinAcceptance",particle->GetMass(),particle->Pt());                      
866         if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){
867           fHistograms->FillHistogram("MC_Mother_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());
868           fHistograms->FillHistogram("MC_Mother_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);
869           fHistograms->FillHistogram("MC_Mother_InvMass_vs_Pt_ConvGamma_withinAcceptance",particle->GetMass(),particle->Pt());                  
870                                         
871         }
872                                 
873                                 
874       }  // end passed R and eta cut
875                         
876     } // end if(particle->GetNDaughters() == 2)
877                 
878   }// end for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++)
879
880 } // end ProcessMCData
881
882
883
884 void AliAnalysisTaskGammaConversion::FillNtuple(){
885   //Fills the ntuple with the different values
886         
887   if(fGammaNtuple == NULL){
888     return;
889   }
890   Int_t numberOfV0s = fV0Reader->GetNumberOfV0s();
891   for(Int_t i=0;i<numberOfV0s;i++){
892     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};
893     AliESDv0 * cV0 = fV0Reader->GetV0(i);
894     Double_t negPID=0;
895     Double_t posPID=0;
896     fV0Reader->GetPIDProbability(negPID,posPID);
897     values[0]=cV0->GetOnFlyStatus();
898     values[1]=fV0Reader->CheckForPrimaryVertex();
899     values[2]=negPID;
900     values[3]=posPID;
901     values[4]=fV0Reader->GetX();
902     values[5]=fV0Reader->GetY();
903     values[6]=fV0Reader->GetZ();
904     values[7]=fV0Reader->GetXYRadius();
905     values[8]=fV0Reader->GetMotherCandidateNDF();
906     values[9]=fV0Reader->GetMotherCandidateChi2();
907     values[10]=fV0Reader->GetMotherCandidateEnergy();
908     values[11]=fV0Reader->GetMotherCandidateEta();
909     values[12]=fV0Reader->GetMotherCandidatePt();
910     values[13]=fV0Reader->GetMotherCandidateMass();
911     values[14]=fV0Reader->GetMotherCandidateWidth();
912     //    values[15]=fV0Reader->GetMotherMCParticle()->Pt();   MOVED TO THE END, HAS TO BE CALLED AFTER HasSameMother NB: still has the same entry in the array
913     values[16]=fV0Reader->GetOpeningAngle();
914     values[17]=fV0Reader->GetNegativeTrackEnergy();
915     values[18]=fV0Reader->GetNegativeTrackPt();
916     values[19]=fV0Reader->GetNegativeTrackEta();
917     values[20]=fV0Reader->GetNegativeTrackPhi();
918     values[21]=fV0Reader->GetPositiveTrackEnergy();
919     values[22]=fV0Reader->GetPositiveTrackPt();
920     values[23]=fV0Reader->GetPositiveTrackEta();
921     values[24]=fV0Reader->GetPositiveTrackPhi();
922     values[25]=fV0Reader->HasSameMCMother();
923     if(values[25] != 0){
924       values[26]=fV0Reader->GetMotherMCParticlePDGCode();
925       values[15]=fV0Reader->GetMotherMCParticle()->Pt();
926     }
927     fTotalNumberOfAddedNtupleEntries++;
928     fGammaNtuple->Fill(values);
929   }
930   fV0Reader->ResetV0IndexNumber();
931         
932 }
933
934 void AliAnalysisTaskGammaConversion::ProcessV0sNoCut(){
935   // Process all the V0's without applying any cuts to it
936         
937   Int_t numberOfV0s = fV0Reader->GetNumberOfV0s();
938   for(Int_t i=0;i<numberOfV0s;i++){
939     /*AliESDv0 * cV0 = */fV0Reader->GetV0(i);
940                 
941     if(fV0Reader->CheckForPrimaryVertex() == kFALSE){
942       continue;
943     }
944
945     //    if( !fV0Reader->GetV0(i)->GetOnFlyStatus()){
946     if( !fV0Reader->CheckV0FinderStatus(i)){
947       continue;
948     }
949
950
951     if( !((fV0Reader->GetNegativeESDTrack())->GetStatus() & AliESDtrack::kTPCrefit) || 
952         !((fV0Reader->GetPositiveESDTrack())->GetStatus() & AliESDtrack::kTPCrefit) ){
953       continue;
954     }
955
956     if(fDoMCTruth){
957                         
958       if(fV0Reader->HasSameMCMother() == kFALSE){
959         continue;
960       }
961                         
962       TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
963       TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
964                         
965       if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){
966         continue;
967       }
968       if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){
969         continue;
970       }
971
972       if(negativeMC->GetUniqueID() != 5 || positiveMC->GetUniqueID() !=5){
973         continue;
974       }
975                         
976       if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
977                                 
978         fHistograms->FillHistogram("ESD_NoCutConvGamma_Pt", fV0Reader->GetMotherCandidatePt());
979         fHistograms->FillHistogram("ESD_NoCutConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());
980         fHistograms->FillHistogram("ESD_NoCutConvGamma_Eta", fV0Reader->GetMotherCandidateEta());                               
981         fHistograms->FillHistogram("ESD_NoCutConvGamma_Phi", fV0Reader->GetMotherCandidatePhi());
982         fHistograms->FillHistogram("ESD_NoCutConvGamma_Mass", fV0Reader->GetMotherCandidateMass());
983         fHistograms->FillHistogram("ESD_NoCutConvGamma_Width", fV0Reader->GetMotherCandidateWidth());
984         fHistograms->FillHistogram("ESD_NoCutConvGamma_Chi2", fV0Reader->GetMotherCandidateChi2());
985         fHistograms->FillHistogram("ESD_NoCutConvGamma_NDF", fV0Reader->GetMotherCandidateNDF());
986         fHistograms->FillHistogram("ESD_NoCutConvGamma_Rapid", fV0Reader->GetMotherCandidateRapidity());
987         fHistograms->FillHistogram("ESD_NoCutConvGamma_Pt_Eta", fV0Reader->GetMotherCandidatePt(),fV0Reader->GetMotherCandidateEta());
988                                 
989         fHistograms->FillHistogram("ESD_NoCutConvGamma_Pt_Chi2", fV0Reader->GetMotherCandidatePt(), fV0Reader->GetMotherCandidateChi2());
990         fHistograms->FillHistogram("ESD_NoCutConvGamma_Eta_Chi2", fV0Reader->GetMotherCandidateEta(), fV0Reader->GetMotherCandidateChi2());
991                                 
992         fHistograms->FillHistogram("ESD_NoCutConversion_XY", fV0Reader->GetX(),fV0Reader->GetY());
993         fHistograms->FillHistogram("ESD_NoCutConversion_R", fV0Reader->GetXYRadius());
994         fHistograms->FillHistogram("ESD_NoCutConversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());
995         fHistograms->FillHistogram("ESD_NoCutConversion_OpeningAngle", fV0Reader->GetOpeningAngle());
996         fHistograms->FillHistogram("ESD_NoCutConvGamma_CosPointingAngle", fV0Reader->GetCosPointingAngle());
997         fHistograms->FillHistogram("ESD_NoCutConvGamma_DcaDaughters", fV0Reader->GetDcaDaughters());
998         fHistograms->FillHistogram("ESD_NoCutConvGamma_NormDcaDistDaughters", fV0Reader->GetNormDcaDistDaughters());
999         fHistograms->FillHistogram("ESD_NoCutConvGamma_LikelihoodAP", fV0Reader->GetLikelihoodAP());
1000
1001         fHistograms->FillHistogram("ESD_NoCutConvGamma_E_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetNegativeTrackP()/fV0Reader->GetMotherCandidateP());
1002         fHistograms->FillHistogram("ESD_NoCutConvGamma_P_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetPositiveTrackP()/fV0Reader->GetMotherCandidateP());
1003         fHistograms->FillHistogram("ESD_NoCutConvGamma_E_dEdxP",fV0Reader->GetNegativeTrackP(),fV0Reader->GetNegativeTrackTPCdEdx());
1004         fHistograms->FillHistogram("ESD_NoCutConvGamma_P_dEdxP",fV0Reader->GetPositiveTrackP(),fV0Reader->GetPositiveTrackTPCdEdx());
1005
1006         //store MCTruth properties
1007         fHistograms->FillHistogram("ESD_NoCutConvGamma_MC_Pt_Eta", fV0Reader->GetMotherMCParticle()->Pt(),fV0Reader->GetMotherMCParticle()->Eta());
1008         fHistograms->FillHistogram("ESD_NoCutConversion_MC_ZR", negativeMC->Vz(),negativeMC->R());
1009         fHistograms->FillHistogram("ESD_NoCutConversion_MC_XY", negativeMC->Vx(),negativeMC->Vy());
1010       }
1011     }
1012   }
1013   fV0Reader->ResetV0IndexNumber();
1014 }
1015
1016 void AliAnalysisTaskGammaConversion::ProcessV0s(){
1017   // see header file for documentation
1018         
1019   if(fWriteNtuple == kTRUE){
1020     FillNtuple();
1021   }
1022         
1023   Int_t nSurvivingV0s=0;
1024   while(fV0Reader->NextV0()){
1025     nSurvivingV0s++;
1026                 
1027                 
1028     //-------------------------- filling v0 information -------------------------------------
1029     fHistograms->FillHistogram("ESD_Conversion_R", fV0Reader->GetXYRadius());
1030     fHistograms->FillHistogram("ESD_Conversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());
1031     fHistograms->FillHistogram("ESD_Conversion_XY", fV0Reader->GetX(),fV0Reader->GetY());
1032     fHistograms->FillHistogram("ESD_Conversion_OpeningAngle", fV0Reader->GetOpeningAngle());    
1033                 
1034     fHistograms->FillHistogram("ESD_E_Energy", fV0Reader->GetNegativeTrackEnergy());
1035     fHistograms->FillHistogram("ESD_E_Pt", fV0Reader->GetNegativeTrackPt());
1036     fHistograms->FillHistogram("ESD_E_Eta", fV0Reader->GetNegativeTrackEta());
1037     fHistograms->FillHistogram("ESD_E_Phi", fV0Reader->GetNegativeTrackPhi());
1038                 
1039     fHistograms->FillHistogram("ESD_P_Energy", fV0Reader->GetPositiveTrackEnergy());
1040     fHistograms->FillHistogram("ESD_P_Pt", fV0Reader->GetPositiveTrackPt());
1041     fHistograms->FillHistogram("ESD_P_Eta", fV0Reader->GetPositiveTrackEta());
1042     fHistograms->FillHistogram("ESD_P_Phi", fV0Reader->GetPositiveTrackPhi());
1043                 
1044     fHistograms->FillHistogram("ESD_ConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());
1045     fHistograms->FillHistogram("ESD_ConvGamma_Pt", fV0Reader->GetMotherCandidatePt());
1046     fHistograms->FillHistogram("ESD_ConvGamma_Eta", fV0Reader->GetMotherCandidateEta());
1047     fHistograms->FillHistogram("ESD_ConvGamma_Phi", fV0Reader->GetMotherCandidatePhi());
1048     fHistograms->FillHistogram("ESD_ConvGamma_Mass", fV0Reader->GetMotherCandidateMass());
1049     fHistograms->FillHistogram("ESD_ConvGamma_Width", fV0Reader->GetMotherCandidateWidth());
1050     fHistograms->FillHistogram("ESD_ConvGamma_Chi2", fV0Reader->GetMotherCandidateChi2());
1051     fHistograms->FillHistogram("ESD_ConvGamma_NDF", fV0Reader->GetMotherCandidateNDF());
1052     fHistograms->FillHistogram("ESD_ConvGamma_Rapid", fV0Reader->GetMotherCandidateRapidity());
1053     fHistograms->FillHistogram("ESD_ConvGamma_Pt_Eta", fV0Reader->GetMotherCandidatePt(),fV0Reader->GetMotherCandidateEta());
1054                 
1055     fHistograms->FillHistogram("ESD_ConvGamma_Pt_Chi2", fV0Reader->GetMotherCandidatePt(), fV0Reader->GetMotherCandidateChi2());
1056     fHistograms->FillHistogram("ESD_ConvGamma_Eta_Chi2", fV0Reader->GetMotherCandidateEta(), fV0Reader->GetMotherCandidateChi2());
1057     
1058     fHistograms->FillHistogram("ESD_ConvGamma_CosPointingAngle", fV0Reader->GetCosPointingAngle());
1059     fHistograms->FillHistogram("ESD_ConvGamma_DcaDaughters", fV0Reader->GetDcaDaughters());
1060     fHistograms->FillHistogram("ESD_ConvGamma_NormDcaDistDaughters", fV0Reader->GetNormDcaDistDaughters());
1061     fHistograms->FillHistogram("ESD_ConvGamma_LikelihoodAP", fV0Reader->GetLikelihoodAP());
1062     
1063     fHistograms->FillHistogram("ESD_ConvGamma_E_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetNegativeTrackP()/fV0Reader->GetMotherCandidateP());
1064     fHistograms->FillHistogram("ESD_ConvGamma_P_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetPositiveTrackP()/fV0Reader->GetMotherCandidateP());
1065     fHistograms->FillHistogram("ESD_ConvGamma_E_dEdxP",fV0Reader->GetNegativeTrackP(),fV0Reader->GetNegativeTrackTPCdEdx());
1066     fHistograms->FillHistogram("ESD_ConvGamma_P_dEdxP",fV0Reader->GetPositiveTrackP(),fV0Reader->GetPositiveTrackTPCdEdx());
1067     
1068                 
1069                 
1070     // begin mapping
1071     Int_t rBin    = fHistograms->GetRBin(fV0Reader->GetXYRadius());
1072     Int_t zBin    = fHistograms->GetZBin(fV0Reader->GetZ());
1073     Int_t phiBin  = fHistograms->GetPhiBin(fV0Reader->GetNegativeTrackPhi());
1074     TVector3 vtxConv(fV0Reader->GetX(),fV0Reader->GetY(), fV0Reader->GetZ());
1075
1076     Double_t motherCandidateEta= fV0Reader->GetMotherCandidateEta();
1077                 
1078     TString nameESDMappingPhiR="";
1079     nameESDMappingPhiR.Form("ESD_Conversion_Mapping_Phi%02d_R%02d",phiBin,rBin);
1080     fHistograms->FillHistogram(nameESDMappingPhiR, fV0Reader->GetZ(), motherCandidateEta);
1081                 
1082     TString nameESDMappingPhi="";
1083     nameESDMappingPhi.Form("ESD_Conversion_Mapping_Phi%02d",phiBin);
1084     fHistograms->FillHistogram(nameESDMappingPhi, fV0Reader->GetZ(), motherCandidateEta);
1085                 
1086     TString nameESDMappingR="";
1087     nameESDMappingR.Form("ESD_Conversion_Mapping_R%02d",rBin);
1088     fHistograms->FillHistogram(nameESDMappingR, fV0Reader->GetZ(), motherCandidateEta);  
1089                 
1090     TString nameESDMappingPhiInR="";
1091     nameESDMappingPhiInR.Form("ESD_Conversion_Mapping_Phi_in_R_%02d",rBin);
1092     //    fHistograms->FillHistogram(nameESDMappingPhiInR, fV0Reader->GetMotherCandidatePhi());
1093     fHistograms->FillHistogram(nameESDMappingPhiInR, vtxConv.Phi());
1094
1095     TString nameESDMappingZInR="";
1096     nameESDMappingZInR.Form("ESD_Conversion_Mapping_Z_in_R_%02d",rBin);
1097     fHistograms->FillHistogram(nameESDMappingZInR, fV0Reader->GetZ());
1098
1099     TString nameESDMappingPhiInZ="";
1100     nameESDMappingPhiInZ.Form("ESD_Conversion_Mapping_Phi_in_Z_%02d",zBin);
1101     //    fHistograms->FillHistogram(nameESDMappingPhiInR, fV0Reader->GetMotherCandidatePhi());
1102     fHistograms->FillHistogram(nameESDMappingPhiInZ, vtxConv.Phi());
1103
1104     TString nameESDMappingRInZ="";
1105     nameESDMappingRInZ.Form("ESD_Conversion_Mapping_R_in_Z_%02d",zBin);
1106     fHistograms->FillHistogram(nameESDMappingRInZ, fV0Reader->GetXYRadius());
1107
1108     if(fV0Reader->GetMotherCandidatePt() > fLowPtMapping && fV0Reader->GetMotherCandidatePt()< fHighPtMapping){
1109       TString nameESDMappingMidPtPhiInR="";
1110       nameESDMappingMidPtPhiInR.Form("ESD_Conversion_Mapping_MidPt_Phi_in_R_%02d",rBin);
1111       fHistograms->FillHistogram(nameESDMappingMidPtPhiInR, vtxConv.Phi());
1112       
1113       TString nameESDMappingMidPtZInR="";
1114       nameESDMappingMidPtZInR.Form("ESD_Conversion_Mapping_MidPt_Z_in_R_%02d",rBin);
1115       fHistograms->FillHistogram(nameESDMappingMidPtZInR, fV0Reader->GetZ());
1116       
1117       TString nameESDMappingMidPtPhiInZ="";
1118       nameESDMappingMidPtPhiInZ.Form("ESD_Conversion_Mapping_MidPt_Phi_in_Z_%02d",zBin);
1119       fHistograms->FillHistogram(nameESDMappingMidPtPhiInZ, vtxConv.Phi());
1120       
1121       TString nameESDMappingMidPtRInZ="";
1122       nameESDMappingMidPtRInZ.Form("ESD_Conversion_Mapping_MidPt_R_in_Z_%02d",zBin);
1123       fHistograms->FillHistogram(nameESDMappingMidPtRInZ, fV0Reader->GetXYRadius());
1124
1125     }
1126
1127
1128     // end mapping
1129                 
1130     new((*fKFReconstructedGammasTClone)[fKFReconstructedGammasTClone->GetEntriesFast()])  AliKFParticle(*fV0Reader->GetMotherCandidateKFCombination());
1131                 
1132     //    fKFReconstructedGammas.push_back(*fV0Reader->GetMotherCandidateKFCombination());
1133     fElectronv1.push_back(fV0Reader->GetCurrentV0()->GetPindex());
1134     fElectronv2.push_back(fV0Reader->GetCurrentV0()->GetNindex());
1135                 
1136     //----------------------------------- checking for "real" conversions (MC match) --------------------------------------
1137     if(fDoMCTruth){
1138                         
1139       if(fV0Reader->HasSameMCMother() == kFALSE){
1140         continue;
1141       }
1142       TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
1143       TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
1144
1145       if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){
1146         continue;
1147       }
1148
1149       if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){
1150         continue;
1151       }
1152
1153       if(negativeMC->GetUniqueID() != 5 || positiveMC->GetUniqueID() !=5){// check if the daughters come from a conversion
1154         continue;
1155       }
1156                         
1157       if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
1158
1159         fHistograms->FillHistogram("ESD_TrueConvGamma_Pt", fV0Reader->GetMotherCandidatePt());
1160         fHistograms->FillHistogram("ESD_TrueConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());
1161         fHistograms->FillHistogram("ESD_TrueConvGamma_Eta", fV0Reader->GetMotherCandidateEta());                                
1162         fHistograms->FillHistogram("ESD_TrueConvGamma_Phi", fV0Reader->GetMotherCandidatePhi());
1163         fHistograms->FillHistogram("ESD_TrueConvGamma_Mass", fV0Reader->GetMotherCandidateMass());
1164         fHistograms->FillHistogram("ESD_TrueConvGamma_Width", fV0Reader->GetMotherCandidateWidth());
1165         fHistograms->FillHistogram("ESD_TrueConvGamma_Chi2", fV0Reader->GetMotherCandidateChi2());
1166         fHistograms->FillHistogram("ESD_TrueConvGamma_NDF", fV0Reader->GetMotherCandidateNDF());
1167         fHistograms->FillHistogram("ESD_TrueConvGamma_Pt_Eta", fV0Reader->GetMotherCandidatePt(),fV0Reader->GetMotherCandidateEta());
1168         fHistograms->FillHistogram("ESD_TrueConvGamma_Rapid", fV0Reader->GetMotherCandidateRapidity());
1169         fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLength", /*fV0Reader->GetNegativeTrackLength()*/fV0Reader->GetNegativeNTPCClusters());
1170         fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLength", /*fV0Reader->GetPositiveTrackLength()*/fV0Reader->GetPositiveNTPCClusters());
1171         fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLengthVSInvMass",/*fV0Reader->GetNegativeTrackLength()*/fV0Reader->GetNegativeNTPCClusters(),fV0Reader->GetMotherCandidateMass());
1172         fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLengthVSInvMass",/*fV0Reader->GetPositiveTrackLength()*/fV0Reader->GetPositiveNTPCClusters(),fV0Reader->GetMotherCandidateMass());
1173                                 
1174         fHistograms->FillHistogram("ESD_TrueConvGamma_Pt_Chi2", fV0Reader->GetMotherCandidatePt(), fV0Reader->GetMotherCandidateChi2());
1175         fHistograms->FillHistogram("ESD_TrueConvGamma_Eta_Chi2", fV0Reader->GetMotherCandidateEta(), fV0Reader->GetMotherCandidateChi2());
1176                                 
1177                                 
1178         fHistograms->FillHistogram("ESD_TrueConversion_XY", fV0Reader->GetX(),fV0Reader->GetY());
1179         fHistograms->FillHistogram("ESD_TrueConversion_R", fV0Reader->GetXYRadius());
1180         fHistograms->FillHistogram("ESD_TrueConversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());
1181         fHistograms->FillHistogram("ESD_TrueConversion_OpeningAngle", fV0Reader->GetOpeningAngle());
1182
1183         fHistograms->FillHistogram("ESD_TrueConvGamma_CosPointingAngle", fV0Reader->GetCosPointingAngle());
1184         fHistograms->FillHistogram("ESD_TrueConvGamma_DcaDaughters", fV0Reader->GetDcaDaughters());
1185         fHistograms->FillHistogram("ESD_TrueConvGamma_NormDcaDistDaughters", fV0Reader->GetNormDcaDistDaughters());
1186         fHistograms->FillHistogram("ESD_TrueConvGamma_LikelihoodAP", fV0Reader->GetLikelihoodAP());
1187
1188         fHistograms->FillHistogram("ESD_TrueConvGamma_E_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetNegativeTrackP()/fV0Reader->GetMotherCandidateP());
1189         fHistograms->FillHistogram("ESD_TrueConvGamma_P_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetPositiveTrackP()/fV0Reader->GetMotherCandidateP());
1190         fHistograms->FillHistogram("ESD_TrueConvGamma_E_dEdxP",fV0Reader->GetNegativeTrackP(),fV0Reader->GetNegativeTrackTPCdEdx());
1191         fHistograms->FillHistogram("ESD_TrueConvGamma_P_dEdxP",fV0Reader->GetPositiveTrackP(),fV0Reader->GetPositiveTrackTPCdEdx());
1192
1193
1194                                 
1195         //store MCTruth properties
1196         fHistograms->FillHistogram("ESD_TrueConvGamma_MC_Pt_Eta", fV0Reader->GetMotherMCParticle()->Pt(),fV0Reader->GetMotherMCParticle()->Eta());
1197         fHistograms->FillHistogram("ESD_TrueConversion_MC_ZR", negativeMC->Vz(),negativeMC->R());
1198         fHistograms->FillHistogram("ESD_TrueConversion_MC_XY", negativeMC->Vx(),negativeMC->Vy());
1199                                 
1200         //resolution
1201         Double_t mcpt   = fV0Reader->GetMotherMCParticle()->Pt();
1202         Double_t esdpt  = fV0Reader->GetMotherCandidatePt();
1203         Double_t resdPt = 0;
1204         if(mcpt > 0){
1205           resdPt = ((esdpt - mcpt)/mcpt)*100;
1206         }
1207         else if(mcpt < 0){
1208           cout<<"Pt of MC particle is negative, this will cause wrong calculation of resPt"<<endl; 
1209         }
1210                                 
1211         fHistograms->FillHistogram("Resolution_dPt", mcpt, resdPt);
1212         fHistograms->FillHistogram("Resolution_MC_Pt", mcpt);
1213         fHistograms->FillHistogram("Resolution_ESD_Pt", esdpt);
1214                                 
1215         Double_t resdZ = 0;
1216         if(fV0Reader->GetNegativeMCParticle()->Vz() != 0){
1217           resdZ = ((fV0Reader->GetZ() -fV0Reader->GetNegativeMCParticle()->Vz())/fV0Reader->GetNegativeMCParticle()->Vz())*100;
1218         }
1219                                 
1220         fHistograms->FillHistogram("Resolution_dZ", fV0Reader->GetNegativeMCParticle()->Vz(), resdZ);
1221         fHistograms->FillHistogram("Resolution_MC_Z", fV0Reader->GetNegativeMCParticle()->Vz());
1222         fHistograms->FillHistogram("Resolution_ESD_Z", fV0Reader->GetZ());
1223                                 
1224         Double_t resdR = 0;
1225         if(fV0Reader->GetNegativeMCParticle()->R() != 0){
1226           resdR = ((fV0Reader->GetXYRadius() - fV0Reader->GetNegativeMCParticle()->R())/fV0Reader->GetNegativeMCParticle()->R())*100;
1227         }
1228                                 
1229         fHistograms->FillHistogram("Resolution_dR", fV0Reader->GetNegativeMCParticle()->R(), resdR);
1230         fHistograms->FillHistogram("Resolution_MC_R", fV0Reader->GetNegativeMCParticle()->R());
1231         fHistograms->FillHistogram("Resolution_ESD_R", fV0Reader->GetXYRadius());
1232         fHistograms->FillHistogram("Resolution_dR_dPt", resdR, resdPt);
1233       }//if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22)
1234     }//if(fDoMCTruth)
1235   }//while(fV0Reader->NextV0)
1236   fHistograms->FillHistogram("ESD_NumberOfSurvivingV0s", nSurvivingV0s);
1237   fHistograms->FillHistogram("ESD_NumberOfV0s", fV0Reader->GetNumberOfV0s());
1238   fHistograms->FillHistogram("ESD_NumberOfContributorsVtx", fV0Reader->GetNumberOfContributorsVtx());
1239
1240   fV0Reader->ResetV0IndexNumber();
1241
1242 }
1243
1244 void AliAnalysisTaskGammaConversion::FillAODWithConversionGammas(){
1245   // Fill AOD with reconstructed Gamma
1246         
1247   for(Int_t gammaIndex=0;gammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();gammaIndex++){
1248     //  for(UInt_t gammaIndex=0;gammaIndex<fKFReconstructedGammas.size();gammaIndex++){
1249     //Create AOD particle object from AliKFParticle
1250                 
1251     /*    AliKFParticle * gammakf = &fKFReconstructedGammas[gammaIndex];
1252     //You could add directly AliKFParticle objects to the AOD, avoiding dependences with PartCorr
1253     //but this means that I have to work a little bit more in my side.
1254     //AODPWG4Particle objects are simpler and lighter, I think
1255     AliAODPWG4Particle gamma = AliAODPWG4Particle(gammakf->Px(),gammakf->Py(),gammakf->Pz(), gammakf->E());
1256     gamma.SetLabel(-1);//How to get the MC label of the reconstructed gamma?
1257     gamma.SetCaloLabel(-1,-1); //How to get the MC label of the 2 electrons that form the gamma?
1258     gamma.SetDetector("CTS"); //tag the gamma as reconstructed in the central barrel
1259     gamma.SetPdg(AliCaloPID::kPhotonConv); //photon id
1260     gamma.SetTag(-1); //Here I usually put a flag saying that montecarlo says it is prompt, decay fragmentation photon, or hadrons or whatever
1261                  
1262     //Add it to the aod list
1263     Int_t i = fAODBranch->GetEntriesFast();
1264     new((*fAODBranch)[i])  AliAODPWG4Particle(gamma);
1265     */
1266     //    AliKFParticle * gammakf = &fKFReconstructedGammas[gammaIndex];
1267     AliKFParticle * gammakf = (AliKFParticle *)fKFReconstructedGammasTClone->At(gammaIndex);
1268     AliGammaConversionAODObject aodObject;
1269     aodObject.SetPx(gammakf->GetPx());
1270     aodObject.SetPy(gammakf->GetPy());
1271     aodObject.SetPz(gammakf->GetPz());
1272     aodObject.SetLabel1(fElectronv1[gammaIndex]);
1273     aodObject.SetLabel2(fElectronv2[gammaIndex]);
1274     Int_t i = fAODBranch->GetEntriesFast();
1275     new((*fAODBranch)[i])  AliGammaConversionAODObject(aodObject);
1276   }
1277         
1278 }
1279
1280
1281 void AliAnalysisTaskGammaConversion::ProcessGammasForNeutralMesonAnalysis(){
1282   // see header file for documentation
1283         
1284   //  for(UInt_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammas.size();firstGammaIndex++){
1285   //    for(UInt_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fKFReconstructedGammas.size();secondGammaIndex++){
1286   for(Int_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();firstGammaIndex++){
1287     for(Int_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();secondGammaIndex++){
1288                         
1289       //      AliKFParticle * twoGammaDecayCandidateDaughter0 = &fKFReconstructedGammas[firstGammaIndex];
1290       //      AliKFParticle * twoGammaDecayCandidateDaughter1 = &fKFReconstructedGammas[secondGammaIndex];
1291                         
1292       AliKFParticle * twoGammaDecayCandidateDaughter0 = (AliKFParticle *)fKFReconstructedGammasTClone->At(firstGammaIndex);
1293       AliKFParticle * twoGammaDecayCandidateDaughter1 = (AliKFParticle *)fKFReconstructedGammasTClone->At(secondGammaIndex);
1294                         
1295       if(fElectronv1[firstGammaIndex]==fElectronv1[secondGammaIndex] || fElectronv1[firstGammaIndex]==fElectronv2[secondGammaIndex]){
1296         continue;
1297       }
1298       if(fElectronv2[firstGammaIndex]==fElectronv1[secondGammaIndex] || fElectronv2[firstGammaIndex]==fElectronv2[secondGammaIndex]){
1299         continue;
1300       }
1301                         
1302       AliKFParticle *twoGammaCandidate = new AliKFParticle(*twoGammaDecayCandidateDaughter0,*twoGammaDecayCandidateDaughter1);
1303                         
1304       Double_t massTwoGammaCandidate = 0.;
1305       Double_t widthTwoGammaCandidate = 0.;
1306       Double_t chi2TwoGammaCandidate =10000.;   
1307       twoGammaCandidate->GetMass(massTwoGammaCandidate,widthTwoGammaCandidate);
1308       if(twoGammaCandidate->GetNDF()>0){
1309         chi2TwoGammaCandidate = twoGammaCandidate->GetChi2()/twoGammaCandidate->GetNDF();
1310                                 
1311         if(chi2TwoGammaCandidate>0 && chi2TwoGammaCandidate<fV0Reader->GetChi2CutMeson()){
1312                                         
1313           TVector3 momentumVectorTwoGammaCandidate(twoGammaCandidate->GetPx(),twoGammaCandidate->GetPy(),twoGammaCandidate->GetPz());
1314           TVector3 spaceVectorTwoGammaCandidate(twoGammaCandidate->GetX(),twoGammaCandidate->GetY(),twoGammaCandidate->GetZ());
1315                                         
1316           Double_t openingAngleTwoGammaCandidate = twoGammaDecayCandidateDaughter0->GetAngle(*twoGammaDecayCandidateDaughter1);                                 
1317           Double_t rapidity;
1318           if(twoGammaCandidate->GetE() - twoGammaCandidate->GetPz() == 0 || twoGammaCandidate->GetE() + twoGammaCandidate->GetPz() == 0){
1319             rapidity=0;
1320           }
1321           else{
1322             rapidity = 0.5*(TMath::Log((twoGammaCandidate->GetE() +twoGammaCandidate->GetPz()) / (twoGammaCandidate->GetE()-twoGammaCandidate->GetPz())));
1323           }
1324                                         
1325           if(openingAngleTwoGammaCandidate < fMinOpeningAngleGhostCut){
1326             delete twoGammaCandidate;
1327             continue;   // minimum opening angle to avoid using ghosttracks
1328           }
1329                         
1330           fHistograms->FillHistogram("ESD_Mother_GammaDaughter_OpeningAngle", openingAngleTwoGammaCandidate);
1331           fHistograms->FillHistogram("ESD_Mother_Energy", twoGammaCandidate->GetE());
1332           fHistograms->FillHistogram("ESD_Mother_Pt", momentumVectorTwoGammaCandidate.Pt());
1333           fHistograms->FillHistogram("ESD_Mother_Eta", momentumVectorTwoGammaCandidate.Eta());
1334           fHistograms->FillHistogram("ESD_Mother_Rapid", rapidity);                                     
1335           fHistograms->FillHistogram("ESD_Mother_Phi", spaceVectorTwoGammaCandidate.Phi());
1336           fHistograms->FillHistogram("ESD_Mother_Mass", massTwoGammaCandidate);
1337           fHistograms->FillHistogram("ESD_Mother_R", spaceVectorTwoGammaCandidate.Pt());    // Pt in Space == R!!!
1338           fHistograms->FillHistogram("ESD_Mother_ZR", twoGammaCandidate->GetZ(), spaceVectorTwoGammaCandidate.Pt());
1339           fHistograms->FillHistogram("ESD_Mother_XY", twoGammaCandidate->GetX(), twoGammaCandidate->GetY());
1340           fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
1341           fHistograms->FillHistogram("ESD_Mother_InvMass",massTwoGammaCandidate);
1342
1343           if ( TMath::Abs(twoGammaDecayCandidateDaughter0->GetEta())<0.9 &&  TMath::Abs(twoGammaDecayCandidateDaughter1->GetEta())<0.9 ){
1344             fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_Fiducial",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
1345             fHistograms->FillHistogram("ESD_Mother_InvMass_Fiducial",massTwoGammaCandidate);
1346           }
1347         }
1348       }
1349       delete twoGammaCandidate;
1350     }
1351   }
1352 }
1353
1354 void AliAnalysisTaskGammaConversion::CalculateBackground(){
1355   // see header file for documentation
1356
1357
1358   TClonesArray * currentEventV0s = fV0Reader->GetCurrentEventGoodV0s();
1359
1360   for(Int_t nEventsInBG=0;nEventsInBG <fV0Reader->GetNBGEvents();nEventsInBG++){
1361     AliGammaConversionKFVector * previousEventV0s = fV0Reader->GetBGGoodV0s(nEventsInBG);
1362     for(Int_t iCurrent=0;iCurrent<currentEventV0s->GetEntriesFast();iCurrent++){
1363       AliKFParticle currentEventGoodV0 = *(AliKFParticle *)(currentEventV0s->At(iCurrent)); 
1364       for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
1365         AliKFParticle previousGoodV0 = (AliKFParticle)(*(previousEventV0s->at(iPrevious)));
1366
1367         AliKFParticle *backgroundCandidate = new AliKFParticle(currentEventGoodV0,previousGoodV0);
1368         
1369         Double_t massBG =0.;
1370         Double_t widthBG = 0.;
1371         Double_t chi2BG =10000.;        
1372         backgroundCandidate->GetMass(massBG,widthBG);
1373         if(backgroundCandidate->GetNDF()>0){
1374           chi2BG = backgroundCandidate->GetChi2()/backgroundCandidate->GetNDF();
1375           if(chi2BG>0 && chi2BG<fV0Reader->GetChi2CutMeson()){
1376                                         
1377             TVector3 momentumVectorbackgroundCandidate(backgroundCandidate->GetPx(),backgroundCandidate->GetPy(),backgroundCandidate->GetPz());
1378             TVector3 spaceVectorbackgroundCandidate(backgroundCandidate->GetX(),backgroundCandidate->GetY(),backgroundCandidate->GetZ());
1379                                         
1380             Double_t openingAngleBG = currentEventGoodV0.GetAngle(previousGoodV0);
1381                                         
1382             Double_t rapidity;
1383             if(backgroundCandidate->GetE() - backgroundCandidate->GetPz() == 0 || backgroundCandidate->GetE() + backgroundCandidate->GetPz() == 0) rapidity=0;
1384             else rapidity = 0.5*(TMath::Log((backgroundCandidate->GetE() +backgroundCandidate->GetPz()) / (backgroundCandidate->GetE()-backgroundCandidate->GetPz())));
1385                                         
1386                                         
1387                                         
1388                                         
1389             if(openingAngleBG < fMinOpeningAngleGhostCut ){
1390               delete backgroundCandidate;   
1391               continue;   // minimum opening angle to avoid using ghosttracks
1392             }                   
1393                                         
1394             fHistograms->FillHistogram("ESD_Background_GammaDaughter_OpeningAngle", openingAngleBG);
1395             fHistograms->FillHistogram("ESD_Background_Energy", backgroundCandidate->GetE());
1396             fHistograms->FillHistogram("ESD_Background_Pt",  momentumVectorbackgroundCandidate.Pt());
1397             fHistograms->FillHistogram("ESD_Background_Eta", momentumVectorbackgroundCandidate.Eta());
1398             fHistograms->FillHistogram("ESD_Background_Rapidity", rapidity);
1399             fHistograms->FillHistogram("ESD_Background_Phi", spaceVectorbackgroundCandidate.Phi());
1400             fHistograms->FillHistogram("ESD_Background_Mass", massBG);
1401             fHistograms->FillHistogram("ESD_Background_R", spaceVectorbackgroundCandidate.Pt());  // Pt in Space == R!!!!
1402             fHistograms->FillHistogram("ESD_Background_ZR", backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
1403             fHistograms->FillHistogram("ESD_Background_XY", backgroundCandidate->GetX(), backgroundCandidate->GetY());
1404             fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt",massBG,momentumVectorbackgroundCandidate.Pt());
1405             fHistograms->FillHistogram("ESD_Background_InvMass",massBG);
1406
1407             if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 &&  TMath::Abs(previousGoodV0.GetEta())<0.9 ){
1408               fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_Fiducial",massBG,momentumVectorbackgroundCandidate.Pt());
1409               fHistograms->FillHistogram("ESD_Background_InvMass_Fiducial",massBG);
1410             }
1411
1412           }
1413         }
1414         delete backgroundCandidate;      
1415       }
1416     }
1417   }
1418 }
1419
1420
1421 void AliAnalysisTaskGammaConversion::ProcessGammasForGammaJetAnalysis(){
1422   //ProcessGammasForGammaJetAnalysis
1423         
1424   Double_t distIsoMin;
1425         
1426   CreateListOfChargedParticles();
1427         
1428         
1429   //  for(UInt_t gammaIndex=0;gammaIndex<fKFReconstructedGammas.size();gammaIndex++){
1430   for(Int_t gammaIndex=0;gammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();gammaIndex++){
1431     AliKFParticle * currentGamma = (AliKFParticle*)fKFReconstructedGammasTClone->At(gammaIndex);
1432     TVector3 momentumVectorCurrentGamma(currentGamma->GetPx(),currentGamma->GetPy(),currentGamma->GetPz());
1433                 
1434     if( momentumVectorCurrentGamma.Pt()> fMinPtForGammaJet){
1435       distIsoMin=GetMinimumDistanceToCharge(gammaIndex);
1436                         
1437       if (distIsoMin > fMinIsoConeSize && fLeadingChargedIndex>=0){
1438         CalculateJetCone(gammaIndex);
1439       }
1440     }
1441   }
1442 }
1443
1444 void AliAnalysisTaskGammaConversion::CreateListOfChargedParticles(){
1445   // CreateListOfChargedParticles
1446         
1447   fESDEvent = fV0Reader->GetESDEvent();
1448   for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
1449     AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);
1450                 
1451     if(!curTrack){
1452       continue;
1453     }
1454                 
1455     if(fEsdTrackCuts->AcceptTrack(curTrack) ){
1456       new((*fChargedParticles)[fChargedParticles->GetEntriesFast()])  AliESDtrack(*curTrack);
1457       //      fChargedParticles.push_back(curTrack);
1458       fChargedParticlesId.push_back(iTracks);
1459     }
1460   }
1461 }
1462 void AliAnalysisTaskGammaConversion::CalculateJetCone(Int_t gammaIndex){
1463   // CaculateJetCone
1464         
1465   Double_t cone;
1466   Double_t coneSize=0.3;
1467   Double_t ptJet=0;
1468         
1469   //  AliKFParticle * currentGamma = &fKFReconstructedGammas[gammaIndex];
1470   AliKFParticle * currentGamma = (AliKFParticle*)fKFReconstructedGammasTClone->At(gammaIndex);
1471
1472   TVector3 momentumVectorCurrentGamma(currentGamma->GetPx(),currentGamma->GetPy(),currentGamma->GetPz());
1473         
1474   AliESDtrack* leadingCharged = (AliESDtrack*)(fChargedParticles->At(fLeadingChargedIndex));
1475
1476   Double_t momLeadingCharged[3];
1477   leadingCharged->GetConstrainedPxPyPz(momLeadingCharged);
1478         
1479   TVector3 momentumVectorLeadingCharged(momLeadingCharged[0],momLeadingCharged[1],momLeadingCharged[2]);
1480         
1481   Double_t phi1=momentumVectorLeadingCharged.Phi();
1482   Double_t eta1=momentumVectorLeadingCharged.Eta();
1483   Double_t phi3=momentumVectorCurrentGamma.Phi();
1484         
1485   for(Int_t iCh=0;iCh<fChargedParticles->GetEntriesFast();iCh++){
1486     AliESDtrack* curTrack = (AliESDtrack*)(fChargedParticles->At(iCh));
1487     Int_t chId = fChargedParticlesId[iCh];
1488     if(fLeadingChargedIndex==chId || fLeadingChargedIndex==chId) continue;
1489     Double_t mom[3];
1490     curTrack->GetConstrainedPxPyPz(mom);
1491     TVector3 momentumVectorChargedParticle(mom[0],mom[1],mom[2]);
1492     Double_t phi2=momentumVectorChargedParticle.Phi();
1493     Double_t eta2=momentumVectorChargedParticle.Eta();
1494                 
1495                 
1496     cone=100.;
1497     if( TMath::Abs(phi2 - phi1) <= ( TMath::TwoPi()-coneSize) ){
1498       cone = TMath::Sqrt(  TMath::Power((eta2-eta1),2)+ TMath::Power((phi2-phi1),2) );
1499     }else{
1500       if( (phi2 - phi1)> TMath::TwoPi()-coneSize ){
1501         cone = TMath::Sqrt(  TMath::Power((eta2-eta1),2)+ TMath::Power((phi2-TMath::TwoPi()-phi1),2) );
1502       }
1503       if( (phi2 - phi1)< -(TMath::TwoPi()-coneSize) ){
1504         cone = TMath::Sqrt(  TMath::Power((eta2-eta1),2)+ TMath::Power((phi2+TMath::TwoPi()-phi1),2) );
1505       }
1506     }
1507                 
1508     if(cone <coneSize&& momentumVectorChargedParticle.Pt()>fMinPtJetCone ){
1509       ptJet+= momentumVectorChargedParticle.Pt();
1510       Double_t ffzHdrGam = momentumVectorChargedParticle.Pt()/momentumVectorCurrentGamma.Pt();
1511       Double_t imbalanceHdrGam=-momentumVectorChargedParticle.Dot(momentumVectorCurrentGamma)/momentumVectorCurrentGamma.Mag2();
1512       fHistograms->FillHistogram("ESD_FFzHdrGam",ffzHdrGam);
1513       fHistograms->FillHistogram("ESD_ImbalanceHdrGam",imbalanceHdrGam);
1514                         
1515     }
1516                 
1517     Double_t dphiHdrGam=phi3-phi2;
1518     if ( dphiHdrGam < (-TMath::PiOver2())){
1519       dphiHdrGam+=(TMath::TwoPi());
1520     }
1521                 
1522     if ( dphiHdrGam > (3.*TMath::PiOver2()) ){
1523       dphiHdrGam-=(TMath::TwoPi());
1524     }
1525                 
1526     if (momentumVectorChargedParticle.Pt()>fMinPtGamChargedCorr){
1527       fHistograms->FillHistogram("ESD_dphiHdrGamIsolated",dphiHdrGam);
1528     }
1529   }//track loop
1530         
1531         
1532 }
1533
1534 Double_t AliAnalysisTaskGammaConversion::GetMinimumDistanceToCharge(Int_t indexHighestPtGamma){
1535   // GetMinimumDistanceToCharge
1536         
1537   Double_t fIsoMin=100.;
1538   Double_t ptLeadingCharged=-1.;
1539
1540   fLeadingChargedIndex=-1;
1541         
1542   AliKFParticle * gammaHighestPt = (AliKFParticle*)fKFReconstructedGammasTClone->At(indexHighestPtGamma);
1543   TVector3 momentumVectorgammaHighestPt(gammaHighestPt->GetPx(),gammaHighestPt->GetPy(),gammaHighestPt->GetPz());
1544         
1545   Double_t phi1=momentumVectorgammaHighestPt.Phi();
1546   Double_t eta1=momentumVectorgammaHighestPt.Eta();
1547         
1548   for(Int_t iCh=0;iCh<fChargedParticles->GetEntriesFast();iCh++){
1549     AliESDtrack* curTrack = (AliESDtrack*)(fChargedParticles->At(iCh));
1550     Int_t chId = fChargedParticlesId[iCh];
1551     if(fElectronv1[indexHighestPtGamma]==chId || fElectronv2[indexHighestPtGamma]==chId) continue;
1552     Double_t mom[3];
1553     curTrack->GetConstrainedPxPyPz(mom);
1554     TVector3 momentumVectorChargedParticle(mom[0],mom[1],mom[2]);
1555     Double_t phi2=momentumVectorChargedParticle.Phi();
1556     Double_t eta2=momentumVectorChargedParticle.Eta();
1557     Double_t iso=pow(  (pow( (eta1-eta2),2)+ pow((phi1-phi2),2)),0.5 );
1558                 
1559     if(momentumVectorChargedParticle.Pt()>fMinPtIsoCone ){
1560       if (iso<fIsoMin){
1561         fIsoMin=iso;
1562       }
1563     }
1564                 
1565     Double_t dphiHdrGam=phi1-phi2;
1566     if ( dphiHdrGam < (-TMath::PiOver2())){
1567       dphiHdrGam+=(TMath::TwoPi());
1568     }
1569                 
1570     if ( dphiHdrGam > (3.*TMath::PiOver2()) ){
1571       dphiHdrGam-=(TMath::TwoPi());
1572     }
1573     if (momentumVectorChargedParticle.Pt()>fMinPtGamChargedCorr){
1574       fHistograms->FillHistogram("ESD_dphiHdrGam",dphiHdrGam);
1575     }
1576                 
1577     if (dphiHdrGam>0.9*TMath::Pi() && dphiHdrGam<1.1*TMath::Pi()){
1578       if (momentumVectorChargedParticle.Pt()> ptLeadingCharged && momentumVectorChargedParticle.Pt()>0.1*momentumVectorgammaHighestPt.Pt()){
1579         ptLeadingCharged=momentumVectorChargedParticle.Pt();
1580         fLeadingChargedIndex=iCh;
1581       }
1582     }
1583                 
1584   }//track loop
1585   fHistograms->FillHistogram("ESD_MinimumIsoDistance",fIsoMin);
1586   return fIsoMin;
1587         
1588 }
1589
1590 Int_t  AliAnalysisTaskGammaConversion::GetIndexHighestPtGamma(){
1591   //GetIndexHighestPtGamma
1592         
1593   Int_t indexHighestPtGamma=-1;
1594   //Double_t 
1595   fGammaPtHighest = -100.;
1596         
1597   for(Int_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();firstGammaIndex++){
1598     AliKFParticle * gammaHighestPtCandidate = (AliKFParticle*)fKFReconstructedGammasTClone->At(firstGammaIndex);
1599     TVector3 momentumVectorgammaHighestPtCandidate(gammaHighestPtCandidate->GetPx(),gammaHighestPtCandidate->GetPy(),gammaHighestPtCandidate->GetPz());
1600     if (momentumVectorgammaHighestPtCandidate.Pt() > fGammaPtHighest){
1601       fGammaPtHighest=momentumVectorgammaHighestPtCandidate.Pt();
1602       //gammaHighestPt = gammaHighestPtCandidate;
1603       indexHighestPtGamma=firstGammaIndex;
1604     }
1605   }
1606         
1607   return indexHighestPtGamma;
1608         
1609 }
1610
1611
1612 void AliAnalysisTaskGammaConversion::Terminate(Option_t */*option*/)
1613 {
1614   // Terminate analysis
1615   //
1616   AliDebug(1,"Do nothing in Terminate");
1617 }
1618
1619 void AliAnalysisTaskGammaConversion::UserCreateOutputObjects()
1620 {
1621   //AOD
1622   if(fAODBranch==NULL){
1623     fAODBranch = new TClonesArray("AliGammaConversionAODObject", 0);
1624   }
1625   fAODBranch->Delete();
1626   fAODBranch->SetName(fAODBranchName); 
1627   AddAODBranch("TClonesArray", &fAODBranch);
1628         
1629   // Create the output container
1630   if(fOutputContainer != NULL){
1631     delete fOutputContainer;
1632     fOutputContainer = NULL;
1633   }
1634   if(fOutputContainer == NULL){
1635     fOutputContainer = new TList();
1636     fOutputContainer->SetOwner(kTRUE);
1637   }
1638         
1639   //Adding the histograms to the output container
1640   fHistograms->GetOutputContainer(fOutputContainer);
1641         
1642         
1643   if(fWriteNtuple){
1644     if(fGammaNtuple == NULL){
1645       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);
1646     }
1647     if(fNeutralMesonNtuple == NULL){
1648       fNeutralMesonNtuple = new TNtuple("NeutralMesonNtuple","NeutralMesonNtuple","test");
1649     }
1650     TList * ntupleTList = new TList();
1651     ntupleTList->SetOwner(kTRUE);
1652     ntupleTList->SetName("Ntuple");
1653     ntupleTList->Add((TNtuple*)fGammaNtuple);
1654     fOutputContainer->Add(ntupleTList);
1655   }
1656         
1657   fOutputContainer->SetName(GetName());
1658 }
1659
1660 Double_t AliAnalysisTaskGammaConversion::GetMCOpeningAngle(TParticle* const daughter0, TParticle* const daughter1) const{
1661   //helper function
1662   TVector3 v3D0(daughter0->Px(),daughter0->Py(),daughter0->Pz());
1663   TVector3 v3D1(daughter1->Px(),daughter1->Py(),daughter1->Pz());
1664   return v3D0.Angle(v3D1);
1665 }
1666
1667 void AliAnalysisTaskGammaConversion::CheckV0Efficiency(){
1668   // see header file for documentation
1669
1670   vector<Int_t> indexOfGammaParticle;
1671         
1672   fStack = fV0Reader->GetMCStack();
1673         
1674   if(fV0Reader->CheckForPrimaryVertex() == kFALSE){
1675     return; // aborts if the primary vertex does not have contributors.
1676   }
1677         
1678   for (Int_t iTracks = 0; iTracks < fStack->GetNprimary(); iTracks++) {
1679     TParticle* particle = (TParticle *)fStack->Particle(iTracks);
1680     if(particle->GetPdgCode()==22){     //Gamma
1681       if(particle->GetNDaughters() >= 2){
1682         TParticle* electron=NULL;
1683         TParticle* positron=NULL; 
1684         for(Int_t daughterIndex=particle->GetFirstDaughter();daughterIndex<=particle->GetLastDaughter();daughterIndex++){
1685           TParticle *tmpDaughter = fStack->Particle(daughterIndex);
1686           if(tmpDaughter->GetUniqueID() == 5){
1687             if(tmpDaughter->GetPdgCode() == 11){
1688               electron = tmpDaughter;
1689             }
1690             else if(tmpDaughter->GetPdgCode() == -11){
1691               positron = tmpDaughter;
1692             }
1693           }
1694         }
1695         if(electron!=NULL && positron!=0){
1696           if(electron->R()<160){
1697             indexOfGammaParticle.push_back(iTracks);
1698           }
1699         }
1700       }
1701     }
1702   }
1703         
1704   Int_t nFoundGammas=0;
1705   Int_t nNotFoundGammas=0;
1706         
1707   Int_t numberOfV0s = fV0Reader->GetNumberOfV0s();
1708   for(Int_t i=0;i<numberOfV0s;i++){
1709     fV0Reader->GetV0(i);
1710                 
1711     if(fV0Reader->HasSameMCMother() == kFALSE){
1712       continue;
1713     }
1714                 
1715     TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
1716     TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
1717                 
1718     if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){
1719       continue;
1720     }
1721     if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){
1722       continue;
1723     }
1724                 
1725     if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
1726       //TParticle * v0Gamma = fV0Reader->GetMotherMCParticle();
1727       for(UInt_t mcIndex=0;mcIndex<indexOfGammaParticle.size();mcIndex++){
1728         if(negativeMC->GetFirstMother()==indexOfGammaParticle[mcIndex]){
1729           nFoundGammas++;
1730         }
1731         else{
1732           nNotFoundGammas++;
1733         }
1734       }
1735     }
1736   }
1737 }
1738
1739
1740 void AliAnalysisTaskGammaConversion::ProcessGammaElectronsForChicAnalysis(){
1741   // see header file for documantation
1742         
1743   fESDEvent = fV0Reader->GetESDEvent();
1744         
1745         
1746   TClonesArray * vESDeNegTemp = new TClonesArray("AliESDtrack",0);
1747   TClonesArray * vESDePosTemp = new TClonesArray("AliESDtrack",0);
1748   TClonesArray * vESDxNegTemp = new TClonesArray("AliESDtrack",0);
1749   TClonesArray * vESDxPosTemp = new TClonesArray("AliESDtrack",0);
1750   TClonesArray * vESDeNegNoJPsi = new TClonesArray("AliESDtrack",0);
1751   TClonesArray * vESDePosNoJPsi = new TClonesArray("AliESDtrack",0);
1752         
1753   /*
1754     vector <AliESDtrack*> vESDeNegTemp(0);
1755     vector <AliESDtrack*> vESDePosTemp(0);
1756     vector <AliESDtrack*> vESDxNegTemp(0);
1757     vector <AliESDtrack*> vESDxPosTemp(0);
1758     vector <AliESDtrack*> vESDeNegNoJPsi(0);
1759     vector <AliESDtrack*> vESDePosNoJPsi(0); 
1760   */
1761         
1762         
1763   fHistograms->FillTable("Table_Electrons",0);//Count number of Events
1764         
1765   for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
1766     AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);
1767                 
1768     if(!curTrack){
1769       //print warning here
1770       continue;
1771     }
1772                 
1773     double p[3];if(!curTrack->GetConstrainedPxPyPz(p))continue;
1774     double r[3];curTrack->GetConstrainedXYZ(r);
1775                 
1776     TVector3 rXYZ(r);
1777                 
1778     fHistograms->FillTable("Table_Electrons",4);//Count number of ESD tracks
1779                 
1780     Bool_t flagKink       =  kTRUE;
1781     Bool_t flagTPCrefit   =  kTRUE;
1782     Bool_t flagTRDrefit   =  kTRUE;
1783     Bool_t flagITSrefit   =  kTRUE;
1784     Bool_t flagTRDout     =  kTRUE;
1785     Bool_t flagVertex     =  kTRUE;
1786                 
1787                 
1788     //Cuts ---------------------------------------------------------------
1789                 
1790     if(curTrack->GetKinkIndex(0) > 0){
1791       fHistograms->FillHistogram("Table_Electrons",5);//Count kink
1792       flagKink = kFALSE;
1793     }
1794                 
1795     ULong_t trkStatus = curTrack->GetStatus();
1796                 
1797     ULong_t tpcRefit = (trkStatus & AliESDtrack::kTPCrefit);
1798                 
1799     if(!tpcRefit){
1800       fHistograms->FillHistogram("Table_Electrons",9);//Count not TPCrefit
1801       flagTPCrefit = kFALSE;
1802     }
1803                 
1804     ULong_t itsRefit = (trkStatus & AliESDtrack::kITSrefit);
1805     if(!itsRefit){
1806       fHistograms->FillHistogram("Table_Electrons",10);//Count not ITSrefit
1807       flagITSrefit = kFALSE;
1808     }
1809                 
1810     ULong_t trdRefit = (trkStatus & AliESDtrack::kTRDrefit);
1811                 
1812     if(!trdRefit){
1813       fHistograms->FillHistogram("Table_Electrons",8); //Count not TRDrefit
1814       flagTRDrefit = kFALSE;
1815     }
1816                 
1817     ULong_t trdOut = (trkStatus & AliESDtrack::kTRDout);
1818                 
1819     if(!trdOut) {
1820       fHistograms->FillHistogram("Table_Electrons",7); //Count not TRDout
1821       flagTRDout = kFALSE;
1822     }
1823                 
1824     double nsigmaToVxt = GetSigmaToVertex(curTrack);
1825                 
1826     if(nsigmaToVxt > 3){
1827       fHistograms->FillHistogram("Table_Electrons",6); //Count Tracks with number of sigmas > 3
1828       flagVertex = kFALSE;
1829     }
1830                 
1831     if(! (flagKink && flagTPCrefit && flagITSrefit && flagTRDrefit && flagTRDout && flagVertex ) ) continue;
1832     fHistograms->FillHistogram("Table_Electrons",11);//Count Tracks passed Cuts
1833                 
1834                 
1835     Stat_t pid, weight;
1836     GetPID(curTrack, pid, weight);
1837                 
1838     if(pid!=0){
1839       fHistograms->FillHistogram("Table_Electrons",12); //Count Tracks with pid != 0
1840     }
1841                 
1842     if(pid == 0){
1843       fHistograms->FillHistogram("Table_Electrons",13); //Count Tracks with pid != 0
1844     }
1845                 
1846                 
1847                 
1848                 
1849                 
1850                 
1851     TLorentzVector curElec;
1852     curElec.SetXYZM(p[0],p[1],p[2],fElectronMass);
1853                 
1854                 
1855     if(fDoMCTruth){             
1856       Int_t labelMC = TMath::Abs(curTrack->GetLabel());
1857       TParticle* curParticle = fStack->Particle(labelMC);
1858       if(curTrack->GetSign() > 0){
1859         if( pid == 0){
1860           fHistograms->FillHistogram("MC_ElectronPosNegPt",curParticle->Pt());
1861           fHistograms->FillHistogram("MC_ElectronPosNegEta",curParticle->Eta());
1862         }
1863         else{
1864           fHistograms->FillHistogram("MC_ElectronPosNegPt",curParticle->Pt());
1865           fHistograms->FillHistogram("MC_ElectronPosNegEta",curParticle->Eta());
1866         }
1867       }
1868     }
1869                 
1870                 
1871     if(curTrack->GetSign() > 0){
1872                         
1873       //     vESDxPosTemp.push_back(curTrack);
1874       new((*vESDxPosTemp)[vESDxPosTemp->GetEntriesFast()])  AliESDtrack(*curTrack);
1875                         
1876       if( pid == 0){
1877                                 
1878         fHistograms->FillHistogram("ESD_ElectronPosNegPt",curElec.Pt());
1879         fHistograms->FillHistogram("ESD_ElectronPosPt",curElec.Pt());
1880         //      fHistograms->FillHistogram("MC_ElectronPosNegPt",curParticle->Pt());
1881         fHistograms->FillHistogram("ESD_ElectronPosNegEta",curElec.Eta());
1882         //      fHistograms->FillHistogram("MC_ElectronPosNegEta",curParticle->Eta());
1883         //      vESDePosTemp.push_back(curTrack);
1884         new((*vESDePosTemp)[vESDePosTemp->GetEntriesFast()])  AliESDtrack(*curTrack);
1885                                 
1886       }
1887                         
1888     }
1889     else {
1890
1891       new((*vESDxNegTemp)[vESDxNegTemp->GetEntriesFast()])  AliESDtrack(*curTrack);
1892                         
1893       if( pid == 0){
1894                                 
1895         fHistograms->FillHistogram("ESD_ElectronPosNegPt",curElec.Pt());
1896         fHistograms->FillHistogram("ESD_ElectronNegPt",curElec.Pt());
1897         fHistograms->FillHistogram("ESD_ElectronPosNegEta",curElec.Eta());
1898         new((*vESDeNegTemp)[vESDeNegTemp->GetEntriesFast()])  AliESDtrack(*curTrack);
1899                                 
1900       }
1901                         
1902     }
1903                 
1904   }
1905         
1906         
1907   Bool_t ePosJPsi = kFALSE;
1908   Bool_t eNegJPsi = kFALSE;             
1909   Bool_t ePosPi0  = kFALSE;
1910   Bool_t eNegPi0  = kFALSE;
1911         
1912   UInt_t iePosJPsi=0,ieNegJPsi=0,iePosPi0=0,ieNegPi0=0;
1913         
1914   for(Int_t iNeg=0; iNeg < vESDeNegTemp->GetEntriesFast(); iNeg++){
1915     if(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDeNegTemp->At(iNeg)))->GetLabel()))->GetPdgCode() == 11)
1916       if(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDeNegTemp->At(iNeg)))->GetLabel()))->GetMother(0) > -1){
1917         Int_t labelMother = fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDeNegTemp->At(iNeg)))->GetLabel()))->GetMother(0);
1918         TParticle* partMother = fStack ->Particle(labelMother);
1919         if (partMother->GetPdgCode() == 111){
1920           ieNegPi0 = iNeg;
1921           eNegPi0 = kTRUE;
1922         }
1923         if(partMother->GetPdgCode() == 443){ //Mother JPsi
1924           fHistograms->FillTable("Table_Electrons",14);
1925           ieNegJPsi = iNeg;
1926           eNegJPsi = kTRUE;
1927         }
1928         else{   
1929           //      vESDeNegNoJPsi.push_back(vESDeNegTemp[iNeg]);
1930           new((*vESDeNegNoJPsi)[vESDeNegNoJPsi->GetEntriesFast()])  AliESDtrack(*(AliESDtrack*)(vESDeNegTemp->At(iNeg)));
1931           //            cout<<"ESD No Positivo JPsi "<<endl;
1932         }
1933                                 
1934       }
1935   }     
1936         
1937   for(Int_t iPos=0; iPos < vESDePosTemp->GetEntriesFast(); iPos++){
1938     if(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDePosTemp->At(iPos)))->GetLabel()))->GetPdgCode() == -11)
1939       if(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDePosTemp->At(iPos)))->GetLabel()))->GetMother(0) > -1){
1940         Int_t labelMother = fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDePosTemp->At(iPos)))->GetLabel()))->GetMother(0);
1941         TParticle* partMother = fStack ->Particle(labelMother);
1942         if (partMother->GetPdgCode() == 111){
1943           iePosPi0 = iPos;
1944           ePosPi0 = kTRUE;
1945         }
1946         if(partMother->GetPdgCode() == 443){ //Mother JPsi
1947           fHistograms->FillTable("Table_Electrons",15);
1948           iePosJPsi = iPos;
1949           ePosJPsi = kTRUE;
1950         }
1951         else{
1952           //      vESDePosNoJPsi.push_back(vESDePosTemp[iPos]);
1953           new((*vESDePosNoJPsi)[vESDePosNoJPsi->GetEntriesFast()])  AliESDtrack(*(AliESDtrack*)(vESDePosTemp->At(iPos)));         
1954           //            cout<<"ESD No Negativo JPsi "<<endl;
1955         }
1956                                 
1957       }
1958   }
1959         
1960   if( eNegJPsi && ePosJPsi ){
1961     TVector3 tempeNegV,tempePosV;
1962     tempeNegV.SetXYZ(((AliESDtrack*)(vESDeNegTemp->At(ieNegJPsi)))->Px(),((AliESDtrack*)(vESDeNegTemp->At(ieNegJPsi)))->Py(),((AliESDtrack*)(vESDeNegTemp->At(ieNegJPsi)))->Pz());                      
1963     tempePosV.SetXYZ(((AliESDtrack*)(vESDePosTemp->At(iePosJPsi)))->Px(),((AliESDtrack*)(vESDePosTemp->At(iePosJPsi)))->Py(),((AliESDtrack*)(vESDePosTemp->At(iePosJPsi)))->Pz());
1964     fHistograms->FillTable("Table_Electrons",16);
1965     fHistograms->FillHistogram("ESD_ElectronPosNegJPsiAngle",tempeNegV.Angle(tempePosV));       
1966     fHistograms->FillHistogram("MC_ElectronPosNegJPsiAngle",GetMCOpeningAngle(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDeNegTemp->At(ieNegJPsi)))->GetLabel())),
1967                                                                               fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDePosTemp->At(iePosJPsi)))->GetLabel()))));        
1968   }
1969         
1970   if( eNegPi0 && ePosPi0 ){
1971     TVector3 tempeNegV,tempePosV;
1972     tempeNegV.SetXYZ(((AliESDtrack*)(vESDeNegTemp->At(ieNegPi0)))->Px(),((AliESDtrack*)(vESDeNegTemp->At(ieNegPi0)))->Py(),((AliESDtrack*)(vESDeNegTemp->At(ieNegPi0)))->Pz());
1973     tempePosV.SetXYZ(((AliESDtrack*)(vESDePosTemp->At(iePosPi0)))->Px(),((AliESDtrack*)(vESDePosTemp->At(iePosPi0)))->Py(),((AliESDtrack*)(vESDePosTemp->At(iePosPi0)))->Pz());
1974     fHistograms->FillHistogram("ESD_ElectronPosNegPi0Angle",tempeNegV.Angle(tempePosV));
1975     fHistograms->FillHistogram("MC_ElectronPosNegPi0Angle",GetMCOpeningAngle(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDeNegTemp->At(ieNegPi0)))->GetLabel())),
1976                                                                              fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDePosTemp->At(iePosPi0)))->GetLabel()))));   
1977   }
1978         
1979         
1980   FillAngle("ESD_eNegePosAngleBeforeCut",GetTLorentzVector(vESDeNegTemp),GetTLorentzVector(vESDePosTemp));
1981         
1982   CleanWithAngleCuts(*vESDeNegTemp,*vESDePosTemp,*fKFReconstructedGammasTClone);
1983         
1984   //  vector <TLorentzVector> vCurrentTLVeNeg = GetTLorentzVector(fCurrentEventNegElectron);
1985   //  vector <TLorentzVector> vCurrentTLVePos = GetTLorentzVector(fCurrentEventPosElectron);
1986         
1987   TClonesArray vCurrentTLVeNeg = GetTLorentzVector(fCurrentEventNegElectronTClone);
1988   TClonesArray vCurrentTLVePos = GetTLorentzVector(fCurrentEventPosElectronTClone);
1989         
1990         
1991   FillAngle("ESD_eNegePosAngleAfterCut",vCurrentTLVeNeg,vCurrentTLVePos);
1992         
1993         
1994         
1995         
1996   //FillAngle("ESD_eNegePosAngleAfterCut",CurrentTLVeNeg,CurrentTLVePos);
1997         
1998         
1999   FillElectronInvMass("ESD_InvMass_ePluseMinus",vCurrentTLVeNeg,vCurrentTLVePos);
2000   FillElectronInvMass("ESD_InvMass_xPlusxMinus",GetTLorentzVector(vESDxNegTemp),GetTLorentzVector(vESDxPosTemp));
2001         
2002         
2003         
2004   FillGammaElectronInvMass("ESD_InvMass_GammaePluseMinusChiC","ESD_InvMass_GammaePluseMinusChiCDiff",*fKFReconstructedGammasCutTClone,vCurrentTLVeNeg,vCurrentTLVePos);
2005         
2006   FillGammaElectronInvMass("ESD_InvMass_GammaePluseMinusPi0","ESD_InvMass_GammaePluseMinusPi0Diff",
2007                            *fKFReconstructedGammasCutTClone,vCurrentTLVeNeg,vCurrentTLVePos);
2008         
2009   //BackGround
2010         
2011   //Like Sign e+e-
2012   ElectronBackground("ESD_ENegBackground",vCurrentTLVeNeg);
2013   ElectronBackground("ESD_EPosBackground",vCurrentTLVePos);
2014   ElectronBackground("ESD_EPosENegBackground",vCurrentTLVeNeg);
2015   ElectronBackground("ESD_EPosENegBackground",vCurrentTLVePos);
2016         
2017   //        Like Sign e+e- no JPsi
2018   ElectronBackground("ESD_EPosENegNoJPsiBG",GetTLorentzVector(vESDeNegNoJPsi));
2019   ElectronBackground("ESD_EPosENegNoJPsiBG",GetTLorentzVector(vESDePosNoJPsi));
2020         
2021   //Mixed Event
2022         
2023   if( fCurrentEventPosElectronTClone->GetEntriesFast() > 0 && fCurrentEventNegElectronTClone->GetEntriesFast() > 0 && fKFReconstructedGammasCutTClone->GetEntriesFast() > 0 ){
2024     FillGammaElectronInvMass("ESD_EPosENegGammaBackgroundMX","ESD_EPosENegGammaBackgroundMXDiff",
2025                              *fKFReconstructedGammasCutTClone,*fPreviousEventTLVNegElectronTClone,*fPreviousEventTLVPosElectronTClone);
2026     *fPreviousEventTLVNegElectronTClone = vCurrentTLVeNeg;
2027     *fPreviousEventTLVPosElectronTClone = vCurrentTLVePos;
2028                 
2029   }
2030         
2031   /*
2032   //Photons P
2033   Double_t vtx[3];
2034   vtx[0]=0;vtx[1]=0;vtx[2]=0;
2035   for(UInt_t i=0;i<fKFReconstructedGammasChic.size();i++){
2036          
2037   //      if(fMCGammaChicTempCut[i]->GetMother(0) < 0) continue;
2038          
2039          
2040          
2041   Int_t tempLabel = fStack->Particle(fMCGammaChicTempCut[i]->GetMother(0))->GetPdgCode();
2042   //      cout<<" Label Pedro Gonzalez " <<tempLabel <<endl;
2043          
2044   //      cout<<" Label Distance"<<fKFReconstructedGammasChic[i].GetDistanceFromVertex(vtx)<<endl;
2045          
2046   if( tempLabel == 10441 || tempLabel == 20443 || tempLabel == 445 )
2047          
2048   fHistograms->FillHistogram("ESD_PhotonsMomentum",fKFReconstructedGammasChic[i].GetMomentum());
2049          
2050          
2051   }
2052          
2053          
2054   */
2055
2056
2057   vESDeNegTemp->Delete();
2058   vESDePosTemp->Delete();
2059   vESDxNegTemp->Delete();
2060   vESDxPosTemp->Delete();
2061   vESDeNegNoJPsi->Delete();
2062   vESDePosNoJPsi->Delete();
2063
2064   delete vESDeNegTemp;
2065   delete vESDePosTemp;
2066   delete vESDxNegTemp;
2067   delete vESDxPosTemp;
2068   delete vESDeNegNoJPsi;
2069   delete vESDePosNoJPsi;        
2070 }
2071
2072 /*
2073   void AliAnalysisTaskGammaConversion::FillAngle(TString histoName,vector <TLorentzVector> tlVeNeg, vector <TLorentzVector> tlVePos){
2074   //see header file for documentation
2075   for( UInt_t iNeg=0; iNeg < tlVeNeg.size(); iNeg++){
2076   for (UInt_t iPos=0; iPos < tlVePos.size(); iPos++){
2077   fHistograms->FillHistogram(histoName.Data(),tlVeNeg[iNeg].Vect().Angle(tlVePos[iPos].Vect()));
2078   }
2079   }
2080   }
2081 */
2082 void AliAnalysisTaskGammaConversion::FillAngle(TString histoName,TClonesArray const tlVeNeg, TClonesArray const tlVePos){
2083   //see header file for documentation
2084   for( Int_t iNeg=0; iNeg < tlVeNeg.GetEntriesFast(); iNeg++){
2085     for (Int_t iPos=0; iPos < tlVePos.GetEntriesFast(); iPos++){
2086       fHistograms->FillHistogram(histoName.Data(),((TLorentzVector*)(tlVeNeg.At(iNeg)))->Vect().Angle(((TLorentzVector*)(tlVePos.At(iPos)))->Vect()));
2087     }
2088   }
2089 }
2090 void AliAnalysisTaskGammaConversion::FillElectronInvMass(TString histoName, TClonesArray const eNeg, TClonesArray const ePos){
2091   //see header file for documentation
2092   for( Int_t n=0; n < eNeg.GetEntriesFast(); n++){
2093     TLorentzVector en = (*(TLorentzVector*)(eNeg.At(n)));
2094     for (Int_t p=0; p < ePos.GetEntriesFast(); p++){
2095       TLorentzVector ep = (*(TLorentzVector*)(ePos.At(p)));
2096       TLorentzVector np = ep + en;
2097       fHistograms->FillHistogram(histoName.Data(),np.M());
2098     }
2099   }
2100 }
2101
2102 void AliAnalysisTaskGammaConversion::FillGammaElectronInvMass(TString histoMass,TString histoDiff,TClonesArray const fKFGammas,
2103                                                               TClonesArray const tlVeNeg,TClonesArray const tlVePos)
2104 {
2105   //see header file for documentation
2106         
2107   for( Int_t iNeg=0; iNeg < tlVeNeg.GetEntriesFast(); iNeg++ ){
2108                 
2109     for (Int_t iPos=0; iPos < tlVePos.GetEntriesFast(); iPos++){
2110                         
2111       TLorentzVector xy = *((TLorentzVector *)(tlVePos.At(iPos))) + *((TLorentzVector *)(tlVeNeg.At(iNeg)));
2112                         
2113       for (Int_t iGam=0; iGam < fKFGammas.GetEntriesFast(); iGam++){
2114                                 
2115         //      AliKFParticle * gammaCandidate = &fKFGammas[iGam];
2116         AliKFParticle * gammaCandidate = (AliKFParticle *)(fKFGammas.At(iGam));
2117         TLorentzVector g;
2118                                 
2119         g.SetXYZM(gammaCandidate->GetPx(),gammaCandidate->GetPy(),gammaCandidate->GetPz(),fGammaMass);
2120         TLorentzVector xyg = xy + g;
2121         fHistograms->FillHistogram(histoMass.Data(),xyg.M());
2122         fHistograms->FillHistogram(histoDiff.Data(),(xyg.M()-xy.M()));
2123       }
2124     }
2125   }
2126         
2127 }
2128 void AliAnalysisTaskGammaConversion::ElectronBackground(TString hBg, TClonesArray e)
2129 {
2130   // see header file for documentation
2131   for(Int_t i=0; i < e.GetEntriesFast(); i++)
2132     {
2133       for (Int_t j=i+1; j < e.GetEntriesFast(); j++)
2134         {
2135           TLorentzVector ee = (*(TLorentzVector*)(e.At(i))) + (*(TLorentzVector*)(e.At(j)));
2136                         
2137           fHistograms->FillHistogram(hBg.Data(),ee.M());
2138         }
2139     }
2140 }
2141
2142
2143 void AliAnalysisTaskGammaConversion::CleanWithAngleCuts(TClonesArray const negativeElectrons,
2144                                                         TClonesArray const positiveElectrons, 
2145                                                         TClonesArray const gammas){
2146   // see header file for documentation
2147         
2148   UInt_t  sizeN = negativeElectrons.GetEntriesFast();
2149   UInt_t  sizeP = positiveElectrons.GetEntriesFast();
2150   UInt_t  sizeG = gammas.GetEntriesFast();
2151         
2152         
2153         
2154   vector <Bool_t> xNegBand(sizeN);
2155   vector <Bool_t> xPosBand(sizeP);
2156   vector <Bool_t> gammaBand(sizeG);
2157         
2158         
2159   for(UInt_t iNeg=0; iNeg < sizeN; iNeg++) xNegBand[iNeg]=kTRUE;
2160   for(UInt_t iPos=0; iPos < sizeP; iPos++) xPosBand[iPos]=kTRUE;
2161   for(UInt_t iGam=0; iGam < sizeG; iGam++) gammaBand[iGam]=kTRUE;
2162         
2163         
2164   for(UInt_t iPos=0; iPos < sizeP; iPos++){
2165                 
2166     Double_t aP[3]; 
2167     ((AliESDtrack*)(positiveElectrons.At(iPos)))->GetConstrainedPxPyPz(aP); 
2168                 
2169     TVector3 ePosV(aP[0],aP[1],aP[2]);
2170                 
2171     for(UInt_t iNeg=0; iNeg < sizeN; iNeg++){
2172                         
2173       Double_t aN[3]; 
2174       ((AliESDtrack*)(negativeElectrons.At(iNeg)))->GetConstrainedPxPyPz(aN); 
2175       TVector3 eNegV(aN[0],aN[1],aN[2]);
2176                         
2177       if(ePosV.Angle(eNegV) < 0.05){ //e+e- from gamma
2178         xPosBand[iPos]=kFALSE;
2179         xNegBand[iNeg]=kFALSE;
2180       }
2181                         
2182       for(UInt_t iGam=0; iGam < sizeG; iGam++){
2183         AliKFParticle* gammaCandidate = (AliKFParticle*)gammas.At(iGam);
2184         TVector3 gammaCandidateVector(gammaCandidate->Px(),gammaCandidate->Py(),gammaCandidate->Pz());
2185         if(ePosV.Angle(gammaCandidateVector) < 0.05 || eNegV.Angle(gammaCandidateVector) < 0.05)
2186           gammaBand[iGam]=kFALSE;
2187       }
2188     }
2189   }
2190         
2191         
2192         
2193         
2194   for(UInt_t iPos=0; iPos < sizeP; iPos++){
2195     if(xPosBand[iPos]){
2196       new((*fCurrentEventPosElectronTClone)[fCurrentEventPosElectronTClone->GetEntriesFast()]) AliESDtrack((*(AliESDtrack*)(positiveElectrons.At(iPos))));
2197       //      fCurrentEventPosElectron.push_back(positiveElectrons[iPos]);
2198     }
2199   }
2200   for(UInt_t iNeg=0;iNeg < sizeN; iNeg++){
2201     if(xNegBand[iNeg]){
2202       new((*fCurrentEventNegElectronTClone)[fCurrentEventNegElectronTClone->GetEntriesFast()]) AliESDtrack((*(AliESDtrack*)(negativeElectrons.At(iNeg))));
2203       //      fCurrentEventNegElectron.push_back(negativeElectrons[iNeg]);
2204     }
2205   }
2206   for(UInt_t iGam=0; iGam < sizeG; iGam++){
2207     if(gammaBand[iGam]){
2208       new((*fKFReconstructedGammasCutTClone)[fKFReconstructedGammasCutTClone->GetEntriesFast()]) AliKFParticle((*(AliKFParticle*)(gammas.At(iGam))));
2209       //fKFReconstructedGammasCut.push_back(*(AliKFParticle*)gammas->At(iGam));
2210     }
2211   }
2212 }
2213
2214
2215 void  AliAnalysisTaskGammaConversion::GetPID(AliESDtrack *track, Stat_t &pid, Stat_t &weight)
2216 {
2217   // see header file for documentation
2218   pid = -1;
2219   weight = -1;
2220         
2221   double wpart[5];
2222   double wpartbayes[5];
2223         
2224   //get probability of the diffenrent particle types
2225   track->GetESDpid(wpart);
2226         
2227   // Tentative particle type "concentrations"
2228   double c[5]={0.01, 0.01, 0.85, 0.10, 0.05};
2229         
2230   //Bayes' formula
2231   double rcc = 0.;
2232   for (int i = 0; i < 5; i++)
2233     {
2234       rcc+=(c[i] * wpart[i]);
2235     }
2236         
2237         
2238         
2239   for (int i=0; i<5; i++) {
2240     if( rcc>0 || rcc<0){//Kenneth: not sure if the rcc<0 should be there, this is from fixing a coding violation where we are not allowed to say: rcc!=0 (RC19)  
2241       wpartbayes[i] = c[i] * wpart[i] / rcc;
2242     }
2243   }
2244         
2245         
2246         
2247   Float_t max=0.;
2248   int ipid=-1;
2249   //find most probable particle in ESD pid
2250   //0:Electron - 1:Muon - 2:Pion - 3:Kaon - 4:Proton
2251   for (int i = 0; i < 5; i++)
2252     {
2253       if (wpartbayes[i] > max)
2254         {
2255           ipid = i;
2256           max = wpartbayes[i];
2257         }
2258     }
2259         
2260   pid = ipid;
2261   weight = max;
2262 }
2263 double AliAnalysisTaskGammaConversion::GetSigmaToVertex(AliESDtrack* t)
2264 {
2265   // Calculates the number of sigma to the vertex.
2266         
2267   Float_t b[2];
2268   Float_t bRes[2];
2269   Float_t bCov[3];
2270   t->GetImpactParameters(b,bCov);
2271   if (bCov[0]<=0 || bCov[2]<=0) {
2272     AliDebug(1, "Estimated b resolution lower or equal zero!");
2273     bCov[0]=0; bCov[2]=0;
2274   }
2275   bRes[0] = TMath::Sqrt(bCov[0]);
2276   bRes[1] = TMath::Sqrt(bCov[2]);
2277         
2278   // -----------------------------------
2279   // How to get to a n-sigma cut?
2280   //
2281   // The accumulated statistics from 0 to d is
2282   //
2283   // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
2284   // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
2285   //
2286   // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
2287   // Can this be expressed in a different way?
2288         
2289   if (bRes[0] == 0 || bRes[1] ==0)
2290     return -1;
2291         
2292   double d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
2293         
2294   // stupid rounding problem screws up everything:
2295   // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
2296   if (TMath::Exp(-d * d / 2) < 1e-10)
2297     return 1000;
2298         
2299         
2300   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
2301   return d;
2302 }
2303
2304 //vector <TLorentzVector> AliAnalysisTaskGammaConversion::GetTLorentzVector(vector <AliESDtrack*> esdTrack){
2305 TClonesArray AliAnalysisTaskGammaConversion::GetTLorentzVector(TClonesArray *const esdTrack){
2306   //Return TLoresntz vector of track?
2307   //  vector <TLorentzVector> tlVtrack(0);
2308   TClonesArray array("TLorentzVector",0); 
2309         
2310   for(Int_t itrack=0; itrack < esdTrack->GetEntriesFast(); itrack++){
2311     double p[3]; 
2312     //esdTrack[itrack]->GetConstrainedPxPyPz(p);
2313     ((AliESDtrack*)(esdTrack->At(itrack)))->GetConstrainedPxPyPz(p);
2314     TLorentzVector currentTrack;
2315     currentTrack.SetXYZM(p[0],p[1],p[2],fElectronMass);
2316     new((array)[array.GetEntriesFast()])  TLorentzVector(currentTrack);
2317     //    tlVtrack.push_back(currentTrack);
2318   }
2319         
2320   return array;
2321         
2322   //  return tlVtrack;
2323 }
2324