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