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