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