Jet and Particle identification tasks moved to different directories
[u/mrichter/AliRoot.git] / PWG4 / PartCorr / 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.0                                                            *
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 // root
17 #include <TROOT.h>
18 #include <TSystem.h>
19 #include <TInterpreter.h>
20 #include <TChain.h>
21 #include <TFile.h>
22 #include <Riostream.h>
23
24 // analysis
25 //#include "AliAnalysisTaskSE.h"
26 #include "AliAnalysisTaskGammaConversion.h"
27 #include "AliAnalysisManager.h"
28 #include "AliESDInputHandler.h"
29 #include "AliMCEventHandler.h"
30 #include "AliMCEvent.h"
31 #include "AliESDEvent.h"
32 #include "AliAODEvent.h"
33 #include "AliAODHandler.h"
34 #include "AliStack.h"
35 #include "AliLog.h"
36 //#include "TLorentzVector.h"
37 #include "AliKFVertex.h"
38
39 ClassImp(AliAnalysisTaskGammaConversion)
40
41
42 AliAnalysisTaskGammaConversion::AliAnalysisTaskGammaConversion():
43   AliAnalysisTaskSE(),
44   fV0Reader(NULL),
45   fStack(NULL),
46   fOutputContainer(NULL),
47   fHistograms(NULL),
48   fDoMCTruth(kFALSE),
49   fMCAllGammas(),
50   fMCPi0s(),
51   fMCEtas(),
52   fMCGammaChi_c(),
53   fKFReconstructedGammas(),
54   fElectronMass(-1),
55   fGammaMass(-1),
56   fPi0Mass(-1),
57   fEtaMass(-1),
58   fGammaWidth(-1),
59   fPi0Width(-1),
60   fEtaWidth(-1),
61   fCalculateBackground(kFALSE)
62 {
63   // Default constructor
64   // Common I/O in slot 0
65   DefineInput (0, TChain::Class());
66   DefineOutput(0, TTree::Class());
67
68   // Your private output
69   DefineOutput(1, TList::Class());
70 }
71
72 AliAnalysisTaskGammaConversion::AliAnalysisTaskGammaConversion(const char* name):
73   AliAnalysisTaskSE(name),
74   fV0Reader(NULL),
75   fStack(NULL),
76   fOutputContainer(0x0),
77   fHistograms(NULL),
78   fDoMCTruth(kFALSE),
79   fMCAllGammas(),
80   fMCPi0s(),
81   fMCEtas(),
82   fMCGammaChi_c(),
83   fKFReconstructedGammas(),
84   fElectronMass(-1),
85   fGammaMass(-1),
86   fPi0Mass(-1),
87   fEtaMass(-1),
88   fGammaWidth(-1),
89   fPi0Width(-1),
90   fEtaWidth(-1),
91   fCalculateBackground(kFALSE)
92 {
93   // Common I/O in slot 0
94   DefineInput (0, TChain::Class());
95   DefineOutput(0, TTree::Class());
96   
97   // Your private output
98   DefineOutput(1, TList::Class());
99 }
100
101 AliAnalysisTaskGammaConversion::~AliAnalysisTaskGammaConversion() 
102 {
103   // Remove all pointers
104  
105   if(fOutputContainer){
106     fOutputContainer->Clear() ; 
107     delete fOutputContainer ;
108   }
109   if(fHistograms){
110     delete fHistograms;
111   }
112   if(fV0Reader){
113     delete fV0Reader;
114   }
115 }
116
117
118 void AliAnalysisTaskGammaConversion::Init()
119 {
120   // Initialization
121 }
122
123
124 void AliAnalysisTaskGammaConversion::Exec(Option_t */*option*/)
125 {
126   // Execute analysis for current event
127   //
128   
129   ConnectInputData("");
130   
131   //clear vectors
132   fMCAllGammas.clear();
133   fMCPi0s.clear();
134   fMCEtas.clear();
135   fMCGammaChi_c.clear();
136   /*
137   for(UInt_t ikf=0;ikf<fKFReconstructedGammas.size();ikf++){
138     delete fKFReconstructedGammas[ikf];
139     fKFReconstructedGammas[ikf]=NULL;
140   }
141   */
142   fKFReconstructedGammas.clear();
143
144   //Clear the data in the v0Reader
145   fV0Reader->UpdateEventByEventData();
146   // Process the v0 information
147   if(fDoMCTruth){
148     ProcessMCData();
149   }
150
151   // Process the v0 information
152   ProcessV0s();
153
154   if(fCalculateBackground){//calculate background if flag is set
155     CalculateBackground();
156   }
157
158   // Process reconstructed gammas
159   ProcessGammasForNeutralMesonAnalysis();
160
161   PostData(1, fOutputContainer);
162   
163 }
164
165 void AliAnalysisTaskGammaConversion::ConnectInputData(Option_t */*option*/){
166
167   if(fV0Reader == NULL){
168     // Write warning here cuts and so on are default if this ever happens
169   }
170   fV0Reader->Initialize();
171 }
172
173 void AliAnalysisTaskGammaConversion::ProcessMCData(){
174   
175   fStack = fV0Reader->GetMCStack();
176
177   for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++) {
178     TParticle* particle = (TParticle *)fStack->Particle(iTracks);
179
180     if (!particle) {
181       //print warning here
182       continue;
183     }
184     
185     if(particle->Pt()<fV0Reader->GetPtCut()){
186       continue;
187     }
188     
189     if(TMath::Abs(particle->Eta())> fV0Reader->GetEtaCut()){
190       continue;
191     }
192
193     if(particle->R()>fV0Reader->GetMaxRCut()){ // cuts on distance from collision point
194       continue;
195     }
196
197     Double_t tmpPhi=particle->Phi();
198     if(particle->Phi()> TMath::Pi()){
199       tmpPhi = particle->Phi()-(2*TMath::Pi());
200     }
201
202     
203     //process the gammas
204     if (particle->GetPdgCode()== 22){
205       fMCAllGammas.push_back(particle);
206       if(particle->GetMother(0)>-1){ //Means we have a mother
207         if( fStack->Particle(particle->GetMother(0))->GetPdgCode() != 22 ){//Checks for a non gamma mother.
208           if(fHistograms->fMC_Gamma_Energy){fHistograms->fMC_Gamma_Energy->Fill(particle->Energy());}
209           if(fHistograms->fMC_Gamma_Pt){fHistograms->fMC_Gamma_Pt->Fill(particle->Pt());}
210           if(fHistograms->fMC_Gamma_Eta){fHistograms->fMC_Gamma_Eta->Fill(particle->Eta());}
211           
212           /*      Double_t tmpPhi=particle->Phi();
213           if(particle->Phi()> TMath::Pi()){
214             tmpPhi = particle->Phi()-(2*TMath::Pi());
215             }*/
216           if(fHistograms->fMC_Gamma_Phi){fHistograms->fMC_Gamma_Phi->Fill(tmpPhi);}
217
218           //adding the conversion points from all gammas with e+e- daughters
219           if(particle->GetNDaughters() == 2){
220             TParticle* daughter0 = (TParticle*)fStack->Particle(particle->GetFirstDaughter());
221             TParticle* daughter1 = (TParticle*)fStack->Particle(particle->GetLastDaughter());
222             
223             if(daughter0->R()>fV0Reader->GetMaxRCut() || daughter1->R()>fV0Reader->GetMaxRCut()){
224               continue;
225             }
226             
227             if(daughter0->GetPdgCode() == -11 && daughter1->GetPdgCode() == 11 ||
228                daughter0->GetPdgCode() == 11 && daughter1->GetPdgCode() == -11){
229
230               // begin Mapping 
231               Int_t rBin    = fHistograms->GetRBin(daughter0->R());
232               Int_t phiBin  = fHistograms->GetPhiBin(daughter0->Phi());
233               
234               if(fHistograms->fMC_Mapping[phiBin][rBin] != NULL){fHistograms->fMC_Mapping[phiBin][rBin]->Fill(daughter0->Vz(), particle->Eta());}
235               if(fHistograms->fMC_Mapping_Phi[phiBin] != NULL){fHistograms->fMC_Mapping_Phi[phiBin]->Fill(daughter0->Vz(), particle->Eta());}
236               if(fHistograms->fMC_Mapping_R[rBin] != NULL){fHistograms->fMC_Mapping_R[rBin]->Fill(daughter0->Vz(), particle->Eta());}
237               //end mapping
238
239
240               if(fHistograms->fMC_EP_R){fHistograms->fMC_EP_R->Fill(daughter0->R());}
241               if(fHistograms->fMC_EP_Z_R){fHistograms->fMC_EP_Z_R->Fill(daughter0->Vz(),daughter0->R());}
242               if(fHistograms->fMC_EP_X_Y){fHistograms->fMC_EP_X_Y->Fill(daughter0->Vx(),daughter0->Vy());}
243               if(fHistograms->fMC_EP_OpeningAngle){fHistograms->fMC_EP_OpeningAngle->Fill(GetMCOpeningAngle(daughter0, daughter1));}
244
245             }
246           }
247         }
248         if( fStack->Particle(particle->GetMother(0))->GetPdgCode()==10441 ||//chi_c0 
249             fStack->Particle(particle->GetMother(0))->GetPdgCode()==20443 ||//psi_2S
250             fStack->Particle(particle->GetMother(0))->GetPdgCode()==445  //chi_c2
251         ){ 
252           fMCGammaChi_c.push_back(particle);
253          }
254       }
255       else{//means we have a primary particle
256         if(fHistograms->fMC_DirectGamma_Energy){fHistograms->fMC_DirectGamma_Energy->Fill(particle->Energy());}
257         if(fHistograms->fMC_DirectGamma_Pt){fHistograms->fMC_DirectGamma_Pt->Fill(particle->Pt());}
258         if(fHistograms->fMC_DirectGamma_Eta){fHistograms->fMC_DirectGamma_Eta->Fill(particle->Eta());}
259         /*
260         Double_t tmpPhi=particle->Phi();
261         if(particle->Phi()> TMath::Pi()){
262           tmpPhi = particle->Phi()-(2*TMath::Pi());
263         }
264         */
265         if(fHistograms->fMC_DirectGamma_Phi){fHistograms->fMC_DirectGamma_Phi->Fill(tmpPhi);}
266
267         //adding the conversion points from all gammas with e+e- daughters
268         if(particle->GetNDaughters() == 2){
269           TParticle* daughter0 = (TParticle*)fStack->Particle(particle->GetFirstDaughter());
270           TParticle* daughter1 = (TParticle*)fStack->Particle(particle->GetLastDaughter());
271           if(daughter0->GetPdgCode() == -11 && daughter1->GetPdgCode() == 11 ||
272              daughter0->GetPdgCode() == 11 && daughter1->GetPdgCode() == -11){
273             
274             if(fHistograms->fMC_EP_R){fHistograms->fMC_EP_R->Fill(daughter0->R());}
275             if(fHistograms->fMC_EP_Z_R){fHistograms->fMC_EP_Z_R->Fill(daughter0->Vz(),daughter0->R());}
276             if(fHistograms->fMC_EP_X_Y){fHistograms->fMC_EP_X_Y->Fill(daughter0->Vx(),daughter0->Vy());}
277             if(fHistograms->fMC_EP_OpeningAngle){fHistograms->fMC_EP_OpeningAngle->Fill(GetMCOpeningAngle(daughter0,daughter1));}
278
279           }
280         }
281
282       }
283     }
284     else if (TMath::Abs(particle->GetPdgCode())== 11){ // Means we have an electron or a positron
285       if(particle->GetMother(0)>-1){ // means we have a mother
286         if( fStack->Particle(particle->GetMother(0))->GetPdgCode()==22 ){ // Means we have a gamma mother
287           if(particle->GetPdgCode() == 11){//electron 
288             if(fHistograms->fMC_E_Energy){fHistograms->fMC_E_Energy->Fill(particle->Energy());}
289             if(fHistograms->fMC_E_Pt){fHistograms->fMC_E_Pt->Fill(particle->Pt());}
290             if(fHistograms->fMC_E_Eta){fHistograms->fMC_E_Eta->Fill(particle->Eta());}
291             /*      Double_t tmpPhi=particle->Phi();
292             if(particle->Phi()> TMath::Pi()){
293               tmpPhi = particle->Phi()-(2*TMath::Pi());
294               }*/
295             if(fHistograms->fMC_E_Phi){fHistograms->fMC_E_Phi->Fill(tmpPhi);}
296           }
297           if(particle->GetPdgCode() == -11){//positron 
298             if(fHistograms->fMC_P_Energy){fHistograms->fMC_P_Energy->Fill(particle->Energy());}
299             if(fHistograms->fMC_P_Pt){fHistograms->fMC_P_Pt->Fill(particle->Pt());}
300             if(fHistograms->fMC_P_Eta){fHistograms->fMC_P_Eta->Fill(particle->Eta());}
301             /*
302             Double_t tmpPhi=particle->Phi();
303             if(particle->Phi()> TMath::Pi()){
304               tmpPhi = particle->Phi()-(2*TMath::Pi());
305             }
306             */
307             if(fHistograms->fMC_P_Phi){fHistograms->fMC_P_Phi->Fill(tmpPhi);}
308           }
309         }
310       }
311     }
312     else if(particle->GetNDaughters() == 2){
313
314       TParticle* daughter0 = (TParticle*)fStack->Particle(particle->GetFirstDaughter());
315       TParticle* daughter1 = (TParticle*)fStack->Particle(particle->GetLastDaughter());
316       if(daughter0->GetPdgCode() == 22 && daughter1->GetPdgCode() == 22){//check for gamma gamma daughters
317         
318         if(particle->GetPdgCode()==111){//Pi0
319           if( iTracks >= fStack->GetNprimary()){
320             
321             if(fHistograms->fMC_Pi0Secondaries_Eta){fHistograms->fMC_Pi0Secondaries_Eta->Fill(particle->Eta());}
322             /*
323               Double_t tmpPhi=particle->Phi();
324               if(particle->Phi()> TMath::Pi()){
325               tmpPhi = particle->Phi()-(2*TMath::Pi());
326               }
327             */
328             if(fHistograms->fMC_Pi0Secondaries_Phi){fHistograms->fMC_Pi0Secondaries_Phi->Fill(tmpPhi);}
329             if(fHistograms->fMC_Pi0Secondaries_Pt){fHistograms->fMC_Pi0Secondaries_Pt->Fill(particle->Pt());}
330             if(fHistograms->fMC_Pi0Secondaries_Energy){fHistograms->fMC_Pi0Secondaries_Energy->Fill(particle->Energy());}
331             if(fHistograms->fMC_Pi0Secondaries_R){fHistograms->fMC_Pi0Secondaries_R->Fill(particle->R());}
332             if(fHistograms->fMC_Pi0Secondaries_Z_R){fHistograms->fMC_Pi0Secondaries_Z_R->Fill(particle->Vz(),particle->R());}
333             if(fHistograms->fMC_Pi0Secondaries_OpeningAngleGamma){fHistograms->fMC_Pi0Secondaries_OpeningAngleGamma->Fill(GetMCOpeningAngle(daughter0,daughter1));}
334             if(fHistograms->fMC_Pi0Secondaries_X_Y){fHistograms->fMC_Pi0Secondaries_X_Y->Fill(particle->Vx(),particle->Vy());}//only fill from one daughter to avoid multiple filling
335           }
336           else{
337             if(fHistograms->fMC_Pi0_Eta){fHistograms->fMC_Pi0_Eta->Fill(particle->Eta());}
338             /*
339             Double_t tmpPhi=particle->Phi();
340             if(particle->Phi()> TMath::Pi()){
341               tmpPhi = particle->Phi()-(2*TMath::Pi());
342             }
343             */
344             if(fHistograms->fMC_Pi0_Phi){fHistograms->fMC_Pi0_Phi->Fill(tmpPhi);}
345             if(fHistograms->fMC_Pi0_Pt){fHistograms->fMC_Pi0_Pt->Fill(particle->Pt());}
346             if(fHistograms->fMC_Pi0_Energy){fHistograms->fMC_Pi0_Energy->Fill(particle->Energy());}
347             if(fHistograms->fMC_Pi0_R){fHistograms->fMC_Pi0_R->Fill(particle->R());}
348             if(fHistograms->fMC_Pi0_Z_R){fHistograms->fMC_Pi0_Z_R->Fill(particle->Vz(),particle->R());}
349             if(fHistograms->fMC_Pi0_OpeningAngleGamma){fHistograms->fMC_Pi0_OpeningAngleGamma->Fill(GetMCOpeningAngle(daughter0,daughter1));}
350             if(fHistograms->fMC_Pi0_X_Y){fHistograms->fMC_Pi0_X_Y->Fill(particle->Vx(),particle->Vy());}//only fill from one daughter to avoid multiple filling
351           }
352         }
353         else if(particle->GetPdgCode()==221){//Eta
354           if(fHistograms->fMC_Eta_Eta){fHistograms->fMC_Eta_Eta->Fill(particle->Eta());}
355           /*
356             Double_t tmpPhi=particle->Phi();
357           if(particle->Phi()> TMath::Pi()){
358             tmpPhi = particle->Phi()-(2*TMath::Pi());
359           }
360           */
361           if(fHistograms->fMC_Eta_Phi){fHistograms->fMC_Eta_Phi->Fill(tmpPhi);}
362           if(fHistograms->fMC_Eta_Pt){fHistograms->fMC_Eta_Pt->Fill(particle->Pt());}
363           if(fHistograms->fMC_Eta_Energy){fHistograms->fMC_Eta_Energy->Fill(particle->Energy());}
364           if(fHistograms->fMC_Eta_R){fHistograms->fMC_Eta_R->Fill(particle->R());}
365           if(fHistograms->fMC_Eta_Z_R){fHistograms->fMC_Eta_Z_R->Fill(particle->Vz(),particle->R());}
366           if(fHistograms->fMC_Eta_OpeningAngleGamma){fHistograms->fMC_Eta_OpeningAngleGamma->Fill(GetMCOpeningAngle(daughter0,daughter1));}
367           if(fHistograms->fMC_Eta_X_Y){fHistograms->fMC_Eta_X_Y->Fill(particle->Vx(),particle->Vy());}//only fill from one daughter to avoid multiple filling
368         }
369         
370         //the match data should be filled no matter which mother the gamma-gamma comes from
371         if(fHistograms->fMC_Match_Gamma_R){fHistograms->fMC_Match_Gamma_R->Fill(particle->R());}
372         if(fHistograms->fMC_Match_Gamma_Z_R){fHistograms->fMC_Match_Gamma_Z_R->Fill(particle->Vz(),particle->R());}
373         if(fHistograms->fMC_Match_Gamma_X_Y){fHistograms->fMC_Match_Gamma_X_Y->Fill(particle->Vx(),particle->Vy());}
374         if(fHistograms->fMC_Match_Gamma_Mass){fHistograms->fMC_Match_Gamma_Mass->Fill(particle->GetCalcMass());}
375         if(fHistograms->fMC_Match_Gamma_OpeningAngle){fHistograms->fMC_Match_Gamma_OpeningAngle->Fill(GetMCOpeningAngle(daughter0,daughter1));}
376         if(fHistograms->fMC_Match_Gamma_Energy){fHistograms->fMC_Match_Gamma_Energy->Fill(particle->Energy());}
377         if(fHistograms->fMC_Match_Gamma_Pt){fHistograms->fMC_Match_Gamma_Pt->Fill(particle->Pt());}
378         if(fHistograms->fMC_Match_Gamma_Eta){fHistograms->fMC_Match_Gamma_Eta->Fill(particle->Eta());}
379         /*
380         Double_t tmpPhi=particle->Phi();
381         if(particle->Phi()> TMath::Pi()){
382           tmpPhi = particle->Phi()-(2*TMath::Pi());
383         }
384         */
385         if(fHistograms->fMC_Match_Gamma_Phi){fHistograms->fMC_Match_Gamma_Phi->Fill(tmpPhi);}
386       }
387     }
388   }
389 }
390
391 void AliAnalysisTaskGammaConversion::ProcessV0s(){
392   Int_t nSurvivingV0s=0;
393   while(fV0Reader->NextV0()){
394     nSurvivingV0s++;
395     //-------------------------- filling v0 information -------------------------------------
396     if(fHistograms->fESD_EP_OpeningAngle){fHistograms->fESD_EP_OpeningAngle->Fill(fV0Reader->GetOpeningAngle());}    
397     if(fHistograms->fESD_EP_R){fHistograms->fESD_EP_R->Fill(fV0Reader->GetXYRadius());}
398     if(fHistograms->fESD_EP_Z_R){fHistograms->fESD_EP_Z_R->Fill(fV0Reader->GetZ(),fV0Reader->GetXYRadius());}
399     if(fHistograms->fESD_EP_X_Y){fHistograms->fESD_EP_X_Y->Fill(fV0Reader->GetX(),fV0Reader->GetY());}
400     
401     
402     if(fHistograms->fESD_E_Energy){fHistograms->fESD_E_Energy->Fill(fV0Reader->GetNegativeTrackEnergy());}
403     if(fHistograms->fESD_E_Pt){fHistograms->fESD_E_Pt->Fill(fV0Reader->GetNegativeTrackPt());}
404     if(fHistograms->fESD_E_Eta){fHistograms->fESD_E_Eta->Fill(fV0Reader->GetNegativeTrackEta());}
405     if(fHistograms->fESD_E_Phi){fHistograms->fESD_E_Phi->Fill(fV0Reader->GetNegativeTrackPhi());}
406     
407     if(fHistograms->fESD_P_Energy){fHistograms->fESD_P_Energy->Fill(fV0Reader->GetPositiveTrackEnergy());}
408     if(fHistograms->fESD_P_Pt){fHistograms->fESD_P_Pt->Fill(fV0Reader->GetPositiveTrackPt());}
409     if(fHistograms->fESD_P_Eta){fHistograms->fESD_P_Eta->Fill(fV0Reader->GetPositiveTrackEta());}
410     if(fHistograms->fESD_P_Phi){fHistograms->fESD_P_Phi->Fill(fV0Reader->GetPositiveTrackPhi());}
411     
412     if(fHistograms->fESD_Gamma_Energy){fHistograms->fESD_Gamma_Energy->Fill(fV0Reader->GetMotherCandidateEnergy());}
413     if(fHistograms->fESD_Gamma_Pt){fHistograms->fESD_Gamma_Pt->Fill(fV0Reader->GetMotherCandidatePt());}
414     if(fHistograms->fESD_Gamma_Eta){fHistograms->fESD_Gamma_Eta->Fill(fV0Reader->GetMotherCandidateEta());}
415     if(fHistograms->fESD_Gamma_Phi){fHistograms->fESD_Gamma_Phi->Fill(fV0Reader->GetMotherCandidatePhi());}
416
417
418     // begin mapping
419     Int_t rBin    = fHistograms->GetRBin(fV0Reader->GetXYRadius());
420     Int_t phiBin  = fHistograms->GetPhiBin(fV0Reader->GetNegativeTrackPhi());
421     Double_t motherCandidateEta= fV0Reader->GetMotherCandidateEta();
422     if(fHistograms->fESD_Mapping[phiBin][rBin] != NULL){fHistograms->fESD_Mapping[phiBin][rBin]->Fill(fV0Reader->GetZ(),motherCandidateEta);}
423     if(fHistograms->fESD_Mapping_Phi[phiBin] != NULL){fHistograms->fESD_Mapping_Phi[phiBin]->Fill(fV0Reader->GetZ(),motherCandidateEta);}
424     if(fHistograms->fESD_Mapping_R[rBin] != NULL){fHistograms->fESD_Mapping_R[rBin]->Fill(fV0Reader->GetZ(),motherCandidateEta);}
425     // end mapping
426     
427
428     fKFReconstructedGammas.push_back(*fV0Reader->GetMotherCandidateKFCombination());
429
430     //----------------------------------- checking for "real" conversions (MC match) --------------------------------------
431     if(fDoMCTruth){
432       if(fV0Reader->HasSameMCMother() == kFALSE){
433         continue;
434       }
435
436       TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
437       TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
438       if(negativeMC->GetPdgCode()!=11 || positiveMC->GetPdgCode()!=-11){
439         continue;
440       }
441
442       if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
443         if(fHistograms->fESD_Match_Gamma_X_Y){fHistograms->fESD_Match_Gamma_X_Y->Fill(fV0Reader->GetX(),fV0Reader->GetY());}
444         if(fHistograms->fESD_Match_Gamma_OpeningAngle){fHistograms->fESD_Match_Gamma_OpeningAngle->Fill(fV0Reader->GetOpeningAngle());}
445         if(fHistograms->fESD_Match_Gamma_Pt){fHistograms->fESD_Match_Gamma_Pt->Fill(fV0Reader->GetMotherCandidatePt());}
446         if(fHistograms->fESD_Match_Gamma_Energy){fHistograms->fESD_Match_Gamma_Energy->Fill(fV0Reader->GetMotherCandidateEnergy());}
447         if(fHistograms->fESD_Match_Gamma_Eta){fHistograms->fESD_Match_Gamma_Eta->Fill(fV0Reader->GetMotherCandidateEta());}
448
449         if(fHistograms->fESD_Match_Gamma_Phi){fHistograms->fESD_Match_Gamma_Phi->Fill(fV0Reader->GetMotherCandidatePhi());}
450         if(fHistograms->fESD_Match_Gamma_Mass){fHistograms->fESD_Match_Gamma_Mass->Fill(fV0Reader->GetMotherCandidateMass());}
451         if(fHistograms->fESD_Match_Gamma_Width){fHistograms->fESD_Match_Gamma_Width->Fill(fV0Reader->GetMotherCandidateWidth());}
452         if(fHistograms->fESD_Match_Gamma_Chi2){fHistograms->fESD_Match_Gamma_Chi2->Fill(fV0Reader->GetMotherCandidateChi2());}
453         if(fHistograms->fESD_Match_Gamma_NDF){fHistograms->fESD_Match_Gamma_NDF->Fill(fV0Reader->GetMotherCandidateNDF());}
454         if(fHistograms->fESD_Match_Gamma_R){fHistograms->fESD_Match_Gamma_R->Fill(fV0Reader->GetXYRadius());}
455         if(fHistograms->fESD_Match_Gamma_Z_R){fHistograms->fESD_Match_Gamma_Z_R->Fill(fV0Reader->GetZ(),fV0Reader->GetXYRadius());}
456
457         //resolution
458         Double_t mc_pt   = fV0Reader->GetMotherMCParticle()->Pt();
459         Double_t esd_pt  = fV0Reader->GetMotherCandidatePt();
460         Double_t res_dPt = ((esd_pt - mc_pt)/mc_pt)*100;
461         if(fHistograms->fResolution_dPt != NULL){fHistograms->fResolution_dPt->Fill(mc_pt,res_dPt);}
462         if(fHistograms->fResolution_MC_Pt != NULL){fHistograms->fResolution_MC_Pt->Fill(mc_pt);}
463         if(fHistograms->fResolution_ESD_Pt != NULL){fHistograms->fResolution_ESD_Pt->Fill(esd_pt);}
464         
465         Double_t res_dZ = ((fV0Reader->GetZ() -fV0Reader->GetNegativeMCParticle()->Vz())/fV0Reader->GetNegativeMCParticle()->Vz())*100;
466         if(fHistograms->fResolution_dZ != NULL){fHistograms->fResolution_dZ->Fill(fV0Reader->GetNegativeMCParticle()->Vz(),res_dZ);}
467         if(fHistograms->fResolution_MC_Z != NULL){fHistograms->fResolution_MC_Z->Fill(fV0Reader->GetNegativeMCParticle()->Vz());}
468         if(fHistograms->fResolution_ESD_Z != NULL){fHistograms->fResolution_ESD_Z->Fill(fV0Reader->GetZ());}
469         
470         Double_t res_dR = ((fV0Reader->GetXYRadius() - fV0Reader->GetNegativeMCParticle()->R())/fV0Reader->GetNegativeMCParticle()->R())*100;
471         if(fHistograms->fResolution_dR != NULL){fHistograms->fResolution_dR->Fill(fV0Reader->GetNegativeMCParticle()->R(),res_dR);}
472         if(fHistograms->fResolution_MC_R != NULL){fHistograms->fResolution_MC_R->Fill(fV0Reader->GetNegativeMCParticle()->R());}
473         if(fHistograms->fResolution_ESD_R != NULL){fHistograms->fResolution_ESD_R->Fill(fV0Reader->GetXYRadius());}
474         if(fHistograms->fResolution_dR_dPt != NULL){fHistograms->fResolution_dR_dPt->Fill(res_dR,res_dPt);}
475       }
476     }
477   }
478   if(fHistograms->fNumberOfSurvivingV0s != NULL){fHistograms->fNumberOfSurvivingV0s->Fill(nSurvivingV0s);}
479   if(fHistograms->fNumberOfV0s != NULL){fHistograms->fNumberOfV0s->Fill(fV0Reader->GetNumberOfV0s());}
480 }
481
482 void AliAnalysisTaskGammaConversion::ProcessGammasForNeutralMesonAnalysis(){
483
484   for(UInt_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammas.size();firstGammaIndex++){
485     for(UInt_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fKFReconstructedGammas.size();secondGammaIndex++){
486       AliKFParticle * twoGammaDecayCandidateDaughter0 = &fKFReconstructedGammas[firstGammaIndex];
487       AliKFParticle * twoGammaDecayCandidateDaughter1 = &fKFReconstructedGammas[secondGammaIndex];
488
489       // Pi0's
490       AliKFParticle *pi0Candidate = new AliKFParticle(*twoGammaDecayCandidateDaughter0,*twoGammaDecayCandidateDaughter1);
491
492       //      pi0Candidate->SetMassConstraint(fPi0Mass,fPi0Width);
493
494       Double_t massPi0 =0.;
495       Double_t widthPi0 = 0.;
496       Double_t chi2Pi0 =10000.; 
497       pi0Candidate->GetMass(massPi0,widthPi0);
498       if(pi0Candidate->GetNDF()>0){
499         chi2Pi0 = pi0Candidate->GetChi2()/pi0Candidate->GetNDF();
500         //      if(chi2Pi0>0 && chi2Pi0<fChi2Cut){//TODO  find this out
501         if(chi2Pi0>0 && chi2Pi0<10000){//TODO  find this out see line above
502
503           /*
504           TVector3 vertexDaughter0(twoGammaDecayCandidateDaughter0->Px(),twoGammaDecayCandidateDaughter0->Py(),twoGammaDecayCandidateDaughter0->Pz());
505           TVector3 vertexDaughter1(twoGammaDecayCandidateDaughter1->Px(),twoGammaDecayCandidateDaughter1->Py(),twoGammaDecayCandidateDaughter1->Pz());
506           TVector3 vertexPi0Candidate(pi0Candidate->Px(),pi0Candidate->Py(),pi0Candidate->Pz());
507           
508           Double_t openingAnglePi0 = vertexDaughter0.Angle(vertexDaughter1);
509           
510           Double_t radiusPi0 = sqrt( vertexPi0Candidate.x()*vertexPi0Candidate.x() + vertexPi0Candidate.y()*vertexPi0Candidate.y() );
511           */
512
513           TVector3 vectorPi0Candidate(pi0Candidate->Px(),pi0Candidate->Py(),pi0Candidate->Pz());
514
515           Double_t openingAnglePi0 = twoGammaDecayCandidateDaughter0->GetAngle(*twoGammaDecayCandidateDaughter1);
516
517           //Double_t vtx000[3] = {0,0,0};
518           //   NOT SURE IF THIS DOES WHAT WE WANT: remember to ask Sergey  Double_t radiusPi0 = pi0Candidate->GetDistanceFromVertex(vtx000);
519           //Calculating by hand the radius
520           Double_t tmpX= pi0Candidate->GetX();
521           Double_t tmpY= pi0Candidate->GetY();
522           
523           Double_t radiusPi0 = TMath::Sqrt(tmpX*tmpX + tmpY*tmpY);
524
525           if(fHistograms->fESD_Pi0_OpeningAngleGamma){fHistograms->fESD_Pi0_OpeningAngleGamma->Fill(openingAnglePi0);}
526           if(fHistograms->fESD_Pi0_Energy){fHistograms->fESD_Pi0_Energy->Fill(pi0Candidate->GetE());}
527           if(fHistograms->fESD_Pi0_Pt){fHistograms->fESD_Pi0_Pt->Fill(sqrt(pi0Candidate->GetPx()*pi0Candidate->GetPx()+pi0Candidate->GetPy()*pi0Candidate->GetPy()));}
528           if(fHistograms->fESD_Pi0_Eta){fHistograms->fESD_Pi0_Eta->Fill(vectorPi0Candidate.Eta());}
529           if(fHistograms->fESD_Pi0_Phi){fHistograms->fESD_Pi0_Phi->Fill(vectorPi0Candidate.Phi());}
530           if(fHistograms->fESD_Pi0_Mass){fHistograms->fESD_Pi0_Mass->Fill(massPi0);}
531           if(fHistograms->fESD_Pi0_R){fHistograms->fESD_Pi0_R->Fill(radiusPi0);}
532           if(fHistograms->fESD_Pi0_Z_R){fHistograms->fESD_Pi0_Z_R->Fill(tmpY,radiusPi0);}
533           if(fHistograms->fESD_Pi0_X_Y){fHistograms->fESD_Pi0_X_Y->Fill(tmpX,tmpY);}      
534         }
535       }
536       delete pi0Candidate;
537
538       //Eta's
539       AliKFParticle *etaCandidate = new AliKFParticle(*twoGammaDecayCandidateDaughter0,*twoGammaDecayCandidateDaughter1);
540
541       //      etaCandidate->SetMassConstraint(fEtaMass,fEtaWidth);
542       Double_t massEta =0.;
543       Double_t widthEta = 0.;
544       Double_t chi2Eta =10000.; 
545       etaCandidate->GetMass(massEta,widthEta);
546       if(etaCandidate->GetNDF()>0){
547         chi2Eta = etaCandidate->GetChi2()/etaCandidate->GetNDF();
548         //      if(chi2Eta>0 && chi2Eta<fChi2Cut){
549         if(chi2Eta>0 && chi2Eta<100000){// TODO find this out se line above
550           /*
551           TVector3 vertexDaughter0(twoGammaDecayCandidateDaughter0->Px(),twoGammaDecayCandidateDaughter0->Py(),twoGammaDecayCandidateDaughter0->Pz());
552           TVector3 vertexDaughter1(twoGammaDecayCandidateDaughter1->Px(),twoGammaDecayCandidateDaughter1->Py(),twoGammaDecayCandidateDaughter1->Pz());
553           TVector3 vertexEtaCandidate(etaCandidate->Px(),etaCandidate->Py(),etaCandidate->Pz());
554           Double_t openingAngleEta = vertexDaughter0.Angle(vertexDaughter1);
555           
556           Double_t radiusEta = sqrt( etaCandidate->GetX()*etaCandidate->GetX() + vertexEtaCandidate.y()*vertexEtaCandidate.y() );
557           */
558
559           //Calculating by hand the radius
560           Double_t tmpX= etaCandidate->GetX();
561           Double_t tmpY= etaCandidate->GetY();
562           
563           Double_t radiusEta = TMath::Sqrt(tmpX*tmpX + tmpY*tmpY);
564           
565           Double_t openingAngleEta = twoGammaDecayCandidateDaughter0->GetAngle(*twoGammaDecayCandidateDaughter1);
566
567           TVector3 vectorEtaCandidate(etaCandidate->Px(),etaCandidate->Py(),etaCandidate->Pz());
568
569           if(fHistograms->fESD_Eta_OpeningAngleGamma){fHistograms->fESD_Eta_OpeningAngleGamma->Fill(openingAngleEta);}
570           if(fHistograms->fESD_Eta_Energy){fHistograms->fESD_Eta_Energy->Fill(etaCandidate->GetE());}
571           if(fHistograms->fESD_Eta_Pt){fHistograms->fESD_Eta_Pt->Fill(sqrt(etaCandidate->GetPx()*etaCandidate->GetPx()+etaCandidate->GetPy()*etaCandidate->GetPy()));}
572           if(fHistograms->fESD_Eta_Eta){fHistograms->fESD_Eta_Eta->Fill(vectorEtaCandidate.Eta());}
573           if(fHistograms->fESD_Eta_Phi){fHistograms->fESD_Eta_Phi->Fill(vectorEtaCandidate.Phi());}
574           if(fHistograms->fESD_Eta_Mass){fHistograms->fESD_Eta_Mass->Fill(massEta);}
575           if(fHistograms->fESD_Eta_R){fHistograms->fESD_Eta_R->Fill(radiusEta);}
576           if(fHistograms->fESD_Eta_Z_R){fHistograms->fESD_Eta_Z_R->Fill(etaCandidate->GetZ(),radiusEta);}
577           if(fHistograms->fESD_Eta_X_Y){fHistograms->fESD_Eta_X_Y->Fill(tmpX,tmpY);}      
578         }
579       }
580       delete etaCandidate;
581     }
582   }
583
584 }
585
586 void AliAnalysisTaskGammaConversion::CalculateBackground(){
587
588   for(UInt_t iCurrent=0;iCurrent<fV0Reader->fCurrentEventGoodV0s.size();iCurrent++){
589     AliKFParticle * currentEventGoodV0 = &fV0Reader->fCurrentEventGoodV0s.at(iCurrent);
590     for(UInt_t iPrevious=0;iPrevious<fV0Reader->fCurrentEventGoodV0s.size();iPrevious++){
591       AliKFParticle * previousEventGoodV0 = &fV0Reader->fPreviousEventGoodV0s.at(iPrevious);
592
593       AliKFParticle *backgroundCandidate = new AliKFParticle(*currentEventGoodV0,*previousEventGoodV0);
594
595       Double_t massBG =0.;
596       Double_t widthBG = 0.;
597       Double_t chi2BG =10000.;  
598       backgroundCandidate->GetMass(massBG,widthBG);
599       if(backgroundCandidate->GetNDF()>0){
600         chi2BG = backgroundCandidate->GetChi2()/backgroundCandidate->GetNDF();
601         //      if(chi2Pi0>0 && chi2Pi0<fChi2Cut){//TODO  find this out
602         if(chi2BG>0 && chi2BG<fV0Reader->GetChi2Cut()){//TODO  find this out see line above
603
604           TVector3 vectorBGCandidate(backgroundCandidate->Px(),backgroundCandidate->Py(),backgroundCandidate->Pz());
605
606           Double_t openingAngleBG = currentEventGoodV0->GetAngle(*previousEventGoodV0);
607
608           //Calculating by hand the radius (find a better way)
609           Double_t tmpX= backgroundCandidate->GetX();
610           Double_t tmpY= backgroundCandidate->GetY();
611           
612           Double_t radiusBG = TMath::Sqrt(tmpX*tmpX + tmpY*tmpY);
613
614           if(fHistograms->fESD_Background_OpeningAngleGamma){fHistograms->fESD_Background_OpeningAngleGamma->Fill(openingAngleBG);}
615           if(fHistograms->fESD_Background_Energy){fHistograms->fESD_Background_Energy->Fill(backgroundCandidate->GetE());}
616           if(fHistograms->fESD_Background_Pt){fHistograms->fESD_Background_Pt->Fill(sqrt(backgroundCandidate->GetPx()*backgroundCandidate->GetPx()+backgroundCandidate->GetPy()*backgroundCandidate->GetPy()));}
617           if(fHistograms->fESD_Background_Eta){fHistograms->fESD_Background_Eta->Fill(vectorBGCandidate.Eta());}
618           if(fHistograms->fESD_Background_Phi){fHistograms->fESD_Background_Phi->Fill(vectorBGCandidate.Phi());}
619           if(fHistograms->fESD_Background_Mass){fHistograms->fESD_Background_Mass->Fill(massBG);}
620           if(fHistograms->fESD_Background_R){fHistograms->fESD_Background_R->Fill(radiusBG);}
621           if(fHistograms->fESD_Background_Z_R){fHistograms->fESD_Background_Z_R->Fill(tmpY,radiusBG);}
622           if(fHistograms->fESD_Background_X_Y){fHistograms->fESD_Background_X_Y->Fill(tmpX,tmpY);}
623         }
624       }
625       delete backgroundCandidate;   
626     }
627   }
628 }
629
630 void AliAnalysisTaskGammaConversion::Terminate(Option_t */*option*/)
631 {
632   // Terminate analysis
633   //
634   AliDebug(1,"Do nothing in Terminate");
635 }
636
637 void AliAnalysisTaskGammaConversion::UserCreateOutputObjects()
638 {
639   // Create the output container
640
641   fOutputContainer = fHistograms->GetOutputContainer();
642   fOutputContainer->SetName(GetName()) ;  
643 }
644
645 Double_t AliAnalysisTaskGammaConversion::GetMCOpeningAngle(TParticle* daughter0,TParticle* daughter1){
646   //helper function
647   TVector3 V3D0(daughter0->Px(),daughter0->Py(),daughter0->Pz());
648   TVector3 V3D1(daughter1->Px(),daughter1->Py(),daughter1->Pz());
649   return V3D0.Angle(V3D1);
650 }
651