AliPhysics  bba8f44 (bba8f44)
AliV0ReaderV1.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Developers: Friederike Bock *
5  * Original Authors: Svein Lindal, Daniel Lohner *
6  * *
7  * Version 1.0 *
8  * *
9  * *
10  * based on: on older version (see aliroot up to v5-04-42-AN) *
11  * AliV0Reader.cxx *
12  * Authors: Kathrin Koch, Kenneth Aamodt, Ana Marin *
13  * *
14  * Permission to use, copy, modify and distribute this software and its *
15  * documentation strictly for non-commercial purposes is hereby granted *
16  * without fee, provided that the above copyright notice appears in all *
17  * copies and that both the copyright notice and this permission notic *
18  * appear in the supporting documentation. The authors make no claims *
19  * about the suitability of this software for any purpose. It is *
20  * provided "as is" without express or implied warranty. *
21  **************************************************************************/
22 
24 //---------------------------------------------
25 // Class reconstructing conversion photons from V0s
26 //---------------------------------------------
28 
29 // The AliAODConversionPhotons will return the Position (== ESDTrack->GetID())of the positive (negative) ESDTrack
30 // in the Array of ESDTracks stored in the ESDEvent via GetTrackLabelPositive() (GetTrackLabelNagative()).
31 // Since it is the Position in the Array it is always positive.
32 // After the conversion ESD --> AOD each AOD track will give you the former Position in the ESDArray via the
33 // Function AODTrack->GetID(). AODTracks are stored for different TrackParameter (e.g. HybridTracks, TPC only ...). No Standard
34 // AODTracks (not copies of AliExternalTrackParam) will be stored with negative values of AODTrack->GetID() (Formula: (ESDTrack->GetID()+1)*-1)
35 // This leads to the possibility of more than one AOD track with the same negative ID. For that reason all AOD tracks are additionally flaged.
36 // In this analysis we should therefore only use AODTracks with positive values of AODTrack->GetID().
37 
38 // If you want to find the AODTrack corresponding to the daugher track of a AliAODConversionPhoton you have to find the AODTrack with
39 // AODTrack->GetID() == GetTrackLabelPositive() (GetTrackLabelNagative()).
40 
41 #include <vector>
42 #include <TGeoGlobalMagField.h>
43 
44 #include "AliV0ReaderV1.h"
45 #include "AliKFParticle.h"
46 #include "AliAODv0.h"
47 #include "AliESDv0.h"
48 #include "AliAODEvent.h"
49 #include "AliESDEvent.h"
50 #include "AliPID.h"
51 #include "AliMCEvent.h"
52 #include "AliMCEventHandler.h"
53 #include "AliESDpid.h"
54 #include "AliESDtrackCuts.h"
55 #include "TRandom3.h"
56 #include "AliGenCocktailEventHeader.h"
57 #include "TList.h"
58 #include "AliKFConversionPhoton.h"
59 #include "AliAODConversionPhoton.h"
60 #include "AliConversionPhotonBase.h"
61 #include "TVector.h"
62 #include "AliKFVertex.h"
63 #include "AliAODTrack.h"
64 #include "AliESDtrack.h"
65 #include "AliAnalysisManager.h"
66 #include "AliInputEventHandler.h"
67 #include "AliAODHandler.h"
68 #include "AliPIDResponse.h"
69 #include "TChain.h"
70 #include "TFile.h"
71 #include "TString.h"
72 #include "TObjArray.h"
73 
74 class iostream;
75 
76 using namespace std;
77 
78 ClassImp(AliV0ReaderV1)
79 
80 //________________________________________________________________________
81 AliV0ReaderV1::AliV0ReaderV1(const char *name) : AliAnalysisTaskSE(name),
82  kAddv0sInESDFilter(kFALSE),
83  fPCMv0BitField(NULL),
84  fConversionCuts(NULL),
85  fEventCuts(NULL),
86  fConversionGammas(NULL),
87  fUseImprovedVertex(kTRUE),
88  fUseOwnXYZCalculation(kTRUE),
89  fUseConstructGamma(kFALSE),
90  kUseAODConversionPhoton(kTRUE),
91  fCreateAOD(kFALSE),
92  fDeltaAODBranchName("GammaConv"),
93  fDeltaAODFilename("AliAODGammaConversion.root"),
94  fRelabelAODs(kFALSE),
95  fPreviousV0ReaderPerformsAODRelabeling(0),
96  fEventIsSelected(kFALSE),
97  fNumberOfPrimaryTracks(0),
98  fNumberOfTPCoutTracks(0),
99  fPeriodName(""),
100  fPtHardBin(0),
101  fUseMassToZero(kTRUE),
102  fProduceV0findingEffi(kFALSE),
103  fProduceImpactParamHistograms(kFALSE),
104  fCurrentInvMassPair(0),
105  fImprovedPsiPair(3),
106  fHistograms(NULL),
107  fImpactParamHistograms(NULL),
108  fHistoMCGammaPtvsR(NULL),
109  fHistoMCGammaPtvsPhi(NULL),
110  fHistoMCGammaPtvsEta(NULL),
111  fHistoMCGammaRvsPhi(NULL),
112  fHistoMCGammaRvsEta(NULL),
113  fHistoMCGammaPhivsEta(NULL),
114  fHistoRecMCGammaPtvsR(NULL),
115  fHistoRecMCGammaPtvsPhi(NULL),
116  fHistoRecMCGammaPtvsEta(NULL),
117  fHistoRecMCGammaRvsPhi(NULL),
118  fHistoRecMCGammaRvsEta(NULL),
119  fHistoRecMCGammaPhivsEta(NULL),
120  fHistoRecMCGammaMultiPt(NULL),
121  fHistoRecMCGammaMultiPtvsEta(NULL),
122  fHistoRecMCGammaMultiR(NULL),
123  fHistoRecMCGammaMultiPhi(NULL),
124  fHistoPosTrackImpactParamZ(NULL),
125  fHistoPosTrackImpactParamY(NULL),
126  fHistoPosTrackImpactParamX(NULL),
127  fHistoPosTrackImpactParamZvsPt(NULL),
128  fHistoPosTrackImpactParamYvsPt(NULL),
129  fHistoPosTrackImpactParamXvsPt(NULL),
130  fHistoNegTrackImpactParamZ(NULL),
131  fHistoNegTrackImpactParamY(NULL),
132  fHistoNegTrackImpactParamX(NULL),
133  fHistoNegTrackImpactParamZvsPt(NULL),
134  fHistoNegTrackImpactParamYvsPt(NULL),
135  fHistoNegTrackImpactParamXvsPt(NULL),
136  fHistoImpactParamZvsR(NULL),
137  fHistoImpactParamZvsR2(NULL),
138  fHistoPt(NULL),
139  fHistoPt2(NULL),
140  fHistoDCAzPhoton(NULL),
141  fHistoDCAzPhoton2(NULL),
142  fHistoR(NULL),
143  fHistoRrecalc(NULL),
144  fHistoRviaAlpha(NULL),
145  fHistoRviaAlphaRecalc(NULL),
146  fHistoRdiff(NULL),
147  fHistoImpactParameterStudy(NULL),
148  fImpactParamTree(NULL),
149  fVectorFoundGammas(0),
150  fCurrentFileName(""),
151  fMCFileChecked(kFALSE)
152 {
153  // Default constructor
154 
155  DefineInput(0, TChain::Class());
156  DefineOutput(1,TBits::Class());
157 }
158 
159 //________________________________________________________________________
160 AliV0ReaderV1::~AliV0ReaderV1()
161 {
162  // default deconstructor
163 
164  if(fConversionGammas){
165  fConversionGammas->Delete();// Clear Objects
166  delete fConversionGammas;
167  fConversionGammas=0x0;
168  }
169 }
170 
179 AliConversionPhotonBase *AliV0ReaderV1::operator[](int index) const {
180  if(fConversionGammas && (index < fConversionGammas->GetEntriesFast())) {
181  TObject *tmp = fConversionGammas->At(index);
182  if(tmp->IsA() == AliAODConversionPhoton::Class()){
183  AliAODConversionPhoton *photon = static_cast<AliAODConversionPhoton *>(tmp);
184  return photon;
185  } else {
186  AliKFConversionPhoton *photon = static_cast<AliKFConversionPhoton *>(tmp);
187  return photon;
188  }
189  } else {
190  return nullptr;
191  }
192 }
193 
194 /*
195 //________________________________________________________________________
196 AliV0ReaderV1::AliV0ReaderV1(AliV0ReaderV1 &original) : AliAnalysisTaskSE(original),
197 fConversionCuts(NULL),
198 fConversionGammas(NULL),
199 fUseImprovedVertex(original.fUseImprovedVertex),
200 fUseOwnXYZCalculation(original.fUseOwnXYZCalculation),
201 fUseConstructGamma(original.fUseConstructGamma),
202 kUseAODConversionPhoton(original.kUseAODConversionPhoton),
203 fCreateAOD(original.fCreateAOD),
204 fDeltaAODBranchName(original.fDeltaAODBranchName),
205 fDeltaAODFilename(original.fDeltaAODFilename),
206 fEventIsSelected(original.fEventIsSelected)
207 {
208 // Default constructor
209 
210 DefineInput(0, TChain::Class());
211 }
212 
213 //____________________________________________________________
214 AliV0ReaderV1 &AliV0ReaderV1::operator=(const AliV0ReaderV1 &ref){
215 //
216 // Assignment operator
217 // Only copies pointers, object is not the owner of the references
218 //
219 if(this != &ref){
220 AliAnalysisTaskSE::operator=(ref);
221 fUseImprovedVertex=ref.fUseImprovedVertex;
222 fUseOwnXYZCalculation=ref.fUseOwnXYZCalculation;
223 fUseConstructGamma=ref.fUseConstructGamma;
224 kUseAODConversionPhoton=ref.kUseAODConversionPhoton;
225 fCreateAOD=ref.fCreateAOD;
226 fDeltaAODBranchName=ref.fDeltaAODBranchName;
227 fDeltaAODFilename=ref.fDeltaAODFilename;
228 fEventIsSelected=ref.fEventIsSelected;
229 }
230 return *this;
231 }
232 */
233 
234 //________________________________________________________________________
235 void AliV0ReaderV1::Init()
236 {
237  // Initialize function to be called once before analysis
238  if(fConversionCuts==NULL){
239  if(fConversionCuts==NULL)AliError("No Conversion Cut Selection initialized");
240  }
241  if(fEventCuts==NULL){
242  if(fEventCuts==NULL)AliError("No Event Cut Selection initialized");
243  }
244 
245  if(fCreateAOD){kUseAODConversionPhoton=kTRUE;}
246 
247  if(fConversionGammas != NULL){
248  delete fConversionGammas;
249  fConversionGammas=NULL;
250  }
251 
252  if(fConversionGammas == NULL){
253  if(kUseAODConversionPhoton){
254  fConversionGammas = new TClonesArray("AliAODConversionPhoton",100);}
255  else{
256  fConversionGammas = new TClonesArray("AliKFConversionPhoton",100);}
257  }
258  fConversionGammas->Delete();//Reset the TClonesArray
259 }
260 
261 //________________________________________________________________________
262 void AliV0ReaderV1::UserCreateOutputObjects()
263 {
264  // Create AODs
265 
266  if(fCreateAOD){
267  fPCMv0BitField = new TBits();
268 
269  if (fEventCuts){
270  fDeltaAODBranchName.Append("_");
271  fDeltaAODBranchName.Append(fEventCuts->GetCutNumber());
272  }
273  if(fConversionCuts){
274  fDeltaAODBranchName.Append("_");
275  fDeltaAODBranchName.Append(fConversionCuts->GetCutNumber());
276  fDeltaAODBranchName.Append("_gamma");
277  }
278  fConversionGammas->SetName(fDeltaAODBranchName.Data());
279 
280  AddAODBranch("TClonesArray", &fConversionGammas, fDeltaAODFilename.Data());
281  AliAnalysisManager::GetAnalysisManager()->RegisterExtraFile(fDeltaAODFilename.Data());
282  }
283 
284  if(fProduceImpactParamHistograms){
285  if(fImpactParamHistograms != NULL){
286  delete fImpactParamHistograms;
287  fImpactParamHistograms = NULL;
288  }
289  if(fImpactParamHistograms==NULL){
290  fImpactParamHistograms = new TList();
291  fImpactParamHistograms->SetOwner(kTRUE);
292  fImpactParamHistograms->SetName(Form("ImpactParamHistograms_%s_%s",fEventCuts->GetCutNumber().Data(),fConversionCuts->GetCutNumber().Data()));
293  }
294  fHistoPosTrackImpactParamZ = new TH1F("fHistoPosTrackImpactParamZ","",480,-80,80);
295  fHistoPosTrackImpactParamZ->SetXTitle("Z (cm)");
296  fImpactParamHistograms->Add(fHistoPosTrackImpactParamZ);
297 
298  fHistoPosTrackImpactParamY = new TH1F("fHistoPosTrackImpactParamY","",720,-120,120);
299  fHistoPosTrackImpactParamY->SetXTitle("Y (cm)");
300  fImpactParamHistograms->Add(fHistoPosTrackImpactParamY);
301 
302  fHistoPosTrackImpactParamX = new TH1F("fHistoPosTrackImpactParamX","",30,-3,3);
303  fHistoPosTrackImpactParamX->SetXTitle("X (cm)");
304  fImpactParamHistograms->Add(fHistoPosTrackImpactParamX);
305 
306  fHistoNegTrackImpactParamZ = new TH1F("fHistoNegTrackImpactParamZ","",480,-80,80);
307  fHistoNegTrackImpactParamZ->SetXTitle("Z (cm)");
308  fImpactParamHistograms->Add(fHistoNegTrackImpactParamZ);
309 
310  fHistoNegTrackImpactParamY = new TH1F("fHistoNegTrackImpactParamY","",720,-120,120);
311  fHistoNegTrackImpactParamY->SetXTitle("Y (cm)");
312  fImpactParamHistograms->Add(fHistoNegTrackImpactParamY);
313 
314  fHistoNegTrackImpactParamX = new TH1F("fHistoNegTrackImpactParamX","",30,-3,3);
315  fHistoNegTrackImpactParamX->SetXTitle("X (cm)");
316  fImpactParamHistograms->Add(fHistoNegTrackImpactParamX);
317 
318  fHistoPosTrackImpactParamZvsPt = new TH2F("fHistoPosTrackImpactParamZvsPt","",100,0,10,480,-80,80);
319  fHistoPosTrackImpactParamZvsPt->SetYTitle("Z (cm)");
320  fHistoPosTrackImpactParamZvsPt->SetXTitle("Pt (GeV)");
321  fImpactParamHistograms->Add(fHistoPosTrackImpactParamZvsPt);
322 
323  fHistoPosTrackImpactParamYvsPt = new TH2F("fHistoPosTrackImpactParamYvsPt","",100,0,10,720,-120,120);
324  fHistoPosTrackImpactParamYvsPt->SetYTitle("Y (cm)");
325  fHistoPosTrackImpactParamYvsPt->SetXTitle("Pt (GeV)");
326  fImpactParamHistograms->Add(fHistoPosTrackImpactParamYvsPt);
327 
328  fHistoPosTrackImpactParamXvsPt = new TH2F("fHistoPosTrackImpactParamXvsPt","",100,0,10,30,-3,5);
329  fHistoPosTrackImpactParamXvsPt->SetYTitle("X (cm)");
330  fHistoPosTrackImpactParamXvsPt->SetXTitle("Pt (GeV)");
331  fImpactParamHistograms->Add(fHistoPosTrackImpactParamXvsPt);
332 
333  fHistoNegTrackImpactParamZvsPt = new TH2F("fHistoNegTrackImpactParamZvsPt","",100,0,10,480,-80,80);
334  fHistoNegTrackImpactParamZvsPt->SetYTitle("Z (cm)");
335  fHistoNegTrackImpactParamZvsPt->SetXTitle("Pt (GeV)");
336  fImpactParamHistograms->Add(fHistoNegTrackImpactParamZvsPt);
337 
338  fHistoNegTrackImpactParamYvsPt = new TH2F("fHistoNegTrackImpactParamYvsPt","",100,0,10,720,-120,120);
339  fHistoNegTrackImpactParamYvsPt->SetYTitle("Y (cm)");
340  fHistoNegTrackImpactParamYvsPt->SetXTitle("Pt (GeV)");
341  fImpactParamHistograms->Add(fHistoNegTrackImpactParamYvsPt);
342 
343  fHistoNegTrackImpactParamXvsPt = new TH2F("fHistoNegTrackImpactParamXvsPt","",100,0,10,30,-3,3);
344  fHistoNegTrackImpactParamXvsPt->SetYTitle("X (cm)");
345  fHistoNegTrackImpactParamXvsPt->SetXTitle("Pt (GeV)");
346  fImpactParamHistograms->Add(fHistoNegTrackImpactParamXvsPt);
347 
348  fHistoImpactParamZvsR = new TH2F("fHistoImpactParamZvsR","Before cuts",300,-150,150,200,0,200);
349  fHistoImpactParamZvsR->SetXTitle("Z (cm)");
350  fHistoImpactParamZvsR->SetYTitle("R (cm)");
351  fImpactParamHistograms->Add(fHistoImpactParamZvsR);
352 
353  fHistoImpactParamZvsR2 = new TH2F("fHistoImpactParamZvsR2","After cuts",300,-150,150,200,0,200);
354  fHistoImpactParamZvsR2->SetXTitle("Z (cm)");
355  fHistoImpactParamZvsR2->SetYTitle("R (cm)");
356  fImpactParamHistograms->Add(fHistoImpactParamZvsR2);
357 
358  fHistoPt = new TH1F("fHistoPt","Before all cuts",100,0,10);
359  fHistoPt->SetXTitle("Pt (GeV)");
360  fImpactParamHistograms->Add(fHistoPt);
361 
362  fHistoPt2 = new TH1F("fHistoPt2","After all cuts",100,0,10);
363  fHistoPt2->SetXTitle("Pt (GeV)");
364  fImpactParamHistograms->Add(fHistoPt2);
365 
366  fHistoDCAzPhoton = new TH1F("fHistoDCAzPhoton","Before cuts",20,-2,2);
367  fHistoDCAzPhoton->SetXTitle("DCAz photon (cm)");
368  fImpactParamHistograms->Add(fHistoDCAzPhoton);
369 
370  fHistoDCAzPhoton2 = new TH1F("fHistoDCAzPhoton2","After cuts",20,-2,2);
371  fHistoDCAzPhoton2->SetXTitle("DCAz photon (cm)");
372  fImpactParamHistograms->Add(fHistoDCAzPhoton2);
373 
374  fHistoR = new TH1F("fHistoR","",200,0,200);
375  fHistoR->SetXTitle("Conversion radius (cm)");
376  fImpactParamHistograms->Add(fHistoR);
377 
378  fHistoRrecalc = new TH1F("fHistoRrecalc","",200,0,200);
379  fHistoRrecalc->SetXTitle("conversion radius (cm)");
380  fImpactParamHistograms->Add(fHistoRrecalc);
381 
382  fHistoRviaAlpha = new TH1F("fHistoRviaAlpha","",200,0,200);
383  fHistoRviaAlpha->SetXTitle("Conversion radius (cm)");
384  fImpactParamHistograms->Add(fHistoRviaAlpha);
385 
386  fHistoRviaAlphaRecalc = new TH1F("fHistoRviaAlphaRecalc","",200,0,200);
387  fHistoRviaAlphaRecalc->SetXTitle("conversion radius (cm)");
388  fImpactParamHistograms->Add(fHistoRviaAlphaRecalc);
389 
390  fHistoRdiff = new TH1F("fHistoRdiff","",200,0,200);
391  fHistoRdiff->SetXTitle("R_conv - R_cluster conflict (cm)");
392  fImpactParamHistograms->Add(fHistoRdiff);
393 
394  fHistoImpactParameterStudy = new TH1F("fHistoImpactParameterStudy","",7,-0.5,6.5);
395  fHistoImpactParameterStudy->GetXaxis()->SetBinLabel(1,"# V0s");
396  fHistoImpactParameterStudy->GetXaxis()->SetBinLabel(2,"two TPC-only tracks");
397  fHistoImpactParameterStudy->GetXaxis()->SetBinLabel(3,"Z cut not passed");
398  fHistoImpactParameterStudy->GetXaxis()->SetBinLabel(4,"Y cut not passed");
399  fHistoImpactParameterStudy->GetXaxis()->SetBinLabel(5,"R>80cm");
400  fHistoImpactParameterStudy->GetXaxis()->SetBinLabel(6,"causality cut not p.");
401  fHistoImpactParameterStudy->GetXaxis()->SetBinLabel(7,"# removed V0s"); // of all V0s
402  fImpactParamHistograms->Add(fHistoImpactParameterStudy);
403 
404  fImpactParamTree = new TTree("fImpactParamTree","");
405  fImpactParamHistograms->Add(fImpactParamTree);
406  }
407 
408  if (fProduceV0findingEffi){
409  TH1::AddDirectory(kFALSE);
410  if(fHistograms != NULL){
411  delete fHistograms;
412  fHistograms = NULL;
413  }
414  if(fHistograms==NULL){
415  fHistograms = new TList();
416  fHistograms->SetOwner(kTRUE);
417  fHistograms->SetName(Form("V0FindingEfficiencyInput_%s_%s",fEventCuts->GetCutNumber().Data(),fConversionCuts->GetCutNumber().Data()));
418  }
419 
420  fHistoMCGammaPtvsR = new TH2F("MCconvGamma_Pt_R","MC converted gamma Pt vs R (|eta| < 0.9)",250,0.0,25,400,0,200);
421  fHistoMCGammaPtvsR->SetXTitle("p_{MC,T} (GeV/c)");
422  fHistoMCGammaPtvsR->SetYTitle("R_{MC,conv} (cm)");
423  fHistograms->Add(fHistoMCGammaPtvsR);
424 
425  fHistoMCGammaPtvsEta = new TH2F("MCconvGamma_Pt_Eta","MC converted gamma Pt vs Eta ",250,0.0,25,280,-1.4,1.4);
426  fHistoMCGammaPtvsEta->SetXTitle("p_{MC,T} (GeV/c)");
427  fHistoMCGammaPtvsEta->SetYTitle("#eta_{MC}");
428  fHistograms->Add(fHistoMCGammaPtvsEta);
429 
430  fHistoMCGammaPtvsPhi = new TH2F("MCconvGamma_Pt_Phi","MC converted gamma Pt vs Phi (|eta| < 0.9) ",250,0.0,25,400,0,2*TMath::Pi());
431  fHistoMCGammaPtvsPhi->SetXTitle("p_{MC,T} (GeV/c)");
432  fHistoMCGammaPtvsPhi->SetYTitle("#varphi_{MC} (rad)");
433  fHistograms->Add(fHistoMCGammaPtvsPhi);
434 
435  fHistoMCGammaRvsPhi = new TH2F("MCconvGamma_R_Phi","MC converted gamma R vs Phi (|eta| < 0.9) ",400,0,200,400,0,2*TMath::Pi());
436  fHistoMCGammaRvsPhi->SetXTitle("R_{MC,conv} (cm)");
437  fHistoMCGammaRvsPhi->SetYTitle("#varphi_{MC} (rad)");
438  fHistograms->Add(fHistoMCGammaRvsPhi);
439 
440  fHistoMCGammaRvsEta = new TH2F("MCconvGamma_R_Eta","MC converted gamma R vs Eta ",400,0,200,280,-1.4,1.4);
441  fHistoMCGammaRvsEta->SetXTitle("R_{MC,conv} (cm)");
442  fHistoMCGammaRvsEta->SetYTitle("#eta_{MC}");
443  fHistograms->Add(fHistoMCGammaRvsEta);
444 
445  fHistoMCGammaPhivsEta = new TH2F("MCconvGamma_Phi_Eta","MC converted gamma Phi vs Eta ",400,0,2*TMath::Pi(),280,-1.4,1.4);
446  fHistoMCGammaPhivsEta->SetXTitle("#phi_{MC} (rad)");
447  fHistoMCGammaPhivsEta->SetYTitle("#eta_{MC}");
448  fHistograms->Add(fHistoMCGammaPhivsEta);
449 
450  fHistoRecMCGammaPtvsR = new TH2F("RecMCconvGamma_Pt_R","rec MC converted gamma Pt vs R (|eta| < 0.9)",250,0.0,25,400,0,200);
451  fHistoRecMCGammaPtvsR->SetXTitle("p_{MC,T} (GeV/c)");
452  fHistoRecMCGammaPtvsR->SetYTitle("R_{MC,conv} (cm)");
453  fHistograms->Add(fHistoRecMCGammaPtvsR);
454 
455  fHistoRecMCGammaPtvsEta = new TH2F("RecMCconvGamma_Pt_Eta","rec MC converted gamma Pt vs Eta ",250,0.0,25,280,-1.4,1.4);
456  fHistoRecMCGammaPtvsEta->SetXTitle("p_{MC,T} (GeV/c)");
457  fHistoRecMCGammaPtvsEta->SetYTitle("#eta_{MC}");
458  fHistograms->Add(fHistoRecMCGammaPtvsEta);
459 
460  fHistoRecMCGammaPtvsPhi = new TH2F("RecMCconvGamma_Pt_Phi","rec MC converted gamma Pt vs Phi (|eta| < 0.9) ",250,0.0,25,400,0,2*TMath::Pi());
461  fHistoRecMCGammaPtvsPhi->SetXTitle("p_{MC,T} (GeV/c)");
462  fHistoRecMCGammaPtvsPhi->SetYTitle("#varphi_{MC} (rad)");
463  fHistograms->Add(fHistoRecMCGammaPtvsPhi);
464 
465  fHistoRecMCGammaRvsPhi = new TH2F("RecMCconvGamma_R_Phi","rec MC converted gamma R vs Phi (|eta| < 0.9) ",400,0,200,400,0,2*TMath::Pi());
466  fHistoRecMCGammaRvsPhi->SetXTitle("R_{MC,conv} (cm)");
467  fHistoRecMCGammaRvsPhi->SetYTitle("#varphi_{MC} (rad)");
468  fHistograms->Add(fHistoRecMCGammaRvsPhi);
469 
470  fHistoRecMCGammaRvsEta = new TH2F("RecMCconvGamma_R_Eta","rec MC converted gamma R vs Eta ",400,0,200,280,-1.4,1.4);
471  fHistoRecMCGammaRvsEta->SetXTitle("R_{MC,conv} (cm)");
472  fHistoRecMCGammaRvsEta->SetYTitle("#eta_{MC}");
473  fHistograms->Add(fHistoRecMCGammaRvsEta);
474 
475  fHistoRecMCGammaPhivsEta = new TH2F("RecMCconvGamma_Phi_Eta","rec MC converted gamma Phi vs Eta ",400,0,2*TMath::Pi(),280,-1.4,1.4);
476  fHistoRecMCGammaPhivsEta->SetXTitle("#phi_{MC} (rad)");
477  fHistoRecMCGammaPhivsEta->SetYTitle("#eta_{MC}");
478  fHistograms->Add(fHistoRecMCGammaPhivsEta);
479 
480  fHistoRecMCGammaMultiPtvsEta = new TH2F("RecMCconvGammaMulti_Pt_Eta","rec MC converted gamma (at least double counted) Pt vs Eta ",250,0.0,25,280,-1.4,1.4);
481  fHistoRecMCGammaMultiPtvsEta->SetXTitle("p_{MC,T} (GeV/c)");
482  fHistoRecMCGammaMultiPtvsEta->SetYTitle("#eta_{MC}");
483  fHistograms->Add(fHistoRecMCGammaMultiPtvsEta);
484 
485  fHistoRecMCGammaMultiPt = new TH1F("RecMCconvGammaMulti_Pt","rec MC converted gamma (at least double counted) Pt (|eta| < 0.9)",250,0.0,25);
486  fHistoRecMCGammaMultiPt->SetXTitle("p_{MC,T} (GeV/c)");
487  fHistograms->Add(fHistoRecMCGammaMultiPt);
488 
489  fHistoRecMCGammaMultiR = new TH1F("RecMCconvGammaMulti_R","rec MC converted gamma (at least double counted) R (|eta| < 0.9)",400,0,200);
490  fHistoRecMCGammaMultiR->SetXTitle("R_{MC,conv} (cm)");
491  fHistograms->Add(fHistoRecMCGammaMultiR);
492 
493  fHistoRecMCGammaMultiPhi = new TH1F("RecMCconvGammaMulti_Phi","rec MC converted gamma (at least double counted) Phi (|eta| < 0.9)",400,0,2*TMath::Pi());
494  fHistoRecMCGammaMultiPhi->SetXTitle("#phi_{MC} (rad)");
495  fHistograms->Add(fHistoRecMCGammaMultiPhi);
496 
497  fVectorFoundGammas.clear();
498  }
499 
500 }
501 //________________________________________________________________________
502 Bool_t AliV0ReaderV1::Notify(){
503 
504  // set file name to empty & reset check flag
505  fCurrentFileName = "";
506  fMCFileChecked = kFALSE;
507 
508  // obtain file name from analysis manager
509  AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
510  if(man) {
511  AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
512  if (inputHandler){
513  TTree* tree = (TTree*) inputHandler->GetTree();
514  TFile* file = (TFile*) tree->GetCurrentFile();
515  fCurrentFileName = file->GetName();
516  }
517  }
518 
519  // try to extract period name from file name if not given
520  if (fPeriodName.CompareTo("") == 0){
521  TObjArray *arr = fCurrentFileName.Tokenize("/");
522  fPtHardBin = -1;
523  for (Int_t i = 0; i < arr->GetEntriesFast();i++ ){
524  TObjString* testObjString = (TObjString*)arr->At(i);
525  if (testObjString->GetString().BeginsWith("LHC")){
526  fPeriodName = testObjString->GetString();
527  i = arr->GetEntriesFast();
528  }
529  }
530  if (fPeriodName.CompareTo("")==0){
531  TObjArray *arr2 = fCurrentFileName.Tokenize("__");
532  for (Int_t i = 0; i < arr->GetEntriesFast();i++ ){
533  TObjString* testObjString = (TObjString*)arr2->At(i);
534  if (testObjString->GetString().BeginsWith("LHC")){
535  fPeriodName = testObjString->GetString();
536  i = arr2->GetEntriesFast();
537  }
538  }
539  }
540  if (fPeriodName.CompareTo("") != 0 && fEventCuts->GetPeriodEnum() == AliConvEventCuts::kNoPeriod ){
541  fEventCuts->SetPeriodEnum (fPeriodName);
542  }
543  } else {
544  if(fEventCuts->GetPeriodEnum() == AliConvEventCuts::kNoPeriod ){
545  fEventCuts->SetPeriodEnum (fPeriodName);
546  }
547  }
548  // set pt hard bin from file name
549  if ( fEventCuts->GetPeriodEnum() == AliConvEventCuts::kLHC15a3a || fEventCuts->GetPeriodEnum() == AliConvEventCuts::kLHC15a3a_plus ||
550  fEventCuts->GetPeriodEnum() == AliConvEventCuts::kLHC15a3b ||
551  fEventCuts->GetPeriodEnum() == AliConvEventCuts::kLHC15g1a || fEventCuts->GetPeriodEnum() == AliConvEventCuts::kLHC15g1b ||
552  fEventCuts->GetPeriodEnum() == AliConvEventCuts::kLHC16c2 || fEventCuts->GetPeriodEnum() == AliConvEventCuts::kLHC16c2_plus ||
553  fEventCuts->GetPeriodEnum() == AliConvEventCuts::kLHC16c3a || fEventCuts->GetPeriodEnum() == AliConvEventCuts::kLHC16c3b ||
554  fEventCuts->GetPeriodEnum() == AliConvEventCuts::kLHC16c3c ||
555  fEventCuts->GetPeriodEnum() == AliConvEventCuts::kLHC16P1JJ || fEventCuts->GetPeriodEnum() == AliConvEventCuts::kLHC16P1JJLowB
556  ){
557  TObjArray *arr = fCurrentFileName.Tokenize("/");
558  fPtHardBin = -1;
559  for (Int_t i = 0; i < arr->GetEntriesFast();i++ ){
560  TObjString* testObjString = (TObjString*)arr->At(i);
561  if (testObjString->GetString().BeginsWith("LHC")){
562  TObjString* testObjString2 = (TObjString*)arr->At(i+1);
563  fPtHardBin = testObjString2->GetString().Atoi();
564  i = arr->GetEntriesFast();
565  }
566  }
567  }
568 
569  if(!fEventCuts->GetDoEtaShift()) return kTRUE; // No Eta Shift requested, continue
570  if(fEventCuts->GetEtaShift() == 0.0){ // Eta Shift requested but not set, get shift automatically
571  fEventCuts->GetCorrectEtaShiftFromPeriod();
572  fEventCuts->DoEtaShift(kFALSE); // Eta Shift Set, make sure that it is called only once
573  return kTRUE;
574  } else {
575  printf(" Gamma Conversion Reader %s_%s :: Eta Shift Manually Set to %f \n\n",
576  (fEventCuts->GetCutNumber()).Data(),(fConversionCuts->GetCutNumber()).Data(),fEventCuts->GetEtaShift());
577  fEventCuts->DoEtaShift(kFALSE); // Eta Shift Set, make sure that it is called only once
578  }
579 
580  return kTRUE;
581 }
582 //________________________________________________________________________
583 void AliV0ReaderV1::UserExec(Option_t *option){
584 
585  AliESDEvent * esdEvent = dynamic_cast<AliESDEvent*>(fInputEvent);
586  if(esdEvent) {
587  if (!TGeoGlobalMagField::Instance()->GetField()) esdEvent->InitMagneticField();
588  }
589 
590  // Check if correctly initialized
591  if(!fConversionGammas)Init();
592 
593  // User Exec
594  fEventIsSelected=ProcessEvent(fInputEvent,fMCEvent);
595 
596  // AOD Output
597  FillAODOutput();
598 
599 }
600 
601 //________________________________________________________________________
602 Bool_t AliV0ReaderV1::ProcessEvent(AliVEvent *inputEvent,AliMCEvent *mcEvent)
603 {
604 
605  //Reset the TClonesArray
606  fConversionGammas->Delete();
607 
608  //Clear TBits object with accepted v0s from previous event
609  if (kAddv0sInESDFilter){fPCMv0BitField->Clear();}
610 
611  fInputEvent = inputEvent;
612  fMCEvent = mcEvent;
613 
614  if(!fInputEvent){
615  AliError("No Input event");
616  return kFALSE;
617  }
618  if(!fEventCuts){AliError("No EventCuts");return kFALSE;}
619  if(!fConversionCuts){AliError("No ConversionCuts");return kFALSE;}
620 
621 
622  // Count Primary Tracks Event
623  CountTracks();
624 
625  //Count Tracks with TPCout flag
626  CountTPCoutTracks();
627 
628  // Event Cuts
629  if(!fEventCuts->EventIsSelected(fInputEvent,fMCEvent)){
630  if (fEventCuts->GetEventQuality() == 2 && !fMCFileChecked ){
631  cout << "ERROR with MC reading for: "<< fCurrentFileName.Data() << endl;
632  fMCFileChecked = kTRUE;
633  }
634  return kFALSE;
635  }
636  // Set Magnetic Field
637  AliKFParticle::SetField(fInputEvent->GetMagneticField());
638 
639  if(fInputEvent->IsA()==AliAODEvent::Class() && fProduceV0findingEffi){
640  fProduceV0findingEffi = kFALSE;
641  AliWarning("V0finding effi cannot be run on AODs ");
642  }
643 
644  if(fProduceV0findingEffi){
645  CreatePureMCHistosForV0FinderEffiESD();
646  fVectorFoundGammas.clear();
647  }
648 
649  if(fInputEvent->IsA()==AliESDEvent::Class()){
650  ProcessESDV0s();
651  }
652  if(fInputEvent->IsA()==AliAODEvent::Class()){
653  GetAODConversionGammas();
654  }
655 
656  return kTRUE;
657 }
659 void AliV0ReaderV1::FillAODOutput()
660 {
661  // Fill AOD Output with reconstructed Photons
662 
663  if(fInputEvent->IsA()==AliESDEvent::Class()){
665  if(fCreateAOD) {
666  PostData(1, fPCMv0BitField);
667  AliAODHandler * aodhandler = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
668  if (aodhandler && aodhandler->GetFillAOD()) {
669  AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillExtension(kTRUE);
670  //PostData(0, fConversionGammas);
671 
672  }
673  }
674  }
675 }
676 
678 const AliExternalTrackParam *AliV0ReaderV1::GetExternalTrackParam(AliESDv0 *fCurrentV0,Int_t &tracklabel,Int_t charge){
679 
680  // Get External Track Parameter with given charge
681 
682  if(!(charge==1||charge==-1)){AliError("Charge not defined");return 0x0;}
683 
684  if(fConversionCuts->GetV0FinderSameSign()==1){
685  if(fCurrentV0){
686  if((fConversionCuts->GetTrack(fInputEvent,fCurrentV0->GetPindex()))->Charge()!=(fConversionCuts->GetTrack(fInputEvent,fCurrentV0->GetNindex()))->Charge())return 0x0;
687  if((fConversionCuts->GetTrack(fInputEvent,fCurrentV0->GetPindex()))->Charge()==(fConversionCuts->GetTrack(fInputEvent,fCurrentV0->GetNindex()))->Charge()){
688  if(charge==1){
689  tracklabel=fCurrentV0->GetPindex();
690  return fCurrentV0->GetParamP();
691  }else{
692  tracklabel=fCurrentV0->GetNindex();
693  return fCurrentV0->GetParamN();
694  }
695  }
696  }
697  }else if(fConversionCuts->GetV0FinderSameSign()==2){
698  if(fCurrentV0){
699  if(charge==1){
700  tracklabel=fCurrentV0->GetPindex();
701  return fCurrentV0->GetParamP();
702  }else{
703  tracklabel=fCurrentV0->GetNindex();
704  return fCurrentV0->GetParamN();
705  }
706  }
707  }else{
708  // Check for sign flip
709  if(fCurrentV0){
710  if(!fCurrentV0->GetParamN()||!fCurrentV0->GetParamP())return 0x0;
711  if(!fConversionCuts->GetTrack(fInputEvent,fCurrentV0->GetNindex())||!fConversionCuts->GetTrack(fInputEvent,fCurrentV0->GetPindex()))return 0x0;
712  if((fConversionCuts->GetTrack(fInputEvent,fCurrentV0->GetPindex()))->Charge()==charge){
713  tracklabel=fCurrentV0->GetPindex();
714  return fCurrentV0->GetParamP();}
715  if((fConversionCuts->GetTrack(fInputEvent,fCurrentV0->GetNindex()))->Charge()==charge){
716  tracklabel=fCurrentV0->GetNindex();
717  return fCurrentV0->GetParamN();}
718  }
719  }
720  return 0x0;
721 }
722 
724 Bool_t AliV0ReaderV1::ProcessESDV0s()
725 {
726 
727  // Process ESD V0s for conversion photon reconstruction
728  AliESDEvent *fESDEvent=dynamic_cast<AliESDEvent*>(fInputEvent);
729 
730  AliKFConversionPhoton *fCurrentMotherKFCandidate=NULL;
731 
732  if(fESDEvent){
733  for(Int_t currentV0Index=0;currentV0Index<fESDEvent->GetNumberOfV0s();currentV0Index++){
734  AliESDv0 *fCurrentV0=(AliESDv0*)(fESDEvent->GetV0(currentV0Index));
735  if(!fCurrentV0){
736  printf("Requested V0 does not exist");
737  continue;
738  }
739 
740  fCurrentMotherKFCandidate=ReconstructV0(fCurrentV0,currentV0Index);
741 
742  if(fCurrentMotherKFCandidate){
743  // Add Gamma to the TClonesArray
744 
745  if(kUseAODConversionPhoton){
746  new((*fConversionGammas)[fConversionGammas->GetEntriesFast()]) AliAODConversionPhoton(fCurrentMotherKFCandidate);
747  AliAODConversionPhoton * currentConversionPhoton = (AliAODConversionPhoton*)(fConversionGammas->At(fConversionGammas->GetEntriesFast()-1));
748  currentConversionPhoton->SetMass(fCurrentMotherKFCandidate->M());
749  if (fUseMassToZero) currentConversionPhoton->SetMassToZero();
750  currentConversionPhoton->SetInvMassPair(fCurrentInvMassPair);
751  if(kAddv0sInESDFilter){fPCMv0BitField->SetBitNumber(currentV0Index, kTRUE);}
752  } else {
753  new((*fConversionGammas)[fConversionGammas->GetEntriesFast()]) AliKFConversionPhoton(*fCurrentMotherKFCandidate);
754  }
755 
756  delete fCurrentMotherKFCandidate;
757  fCurrentMotherKFCandidate=NULL;
758  }
759  }
760  if(kAddv0sInESDFilter){fPCMv0BitField->Compact();}
761  }
762  return kTRUE;
763 }
764 
766 AliKFConversionPhoton *AliV0ReaderV1::ReconstructV0(AliESDv0 *fCurrentV0,Int_t currentV0Index)
767 {
768  // cout << currentV0Index << endl;
769  // Reconstruct conversion photon from ESD v0
770  fConversionCuts->FillPhotonCutIndex(AliConversionPhotonCuts::kPhotonIn);
771 
772  //checks if on the fly mode is set
773  if(!fConversionCuts->SelectV0Finder(fCurrentV0->GetOnFlyStatus())){
774  fConversionCuts->FillPhotonCutIndex(AliConversionPhotonCuts::kOnFly);
775  return 0x0;
776  }
777 
778  if (fMCEvent && fProduceV0findingEffi ) FillRecMCHistosForV0FinderEffiESD(fCurrentV0);
779 
780  // TrackLabels
781  Int_t currentTrackLabels[2]={-1,-1};
782 
783  // Get Daughter KF Particles
784 
785  const AliExternalTrackParam *fCurrentExternalTrackParamPositive=GetExternalTrackParamP(fCurrentV0,currentTrackLabels[0]);
786  // cout << fCurrentExternalTrackParamPositive << "\t" << currentTrackLabels[0] << endl;
787  const AliExternalTrackParam *fCurrentExternalTrackParamNegative=GetExternalTrackParamN(fCurrentV0,currentTrackLabels[1]);
788  // cout << fCurrentExternalTrackParamNegative << "\t" << currentTrackLabels[1] << endl;
789  if(!fCurrentExternalTrackParamPositive||!fCurrentExternalTrackParamNegative)return 0x0;
790 
791  // Apply some Cuts before Reconstruction
792 
793  AliVTrack * posTrack = fConversionCuts->GetTrack(fInputEvent,currentTrackLabels[0]);
794  AliVTrack * negTrack = fConversionCuts->GetTrack(fInputEvent,currentTrackLabels[1]);
795  if(!negTrack || !posTrack) {
796  fConversionCuts->FillPhotonCutIndex(AliConversionPhotonCuts::kNoTracks);
797  return 0x0;
798  }
799  // Track Cuts
800  if(!fConversionCuts->TracksAreSelected(negTrack, posTrack)){
801  fConversionCuts->FillPhotonCutIndex(AliConversionPhotonCuts::kTrackCuts);
802  return 0x0;
803  }
804 
805  fConversionCuts->FillV0EtaBeforedEdxCuts(fCurrentV0->Eta());
806  if (!fConversionCuts->dEdxCuts(posTrack)) {
807  fConversionCuts->FillPhotonCutIndex(AliConversionPhotonCuts::kdEdxCuts);
808  return 0x0;
809  }
810  // PID Cuts
811  if(!fConversionCuts->dEdxCuts(negTrack)) {
812  fConversionCuts->FillPhotonCutIndex(AliConversionPhotonCuts::kdEdxCuts);
813  return 0x0;
814  }
815  fConversionCuts->FillV0EtaAfterdEdxCuts(fCurrentV0->Eta());
816  // Reconstruct Photon
817  AliKFConversionPhoton *fCurrentMotherKF=NULL;
818  // fUseConstructGamma = kFALSE;
819  // cout << "construct gamma " << endl;
820  AliKFParticle fCurrentNegativeKFParticle(*(fCurrentExternalTrackParamNegative),11);
821  // cout << fCurrentExternalTrackParamNegative << "\t" << endl;
822  AliKFParticle fCurrentPositiveKFParticle(*(fCurrentExternalTrackParamPositive),-11);
823  // cout << fCurrentExternalTrackParamPositive << "\t" << endl;
824  // cout << currentTrackLabels[0] << "\t" << currentTrackLabels[1] << endl;
825  // cout << "construct gamma " <<fUseConstructGamma << endl;
826 
827  // Reconstruct Gamma
828  if(fUseConstructGamma){
829  fCurrentMotherKF = new AliKFConversionPhoton();
830  fCurrentMotherKF->ConstructGamma(fCurrentNegativeKFParticle,fCurrentPositiveKFParticle);
831  }else{
832  fCurrentMotherKF = new AliKFConversionPhoton(fCurrentNegativeKFParticle,fCurrentPositiveKFParticle);
833  fCurrentMotherKF->SetMassConstraint(0,0.0001);
834  }
835 
836  // Set Track Labels
837 
838  fCurrentMotherKF->SetTrackLabels(currentTrackLabels[0],currentTrackLabels[1]);
839 
840  // Set V0 index
841 
842  fCurrentMotherKF->SetV0Index(currentV0Index);
843 
844  //Set MC Label
845  if(fMCEvent){
846 
847  Int_t labelp=TMath::Abs(fConversionCuts->GetTrack(fInputEvent,fCurrentMotherKF->GetTrackLabelPositive())->GetLabel());
848  Int_t labeln=TMath::Abs(fConversionCuts->GetTrack(fInputEvent,fCurrentMotherKF->GetTrackLabelNegative())->GetLabel());
849 
850 // cout << "rec: " << currentTrackLabels[0] << "\t" << currentTrackLabels[1] << endl;
851 // cout << "recProp: " << fCurrentMotherKF->GetTrackLabelPositive() << "\t" << fCurrentMotherKF->GetTrackLabelNegative() << endl;
852 // cout << "MC: " << labeln << "\t" << labelp << endl;
853 
854  TParticle *fNegativeMCParticle = 0x0;
855  if(labeln>-1) fNegativeMCParticle = fMCEvent->Particle(labeln);
856  TParticle *fPositiveMCParticle = 0x0;
857  if(labelp>-1) fPositiveMCParticle = fMCEvent->Particle(labelp);
858 
859  if(fPositiveMCParticle&&fNegativeMCParticle){
860  fCurrentMotherKF->SetMCLabelPositive(labelp);
861  fCurrentMotherKF->SetMCLabelNegative(labeln);
862  }
863  }
864 
865  // Update Vertex (moved for same eta compared to old)
866  // cout << currentV0Index <<" \t before: \t" << fCurrentMotherKF->GetPx() << "\t" << fCurrentMotherKF->GetPy() << "\t" << fCurrentMotherKF->GetPz() << endl;
867  if(fUseImprovedVertex == kTRUE){
868  AliKFVertex primaryVertexImproved(*fInputEvent->GetPrimaryVertex());
869  // cout << "Prim Vtx: " << primaryVertexImproved.GetX() << "\t" << primaryVertexImproved.GetY() << "\t" << primaryVertexImproved.GetZ() << endl;
870  primaryVertexImproved+=*fCurrentMotherKF;
871  fCurrentMotherKF->SetProductionVertex(primaryVertexImproved);
872  }
873  // SetPsiPair
874  Double_t convpos[3]={0,0,0};
875  if (fImprovedPsiPair == 0){
876  Double_t PsiPair=GetPsiPair(fCurrentV0,fCurrentExternalTrackParamPositive,fCurrentExternalTrackParamNegative, convpos);
877  fCurrentMotherKF->SetPsiPair(PsiPair);
878  }
879 
880  // Recalculate ConversionPoint
881  Double_t dca[2]={0,0};
882  if(fUseOwnXYZCalculation){
883  // Double_t convpos[3]={0,0,0};
884  if(!GetConversionPoint(fCurrentExternalTrackParamPositive,fCurrentExternalTrackParamNegative,convpos,dca)){
885  fConversionCuts->FillPhotonCutIndex(AliConversionPhotonCuts::kConvPointFail);
886  delete fCurrentMotherKF;
887  fCurrentMotherKF=NULL;
888  return 0x0;
889  }
890 
891  fCurrentMotherKF->SetConversionPoint(convpos);
892  }
893 
894 
895  // SetPsiPair
896  if (fImprovedPsiPair >= 1){
897  // the propagation can be more precise after the precise conversion point calculation
898  Double_t PsiPair=GetPsiPair(fCurrentV0,fCurrentExternalTrackParamPositive,fCurrentExternalTrackParamNegative,convpos);
899  fCurrentMotherKF->SetPsiPair(PsiPair);
900  //cout<<" GetPsiPair::"<<fCurrentMotherKF->GetPsiPair() <<endl;
901  }
902 
903  if(fCurrentMotherKF->GetNDF() > 0.)
904  fCurrentMotherKF->SetChi2perNDF(fCurrentMotherKF->GetChi2()/fCurrentMotherKF->GetNDF()); //->Photon is created before all chi2 relevant changes are performed, set it "by hand"
905 
906 
907  // Set Dilepton Mass (moved down for same eta compared to old)
908  fCurrentMotherKF->SetMass(fCurrentMotherKF->M());
909 
910  // Calculating invariant mass
911  Double_t mass=-99.0, mass_width=-99.0, Pt=-99.0, Pt_width=-99.0;
912  AliKFParticle fCurrentMotherKFForMass(fCurrentNegativeKFParticle,fCurrentPositiveKFParticle);
913  fCurrentMotherKFForMass.GetMass(mass,mass_width);
914  fCurrentMotherKFForMass.GetPt(Pt,Pt_width);
915  fCurrentInvMassPair=mass;
916 
917  // apply possible Kappa cut
918  if (!fConversionCuts->KappaCuts(fCurrentMotherKF,fInputEvent)){
919  fConversionCuts->FillPhotonCutIndex(AliConversionPhotonCuts::kdEdxCuts);
920  delete fCurrentMotherKF;
921  fCurrentMotherKF=NULL;
922  return 0x0;
923  }
924 
925  // Apply Photon Cuts
926  if(!fConversionCuts->PhotonCuts(fCurrentMotherKF,fInputEvent)){
927  fConversionCuts->FillPhotonCutIndex(AliConversionPhotonCuts::kPhotonCuts);
928  delete fCurrentMotherKF;
929  fCurrentMotherKF=NULL;
930  return 0x0;
931  }
932 
933  // cout << currentV0Index <<" \t after: \t" <<fCurrentMotherKF->GetPx() << "\t" << fCurrentMotherKF->GetPy() << "\t" << fCurrentMotherKF->GetPz() << endl;
934 
935  if(fProduceImpactParamHistograms) FillImpactParamHistograms(posTrack, negTrack, fCurrentV0, fCurrentMotherKF);
936 
937  fConversionCuts->FillPhotonCutIndex(AliConversionPhotonCuts::kPhotonOut);
938  return fCurrentMotherKF;
939 }
940 
942 Double_t AliV0ReaderV1::GetPsiPair(const AliESDv0* v0, const AliExternalTrackParam *positiveparam,const AliExternalTrackParam *negativeparam,const Double_t convpos[3]) const {
943  //
944  // Angle between daughter momentum plane and plane
945  //
946 
947  AliExternalTrackParam nt(*negativeparam);
948  AliExternalTrackParam pt(*positiveparam);
949 
950  Float_t magField = fInputEvent->GetMagneticField();
951 
952  Double_t xyz[3] = {0.,0.,0.};
953  if (fImprovedPsiPair==0 ) {
954  v0->GetXYZ(xyz[0],xyz[1],xyz[2]);
955  } else if (fImprovedPsiPair>=1 ) {
956  xyz[0]= convpos[0];
957  xyz[1]= convpos[1];
958  xyz[2]= convpos[2];
959  }
960 
961 
962  // Double_t pPlus[3] = {pt.Px(),pt.Py(),pt.Pz()};
963  // Double_t pMinus[3] = {nt.Px(),nt.Py(),nt.Pz()};
964 
965  // Double_t u[3] = {pPlus[0]+pMinus[0],pPlus[1]+pMinus[1],pPlus[2]+pMinus[2]};
966  // Double_t normu = sqrt( (u[0]*u[0]) + (u[1]*u[1]) + (u[2]*u[2]) );
967 
968  // u[0] = u[0] / normu;
969  // u[1] = u[1] / normu;
970  // u[2] = u[2] / normu;
971 
972  // Double_t normpPlus = sqrt( (pPlus[0]*pPlus[0]) + (pPlus[1]*pPlus[1]) + (pPlus[2]*pPlus[2]) );
973  // Double_t normpMinus = sqrt( (pMinus[0]*pMinus[0]) + (pMinus[1]*pMinus[1]) + (pMinus[2]*pMinus[2]) );
974 
975  // pPlus[0] = pPlus[0] / normpPlus;
976  // pPlus[1] = pPlus[1] / normpPlus;
977  // pPlus[2] = pPlus[2] / normpPlus;
978 
979  // pMinus[0] = pMinus[0] / normpMinus;
980  // pMinus[1] = pMinus[1] / normpMinus;
981  // pMinus[2] = pMinus[2] / normpMinus;
982 
983  // Double_t v[3] = {0,0,0}; // pPlus X pMinus
984  // v[0] = (pPlus[1]*pMinus[2]) - (pPlus[2]*pMinus[1]);
985  // v[1] = (pPlus[2]*pMinus[0]) - (pPlus[0]*pMinus[2]);
986  // v[2] = (pPlus[0]*pMinus[1]) - (pPlus[1]*pMinus[0]);
987 
988  // Double_t w[3] = {0,0,0}; // u X v
989  // w[0] = (u[1]*v[2]) - (u[2]*v[1]);
990  // w[1] = (u[2]*v[0]) - (u[0]*v[2]);
991  // w[2] = (u[0]*v[1]) - (u[1]*v[0]);
992 
993  // Double_t z[3] = {0,0,1};
994  // Double_t wc[3] = {0,0,0}; // u X v
995  // wc[0] = (u[1]*z[2]) - (u[2]*z[1]);
996  // wc[1] = (u[2]*z[0]) - (u[0]*z[2]);
997  // wc[2] = (u[0]*z[1]) - (u[1]*z[0]);
998 
999  // Double_t PhiV = TMath::ACos((w[0]*wc[0]) + (w[1]*wc[1]) + (w[2]*wc[2]));
1000  //return TMath::Abs(PhiV);
1001 
1002 
1003  // TVector3 pPlus(pt.Px(),pt.Py(),pt.Pz());
1004  // TVector3 pMinus(nt.Px(),nt.Py(),nt.Pz());
1005 
1006  // TVector3 u = pMinus + pPlus;
1007  // u = u*(1/u.Mag());
1008 
1009  // TVector3 pHPlus = pPlus*(1/pPlus.Mag());
1010  // TVector3 pHMinus = pMinus*(1/pMinus.Mag());
1011 
1012  // TVector3 v = pHPlus.Cross(pHMinus);
1013  // TVector3 w = u.Cross(v);
1014  // TVector3 z(0,0,1);
1015  // TVector3 wc = u.Cross(z);
1016 
1017  // Double_t PhiV = w * wc;
1018 
1019  Double_t mn[3] = {0,0,0};
1020  Double_t mp[3] = {0,0,0};
1021 
1022  v0->GetNPxPyPz(mn[0],mn[1],mn[2]);//reconstructed cartesian momentum components of negative daughter;
1023  v0->GetPPxPyPz(mp[0],mp[1],mp[2]);//reconstructed cartesian momentum components of positive daughter;
1024 
1025  Double_t deltat = 1.;
1026  deltat = TMath::ATan(mp[2]/(TMath::Sqrt(mp[0]*mp[0] + mp[1]*mp[1])+1.e-13)) - TMath::ATan(mn[2]/(TMath::Sqrt(mn[0]*mn[0] + mn[1]*mn[1])+1.e-13));//difference of angles of the two daughter tracks with z-axis
1027  Double_t radiussum = TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]) + 50;//radius to which tracks shall be propagated
1028 
1029  Double_t momPosProp[3] = {0,0,0};
1030  Double_t momNegProp[3] = {0,0,0};
1031 
1032  Double_t psiPair = 4.;
1033  // cout<< "Momentum before propagation at radius ::"<< TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1])<< endl;
1034  // cout<< mn[0]<< " " <<mn[1]<<" "<< mn[2]<<" "<<TMath::Sqrt(mn[0]*mn[0] + mn[1]*mn[1] +mn[2]*mn[2])<< endl;
1035  // cout<< mp[0]<< " " <<mp[1]<<" "<< mp[2]<<" "<<TMath::Sqrt(mp[0]*mp[0] + mp[1]*mp[1] +mp[2]*mp[2])<< endl;
1036  Double_t pEle,pPos;
1037 
1038  if (fImprovedPsiPair==1 || fImprovedPsiPair==0 ){
1039  if(nt.PropagateTo(radiussum,magField) == 0) return psiPair; //propagate tracks to the outside -> Better Purity and Efficiency
1040  if(pt.PropagateTo(radiussum,magField) == 0) return psiPair; //propagate tracks to the outside -> Better Purity and Efficiency
1041 
1042  pt.GetPxPyPz(momPosProp);//Get momentum vectors of tracks after propagation
1043  nt.GetPxPyPz(momNegProp);
1044  } else if (fImprovedPsiPair>=2) {
1045  momPosProp[0] = pt.GetParameterAtRadius(radiussum,magField,3);
1046  momPosProp[1] = pt.GetParameterAtRadius(radiussum,magField,4);
1047  momPosProp[2] = pt.GetParameterAtRadius(radiussum,magField,5);
1048 
1049  momNegProp[0] = nt.GetParameterAtRadius(radiussum,magField,3);
1050  momNegProp[1] = nt.GetParameterAtRadius(radiussum,magField,4);
1051  momNegProp[2] = nt.GetParameterAtRadius(radiussum,magField,5);
1052  pEle = TMath::Sqrt(momNegProp[0]*momNegProp[0]+momNegProp[1]*momNegProp[1]+momNegProp[2]*momNegProp[2]);//absolute momentum value of negative daughter
1053 
1054  pPos = TMath::Sqrt(momPosProp[0]*momPosProp[0]+momPosProp[1]*momPosProp[1]+momPosProp[2]*momPosProp[2]);//absolute momentum value of positive daughter
1055 
1056  if ( (pEle==0 || pPos==0) && fImprovedPsiPair==3) {
1057  radiussum = TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]) + 30;
1058  momPosProp[0] = pt.GetParameterAtRadius(radiussum,magField,3);
1059  momPosProp[1] = pt.GetParameterAtRadius(radiussum,magField,4);
1060  momPosProp[2] = pt.GetParameterAtRadius(radiussum,magField,5);
1061 
1062  momNegProp[0] = nt.GetParameterAtRadius(radiussum,magField,3);
1063  momNegProp[1] = nt.GetParameterAtRadius(radiussum,magField,4);
1064  momNegProp[2] = nt.GetParameterAtRadius(radiussum,magField,5);
1065 
1066  }
1067  }
1068 
1069  pEle = TMath::Sqrt(momNegProp[0]*momNegProp[0]+momNegProp[1]*momNegProp[1]+momNegProp[2]*momNegProp[2]);//absolute momentum value of negative daughter
1070  pPos = TMath::Sqrt(momPosProp[0]*momPosProp[0]+momPosProp[1]*momPosProp[1]+momPosProp[2]*momPosProp[2]);//absolute momentum value of positive daughter
1071 
1072  Double_t scalarproduct =
1073  momPosProp[0]*momNegProp[0]+momPosProp[1]*momNegProp[1]+momPosProp[2]*momNegProp[2];//scalar product of propagated positive and negative daughters' momenta
1074 
1075  if (pEle==0 || pPos==0) return psiPair;
1076 
1077 
1078  Double_t chipair = TMath::ACos(scalarproduct/(pEle*pPos));//Angle between propagated daughter tracks
1079 
1080  // psiPair = TMath::Abs(TMath::ASin(deltat/chipair));
1081  psiPair = TMath::ASin(deltat/chipair);
1082  return psiPair;
1083 }
1084 
1086 Bool_t AliV0ReaderV1::GetHelixCenter(const AliExternalTrackParam *track,Double_t center[2]){
1087 
1088  // Get Center of the helix track parametrization
1089 
1090  Int_t charge=track->Charge();
1091  Double_t b=fInputEvent->GetMagneticField();
1092 
1093  Double_t helix[6];
1094  track->GetHelixParameters(helix,b);
1095 
1096  Double_t xpos = helix[5];
1097  Double_t ypos = helix[0];
1098  Double_t radius = TMath::Abs(1./helix[4]);
1099  Double_t phi = helix[2];
1100 
1101  if(phi < 0){
1102  phi = phi + 2*TMath::Pi();
1103  }
1104 
1105  phi -= TMath::Pi()/2.;
1106  Double_t xpoint = radius * TMath::Cos(phi);
1107  Double_t ypoint = radius * TMath::Sin(phi);
1108 
1109  if(b<0){
1110  if(charge > 0){
1111  xpoint = - xpoint;
1112  ypoint = - ypoint;
1113  }
1114  }
1115  if(b>0){
1116  if(charge < 0){
1117  xpoint = - xpoint;
1118  ypoint = - ypoint;
1119  }
1120  }
1121  center[0] = xpos + xpoint;
1122  center[1] = ypos + ypoint;
1123 
1124  return 1;
1125 }
1127 Bool_t AliV0ReaderV1::GetConversionPoint(const AliExternalTrackParam *pparam,const AliExternalTrackParam *nparam,Double_t convpos[3],Double_t dca[2]){
1128 
1129  // Recalculate Conversion Point
1130 
1131  if(!pparam||!nparam)return kFALSE;
1132 
1133  Double_t helixcenterpos[2];
1134  GetHelixCenter(pparam,helixcenterpos);
1135 
1136  Double_t helixcenterneg[2];
1137  GetHelixCenter(nparam,helixcenterneg);
1138 
1139  Double_t helixpos[6];
1140  pparam->GetHelixParameters(helixpos,fInputEvent->GetMagneticField());
1141  Double_t posradius = TMath::Abs(1./helixpos[4]);
1142 
1143  Double_t helixneg[6];
1144  nparam->GetHelixParameters(helixneg,fInputEvent->GetMagneticField());
1145  Double_t negradius = TMath::Abs(1./helixneg[4]);
1146 
1147  // Calculate xy-position
1148 
1149  Double_t xpos = helixcenterpos[0];
1150  Double_t ypos = helixcenterpos[1];
1151  Double_t xneg = helixcenterneg[0];
1152  Double_t yneg = helixcenterneg[1];
1153 
1154  convpos[0] = (xpos*negradius + xneg*posradius)/(negradius+posradius);
1155  convpos[1] = (ypos*negradius+ yneg*posradius)/(negradius+posradius);
1156 
1157 
1158  // Calculate Track XY vertex position
1159 
1160  Double_t deltaXPos = convpos[0] - xpos;
1161  Double_t deltaYPos = convpos[1] - ypos;
1162 
1163  Double_t deltaXNeg = convpos[0] - xneg;
1164  Double_t deltaYNeg = convpos[1] - yneg;
1165 
1166  Double_t alphaPos = TMath::Pi() + TMath::ATan2(-deltaYPos,-deltaXPos);
1167  Double_t alphaNeg = TMath::Pi() + TMath::ATan2(-deltaYNeg,-deltaXNeg);
1168 
1169  Double_t vertexXNeg = xneg + TMath::Abs(negradius)*TMath::Cos(alphaNeg);
1170  Double_t vertexYNeg = yneg + TMath::Abs(negradius)*TMath::Sin(alphaNeg);
1171 
1172  Double_t vertexXPos = xpos + TMath::Abs(posradius)*TMath::Cos(alphaPos);
1173  Double_t vertexYPos = ypos + TMath::Abs(posradius)*TMath::Sin(alphaPos);
1174 
1175  AliExternalTrackParam p(*pparam);
1176  AliExternalTrackParam n(*nparam);
1177 
1178  TVector2 vertexPos(vertexXPos,vertexYPos);
1179  TVector2 vertexNeg(vertexXNeg,vertexYNeg);
1180 
1181  // Convert to local coordinate system
1182  TVector2 vertexPosRot=vertexPos.Rotate(-p.GetAlpha());
1183  TVector2 vertexNegRot=vertexNeg.Rotate(-n.GetAlpha());
1184 
1185  // Propagate Track Params to Vertex
1186 
1187  if(!p.PropagateTo(vertexPosRot.X(),fInputEvent->GetMagneticField()))return kFALSE;
1188  if(!n.PropagateTo(vertexNegRot.X(),fInputEvent->GetMagneticField()))return kFALSE;
1189 
1190  // Check whether propagation was sucessful
1191 
1192  if(TMath::Abs(vertexPos.Mod()-TMath::Sqrt(p.GetX()*p.GetX()+p.GetY()*p.GetY()))>0.01)return kFALSE;
1193  if(TMath::Abs(vertexNeg.Mod()-TMath::Sqrt(n.GetX()*n.GetX()+n.GetY()*n.GetY()))>0.01)return kFALSE;
1194 
1195  // Calculate z position
1196 
1197  convpos[2] = (p.GetZ()*negradius+n.GetZ()*posradius)/(negradius+posradius);
1198 
1199  // Calculate DCA
1200  TVector2 vdca=vertexPos-vertexNeg;
1201  dca[0]=vdca.Mod();
1202  dca[1]=TMath::Abs(n.GetZ()-p.GetZ());
1203 
1204  return kTRUE;
1205 }
1206 
1207 //________________________________________________________________________
1208 Bool_t AliV0ReaderV1::GetAODConversionGammas(){
1209 
1210  // Get reconstructed conversion photons from satellite AOD file
1211 
1212  AliAODEvent *fAODEvent=dynamic_cast<AliAODEvent*>(fInputEvent);
1213 
1214  if(fAODEvent){
1215 
1216  if(fConversionGammas == NULL){
1217  fConversionGammas = new TClonesArray("AliAODConversionPhoton",100);
1218  }
1219  fConversionGammas->Delete();//Reset the TClonesArray
1220 
1221  //Get Gammas from satellite AOD gamma branch
1222 
1223  AliAODConversionPhoton *gamma=0x0;
1224 
1225  TClonesArray *fInputGammas=dynamic_cast<TClonesArray*>(fAODEvent->FindListObject(fDeltaAODBranchName.Data()));
1226 
1227  if(!fInputGammas){
1228  FindDeltaAODBranchName();
1229  fInputGammas=dynamic_cast<TClonesArray*>(fAODEvent->FindListObject(fDeltaAODBranchName.Data()));}
1230  if(!fInputGammas){AliError("No Gamma Satellites found");return kFALSE;}
1231  // Apply Selection Cuts to Gammas and create local working copy
1232  if(fInputGammas){
1233  for(Int_t i=0;i<fInputGammas->GetEntriesFast();i++){
1234  gamma=dynamic_cast<AliAODConversionPhoton*>(fInputGammas->At(i));
1235  if(gamma){
1236  if(fRelabelAODs)RelabelAODPhotonCandidates(gamma);
1237  if(fConversionCuts->PhotonIsSelected(gamma,fInputEvent)){
1238  new((*fConversionGammas)[fConversionGammas->GetEntriesFast()]) AliAODConversionPhoton(*gamma);}
1239  }
1240  }
1241  }
1242  }
1243 
1244  if(fConversionGammas->GetEntries()){return kTRUE;}
1245 
1246  return kFALSE;
1247 }
1248 
1249 //________________________________________________________________________
1250 void AliV0ReaderV1::FindDeltaAODBranchName(){
1251 
1252  // Find delta AOD branchname containing reconstructed photons
1253 
1254  TList *list=fInputEvent->GetList();
1255  for(Int_t ii=0;ii<list->GetEntries();ii++){
1256  TString name((list->At(ii))->GetName());
1257  if(name.BeginsWith(fDeltaAODBranchName)&&name.EndsWith("gamma")){
1258  fDeltaAODBranchName=name;
1259  AliInfo(Form("Set DeltaAOD BranchName to: %s",fDeltaAODBranchName.Data()));
1260  return;
1261  }
1262  }
1263 }
1264 
1265 //________________________________________________________________________
1266 void AliV0ReaderV1::RelabelAODPhotonCandidates(AliAODConversionPhoton *PhotonCandidate){
1267 
1268  if(fPreviousV0ReaderPerformsAODRelabeling == 2) return;
1269  else if(fPreviousV0ReaderPerformsAODRelabeling == 0){
1270  printf("Running AODs! Determine if V0Reader '%s' should perform relabeling\n",this->GetName());
1271  TObjArray* obj = (TObjArray*)AliAnalysisManager::GetAnalysisManager()->GetTasks();
1272  Int_t iPosition = obj->IndexOf(this);
1273  Bool_t prevV0ReaderRunningButNotRelabeling = kFALSE;
1274  for(Int_t i=iPosition-1; i>=0; i--){
1275  if( (obj->At(i))->IsA() == AliV0ReaderV1::Class()){
1276  AliV0ReaderV1* tempReader = (AliV0ReaderV1*) obj->At(i);
1277  if( tempReader->AreAODsRelabeled() && tempReader->IsReaderPerformingRelabeling() == 1){
1278  fPreviousV0ReaderPerformsAODRelabeling = 2;
1279  prevV0ReaderRunningButNotRelabeling = kFALSE;
1280  printf("V0Reader '%s' is running before this V0Reader '%s': do _NOT_ relabel AODs by current reader!\n",tempReader->GetName(),this->GetName());
1281  break;
1282  }else prevV0ReaderRunningButNotRelabeling = kTRUE;
1283  }
1284  }
1285  if(prevV0ReaderRunningButNotRelabeling) AliFatal(Form("There are V0Readers before '%s', but none of them is relabeling!",this->GetName()));
1286 
1287  if(fPreviousV0ReaderPerformsAODRelabeling == 2) return;
1288  else{
1289  printf("This V0Reader '%s' is first to be processed: do relabel AODs by current reader!\n",this->GetName());
1290  fPreviousV0ReaderPerformsAODRelabeling = 1;
1291  }
1292  }
1293 
1294  if(fPreviousV0ReaderPerformsAODRelabeling != 1) AliFatal(Form("In %s: fPreviousV0ReaderPerformsAODRelabeling = '%i' - while it should be impossible it is something different than '1'!",this->GetName(),fPreviousV0ReaderPerformsAODRelabeling));
1295 
1296  // Relabeling For AOD Event
1297  // ESDiD -> AODiD
1298  // MCLabel -> AODMCLabel
1299  Bool_t AODLabelPos = kFALSE;
1300  Bool_t AODLabelNeg = kFALSE;
1301 
1302  for(Int_t i = 0; i<fInputEvent->GetNumberOfTracks();i++){
1303  AliAODTrack *tempDaughter = static_cast<AliAODTrack*>(fInputEvent->GetTrack(i));
1304  if(!AODLabelPos){
1305  if( tempDaughter->GetID() == PhotonCandidate->GetTrackLabelPositive() ){
1306  PhotonCandidate->SetMCLabelPositive(TMath::Abs(tempDaughter->GetLabel()));
1307  PhotonCandidate->SetLabelPositive(i);
1308  AODLabelPos = kTRUE;
1309  }
1310  }
1311  if(!AODLabelNeg){
1312  if( tempDaughter->GetID() == PhotonCandidate->GetTrackLabelNegative()){
1313  PhotonCandidate->SetMCLabelNegative(TMath::Abs(tempDaughter->GetLabel()));
1314  PhotonCandidate->SetLabelNegative(i);
1315  AODLabelNeg = kTRUE;
1316  }
1317  }
1318  if(AODLabelNeg && AODLabelPos){
1319  return;
1320  }
1321  }
1322  if(!AODLabelPos || !AODLabelNeg){
1323  AliError(Form("NO AOD Daughters Found Pos: %i %i Neg: %i %i, setting all labels to -999999",AODLabelPos,PhotonCandidate->GetTrackLabelPositive(),AODLabelNeg,PhotonCandidate->GetTrackLabelNegative()));
1324  if(!AODLabelNeg){
1325  PhotonCandidate->SetMCLabelNegative(-999999);
1326  PhotonCandidate->SetLabelNegative(-999999);
1327  }
1328  if(!AODLabelPos){
1329  PhotonCandidate->SetMCLabelPositive(-999999);
1330  PhotonCandidate->SetLabelPositive(-999999);
1331  }
1332  }
1333 
1334 }
1335 
1336 //************************************************************************
1337 // This function counts the number of primary tracks in the event, the
1338 // implementation for ESD and AOD had to be different due to the different
1339 // filters which are already applied on AOD level the so-called filter
1340 // bits, we tried to replicate the conditions in both but an exact match
1341 // could not be reached. The cuts are similar to the ones used by the
1342 // jet-group
1343 //************************************************************************
1344 //________________________________________________________________________
1345 void AliV0ReaderV1::CountTracks(){
1346 
1347  // In Case of MC count only MB tracks
1348  // if(fMCEvent && fInputEvent->IsA()==AliESDEvent::Class()){
1349  // fEventCuts->GetNotRejectedParticles(1,NULL,fMCEvent);
1350  // }
1351  // if(fMCEvent && fInputEvent->IsA()==AliAODEvent::Class()){
1352  // fEventCuts->GetNotRejectedParticles(1,NULL,fInputEvent);
1353  // }
1354 
1355  if(fInputEvent->IsA()==AliESDEvent::Class()){
1356  static AliESDtrackCuts *EsdTrackCuts = 0x0;
1357  static int prevRun = -1;
1358  // Using standard function for setting Cuts
1359  Int_t runNumber = fInputEvent->GetRunNumber();
1360 
1361  if (prevRun!=runNumber) {
1362  delete EsdTrackCuts;
1363  EsdTrackCuts = 0;
1364  prevRun = runNumber;
1365  }
1366  if (!EsdTrackCuts) {
1367  // if LHC11a or earlier or if LHC13g or if LHC12a-i -> use 2010 cuts
1368  if( (runNumber<=146860) || (runNumber>=197470 && runNumber<=197692) || (runNumber>=172440 && runNumber<=193766) ){
1369  EsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010();
1370  // else if run2 data use 2015 PbPb cuts
1371  } else if (runNumber>=209122){
1372  //EsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2015PbPb();
1373  // hard coded track cuts for the moment, because AliESDtrackCuts::GetStandardITSTPCTrackCuts2015PbPb() gives spams warnings
1374  EsdTrackCuts = new AliESDtrackCuts();
1375  // TPC; clusterCut = 1, cutAcceptanceEdges = kTRUE, removeDistortedRegions = kFALSE
1376  EsdTrackCuts->AliESDtrackCuts::SetMinNCrossedRowsTPC(70);
1377  EsdTrackCuts->AliESDtrackCuts::SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);
1378  EsdTrackCuts->SetCutGeoNcrNcl(2., 130., 1.5, 0.0, 0.0); // only dead zone and not clusters per length
1379  //EsdTrackCuts->AliESDtrackCuts::SetCutOutDistortedRegionsTPC(kTRUE);
1380  EsdTrackCuts->AliESDtrackCuts::SetMaxChi2PerClusterTPC(4);
1381  EsdTrackCuts->AliESDtrackCuts::SetAcceptKinkDaughters(kFALSE);
1382  EsdTrackCuts->AliESDtrackCuts::SetRequireTPCRefit(kTRUE);
1383  // ITS; selPrimaries = 1
1384  EsdTrackCuts->AliESDtrackCuts::SetRequireITSRefit(kTRUE);
1385  EsdTrackCuts->AliESDtrackCuts::SetClusterRequirementITS(AliESDtrackCuts::kSPD,
1386  AliESDtrackCuts::kAny);
1387  EsdTrackCuts->AliESDtrackCuts::SetMaxDCAToVertexXYPtDep("0.0105+0.0350/pt^1.1");
1388  EsdTrackCuts->AliESDtrackCuts::SetMaxChi2TPCConstrainedGlobal(36);
1389  EsdTrackCuts->AliESDtrackCuts::SetMaxDCAToVertexZ(2);
1390  EsdTrackCuts->AliESDtrackCuts::SetDCAToVertex2D(kFALSE);
1391  EsdTrackCuts->AliESDtrackCuts::SetRequireSigmaToVertex(kFALSE);
1392  EsdTrackCuts->AliESDtrackCuts::SetMaxChi2PerClusterITS(36);
1393  // else use 2011 version of track cuts
1394  }else{
1395  EsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011();
1396  }
1397  EsdTrackCuts->SetMaxDCAToVertexZ(2);
1398  EsdTrackCuts->SetEtaRange(-0.8, 0.8);
1399  EsdTrackCuts->SetPtRange(0.15);
1400  }
1401  fNumberOfPrimaryTracks = 0;
1402  for(Int_t iTracks = 0; iTracks < fInputEvent->GetNumberOfTracks(); iTracks++){
1403  AliESDtrack* curTrack = (AliESDtrack*) fInputEvent->GetTrack(iTracks);
1404  if(!curTrack) continue;
1405  if(!EsdTrackCuts->AcceptTrack(curTrack)) continue;
1406  //if(fMCEvent && !(fEventCuts->IsParticleFromBGEvent(TMath::Abs(curTrack->GetLabel()),fMCEvent->Stack(),fInputEvent))) continue;
1407  fNumberOfPrimaryTracks++;
1408  }
1409  }
1410  else if(fInputEvent->IsA()==AliAODEvent::Class()){
1411  fNumberOfPrimaryTracks = 0;
1412  for(Int_t iTracks = 0; iTracks<fInputEvent->GetNumberOfTracks(); iTracks++){
1413  AliAODTrack* curTrack = (AliAODTrack*) fInputEvent->GetTrack(iTracks);
1414  if(curTrack->GetID()<0) continue; // Avoid double counting of tracks
1415  if(!curTrack->IsHybridGlobalConstrainedGlobal()) continue;
1416  if(TMath::Abs(curTrack->Eta())>0.8) continue;
1417  if(curTrack->Pt()<0.15) continue;
1418  //if(fMCEvent && !(fEventCuts->IsParticleFromBGEvent(TMath::Abs(curTrack->GetLabel()),NULL,fInputEvent))) continue;
1419  //if(TMath::Abs(curTrack->ZAtDCA())>2) continue; // Only Set For TPCOnly tracks
1420  fNumberOfPrimaryTracks++;
1421  }
1422  }
1423 
1424  return;
1425 }
1426 
1428 void AliV0ReaderV1::CountTPCoutTracks(){
1429  fNumberOfTPCoutTracks = 0;
1430 
1431  for (Int_t itrk = 0; itrk < fInputEvent->GetNumberOfTracks(); itrk++) {
1432  AliVTrack *trk = dynamic_cast<AliVTrack*>(fInputEvent->GetTrack(itrk));
1433 
1434  if (trk != NULL) {
1435  /* the initial method of counting TPC out tracks */
1436  if (!(trk->Pt() < 0.15) && (TMath::Abs(trk->Eta()) < 0.8)) {
1437  if ((trk->GetStatus() & AliVTrack::kTPCout) == AliVTrack::kTPCout) {
1438  fNumberOfTPCoutTracks++;
1439  }
1440  }
1441  }
1442  }
1443  return;
1444 }
1445 
1446 
1448 Bool_t AliV0ReaderV1::ParticleIsConvertedPhoton(AliMCEvent *mcEvent, TParticle *particle, Double_t etaMax, Double_t rMax, Double_t zMax){
1449  // MonteCarlo Photon Selection
1450  if(!mcEvent)return kFALSE;
1451 
1452  if (particle->GetPdgCode() == 22){
1453  // check whether particle is within eta range
1454  if( TMath::Abs(particle->Eta()) > etaMax ) return kFALSE;
1455  // check if particle doesn't have a photon as mother
1456  if(particle->GetMother(0) >-1 && mcEvent->Particle(particle->GetMother(0))->GetPdgCode() == 22){
1457  return kFALSE; // no photon as mothers!
1458  }
1459  // looking for conversion gammas (electron + positron from pairbuilding (= 5) )
1460  TParticle* ePos = NULL;
1461  TParticle* eNeg = NULL;
1462  if(particle->GetNDaughters() >= 2){
1463  for(Int_t daughterIndex=particle->GetFirstDaughter();daughterIndex<=particle->GetLastDaughter();daughterIndex++){
1464  if(daughterIndex<0) continue;
1465  TParticle *tmpDaughter = mcEvent->Particle(daughterIndex);
1466  if(tmpDaughter->GetUniqueID() == 5){
1467  if(tmpDaughter->GetPdgCode() == 11){
1468  eNeg = tmpDaughter;
1469  } else if(tmpDaughter->GetPdgCode() == -11){
1470  ePos = tmpDaughter;
1471  }
1472  }
1473  }
1474  }
1475  if(ePos == NULL || eNeg == NULL){ // means we do not have two daughters from pair production
1476  return kFALSE;
1477  }
1478  // check if electrons are in correct eta window
1479  if( TMath::Abs(ePos->Eta()) > etaMax ||
1480  TMath::Abs(eNeg->Eta()) > etaMax )
1481  return kFALSE;
1482 
1483  // check if photons have converted in reconstructable range
1484  if(ePos->R() > rMax){
1485  return kFALSE; // cuts on distance from collision point
1486  }
1487  if(TMath::Abs(ePos->Vz()) > zMax){
1488  return kFALSE; // outside material
1489  }
1490  if(TMath::Abs(eNeg->Vz()) > zMax){
1491  return kFALSE; // outside material
1492  }
1493 
1494 
1495  Double_t lineCutZRSlope = tan(2*atan(exp(-etaMax)));
1496  Double_t lineCutZValue = 7.;
1497  if( ePos->R() <= ((TMath::Abs(ePos->Vz()) * lineCutZRSlope) - lineCutZValue)){
1498  return kFALSE; // line cut to exclude regions where we do not reconstruct
1499  }
1500  if( eNeg->R() <= ((TMath::Abs(eNeg->Vz()) * lineCutZRSlope) - lineCutZValue)){
1501  return kFALSE; // line cut to exclude regions where we do not reconstruct
1502  }
1503  if (ePos->Pt() < 0.05 || eNeg->Pt() < 0.05){
1504  return kFALSE;
1505  }
1506 
1507  return kTRUE;
1508  }
1509  return kFALSE;
1510 }
1511 
1513 void AliV0ReaderV1::CreatePureMCHistosForV0FinderEffiESD(){
1514  const AliVVertex* primVtxMC = fMCEvent->GetPrimaryVertex();
1515  Double_t mcProdVtxX = primVtxMC->GetX();
1516  Double_t mcProdVtxY = primVtxMC->GetY();
1517  Double_t mcProdVtxZ = primVtxMC->GetZ();
1518  //cout << mcProdVtxX <<"\t" << mcProdVtxY << "\t" << mcProdVtxZ << endl;
1519 
1520  // Loop over all primary MC particle
1521  for(Long_t i = 0; i < fMCEvent->GetNumberOfTracks(); i++) {
1522  if (fEventCuts->IsConversionPrimaryESD( fMCEvent, i, mcProdVtxX, mcProdVtxY, mcProdVtxZ)){
1523  // fill primary histogram
1524  TParticle* particle = (TParticle *)fMCEvent->Particle(i);
1525  if (!particle) continue;
1526  if (ParticleIsConvertedPhoton(fMCEvent, particle, 0.9, 180.,250. )){
1527  if(particle->GetFirstDaughter()<0) continue;
1528  TParticle *tmpDaughter = fMCEvent->Particle(particle->GetFirstDaughter());
1529  if (!tmpDaughter) continue;
1530  fHistoMCGammaPtvsR->Fill(particle->Pt(),tmpDaughter->R());
1531  fHistoMCGammaPtvsPhi->Fill(particle->Pt(),particle->Phi());
1532  fHistoMCGammaRvsPhi->Fill(tmpDaughter->R(),particle->Phi());
1533  }
1534  if (ParticleIsConvertedPhoton(fMCEvent, particle, 1.4, 180.,250. )){
1535  if(particle->GetFirstDaughter()<0) continue;
1536  TParticle *tmpDaughter = fMCEvent->Particle(particle->GetFirstDaughter());
1537  if (!tmpDaughter) continue;
1538  fHistoMCGammaPtvsEta->Fill(particle->Pt(),particle->Eta());
1539  fHistoMCGammaRvsEta->Fill(tmpDaughter->R(),particle->Eta());
1540  fHistoMCGammaPhivsEta->Fill(particle->Phi(),particle->Eta());
1541  }
1542  }
1543  }
1544 }
1545 
1547 void AliV0ReaderV1::FillRecMCHistosForV0FinderEffiESD( AliESDv0* currentV0){
1548 
1549  const AliVVertex* primVtxMC = fMCEvent->GetPrimaryVertex();
1550  Double_t mcProdVtxX = primVtxMC->GetX();
1551  Double_t mcProdVtxY = primVtxMC->GetY();
1552  Double_t mcProdVtxZ = primVtxMC->GetZ();
1553 // cout << mcProdVtxX <<"\t" << mcProdVtxY << "\t" << mcProdVtxZ << endl;
1554 
1555  Int_t tracklabelPos=currentV0->GetPindex();
1556  Int_t tracklabelNeg=currentV0->GetNindex();
1557 
1558  Int_t labelp=TMath::Abs(fConversionCuts->GetTrack(fInputEvent,tracklabelPos)->GetLabel());
1559  Int_t labeln=TMath::Abs(fConversionCuts->GetTrack(fInputEvent,tracklabelNeg)->GetLabel());
1560 
1561  TParticle* negPart = 0x0;
1562  if(labeln>-1) negPart = (TParticle *)fMCEvent->Particle(labeln);
1563  TParticle* posPart = 0x0;
1564  if(labelp>-1) posPart = (TParticle *)fMCEvent->Particle(labelp);
1565 
1566  if ( negPart == NULL || posPart == NULL ) return;
1567 // if (!(negPart->GetPdgCode() == 11)) return;
1568 // if (!(posPart->GetPdgCode() == -11)) return;
1569  Long_t motherlabelNeg = negPart->GetFirstMother();
1570  Long_t motherlabelPos = posPart->GetFirstMother();
1571 
1572 // cout << "mother neg " << motherlabelNeg << " mother pos " << motherlabelPos << endl;
1573  if (motherlabelNeg>-1 && motherlabelNeg == motherlabelPos && negPart->GetFirstMother() != -1){
1574  if (fEventCuts->IsConversionPrimaryESD( fMCEvent, negPart->GetFirstMother(), mcProdVtxX, mcProdVtxY, mcProdVtxZ)){
1575 
1576  TParticle* mother = (TParticle *)fMCEvent->Particle(motherlabelNeg);
1577  if (mother->GetPdgCode() == 22 ){
1578  if (!CheckVectorForDoubleCount(fVectorFoundGammas,motherlabelNeg ) ){
1579  if (ParticleIsConvertedPhoton(fMCEvent, mother, 0.9, 180.,250. )){
1580  fHistoRecMCGammaPtvsR->Fill(mother->Pt(),negPart->R());
1581  fHistoRecMCGammaPtvsPhi->Fill(mother->Pt(),mother->Phi());
1582  fHistoRecMCGammaRvsPhi->Fill(negPart->R(),mother->Phi());
1583  }
1584  if (ParticleIsConvertedPhoton(fMCEvent, mother, 1.4, 180.,250. )){
1585  fHistoRecMCGammaPtvsEta->Fill(mother->Pt(),mother->Eta());
1586  fHistoRecMCGammaRvsEta->Fill(negPart->R(),mother->Eta());
1587  fHistoRecMCGammaPhivsEta->Fill(mother->Phi(),mother->Eta());
1588  }
1589 // cout << "new gamma found" << endl;
1590  } else {
1591  if (ParticleIsConvertedPhoton(fMCEvent, mother, 0.9, 180.,250. )){
1592  fHistoRecMCGammaMultiPt->Fill(mother->Pt());
1593  fHistoRecMCGammaMultiPhi->Fill(mother->Phi());
1594  fHistoRecMCGammaMultiR->Fill(negPart->R());
1595  }
1596  if (ParticleIsConvertedPhoton(fMCEvent, mother, 1.4, 180.,250. )){
1597  fHistoRecMCGammaMultiPtvsEta->Fill(mother->Pt(),mother->Eta());
1598  }
1599 // cout << "this one I had already: " << motherlabelNeg << endl << "-----------------------" << endl;
1600  }
1601 // cout << "event gammas: " << endl;
1602 // for(Int_t iGamma=0; iGamma<fVectorFoundGammas.size(); iGamma++){cout << fVectorFoundGammas.at(iGamma) << ", ";}
1603  }
1604  }
1605  }
1606 }
1607 
1608 //_________________________________________________________________________________
1609 void AliV0ReaderV1::FillImpactParamHistograms( AliVTrack* pTrack, AliVTrack* nTrack, AliESDv0 *fCurrentV0, AliKFConversionPhoton *fCurrentMotherKF){
1610 
1611  // values of cuts to be introduced in 2016 to reduce ESD size
1612  Float_t fZmax = 25; //cm
1613  Float_t fYmax = 25; //cm
1614  Double_t kTPCMargin = 1.0; //cm
1615  Int_t kMaxTPCV0Conflicts = 1; //# conflicting clusters tolerated
1616 
1617  //conversion point
1618  Double_t convX, convY, convZ;
1619  fCurrentV0->GetXYZ(convX,convY,convZ);
1620  Double_t convR = TMath::Sqrt(convX*convX+convY*convY);
1621  //recalculated conversion point
1622  Double_t convXrecalc=fCurrentMotherKF->GetConversionX();
1623  Double_t convYrecalc=fCurrentMotherKF->GetConversionY();
1624  Double_t convZrecalc=fCurrentMotherKF->GetConversionZ();
1625  Double_t convRrecalc = fCurrentMotherKF->GetConversionRadius();
1626 
1627  //count V0s
1628  fHistoImpactParameterStudy->Fill(0);
1629 
1630  //count V0s with two TPC-only tracks
1631  AliESDtrack* positiveTrack = (AliESDtrack*) pTrack;
1632  AliESDtrack* negativeTrack = (AliESDtrack*) nTrack;
1633  Bool_t TPConly = kFALSE;
1634  if (!positiveTrack->IsOn(AliESDtrack::kITSin) && !negativeTrack->IsOn(AliESDtrack::kITSin)){
1635  fHistoImpactParameterStudy->Fill(1);
1636  TPConly = kTRUE;
1637  }
1638 
1639  Bool_t RemovedByZcut = kFALSE;
1640  Bool_t RemovedByYcut = kFALSE;
1641 
1642  //count V0s which have...
1643  if(TPConly){
1644  //not passed z cut: pos.tr. or neg.tr.
1645  if(((TMath::Abs(positiveTrack->GetZ()))>fZmax) || ((TMath::Abs(negativeTrack->GetZ()))>fZmax)){
1646  fHistoImpactParameterStudy->Fill(2);
1647  RemovedByZcut=kTRUE;
1648  }
1649 
1650  //not passed y cut: pos.tr. or neg.tr.
1651  if(((TMath::Abs(positiveTrack->GetY()))>fYmax) || ((TMath::Abs(negativeTrack->GetY()))>fYmax)){
1652  fHistoImpactParameterStudy->Fill(3);
1653  RemovedByYcut=kTRUE;
1654  }
1655  }
1656 
1657  //causality cut
1658  Bool_t RemovedByCausality=kFALSE;
1659  AliESDEvent *fESDEvent=dynamic_cast<AliESDEvent*>(fInputEvent);
1660  const Float_t rTPC[159]={
1661  85.225, 85.975, 86.725, 87.475, 88.225, 88.975, 89.725, 90.475, 91.225, 91.975, 92.725, 93.475, 94.225, 94.975, 95.725,
1662  96.475, 97.225, 97.975, 98.725, 99.475,100.225,100.975,101.725,102.475,103.225,103.975,104.725,105.475,106.225,106.975,
1663  107.725,108.475,109.225,109.975,110.725,111.475,112.225,112.975,113.725,114.475,115.225,115.975,116.725,117.475,118.225,
1664  118.975,119.725,120.475,121.225,121.975,122.725,123.475,124.225,124.975,125.725,126.475,127.225,127.975,128.725,129.475,
1665  130.225,130.975,131.725,135.100,136.100,137.100,138.100,139.100,140.100,141.100,142.100,143.100,144.100,145.100,146.100,
1666  147.100,148.100,149.100,150.100,151.100,152.100,153.100,154.100,155.100,156.100,157.100,158.100,159.100,160.100,161.100,
1667  162.100,163.100,164.100,165.100,166.100,167.100,168.100,169.100,170.100,171.100,172.100,173.100,174.100,175.100,176.100,
1668  177.100,178.100,179.100,180.100,181.100,182.100,183.100,184.100,185.100,186.100,187.100,188.100,189.100,190.100,191.100,
1669  192.100,193.100,194.100,195.100,196.100,197.100,198.100,199.350,200.850,202.350,203.850,205.350,206.850,208.350,209.850,
1670  211.350,212.850,214.350,215.850,217.350,218.850,220.350,221.850,223.350,224.850,226.350,227.850,229.350,230.850,232.350,
1671  233.850,235.350,236.850,238.350,239.850,241.350,242.850,244.350,245.850};
1672 
1673  // fill conversion radius histograms
1674  fHistoR->Fill(convR);
1675  fHistoRrecalc->Fill(convRrecalc);
1676  Double_t alpha = TMath::ATan2(convY,convX);
1677  if (alpha<0) alpha += TMath::Pi()*2;
1678  Int_t sec = alpha/(TMath::Pi()/9);
1679  alpha = (10.+sec*20.)*TMath::DegToRad();
1680  Double_t cs = TMath::Cos(alpha);
1681  Double_t sn = TMath::Sin(alpha);
1682  Double_t xsV0 = convX*cs - convY*sn;
1683  fHistoRviaAlpha->Fill(xsV0);
1684  Double_t alpha_r = TMath::ATan2(convYrecalc,convXrecalc);
1685  if (alpha_r<0) alpha_r += TMath::Pi()*2;
1686  Int_t sec_r = alpha_r/(TMath::Pi()/9);
1687  alpha_r = (10.+sec_r*20.)*TMath::DegToRad();
1688  Double_t cs_r = TMath::Cos(alpha_r);
1689  Double_t sn_r = TMath::Sin(alpha_r);
1690  Double_t xsV0_r = convXrecalc*cs_r - convYrecalc*sn_r;
1691  fHistoRviaAlphaRecalc->Fill(xsV0_r);
1692 
1693  if (convR > 80) { // conversion within TPC <-> TPC-only tracks
1694  fHistoImpactParameterStudy->Fill(4);
1695 
1696  for (Int_t it=2;it--;) {
1697  Int_t trId = fCurrentV0->GetIndex(it);
1698  AliESDtrack* tr = fESDEvent->GetTrack(trId);
1699  const TBits& bits = tr->GetTPCClusterMap();
1700  Int_t nConflict = 0;
1701  for (Int_t ic=0;ic<159;ic++) {
1702  if (rTPC[ic]>(xsV0-kTPCMargin)) break;
1703  if (bits.TestBitNumber(ic)){
1704  nConflict++;
1705  fHistoRdiff->Fill(xsV0-rTPC[ic]);
1706  }
1707  if (nConflict>kMaxTPCV0Conflicts) {
1708  fHistoImpactParameterStudy->Fill(5);
1709  RemovedByCausality=kTRUE;
1710  break;
1711  }
1712  }
1713  }
1714  }
1715 
1716  //not passed y or z o causality cut:
1717  Bool_t RemovedByAnyCut=kFALSE;
1718  if(RemovedByZcut||RemovedByYcut||RemovedByCausality){
1719  fHistoImpactParameterStudy->Fill(6);
1720  RemovedByAnyCut=kTRUE;
1721  }
1722 
1723  //Fill tree for further analysis
1724  Float_t posZ;
1725  Float_t posY;
1726  Float_t posX;
1727  Float_t posPt;
1728  Float_t negZ;
1729  Float_t negY;
1730  Float_t negX;
1731  Float_t negPt;
1732  Float_t R;
1733  fImpactParamTree->Branch("posPt",&posPt,"posPt/F");
1734  fImpactParamTree->Branch("posY",&posY,"posY/F");
1735  fImpactParamTree->Branch("R",&R,"R/F");
1736  posZ = positiveTrack->GetZ();
1737  posY = positiveTrack->GetY();
1738  posX = positiveTrack->GetX();
1739  posPt = positiveTrack->Pt();
1740  negZ = negativeTrack->GetZ();
1741  negY = negativeTrack->GetY();
1742  negX = negativeTrack->GetX();
1743  negPt = negativeTrack->Pt();
1744  R=convRrecalc;
1745  fImpactParamTree->Fill();
1746 
1747  // fill impact parameter histograms
1748  fHistoPosTrackImpactParamZ->Fill(posZ);
1749  fHistoPosTrackImpactParamY->Fill(posY);
1750  fHistoPosTrackImpactParamX->Fill(posX);
1751  fHistoPosTrackImpactParamZvsPt->Fill(posPt, posZ);
1752  fHistoPosTrackImpactParamYvsPt->Fill(posPt, posY);
1753  fHistoPosTrackImpactParamXvsPt->Fill(posPt, posX);
1754  fHistoNegTrackImpactParamZ->Fill(negZ);
1755  fHistoNegTrackImpactParamY->Fill(negY);
1756  fHistoNegTrackImpactParamX->Fill(negX);
1757  fHistoNegTrackImpactParamZvsPt->Fill(negPt, negZ);
1758  fHistoNegTrackImpactParamYvsPt->Fill(negPt, negY);
1759  fHistoNegTrackImpactParamXvsPt->Fill(negPt, negX);
1760 
1761  // fill comparison histos before / after new cuts
1762  fHistoPt->Fill(positiveTrack->Pt());
1763  fHistoImpactParamZvsR->Fill(convZrecalc,convRrecalc);
1764  //Float_t *DCAzPhoton;
1765  Float_t DCAzPhoton;
1766  //const AliVVertex *PrimVertex = fInputEvent->GetPrimaryVertex();
1767  AliAODConversionPhoton* fCurrentMother=(AliAODConversionPhoton*)fCurrentMotherKF;
1768  //AliAODConversionPhoton *fCurrentMother=dynamic_cast<AliAODConversionPhoton*>(fCurrentMotherKF);
1769  //fCurrentMotherKF->GetDistanceOfClossetApproachToPrimVtx(PrimVertex, DCAzPhoton);
1770  DCAzPhoton = fCurrentMother->GetDCAzToPrimVtx();
1771  fHistoDCAzPhoton->Fill(DCAzPhoton);
1772  if(!RemovedByAnyCut) {
1773  fHistoPt2->Fill(positiveTrack->Pt());
1774  fHistoImpactParamZvsR2->Fill(convZrecalc,convRrecalc);
1775  fHistoDCAzPhoton2->Fill(DCAzPhoton);
1776  }
1777 
1778  return;
1779 }
1780 
1781 //_________________________________________________________________________________
1782 Bool_t AliV0ReaderV1::CheckVectorOnly(vector<Int_t> &vec, Int_t tobechecked)
1783 {
1784  if(tobechecked > -1)
1785  {
1786  vector<Int_t>::iterator it;
1787  it = find (vec.begin(), vec.end(), tobechecked);
1788  if (it != vec.end()) return true;
1789  else return false;
1790  }
1791  return false;
1792 }
1793 
1794 //_________________________________________________________________________________
1795 Bool_t AliV0ReaderV1::CheckVectorForDoubleCount(vector<Int_t> &vec, Int_t tobechecked)
1796 {
1797  if(tobechecked > -1)
1798  {
1799  vector<Int_t>::iterator it;
1800  it = find (vec.begin(), vec.end(), tobechecked);
1801  if (it != vec.end()) return true;
1802  else{
1803  vec.push_back(tobechecked);
1804  return false;
1805  }
1806  }
1807  return false;
1808 }
1809 
1810 //________________________________________________________________________
1811 void AliV0ReaderV1::Terminate(Option_t *)
1812 {
1813 
1814 }
1815 
1816 AliV0ReaderV1::iterator::iterator(const AliV0ReaderV1 *reader, AliV0ReaderV1::iterator::Direction_t dir, int position):
1817  std::iterator<std::bidirectional_iterator_tag, AliConversionPhotonBase>(),
1818  fkData(reader),
1819  fCurrentIndex(position),
1820  fDirection(dir)
1821 {
1822 }
1823 
1824 AliV0ReaderV1::iterator::iterator(const AliV0ReaderV1::iterator &other):
1825  std::iterator<std::bidirectional_iterator_tag, AliConversionPhotonBase>(other),
1826  fkData(other.fkData),
1827  fCurrentIndex(other.fCurrentIndex),
1828  fDirection(other.fDirection)
1829 {
1830 }
1831 
1832 AliV0ReaderV1::iterator &AliV0ReaderV1::iterator::operator=(const AliV0ReaderV1::iterator &other){
1833  if(this != &other){
1834  fkData = other.fkData;
1835  fCurrentIndex = other.fCurrentIndex;
1836  fDirection = other.fDirection;
1837  }
1838  return *this;
1839 }
1840 
1841 bool AliV0ReaderV1::iterator::operator!=(AliV0ReaderV1::iterator &other) const{
1842  if(fkData != other.fkData) return true;
1843  return fCurrentIndex != other.fCurrentIndex;
1844 }
1845 
1846 AliV0ReaderV1::iterator AliV0ReaderV1::iterator::operator++(int index){
1847  iterator tmp(*this);
1848  operator++();
1849  return tmp;
1850 }
1851 
1852 AliV0ReaderV1::iterator &AliV0ReaderV1::iterator::operator++(){
1853  if(fDirection == AliV0ReaderV1::iterator::kForwardDirection){
1854  fCurrentIndex++;
1855  } else {
1856  fCurrentIndex--;
1857  }
1858  return *this;
1859 }
1860 
1861 AliV0ReaderV1::iterator AliV0ReaderV1::iterator::operator--(int){
1862  iterator tmp(*this);
1863  operator--();
1864  return tmp;
1865 }
1866 
1867 AliV0ReaderV1::iterator &AliV0ReaderV1::iterator::operator--(){
1868  if(fDirection == AliV0ReaderV1::iterator::kForwardDirection){
1869  fCurrentIndex++;
1870  } else {
1871  fCurrentIndex--;
1872  }
1873  return *this;
1874 }
1875 
1876 AliConversionPhotonBase *AliV0ReaderV1::iterator::operator*(){
1877  return (*fkData)[fCurrentIndex];
1878 }
AliV0ReaderV1::IsReaderPerformingRelabeling
Int_t IsReaderPerformingRelabeling()
Definition: AliV0ReaderV1.h:122
charge
Int_t charge
Definition: MuonTrackingEfficiency.C:54
AliV0ReaderV1::iterator::fDirection
Direction_t fDirection
Iterator in forward direction.
Definition: AliV0ReaderV1.h:69
AliV0ReaderV1::FillAODOutput
void FillAODOutput()
Definition: AliV0ReaderV1.cxx:659
Double_t
double Double_t
Definition: External.C:58
AliV0ReaderV1::ParticleIsConvertedPhoton
Bool_t ParticleIsConvertedPhoton(AliMCEvent *mcEvent, TParticle *particle, Double_t etaMax, Double_t rMax, Double_t zMax)
Definition: AliV0ReaderV1.cxx:1448
AliConvEventCuts::kLHC16c3b
anchored LHC13[d-e] pass 2 - JJ
Definition: AliConvEventCuts.h:156
TH2F
Definition: External.C:236
AliV0ReaderV1::FindDeltaAODBranchName
void FindDeltaAODBranchName()
Definition: AliV0ReaderV1.cxx:1250
AliConversionPhotonBase::SetConversionPoint
void SetConversionPoint(Double_t convpoint[3])
Definition: AliConversionPhotonBase.h:89
AliV0ReaderV1::CheckVectorForDoubleCount
Bool_t CheckVectorForDoubleCount(vector< Int_t > &vec, Int_t tobechecked)
Definition: AliV0ReaderV1.cxx:1795
AliConversionPhotonBase::SetMCLabelNegative
void SetMCLabelNegative(Int_t label)
Definition: AliConversionPhotonBase.h:66
AliV0ReaderV1::AreAODsRelabeled
Bool_t AreAODsRelabeled()
Definition: AliV0ReaderV1.h:121
AliV0ReaderV1::iterator::fkData
const AliV0ReaderV1 * fkData
V0 reader used to iterate over.
Definition: AliV0ReaderV1.h:67
AliV0ReaderV1::CreatePureMCHistosForV0FinderEffiESD
void CreatePureMCHistosForV0FinderEffiESD()
Definition: AliV0ReaderV1.cxx:1513
mass
Double_t mass
Definition: charmCutsOptimization.C:35
AliV0ReaderV1::iterator::operator*
AliConversionPhotonBase * operator*()
Definition: AliV0ReaderV1.cxx:1876
AliAODConversionPhoton::GetDCAzToPrimVtx
Float_t GetDCAzToPrimVtx() const
Definition: AliAODConversionPhoton.h:56
AliV0ReaderV1::Init
virtual void Init()
Definition: AliV0ReaderV1.cxx:235
AliV0ReaderV1::GetPsiPair
Double_t GetPsiPair(const AliESDv0 *v0, const AliExternalTrackParam *positiveparam, const AliExternalTrackParam *negativeparam, const Double_t convpos[3]) const
Definition: AliV0ReaderV1.cxx:942
AliV0ReaderV1::FillRecMCHistosForV0FinderEffiESD
void FillRecMCHistosForV0FinderEffiESD(AliESDv0 *currentV0)
Definition: AliV0ReaderV1.cxx:1547
AliV0ReaderV1::CheckVectorOnly
Bool_t CheckVectorOnly(vector< Int_t > &vec, Int_t tobechecked)
Definition: AliV0ReaderV1.cxx:1782
AliConvEventCuts::kLHC16P1JJLowB
anchored LHC16f pass 1 low B-field - Pythia8 JJ
Definition: AliConvEventCuts.h:232
AliAODConversionPhoton.h
AliConvEventCuts::kLHC15a3a
anchored LHC13g pass 1 - JJ
Definition: AliConvEventCuts.h:159
AliConvEventCuts::kLHC15a3a_plus
anchored LHC13g pass 1 - JJ
Definition: AliConvEventCuts.h:160
AliV0ReaderV1::iterator::fCurrentIndex
int fCurrentIndex
Index of the current element.
Definition: AliV0ReaderV1.h:68
Int_t
int Int_t
Definition: External.C:63
AliV0ReaderV1::~AliV0ReaderV1
virtual ~AliV0ReaderV1()
Definition: AliV0ReaderV1.cxx:160
AliV0ReaderV1::ProcessEvent
Bool_t ProcessEvent(AliVEvent *inputEvent, AliMCEvent *mcEvent=NULL)
Definition: AliV0ReaderV1.cxx:602
Float_t
float Float_t
Definition: External.C:68
AliConversionPhotonBase.h
AliConversionPhotonBase::SetMCLabelPositive
void SetMCLabelPositive(Int_t label)
Definition: AliConversionPhotonBase.h:65
AliKFConversionPhoton::ConstructGamma
void ConstructGamma(const AliKFParticle &fCurrentNegativeKFParticle, const AliKFParticle &fCurrentPositiveKFParticle)
Definition: AliKFConversionPhoton.cxx:74
AliAODConversionPhoton::SetInvMassPair
void SetInvMassPair(Float_t mass)
Definition: AliAODConversionPhoton.h:48
AliV0ReaderV1::iterator::iterator
iterator(const AliV0ReaderV1 *reader, Direction_t dir, int position)
Definition: AliV0ReaderV1.cxx:1816
AliV0ReaderV1::GetHelixCenter
Bool_t GetHelixCenter(const AliExternalTrackParam *track, Double_t center[2])
Definition: AliV0ReaderV1.cxx:1086
AliV0ReaderV1::RelabelAODPhotonCandidates
void RelabelAODPhotonCandidates(AliAODConversionPhoton *PhotonCandidate)
Definition: AliV0ReaderV1.cxx:1266
AliConvEventCuts::kLHC16c2
anchored LHC12[a-h] pass 2 - JJ
Definition: AliConvEventCuts.h:141
AliConversionPhotonBase::SetChi2perNDF
void SetChi2perNDF(Float_t chi2)
Definition: AliConversionPhotonBase.h:47
AliV0ReaderV1::CountTPCoutTracks
void CountTPCoutTracks()
Definition: AliV0ReaderV1.cxx:1428
AliConvEventCuts::kLHC16c3c
anchored LHC13[d-e] pass 2 - GJ
Definition: AliConvEventCuts.h:157
Data
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
Definition: charmCutsOptimization.C:371
AliV0ReaderV1::Notify
virtual Bool_t Notify()
Definition: AliV0ReaderV1.cxx:502
AliV0ReaderV1::iterator::operator=
iterator & operator=(const iterator &other)
Definition: AliV0ReaderV1.cxx:1832
AliConvEventCuts::kLHC16c3a
anchored LHC13[d-e] pass 2 - JJ
Definition: AliConvEventCuts.h:155
AliConversionPhotonBase::SetTrackLabels
void SetTrackLabels(Int_t label1, Int_t label2)
Definition: AliConversionPhotonBase.h:54
AliConvEventCuts::kLHC15g1b
anchored LHC11a pass 4 - JJ
Definition: AliConvEventCuts.h:121
AliV0ReaderV1::Terminate
virtual void Terminate(Option_t *)
Definition: AliV0ReaderV1.cxx:1811
AliV0ReaderV1::GetExternalTrackParam
const AliExternalTrackParam * GetExternalTrackParam(AliESDv0 *fCurrentV0, Int_t &tracklabel, Int_t charge)
Definition: AliV0ReaderV1.cxx:678
AliConvEventCuts::kLHC15a3b
anchored LHC13g pass 1 - JJ
Definition: AliConvEventCuts.h:161
AliV0ReaderV1::ReconstructV0
AliKFConversionPhoton * ReconstructV0(AliESDv0 *fCurrentV0, Int_t currentV0Index)
Definition: AliV0ReaderV1.cxx:766
AliConversionPhotonBase::SetLabelPositive
void SetLabelPositive(Int_t label)
Track labels.
Definition: AliConversionPhotonBase.h:52
AliConversionPhotonBase::SetLabelNegative
void SetLabelNegative(Int_t label)
Definition: AliConversionPhotonBase.h:53
AliV0ReaderV1.h
AliV0ReaderV1::iterator::operator!=
bool operator!=(iterator &other) const
Definition: AliV0ReaderV1.cxx:1841
file
TFile * file
TList with histograms for a given trigger.
Definition: DrawAnaCaloTrackQA.C:89
AliConvEventCuts::kLHC16c2_plus
anchored LHC12[a-h] pass 2 - JJ - additional stat
Definition: AliConvEventCuts.h:142
AliConversionPhotonBase::SetPsiPair
void SetPsiPair(Float_t PsiPair)
Definition: AliConversionPhotonBase.h:120
AliV0ReaderV1::FillImpactParamHistograms
void FillImpactParamHistograms(AliVTrack *ptrack, AliVTrack *ntrack, AliESDv0 *fCurrentV0, AliKFConversionPhoton *fCurrentMotherKF)
Definition: AliV0ReaderV1.cxx:1609
AliV0ReaderV1::operator[]
AliConversionPhotonBase * operator[](int index) const
Array index operator.
Definition: AliV0ReaderV1.cxx:179
Option_t
const char Option_t
Definition: External.C:48
AliConvEventCuts::kLHC15g1a
anchored LHC11a pass 4 - JJ
Definition: AliConvEventCuts.h:120
AliV0ReaderV1::ProcessESDV0s
Bool_t ProcessESDV0s()
Definition: AliV0ReaderV1.cxx:724
Bool_t
bool Bool_t
Definition: External.C:53
AliKFConversionPhoton.h
AliConvEventCuts::kLHC16P1JJ
anchored LHC16x pass 1 nom B-field - Pythia8 JJ
Definition: AliConvEventCuts.h:231
AliV0ReaderV1::GetConversionPoint
Bool_t GetConversionPoint(const AliExternalTrackParam *pparam, const AliExternalTrackParam *nparam, Double_t convpos[3], Double_t dca[2])
Definition: AliV0ReaderV1.cxx:1127
AliV0ReaderV1::GetAODConversionGammas
Bool_t GetAODConversionGammas()
Definition: AliV0ReaderV1.cxx:1208
AliV0ReaderV1::UserCreateOutputObjects
void UserCreateOutputObjects()
Definition: AliV0ReaderV1.cxx:262
AliV0ReaderV1::UserExec
virtual void UserExec(Option_t *option)
Definition: AliV0ReaderV1.cxx:583
AliConversionPhotonBase::GetConversionRadius
Double_t GetConversionRadius() const
Definition: AliConversionPhotonBase.h:91
dir
TDirectoryFile * dir
Definition: DrawAnaCaloTrackQA.C:87