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