]> git.uio.no Git - u/mrichter/AliRoot.git/blob - PWGGA/GammaConv/AliConversionCuts.cxx
- changes for the reweighting part
[u/mrichter/AliRoot.git] / PWGGA / GammaConv / AliConversionCuts.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Authors: Svein Lindal, Daniel Lohner                                   *
5  * Version 1.0                                                            *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 ////////////////////////////////////////////////
17 //---------------------------------------------
18 // Class handling all kinds of selection cuts for
19 // Gamma Conversion analysis
20 //---------------------------------------------
21 ////////////////////////////////////////////////
22
23 #include "AliConversionCuts.h"
24
25 #include "AliKFVertex.h"
26 #include "AliAODTrack.h"
27 #include "AliESDtrack.h"
28 #include "AliAnalysisManager.h"
29 #include "AliInputEventHandler.h"
30 #include "AliMCEventHandler.h"
31 #include "AliAODHandler.h"
32 #include "AliPIDResponse.h"
33 #include "TH1.h"
34 #include "TH2.h"
35 #include "TF1.h"
36 #include "AliStack.h"
37 #include "AliAODConversionMother.h"
38 #include "TObjString.h"
39 #include "AliAODEvent.h"
40 #include "AliESDEvent.h"
41 #include "AliCentrality.h"
42 #include "TList.h"
43 #include "TFile.h"
44 #include "AliLog.h"
45 #include "AliGenCocktailEventHeader.h"
46 #include "AliGenDPMjetEventHeader.h"
47 #include "AliGenPythiaEventHeader.h"
48 #include "AliGenHijingEventHeader.h"
49 #include "AliTriggerAnalysis.h"
50 #include "AliV0ReaderV1.h"
51 #include "AliAODMCParticle.h"
52 #include "AliAODMCHeader.h"
53
54 class iostream;
55
56 using namespace std;
57
58 ClassImp(AliConversionCuts)
59
60
61 const char* AliConversionCuts::fgkCutNames[AliConversionCuts::kNCuts] = {
62    "HeavyIon",//0
63    "CentralityMin",//1
64    "CentralityMax",//2
65    "SelectV0AND",//3
66    "MultiplicityMethod",//4
67    "RemovePileUp",//5
68    "RejectExtraSignals",//6
69    "V0FinderType",//7
70    "EtaCut",//8
71    "MinRCut",//9
72    "SinglePtCut",//10
73    "ClsTPCCut", //11
74    "ededxSigmaCut",//12
75    "pidedxSigmaCut",//13
76    "piMomdedxSigmaCut",//14
77    "piMaxMomdedxSigmaCut",//15
78    "LowPRejectionSigmaCut",//16
79    "TOFelectronPID",//17
80    "QtMaxCut",//18
81    "Chi2GammaCut", //19
82    "PsiPair", //20
83    "DoPhotonAsymmetryCut",//21
84    "CosinePointingAngle", //22
85    "SharedElectronCuts", //23
86    "RejectToCloseV0s", //24
87    "DcaRPrimVtx", //25
88    "DcaZPrimVtx", //26
89    "EvetPlane" //27
90 };
91
92
93 //________________________________________________________________________
94 AliConversionCuts::AliConversionCuts(const char *name,const char *title) :
95    AliAnalysisCuts(name,title),
96    fHistograms(NULL),
97    fHeaderList(NULL),
98    fPIDResponse(NULL),
99    fEventQuality(-1),
100    fMaxR(200),
101    fMinR(0),
102    fEtaCut(0.9),
103    fEtaCutMin(-0.1),
104    fPtCut(0.02),
105    fSinglePtCut(0),
106    fMaxZ(1000),
107    fMinClsTPC(0.),
108    fMinClsTPCToF(0.),
109    fLineCutZRSlope(0.),
110    fLineCutZValue(0),
111    fLineCutZRSlopeMin(0.),
112    fLineCutZValueMin(0),
113    fChi2CutConversion(1000),
114    fPIDProbabilityCutNegativeParticle(0),
115    fPIDProbabilityCutPositiveParticle(0),
116    fDodEdxSigmaCut(kTRUE),
117    fDoTOFsigmaCut(kFALSE),
118    fPIDTRDEfficiency(1),
119    fDoTRDPID(kFALSE),
120    fPIDnSigmaAboveElectronLine(100),
121    fPIDnSigmaBelowElectronLine(-100),
122    fTofPIDnSigmaAboveElectronLine(100),
123    fTofPIDnSigmaBelowElectronLine(-100),
124    fPIDnSigmaAbovePionLine(0),
125    fPIDnSigmaAbovePionLineHighPt(-100),
126    fPIDMinPnSigmaAbovePionLine(0),
127    fPIDMaxPnSigmaAbovePionLine(0),
128    fDoKaonRejectionLowP(kFALSE),
129    fDoProtonRejectionLowP(kFALSE),
130    fDoPionRejectionLowP(kFALSE),
131    fPIDnSigmaAtLowPAroundKaonLine(0),
132    fPIDnSigmaAtLowPAroundProtonLine(0),
133    fPIDnSigmaAtLowPAroundPionLine(0),
134    fPIDMinPKaonRejectionLowP(1.5),
135    fPIDMinPProtonRejectionLowP(2),
136    fPIDMinPPionRejectionLowP(0),
137    fDoQtGammaSelection(kTRUE),
138    fDo2DQt(kFALSE),
139    fQtMax(100),
140    fXVertexCut(0.),
141    fYVertexCut(0.),
142    fZVertexCut(0.),
143    fNSigmaMass(0.),
144    fUseEtaMinCut(kFALSE),
145    fUseOnFlyV0Finder(kTRUE),
146    fDoPhotonAsymmetryCut(kTRUE),
147    fMinPPhotonAsymmetryCut(100.),
148    fMinPhotonAsymmetry(0.),
149    fIsHeavyIon(0),
150    fDetectorCentrality(0),
151    fModCentralityClass(0),
152    fMaxVertexZ(10),
153    fCentralityMin(0),
154    fCentralityMax(0),
155    fUseCorrectedTPCClsInfo(kFALSE),
156    fUseTOFpid(kFALSE),
157    fMultiplicityMethod(0),
158    fSpecialTrigger(0),
159    fRemovePileUp(kFALSE),
160    fOpeningAngle(0.005),
161    fPsiPairCut(10000),
162    fDo2DPsiPairChi2(kFALSE),
163    fCosPAngleCut(10000),
164    fDoToCloseV0sCut(kFALSE),
165    fRejectExtraSignals(0),
166    fminV0Dist(200.),
167    fDoSharedElecCut(kFALSE),
168    fDoPhotonQualitySelectionCut(kFALSE),
169    fPhotonQualityCut(0),
170    fOfflineTriggerMask(0),
171    fHasV0AND(kTRUE),
172    fIsSDDFired(kTRUE),
173    fRandom(0),
174    fElectronArraySize(500),
175    fElectronLabelArray(NULL),
176    fDCAZPrimVtxCut(1000),
177    fDCARPrimVtxCut(1000),
178    fInPlaneOutOfPlane(0),
179    fConversionPointXArray(0.0),
180    fConversionPointYArray(0.0),
181    fConversionPointZArray(0.0),
182    fnHeaders(0),
183    fNotRejectedStart(NULL),
184    fNotRejectedEnd(NULL),
185    fGeneratorNames(NULL),
186    fCutString(NULL),
187    fUtils(NULL),
188    fEtaShift(0.0),
189    fDoEtaShift(kFALSE),
190    fDoReweightHistoMCPi0(kFALSE),
191    fDoReweightHistoMCEta(kFALSE),
192    fDoReweightHistoMCK0s(kFALSE),
193    fPathTrFReweighting(""),
194    fNameHistoReweightingPi0(""),
195    fNameHistoReweightingEta(""),
196    fNameHistoReweightingK0s(""),
197    fNameFitDataPi0(""),
198    fNameFitDataEta(""),
199    fNameFitDataK0s(""),
200    hEtaDistV0s(NULL),
201    hEtaDistV0sAfterdEdxCuts(NULL),
202    hdEdxCuts(NULL),
203    hTPCdEdxbefore(NULL),
204    hTPCdEdxafter(NULL),
205    hTPCdEdxSigbefore(NULL),
206    hTPCdEdxSigafter(NULL),
207    hTOFbefore(NULL),
208    hTOFSigbefore(NULL),
209    hTOFSigafter(NULL),
210    hPsiPairDeltaPhiafter(NULL),
211    hTrackCuts(NULL),
212    hPhotonCuts(NULL),
213    hInvMassbefore(NULL),
214    hArmenterosbefore(NULL),
215    hInvMassafter(NULL),
216    hArmenterosafter(NULL),
217    hAcceptanceCuts(NULL),
218    hCutIndex(NULL),
219    hV0EventCuts(NULL),
220    hCentrality(NULL),
221    hCentralityVsNumberOfPrimaryTracks(NULL),
222    hVertexZ(NULL),
223    hEventPlanePhi(NULL),
224    hTriggerClass(NULL),
225    hTriggerClassSelected(NULL),
226    hReweightMCHistPi0(NULL),
227    hReweightMCHistEta(NULL),
228    hReweightMCHistK0s(NULL),
229    fFitDataPi0(NULL),
230    fFitDataEta(NULL),
231    fFitDataK0s(NULL),
232    fPreSelCut(kFALSE),
233    fTriggerSelectedManually(kFALSE),
234    fSpecialTriggerName("")
235
236 {
237    InitPIDResponse();
238    for(Int_t jj=0;jj<kNCuts;jj++){fCuts[jj]=0;}
239    fCutString=new TObjString((GetCutNumber()).Data());
240
241    fElectronLabelArray = new Int_t[fElectronArraySize];
242    fUtils = new AliAnalysisUtils();
243    //if you do not want to apply the cut on the distance between the SPD and TRK vertex:
244    //fUtils->SetCutOnZVertexSPD(kFALSE);
245
246
247 }
248
249 //________________________________________________________________________
250 AliConversionCuts::AliConversionCuts(const AliConversionCuts &ref) :
251    AliAnalysisCuts(ref),
252    fHistograms(NULL),
253    fHeaderList(ref.fHeaderList),
254    fPIDResponse(NULL),
255    fEventQuality(ref.fEventQuality),
256    fMaxR(ref.fMaxR),
257    fMinR(ref.fMinR),
258    fEtaCut(ref.fEtaCut),
259    fEtaCutMin(ref.fEtaCutMin),
260    fPtCut(ref.fPtCut),
261    fSinglePtCut(ref.fSinglePtCut),
262    fMaxZ(ref.fMaxZ),
263    fMinClsTPC(ref.fMinClsTPC),
264    fMinClsTPCToF(ref.fMinClsTPCToF),
265    fLineCutZRSlope(ref.fLineCutZRSlope),
266    fLineCutZValue(ref.fLineCutZValue),
267    fLineCutZRSlopeMin(ref.fLineCutZRSlopeMin),
268    fLineCutZValueMin(ref.fLineCutZValueMin),
269    fChi2CutConversion(ref.fChi2CutConversion),
270    fPIDProbabilityCutNegativeParticle(ref.fPIDProbabilityCutNegativeParticle),
271    fPIDProbabilityCutPositiveParticle(ref.fPIDProbabilityCutPositiveParticle),
272    fDodEdxSigmaCut(ref. fDodEdxSigmaCut),
273    fDoTOFsigmaCut(ref.fDoTOFsigmaCut),
274    fPIDTRDEfficiency(ref.fPIDTRDEfficiency),
275    fDoTRDPID(ref.fDoTRDPID),
276    fPIDnSigmaAboveElectronLine(ref.fPIDnSigmaAboveElectronLine),
277    fPIDnSigmaBelowElectronLine(ref.fPIDnSigmaBelowElectronLine),
278    fTofPIDnSigmaAboveElectronLine(ref.fTofPIDnSigmaAboveElectronLine),
279    fTofPIDnSigmaBelowElectronLine(ref.fTofPIDnSigmaBelowElectronLine),
280    fPIDnSigmaAbovePionLine(ref.fPIDnSigmaAbovePionLine),
281    fPIDnSigmaAbovePionLineHighPt(ref.fPIDnSigmaAbovePionLineHighPt),
282    fPIDMinPnSigmaAbovePionLine(ref.fPIDMinPnSigmaAbovePionLine),
283    fPIDMaxPnSigmaAbovePionLine(ref.fPIDMaxPnSigmaAbovePionLine),
284    fDoKaonRejectionLowP(ref.fDoKaonRejectionLowP),
285    fDoProtonRejectionLowP(ref.fDoProtonRejectionLowP),
286    fDoPionRejectionLowP(ref.fDoPionRejectionLowP),
287    fPIDnSigmaAtLowPAroundKaonLine(ref.fPIDnSigmaAtLowPAroundKaonLine),
288    fPIDnSigmaAtLowPAroundProtonLine(ref.fPIDnSigmaAtLowPAroundProtonLine),
289    fPIDnSigmaAtLowPAroundPionLine(ref.fPIDnSigmaAtLowPAroundPionLine),
290    fPIDMinPKaonRejectionLowP(ref.fPIDMinPKaonRejectionLowP),
291    fPIDMinPProtonRejectionLowP(ref.fPIDMinPProtonRejectionLowP),
292    fPIDMinPPionRejectionLowP(ref.fPIDMinPPionRejectionLowP),
293    fDoQtGammaSelection(ref.fDoQtGammaSelection),
294    fDo2DQt(ref.fDo2DQt),
295    fQtMax(ref.fQtMax),
296    fXVertexCut(ref.fXVertexCut),
297    fYVertexCut(ref.fYVertexCut),
298    fZVertexCut(ref.fZVertexCut),
299    fNSigmaMass(ref.fNSigmaMass),
300    fUseEtaMinCut(ref.fUseEtaMinCut),
301    fUseOnFlyV0Finder(ref.fUseOnFlyV0Finder),
302    fDoPhotonAsymmetryCut(ref.fDoPhotonAsymmetryCut),
303    fMinPPhotonAsymmetryCut(ref.fMinPPhotonAsymmetryCut),
304    fMinPhotonAsymmetry(ref.fMinPhotonAsymmetry),
305    fIsHeavyIon(ref.fIsHeavyIon),
306    fDetectorCentrality(ref.fDetectorCentrality),
307    fModCentralityClass(ref.fModCentralityClass),
308    fMaxVertexZ(ref.fMaxVertexZ),
309    fCentralityMin(ref.fCentralityMin),
310    fCentralityMax(ref.fCentralityMax),
311    fUseCorrectedTPCClsInfo(ref.fUseCorrectedTPCClsInfo),
312    fUseTOFpid(ref.fUseTOFpid),
313    fMultiplicityMethod(ref.fMultiplicityMethod),
314    fSpecialTrigger(ref.fSpecialTrigger),
315    fRemovePileUp(ref.fRemovePileUp),
316    fOpeningAngle(ref.fOpeningAngle),
317    fPsiPairCut(ref.fPsiPairCut),
318    fDo2DPsiPairChi2(ref.fDo2DPsiPairChi2),
319    fCosPAngleCut(ref.fCosPAngleCut),
320    fDoToCloseV0sCut(ref.fDoToCloseV0sCut),
321    fRejectExtraSignals(ref.fRejectExtraSignals),
322    fminV0Dist(ref.fminV0Dist),
323    fDoSharedElecCut(ref.fDoSharedElecCut),
324    fDoPhotonQualitySelectionCut(ref.fDoPhotonQualitySelectionCut),
325    fPhotonQualityCut(ref.fPhotonQualityCut),
326    fOfflineTriggerMask(ref.fOfflineTriggerMask),
327    fHasV0AND(ref.fHasV0AND),
328    fIsSDDFired(ref.fIsSDDFired),
329    fRandom(ref.fRandom),
330    fElectronArraySize(ref.fElectronArraySize),
331    fElectronLabelArray(NULL),
332    fDCAZPrimVtxCut(ref.fDCAZPrimVtxCut),
333    fDCARPrimVtxCut(ref.fDCAZPrimVtxCut),
334    fInPlaneOutOfPlane(ref.fInPlaneOutOfPlane),
335    fConversionPointXArray(ref.fConversionPointXArray),
336    fConversionPointYArray(ref.fConversionPointYArray),
337    fConversionPointZArray(ref.fConversionPointZArray),
338    fnHeaders(ref.fnHeaders),
339    fNotRejectedStart(NULL),
340    fNotRejectedEnd(NULL),
341    fGeneratorNames(ref.fGeneratorNames),
342    fCutString(NULL),
343    fUtils(NULL),
344    fEtaShift(ref.fEtaShift),
345    fDoEtaShift(ref.fDoEtaShift),
346    fDoReweightHistoMCPi0(ref.fDoReweightHistoMCPi0),
347    fDoReweightHistoMCEta(ref.fDoReweightHistoMCEta),
348    fDoReweightHistoMCK0s(ref.fDoReweightHistoMCK0s),
349    fPathTrFReweighting(ref.fPathTrFReweighting),
350    fNameHistoReweightingPi0(ref.fNameHistoReweightingPi0),
351    fNameHistoReweightingEta(ref.fNameHistoReweightingEta),
352    fNameHistoReweightingK0s(ref.fNameHistoReweightingK0s),
353    fNameFitDataPi0(ref.fNameFitDataPi0),
354    fNameFitDataEta(ref.fNameFitDataEta),
355    fNameFitDataK0s(ref.fNameFitDataK0s),
356    hEtaDistV0s(NULL),
357    hEtaDistV0sAfterdEdxCuts(NULL),
358    hdEdxCuts(NULL),
359    hTPCdEdxbefore(NULL),
360    hTPCdEdxafter(NULL),
361    hTPCdEdxSigbefore(NULL),
362    hTPCdEdxSigafter(NULL),
363    hTOFbefore(NULL),
364    hTOFSigbefore(NULL),
365    hTOFSigafter(NULL),
366    hPsiPairDeltaPhiafter(NULL),
367    hTrackCuts(NULL),
368    hPhotonCuts(NULL),
369    hInvMassbefore(NULL),
370    hArmenterosbefore(NULL),
371    hInvMassafter(NULL),
372    hArmenterosafter(NULL),
373    hAcceptanceCuts(NULL),
374    hCutIndex(NULL),
375    hV0EventCuts(NULL),
376    hCentrality(NULL),
377    hCentralityVsNumberOfPrimaryTracks(NULL),
378    hVertexZ(NULL),
379    hEventPlanePhi(NULL),
380    hTriggerClass(NULL),
381    hTriggerClassSelected(NULL),
382    hReweightMCHistPi0(ref.hReweightMCHistPi0),
383    hReweightMCHistEta(ref.hReweightMCHistEta),
384    hReweightMCHistK0s(ref.hReweightMCHistK0s),
385    fFitDataPi0(ref.fFitDataPi0),
386    fFitDataEta(ref.fFitDataEta),
387    fFitDataK0s(ref.fFitDataK0s),
388    fPreSelCut(ref.fPreSelCut),
389    fTriggerSelectedManually(ref.fTriggerSelectedManually),
390    fSpecialTriggerName(ref.fSpecialTriggerName)
391 {
392    // Copy Constructor
393    for(Int_t jj=0;jj<kNCuts;jj++){fCuts[jj]=ref.fCuts[jj];}
394    fCutString=new TObjString((GetCutNumber()).Data());
395    fElectronLabelArray = new Int_t[fElectronArraySize];
396    fUtils = new AliAnalysisUtils();
397    // dont copy histograms (if you like histograms, call InitCutHistograms())
398
399 }
400
401
402 //________________________________________________________________________
403 AliConversionCuts::~AliConversionCuts() {
404    // Destructor
405    //Deleting fHistograms leads to seg fault it it's added to output collection of a task
406    // if(fHistograms)
407    //    delete fHistograms;
408    // fHistograms = NULL;
409    if(fCutString != NULL){
410       delete fCutString;
411       fCutString = NULL;
412    }
413    if(fElectronLabelArray){
414       delete fElectronLabelArray;
415       fElectronLabelArray = NULL;
416    }
417    if(fNotRejectedStart){
418       delete[] fNotRejectedStart;
419       fNotRejectedStart = NULL;
420    }
421    if(fNotRejectedEnd){
422       delete[] fNotRejectedEnd;
423       fNotRejectedEnd = NULL;
424    }
425    if(fGeneratorNames){
426       delete[] fGeneratorNames;
427       fGeneratorNames = NULL;
428    }
429    if(fUtils){
430      delete fUtils;
431      fUtils = NULL;
432    }
433
434 }
435
436 //________________________________________________________________________
437 void AliConversionCuts::InitCutHistograms(TString name, Bool_t preCut){
438
439    // Initialize Cut Histograms for QA (only initialized and filled if function is called)
440    TH1::AddDirectory(kFALSE);
441
442    if(fHistograms != NULL){
443       delete fHistograms;
444       fHistograms=NULL;
445    }
446    if(fHistograms==NULL){
447       fHistograms=new TList();
448       fHistograms->SetOwner(kTRUE);
449       if(name=="")fHistograms->SetName(Form("ConvCuts_%s",GetCutNumber().Data()));
450       else fHistograms->SetName(Form("%s_%s",name.Data(),GetCutNumber().Data()));
451    }
452
453    if (hReweightMCHistPi0){
454       hReweightMCHistPi0->SetName("MCInputForWeightingPi0");
455       fHistograms->Add(hReweightMCHistPi0);
456    }
457    if (hReweightMCHistEta){
458       hReweightMCHistEta->SetName("MCInputForWeightingEta");
459       fHistograms->Add(hReweightMCHistEta);
460    }
461    if (hReweightMCHistK0s){
462       hReweightMCHistK0s->SetName("MCInputForWeightingK0s");
463       fHistograms->Add(hReweightMCHistK0s);
464    }
465 //    if (fFitDataPi0){
466 //       fFitDataPi0->SetName("DataFitForWeightingPi0");
467 //       fHistograms->Add(fFitDataPi0);
468 //    }
469 //    if (fFitDataEta){
470 //       fFitDataEta->SetName("DataFitForWeightingEta");
471 //       fHistograms->Add(fFitDataEta);
472 //    }
473 //    if (fFitDataK0s){
474 //       fFitDataK0s->SetName("DataFitForWeightingK0s");
475 //       fHistograms->Add(fFitDataK0s);
476 //    }
477    // IsPhotonSelected
478    hCutIndex=new TH1F(Form("IsPhotonSelected %s",GetCutNumber().Data()),"IsPhotonSelected",11,-0.5,10.5);
479    hCutIndex->GetXaxis()->SetBinLabel(kPhotonIn+1,"in");
480    hCutIndex->GetXaxis()->SetBinLabel(kOnFly+1,"onfly");
481    hCutIndex->GetXaxis()->SetBinLabel(kNoTracks+1,"no tracks");
482    hCutIndex->GetXaxis()->SetBinLabel(kdEdxCuts+1,"dEdx");
483    hCutIndex->GetXaxis()->SetBinLabel(kTrackCuts+1,"Track cuts");
484    hCutIndex->GetXaxis()->SetBinLabel(kConvPointFail+1,"ConvPoint fail");
485    hCutIndex->GetXaxis()->SetBinLabel(kPhotonCuts+1,"PhotonCuts");
486    hCutIndex->GetXaxis()->SetBinLabel(kEventPlane+1,"EventPlane");
487    hCutIndex->GetXaxis()->SetBinLabel(kPhotonOut+1,"out");
488    fHistograms->Add(hCutIndex);
489
490    // Track Cuts
491    hTrackCuts=new TH1F(Form("TrackCuts %s",GetCutNumber().Data()),"TrackCuts",9,-0.5,8.5);
492    hTrackCuts->GetXaxis()->SetBinLabel(1,"in");
493    hTrackCuts->GetXaxis()->SetBinLabel(2,"likesign");
494    hTrackCuts->GetXaxis()->SetBinLabel(3,"ntpccl");
495    hTrackCuts->GetXaxis()->SetBinLabel(4,"acceptance");
496    hTrackCuts->GetXaxis()->SetBinLabel(5,"singlept");
497    hTrackCuts->GetXaxis()->SetBinLabel(6,"TPCrefit");
498    hTrackCuts->GetXaxis()->SetBinLabel(7,"kink");
499    hTrackCuts->GetXaxis()->SetBinLabel(8,"out");
500    fHistograms->Add(hTrackCuts);
501
502    // Photon Cuts
503    hPhotonCuts=new TH1F(Form("PhotonCuts %s",GetCutNumber().Data()),"PhotonCuts",15,-0.5,14.5);
504    hPhotonCuts->GetXaxis()->SetBinLabel(1,"in");
505    hPhotonCuts->GetXaxis()->SetBinLabel(2,"qtcut");
506    hPhotonCuts->GetXaxis()->SetBinLabel(3,"chi2");
507    hPhotonCuts->GetXaxis()->SetBinLabel(4,"acceptance");
508    hPhotonCuts->GetXaxis()->SetBinLabel(5,"asymmetry");
509    hPhotonCuts->GetXaxis()->SetBinLabel(6,"pidprob");
510    hPhotonCuts->GetXaxis()->SetBinLabel(7,"cortpcclinfo");
511    hPhotonCuts->GetXaxis()->SetBinLabel(8,"PsiPair");
512    hPhotonCuts->GetXaxis()->SetBinLabel(9,"CosPAngle");
513    hPhotonCuts->GetXaxis()->SetBinLabel(10,"DCA R");
514    hPhotonCuts->GetXaxis()->SetBinLabel(11,"DCA Z");
515    hPhotonCuts->GetXaxis()->SetBinLabel(12,"Photon Quality");
516    hPhotonCuts->GetXaxis()->SetBinLabel(13,"out");
517    fHistograms->Add(hPhotonCuts);
518
519    if(preCut){
520       hInvMassbefore=new TH1F(Form("InvMass_before %s",GetCutNumber().Data()),"InvMass_before",1000,0,0.3);
521       fHistograms->Add(hInvMassbefore);
522       hArmenterosbefore=new TH2F(Form("Armenteros_before %s",GetCutNumber().Data()),"Armenteros_before",200,-1,1,1000,0,1.);
523       fHistograms->Add(hArmenterosbefore);
524       hEtaDistV0s = new TH1F(Form("Eta_before %s",GetCutNumber().Data()),"Eta_before",2000,-2,2);
525       fHistograms->Add(hEtaDistV0s);
526
527    }
528    hInvMassafter=new TH1F(Form("InvMass_after %s",GetCutNumber().Data()),"InvMass_after",1000,0,0.3);
529    fHistograms->Add(hInvMassafter);
530    hArmenterosafter=new TH2F(Form("Armenteros_after %s",GetCutNumber().Data()),"Armenteros_after",200,-1,1,250,0,0.25);
531    fHistograms->Add(hArmenterosafter);
532
533    hAcceptanceCuts=new TH1F(Form("PhotonAcceptanceCuts %s",GetCutNumber().Data()),"PhotonAcceptanceCuts",10,-0.5,9.5);
534    hAcceptanceCuts->GetXaxis()->SetBinLabel(1,"in");
535    hAcceptanceCuts->GetXaxis()->SetBinLabel(2,"maxR");
536    hAcceptanceCuts->GetXaxis()->SetBinLabel(3,"minR");
537    hAcceptanceCuts->GetXaxis()->SetBinLabel(4,"line");
538    hAcceptanceCuts->GetXaxis()->SetBinLabel(5,"maxZ");
539    hAcceptanceCuts->GetXaxis()->SetBinLabel(6,"eta");
540    hAcceptanceCuts->GetXaxis()->SetBinLabel(7,"minpt");
541    hAcceptanceCuts->GetXaxis()->SetBinLabel(8,"out");
542    fHistograms->Add(hAcceptanceCuts);
543
544    // dEdx Cuts
545    hdEdxCuts=new TH1F(Form("dEdxCuts %s",GetCutNumber().Data()),"dEdxCuts",10,-0.5,9.5);
546    hdEdxCuts->GetXaxis()->SetBinLabel(1,"in");
547    hdEdxCuts->GetXaxis()->SetBinLabel(2,"TPCelectron");
548    hdEdxCuts->GetXaxis()->SetBinLabel(3,"TPCpion");
549    hdEdxCuts->GetXaxis()->SetBinLabel(4,"TPCpionhighp");
550    hdEdxCuts->GetXaxis()->SetBinLabel(5,"TPCkaonlowprej");
551    hdEdxCuts->GetXaxis()->SetBinLabel(6,"TPCprotonlowprej");
552    hdEdxCuts->GetXaxis()->SetBinLabel(7,"TPCpionlowprej");
553    hdEdxCuts->GetXaxis()->SetBinLabel(8,"TOFelectron");
554    hdEdxCuts->GetXaxis()->SetBinLabel(9,"TRDelectron");
555    hdEdxCuts->GetXaxis()->SetBinLabel(10,"out");
556    fHistograms->Add(hdEdxCuts);
557
558    TAxis *AxisBeforedEdx = NULL;
559    TAxis *AxisBeforedEdxSig = NULL;
560    TAxis *AxisBeforeTOF = NULL;
561    TAxis *AxisBeforeTOFSig = NULL;
562    if(preCut){
563       hTPCdEdxbefore=new TH2F(Form("Gamma_dEdx_before %s",GetCutNumber().Data()),"dEdx Gamma before" ,150,0.03,20,800,0,200);
564       fHistograms->Add(hTPCdEdxbefore);
565       AxisBeforedEdx = hTPCdEdxbefore->GetXaxis();
566       hTPCdEdxSigbefore=new TH2F(Form("Gamma_dEdxSig_before %s",GetCutNumber().Data()),"dEdx Sigma Gamma before" ,150,0.03,20,400,-10,10);
567       fHistograms->Add(hTPCdEdxSigbefore);
568       AxisBeforedEdxSig = hTPCdEdxSigbefore->GetXaxis();
569
570       hTOFbefore=new TH2F(Form("Gamma_TOF_before %s",GetCutNumber().Data()),"TOF Gamma before" ,150,0.03,20,11000,-1000,10000);
571       fHistograms->Add(hTOFbefore);
572       AxisBeforeTOF = hTOFbefore->GetXaxis();
573       hTOFSigbefore=new TH2F(Form("Gamma_TOFSig_before %s",GetCutNumber().Data()),"TOF Sigma Gamma before" ,150,0.03,20,400,-6,10);
574       fHistograms->Add(hTOFSigbefore);
575       AxisBeforeTOFSig = hTOFSigbefore->GetXaxis();
576
577    }
578    hTPCdEdxSigafter=new TH2F(Form("Gamma_dEdxSig_after %s",GetCutNumber().Data()),"dEdx Sigma Gamma after" ,150,0.03,20,400, -10,10);
579    fHistograms->Add(hTPCdEdxSigafter);
580
581    hTPCdEdxafter=new TH2F(Form("Gamma_dEdx_after %s",GetCutNumber().Data()),"dEdx Gamma after" ,150,0.03,20,800,0,200);
582    fHistograms->Add(hTPCdEdxafter);
583
584    hTOFSigafter=new TH2F(Form("Gamma_TOFSig_after %s",GetCutNumber().Data()),"TOF Sigma Gamma after" ,150,0.03,20,400,-6,10);
585    fHistograms->Add(hTOFSigafter);
586
587    hEtaDistV0sAfterdEdxCuts = new TH1F(Form("Eta_afterdEdx %s",GetCutNumber().Data()),"Eta_afterdEdx",2000,-2,2);
588    fHistograms->Add(hEtaDistV0sAfterdEdxCuts);
589
590    hPsiPairDeltaPhiafter=new TH2F(Form("Gamma_PsiPairDeltaPhi_after %s",GetCutNumber().Data()),"Psi Pair vs Delta Phi Gamma after" ,200,-2,2,200,-2,2);
591    fHistograms->Add(hPsiPairDeltaPhiafter);
592
593    TAxis *AxisAfter = hTPCdEdxSigafter->GetXaxis();
594    Int_t bins = AxisAfter->GetNbins();
595    Double_t from = AxisAfter->GetXmin();
596    Double_t to = AxisAfter->GetXmax();
597    Double_t *newBins = new Double_t[bins+1];
598    newBins[0] = from;
599    Double_t factor = TMath::Power(to/from, 1./bins);
600    for(Int_t i=1; i<=bins; ++i) newBins[i] = factor * newBins[i-1];
601    AxisAfter->Set(bins, newBins);
602    AxisAfter = hTOFSigafter->GetXaxis();
603    AxisAfter->Set(bins, newBins);
604    AxisAfter = hTPCdEdxafter->GetXaxis();
605    AxisAfter->Set(bins, newBins);
606    if(preCut){
607       AxisBeforedEdx->Set(bins, newBins);
608       AxisBeforeTOF->Set(bins, newBins);
609       AxisBeforedEdxSig->Set(bins, newBins);
610       AxisBeforeTOFSig->Set(bins, newBins);
611    }
612    delete [] newBins;
613
614    hCentrality=new TH1F(Form("Centrality %s",GetCutNumber().Data()),"Centrality",100,0,100);
615    fHistograms->Add(hCentrality);
616    hCentralityVsNumberOfPrimaryTracks=new TH2F(Form("Centrality vs Primary Tracks %s",GetCutNumber().Data()),"Centrality vs Primary Tracks ",100,0,100,4000,0,4000);
617    fHistograms->Add(hCentralityVsNumberOfPrimaryTracks);
618
619    // Event Cuts and Info
620    if(preCut){
621       hV0EventCuts=new TH1F(Form("ESD_EventCuts %s",GetCutNumber().Data()),"Event Cuts",7,-0.5,6.5);
622       hV0EventCuts->GetXaxis()->SetBinLabel(1,"in");
623       hV0EventCuts->GetXaxis()->SetBinLabel(2,"OfflineTrigger");
624       hV0EventCuts->GetXaxis()->SetBinLabel(3,"nvtxcontr");
625       hV0EventCuts->GetXaxis()->SetBinLabel(4,"VertexZ");
626       hV0EventCuts->GetXaxis()->SetBinLabel(5,"pileup");
627       hV0EventCuts->GetXaxis()->SetBinLabel(6,"centrsel");
628       hV0EventCuts->GetXaxis()->SetBinLabel(7,"out");
629       fHistograms->Add(hV0EventCuts);
630
631       hVertexZ=new TH1F(Form("VertexZ %s",GetCutNumber().Data()),"VertexZ",1000,-50,50);
632       fHistograms->Add(hVertexZ);
633
634       hTriggerClass= new TH1F(Form("OfflineTrigger %s",GetCutNumber().Data()),"OfflineTrigger",35,-0.5,34.5);
635       hTriggerClass->GetXaxis()->SetBinLabel( 1,"kMB");
636       hTriggerClass->GetXaxis()->SetBinLabel( 2,"kINT7");
637       hTriggerClass->GetXaxis()->SetBinLabel( 3,"kMUON");
638       hTriggerClass->GetXaxis()->SetBinLabel( 4,"kHighMult");
639       hTriggerClass->GetXaxis()->SetBinLabel( 5,"kKEMC1");
640       hTriggerClass->GetXaxis()->SetBinLabel( 6,"kCINT5");
641       hTriggerClass->GetXaxis()->SetBinLabel( 7,"kCMUS5/kMUSPB");
642       hTriggerClass->GetXaxis()->SetBinLabel( 8,"kMUSH7/kMUSHPB");
643       hTriggerClass->GetXaxis()->SetBinLabel( 9,"kMUL7/kMuonLikePB");
644       hTriggerClass->GetXaxis()->SetBinLabel(10,"kMUU7/kMuonUnlikePB");
645       hTriggerClass->GetXaxis()->SetBinLabel(11,"kEMC7/kEMC8");
646       hTriggerClass->GetXaxis()->SetBinLabel(12,"kMUS7");
647       hTriggerClass->GetXaxis()->SetBinLabel(13,"kPHI1");
648       hTriggerClass->GetXaxis()->SetBinLabel(14,"kPHI7/kPHI8/kPHOSPb");
649       hTriggerClass->GetXaxis()->SetBinLabel(15,"kEMCEJE");
650       hTriggerClass->GetXaxis()->SetBinLabel(16,"kEMCEGA");
651       hTriggerClass->GetXaxis()->SetBinLabel(17,"kCentral");
652       hTriggerClass->GetXaxis()->SetBinLabel(18,"kSemiCentral");
653       hTriggerClass->GetXaxis()->SetBinLabel(19,"kDG5");
654       hTriggerClass->GetXaxis()->SetBinLabel(20,"kZED");
655       hTriggerClass->GetXaxis()->SetBinLabel(21,"kSPI7/kSPI");
656       hTriggerClass->GetXaxis()->SetBinLabel(22,"kINT8");
657       hTriggerClass->GetXaxis()->SetBinLabel(23,"kMuonSingleLowPt8");
658       hTriggerClass->GetXaxis()->SetBinLabel(24,"kMuonSingleHighPt8");
659       hTriggerClass->GetXaxis()->SetBinLabel(25,"kMuonLikeLowPt8");
660       hTriggerClass->GetXaxis()->SetBinLabel(26,"kMuonUnlikeLowPt8");
661       hTriggerClass->GetXaxis()->SetBinLabel(27,"kMuonUnlikeLowPt0");
662       hTriggerClass->GetXaxis()->SetBinLabel(28,"kUserDefined");
663       hTriggerClass->GetXaxis()->SetBinLabel(29,"kTRD");
664       hTriggerClass->GetXaxis()->SetBinLabel(30,"kFastOnly");
665       hTriggerClass->GetXaxis()->SetBinLabel(31,"kAnyINT");
666       hTriggerClass->GetXaxis()->SetBinLabel(32,"kAny");
667       hTriggerClass->GetXaxis()->SetBinLabel(33,"V0AND");
668       hTriggerClass->GetXaxis()->SetBinLabel(34,"NOT kFastOnly");
669       hTriggerClass->GetXaxis()->SetBinLabel(35,"failed Physics Selection");
670       fHistograms->Add(hTriggerClass);
671    }
672    if(!preCut){
673       hTriggerClassSelected= new TH1F(Form("OfflineTriggerSelected %s",GetCutNumber().Data()),"OfflineTriggerSelected",34,-0.5,33.5);
674       hTriggerClassSelected->GetXaxis()->SetBinLabel( 1,"kMB");
675       hTriggerClassSelected->GetXaxis()->SetBinLabel( 2,"kINT7");
676       hTriggerClassSelected->GetXaxis()->SetBinLabel( 3,"kMUON");
677       hTriggerClassSelected->GetXaxis()->SetBinLabel( 4,"kHighMult");
678       hTriggerClassSelected->GetXaxis()->SetBinLabel( 5,"kKEMC1");
679       hTriggerClassSelected->GetXaxis()->SetBinLabel( 6,"kCINT5");
680       hTriggerClassSelected->GetXaxis()->SetBinLabel( 7,"kCMUS5/kMUSPB");
681       hTriggerClassSelected->GetXaxis()->SetBinLabel( 8,"kMUSH7/kMUSHPB");
682       hTriggerClassSelected->GetXaxis()->SetBinLabel( 9,"kMUL7/kMuonLikePB");
683       hTriggerClassSelected->GetXaxis()->SetBinLabel(10,"kMUU7/kMuonUnlikePB");
684       hTriggerClassSelected->GetXaxis()->SetBinLabel(11,"kEMC7/kEMC8");
685       hTriggerClassSelected->GetXaxis()->SetBinLabel(12,"kMUS7");
686       hTriggerClassSelected->GetXaxis()->SetBinLabel(13,"kPHI1");
687       hTriggerClassSelected->GetXaxis()->SetBinLabel(14,"kPHI7/kPHI8/kPHOSPb");
688       hTriggerClassSelected->GetXaxis()->SetBinLabel(15,"kEMCEJE");
689       hTriggerClassSelected->GetXaxis()->SetBinLabel(16,"kEMCEGA");
690       hTriggerClassSelected->GetXaxis()->SetBinLabel(17,"kCentral");
691       hTriggerClassSelected->GetXaxis()->SetBinLabel(18,"kSemiCentral");
692       hTriggerClassSelected->GetXaxis()->SetBinLabel(19,"kDG5");
693       hTriggerClassSelected->GetXaxis()->SetBinLabel(20,"kZED");
694       hTriggerClassSelected->GetXaxis()->SetBinLabel(21,"kSPI7/kSPI");
695       hTriggerClassSelected->GetXaxis()->SetBinLabel(22,"kINT8");
696       hTriggerClassSelected->GetXaxis()->SetBinLabel(23,"kMuonSingleLowPt8");
697       hTriggerClassSelected->GetXaxis()->SetBinLabel(24,"kMuonSingleHighPt8");
698       hTriggerClassSelected->GetXaxis()->SetBinLabel(25,"kMuonLikeLowPt8");
699       hTriggerClassSelected->GetXaxis()->SetBinLabel(26,"kMuonUnlikeLowPt8");
700       hTriggerClassSelected->GetXaxis()->SetBinLabel(27,"kMuonUnlikeLowPt0");
701       hTriggerClassSelected->GetXaxis()->SetBinLabel(28,"kUserDefined");
702       hTriggerClassSelected->GetXaxis()->SetBinLabel(29,"kTRD");
703       hTriggerClassSelected->GetXaxis()->SetBinLabel(30,"kFastOnly");
704       hTriggerClassSelected->GetXaxis()->SetBinLabel(31,"kAnyINT");
705       hTriggerClassSelected->GetXaxis()->SetBinLabel(32,"kAny");
706       hTriggerClassSelected->GetXaxis()->SetBinLabel(33,"V0AND");
707       hTriggerClassSelected->GetXaxis()->SetBinLabel(34,"NOT kFastOnly");
708       fHistograms->Add(hTriggerClassSelected);
709       
710       hEventPlanePhi=new TH1F(Form("EventPlaneMinusPhotonAngle %s",GetCutNumber().Data()),"EventPlaneMinusPhotonAngle",360,-TMath::Pi(),TMath::Pi());
711       fHistograms->Add(hEventPlanePhi);
712
713       
714    }
715    TH1::AddDirectory(kTRUE);
716 }
717
718 //________________________________________________________________________
719 Bool_t AliConversionCuts::InitPIDResponse(){
720    // Set Pointer to AliPIDResponse
721
722    AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
723    if(man) {
724       AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
725       fPIDResponse = (AliPIDResponse*)inputHandler->GetPIDResponse();
726       if(fPIDResponse)return kTRUE;
727
728    }
729
730
731    return kFALSE;
732 }
733 ///________________________________________________________________________
734 Bool_t AliConversionCuts::EventIsSelected(AliVEvent *fInputEvent, AliVEvent *fMCEvent){
735    // Process Event Selection
736
737    Int_t cutindex=0;
738    if(hV0EventCuts)hV0EventCuts->Fill(cutindex);
739    cutindex++;
740
741    // Check for MC event
742    if(fMCEvent && fInputEvent->IsA()==AliESDEvent::Class()){
743       // Check if MC event is correctly loaded
744       AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
745       if (!mcHandler){
746          fEventQuality = 2;
747          return kFALSE;
748       }
749       if (!mcHandler->InitOk() ){
750          fEventQuality = 2;
751          return kFALSE;
752       }
753       if (!mcHandler->TreeK() ){
754          fEventQuality = 2;
755          return kFALSE;
756       }
757       if (!mcHandler->TreeTR() ) {
758          fEventQuality = 2;
759          return kFALSE;
760       }
761    }
762
763    // Event Trigger
764 //    cout << "before event trigger" << endl;
765    if(!IsTriggerSelected(fInputEvent)){
766       if(hV0EventCuts)hV0EventCuts->Fill(cutindex);
767       fEventQuality = 3;
768       return kFALSE;
769    }
770    cutindex++;
771
772    if(fInputEvent->IsA()==AliESDEvent::Class()){
773       AliTriggerAnalysis fTriggerAnalysis;// = new AliTriggerAnalysis;
774       fHasV0AND = fTriggerAnalysis.IsOfflineTriggerFired((AliESDEvent*)fInputEvent, AliTriggerAnalysis::kV0AND);
775       if(fHasV0AND&&hTriggerClass)hTriggerClass->Fill(32);
776    }
777 //   cout << "event number " << ((AliESDEvent*)fInputEvent)->GetEventNumberInFile() << " entered"<< endl;
778
779
780    // Number of Contributors Cut
781    if(GetNumberOfContributorsVtx(fInputEvent)<=0) {
782       if(hV0EventCuts)hV0EventCuts->Fill(cutindex);
783       fEventQuality = 5;
784       return kFALSE;
785    }
786    cutindex++;
787
788    // Z Vertex Position Cut
789    if(!VertexZCut(fInputEvent)){
790       if(hV0EventCuts)hV0EventCuts->Fill(cutindex);
791       fEventQuality = 4;
792       return kFALSE;
793    }
794    cutindex++;
795
796    // Pile Up Rejection
797
798    if(fRemovePileUp){
799       if(fInputEvent->IsPileupFromSPD(3,0.8,3.,2.,5.)){
800          if(hV0EventCuts)hV0EventCuts->Fill(cutindex);
801          fEventQuality = 6;
802          return kFALSE;
803       }
804    }
805    cutindex++;
806
807    // Centrality Selection
808    if(!IsCentralitySelected(fInputEvent,fMCEvent)){
809       if(hV0EventCuts)hV0EventCuts->Fill(cutindex);
810       fEventQuality = 1;
811       return kFALSE;
812    }
813    cutindex++;
814
815    // Fill Event Histograms
816    if(hV0EventCuts)hV0EventCuts->Fill(cutindex);
817    if(hVertexZ)hVertexZ->Fill(fInputEvent->GetPrimaryVertex()->GetZ());
818    if(hCentrality)hCentrality->Fill(GetCentrality(fInputEvent));
819    if(hCentralityVsNumberOfPrimaryTracks)
820       hCentralityVsNumberOfPrimaryTracks->Fill(GetCentrality(fInputEvent),
821                                                ((AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()
822                                                 ->GetTask("V0ReaderV1"))->GetNumberOfPrimaryTracks());
823    fEventQuality = 0;
824    return kTRUE;
825 }
826
827 ///________________________________________________________________________
828 Bool_t AliConversionCuts::PhotonIsSelectedMC(TParticle *particle,AliStack *fMCStack,Bool_t checkForConvertedGamma){
829    // MonteCarlo Photon Selection
830
831    if(!fMCStack)return kFALSE;
832
833    if (particle->GetPdgCode() == 22){
834
835
836       if( particle->Eta() > (fEtaCut + fEtaShift) || particle->Eta() < (-fEtaCut + fEtaShift) )
837          return kFALSE;
838       if(fEtaCutMin>-0.1){
839          if( particle->Eta() < (fEtaCutMin + fEtaShift) && particle->Eta() > (-fEtaCutMin + fEtaShift) )
840             return kFALSE;
841       }
842
843       if(particle->GetMother(0) >-1 && fMCStack->Particle(particle->GetMother(0))->GetPdgCode() == 22){
844          return kFALSE; // no photon as mothers!
845       }
846
847       if(particle->GetMother(0) >= fMCStack->GetNprimary()){
848          return kFALSE; // the gamma has a mother, and it is not a primary particle
849       }
850
851       if(!checkForConvertedGamma) return kTRUE; // return in case of accepted gamma
852
853       // looking for conversion gammas (electron + positron from pairbuilding (= 5) )
854       TParticle* ePos = NULL;
855       TParticle* eNeg = NULL;
856
857       if(particle->GetNDaughters() >= 2){
858          for(Int_t daughterIndex=particle->GetFirstDaughter();daughterIndex<=particle->GetLastDaughter();daughterIndex++){
859             TParticle *tmpDaughter = fMCStack->Particle(daughterIndex);
860             if(tmpDaughter->GetUniqueID() == 5){
861                if(tmpDaughter->GetPdgCode() == 11){
862                   eNeg = tmpDaughter;
863                } else if(tmpDaughter->GetPdgCode() == -11){
864                   ePos = tmpDaughter;
865                }
866             }
867          }
868       }
869
870       if(ePos == NULL || eNeg == NULL){ // means we do not have two daughters from pair production
871          return kFALSE;
872       }
873
874       if(ePos->Pt()<fSinglePtCut || eNeg->Pt()<fSinglePtCut){
875          return kFALSE; // no reconstruction below the Pt cut
876       }
877
878       if( ePos->Eta() > (fEtaCut + fEtaShift) || ePos->Eta() < (-fEtaCut + fEtaShift) ||
879           eNeg->Eta() > (fEtaCut + fEtaShift) || eNeg->Eta() < (-fEtaCut + fEtaShift) )
880          return kFALSE;
881
882       if(fEtaCutMin > -0.1){
883          if( (ePos->Eta() < (fEtaCutMin + fEtaShift) && ePos->Eta() > (-fEtaCutMin + fEtaShift)) ||
884              (eNeg->Eta() < (fEtaCutMin + fEtaShift) && eNeg->Eta() > (-fEtaCutMin + fEtaShift)) )
885             return kFALSE;
886       }
887
888       if(ePos->R()>fMaxR){
889          return kFALSE; // cuts on distance from collision point
890       }
891
892       if(abs(ePos->Vz()) > fMaxZ){
893          return kFALSE;  // outside material
894       }
895       if(abs(eNeg->Vz()) > fMaxZ){
896          return kFALSE;  // outside material
897       }
898
899       if( ePos->R() <= ((abs(ePos->Vz()) * fLineCutZRSlope) - fLineCutZValue)){
900          return kFALSE;  // line cut to exclude regions where we do not reconstruct
901       } else if ( fEtaCutMin != -0.1 &&   ePos->R() >= ((abs(ePos->Vz()) * fLineCutZRSlopeMin) - fLineCutZValueMin)){
902          return kFALSE;
903       }
904
905       if( eNeg->R() <= ((abs(eNeg->Vz()) * fLineCutZRSlope) - fLineCutZValue)){
906          return kFALSE; // line cut to exclude regions where we do not reconstruct
907       } else if ( fEtaCutMin != -0.1 &&   eNeg->R() >= ((abs(eNeg->Vz()) * fLineCutZRSlopeMin) - fLineCutZValueMin)){
908          return kFALSE;
909       }
910
911       return kTRUE;
912       //if(AcceptanceCut(particle,ePos,eNeg))return kTRUE;
913    }
914    return kFALSE;
915 }
916 ///________________________________________________________________________
917 Bool_t AliConversionCuts::PhotonIsSelectedAODMC(AliAODMCParticle *particle,TClonesArray *aodmcArray,Bool_t checkForConvertedGamma){
918    // MonteCarlo Photon Selection
919
920    if(!aodmcArray)return kFALSE;
921
922    if (particle->GetPdgCode() == 22){
923       if( particle->Eta() > (fEtaCut + fEtaShift) || particle->Eta() < (-fEtaCut + fEtaShift) )
924          return kFALSE;
925       if(fEtaCutMin>-0.1){
926          if( particle->Eta() < (fEtaCutMin + fEtaShift) && particle->Eta() > (-fEtaCutMin + fEtaShift) )
927             return kFALSE;
928       }
929
930       if(particle->GetMother() > -1){
931          if((static_cast<AliAODMCParticle*>(aodmcArray->At(particle->GetMother())))->GetPdgCode() == 22){
932             return kFALSE; // no photon as mothers!
933          }
934          if(!(static_cast<AliAODMCParticle*>(aodmcArray->At(particle->GetMother()))->IsPrimary())){
935             return kFALSE; // the gamma has a mother, and it is not a primary particle
936          }
937       }
938
939       if(!checkForConvertedGamma) return kTRUE; // return in case of accepted gamma
940
941       // looking for conversion gammas (electron + positron from pairbuilding (= 5) )
942       AliAODMCParticle* ePos = NULL;
943       AliAODMCParticle* eNeg = NULL;
944
945       if(particle->GetNDaughters() >= 2){
946          for(Int_t daughterIndex=particle->GetDaughter(0);daughterIndex<=particle->GetDaughter(1);daughterIndex++){
947             AliAODMCParticle *tmpDaughter = static_cast<AliAODMCParticle*>(aodmcArray->At(daughterIndex));
948             if(!tmpDaughter) continue;
949             if(((tmpDaughter->GetMCProcessCode())) == 5){    // STILL A BUG IN ALIROOT >>8 HAS TPO BE REMOVED AFTER FIX
950                if(tmpDaughter->GetPdgCode() == 11){
951                   eNeg = tmpDaughter;
952                } else if(tmpDaughter->GetPdgCode() == -11){
953                   ePos = tmpDaughter;
954                }
955             }
956          }
957       }
958
959       if(ePos == NULL || eNeg == NULL){ // means we do not have two daughters from pair production
960          return kFALSE;
961       }
962
963       if(ePos->Pt()<fSinglePtCut || eNeg->Pt()<fSinglePtCut){
964          return kFALSE; // no reconstruction below the Pt cut
965       }
966
967       if( ePos->Eta() > (fEtaCut + fEtaShift) || ePos->Eta() < (-fEtaCut + fEtaShift) ||
968           eNeg->Eta() > (fEtaCut + fEtaShift) || eNeg->Eta() < (-fEtaCut + fEtaShift) )
969          return kFALSE;
970
971       if(fEtaCutMin > -0.1){
972          if( (ePos->Eta() < (fEtaCutMin + fEtaShift) && ePos->Eta() > (-fEtaCutMin + fEtaShift)) ||
973              (eNeg->Eta() < (fEtaCutMin + fEtaShift) && eNeg->Eta() > (-fEtaCutMin + fEtaShift)) )
974             return kFALSE;
975       }
976
977       Double_t rPos = sqrt( (ePos->Xv()*ePos->Xv()) + (ePos->Yv()*ePos->Yv()) );
978       Double_t rNeg = sqrt( (eNeg->Xv()*eNeg->Xv()) + (eNeg->Yv()*eNeg->Yv()) );
979
980       if(rPos>fMaxR){
981          return kFALSE; // cuts on distance from collision point
982       }
983       if(abs(ePos->Zv()) > fMaxZ){
984          return kFALSE;  // outside material
985       }
986       if(abs(eNeg->Zv()) > fMaxZ){
987          return kFALSE;  // outside material
988       }
989
990       if( rPos <= ((abs(ePos->Zv()) * fLineCutZRSlope) - fLineCutZValue)){
991          return kFALSE;  // line cut to exclude regions where we do not reconstruct
992       } else if ( fEtaCutMin != -0.1 &&   rPos >= ((abs(ePos->Zv()) * fLineCutZRSlopeMin) - fLineCutZValueMin)){
993          return kFALSE;
994       }
995
996       if( rNeg <= ((abs(eNeg->Zv()) * fLineCutZRSlope) - fLineCutZValue)){
997          return kFALSE; // line cut to exclude regions where we do not reconstruct
998       } else if ( fEtaCutMin != -0.1 &&   rNeg >= ((abs(eNeg->Zv()) * fLineCutZRSlopeMin) - fLineCutZValueMin)){
999          return kFALSE;
1000       }
1001
1002       return kTRUE;
1003       //if(AcceptanceCut(particle,ePos,eNeg))return kTRUE;
1004    }
1005    return kFALSE;
1006 }
1007 ///________________________________________________________________________
1008 Bool_t AliConversionCuts::PhotonCuts(AliConversionPhotonBase *photon,AliVEvent *event)
1009 {   // Specific Photon Cuts
1010
1011    Int_t cutIndex = 0;
1012    if(hPhotonCuts)hPhotonCuts->Fill(cutIndex);
1013    cutIndex++;
1014
1015    // Fill Histos before Cuts
1016    if(hInvMassbefore)hInvMassbefore->Fill(photon->GetMass());
1017    if(hArmenterosbefore)hArmenterosbefore->Fill(photon->GetArmenterosAlpha(),photon->GetArmenterosQt());
1018
1019    // Gamma selection based on QT from Armenteros
1020    if(fDoQtGammaSelection == kTRUE){
1021       if(!ArmenterosQtCut(photon)){
1022          if(hPhotonCuts)hPhotonCuts->Fill(cutIndex); //1
1023          return kFALSE;
1024       }
1025    }
1026    cutIndex++; //2
1027
1028    // Chi Cut
1029    if(photon->GetChi2perNDF() > fChi2CutConversion || photon->GetChi2perNDF() <=0){
1030       {
1031          if(hPhotonCuts)hPhotonCuts->Fill(cutIndex); //2
1032          return kFALSE;
1033       }
1034    }
1035    cutIndex++;//3
1036
1037    // Reconstruction Acceptance Cuts
1038    if(!AcceptanceCuts(photon)){
1039       if(hPhotonCuts)hPhotonCuts->Fill(cutIndex); //3
1040       return kFALSE;
1041    }
1042
1043    cutIndex++; //4
1044    // Asymmetry Cut
1045    if(fDoPhotonAsymmetryCut == kTRUE){
1046       if(!AsymmetryCut(photon,event)){
1047          if(hPhotonCuts)hPhotonCuts->Fill(cutIndex); //4
1048          return kFALSE;
1049       }
1050    }
1051
1052    //Check the pid probability
1053    cutIndex++; //5
1054    if(!PIDProbabilityCut(photon, event)) {
1055       if(hPhotonCuts)hPhotonCuts->Fill(cutIndex); //5
1056       return kFALSE;
1057    }
1058
1059    cutIndex++; //6
1060    if(!CorrectedTPCClusterCut(photon, event)) {
1061       if(hPhotonCuts)hPhotonCuts->Fill(cutIndex); //6
1062       return kFALSE;
1063    }
1064
1065    Double_t magField = event->GetMagneticField();
1066    if( magField  < 0.0 ){
1067       magField =  1.0;
1068    } else {
1069       magField =  -1.0;
1070    }
1071    
1072    AliVTrack * electronCandidate = GetTrack(event,photon->GetTrackLabelNegative() );
1073    AliVTrack * positronCandidate = GetTrack(event,photon->GetTrackLabelPositive() );
1074    Double_t deltaPhi = magField * TVector2::Phi_mpi_pi( electronCandidate->Phi()-positronCandidate->Phi());
1075
1076    cutIndex++; //7
1077    if(!PsiPairCut(photon)) {
1078       if(hPhotonCuts)hPhotonCuts->Fill(cutIndex); //7
1079       return kFALSE;
1080    }
1081
1082    cutIndex++; //8
1083    if(!CosinePAngleCut(photon, event)) {
1084       if(hPhotonCuts)hPhotonCuts->Fill(cutIndex); //8
1085       return kFALSE;
1086    }
1087
1088    AliAODConversionPhoton* photonAOD = dynamic_cast<AliAODConversionPhoton*>(photon);
1089    if (photonAOD){
1090       photonAOD->CalculateDistanceOfClossetApproachToPrimVtx(event->GetPrimaryVertex());
1091
1092       cutIndex++; //9
1093       if(photonAOD->GetDCArToPrimVtx() > fDCARPrimVtxCut) { //DCA R cut of photon to primary vertex
1094          if(hPhotonCuts)hPhotonCuts->Fill(cutIndex); //9
1095          return kFALSE;
1096       }
1097
1098       cutIndex++; //10
1099       if(abs(photonAOD->GetDCAzToPrimVtx()) > fDCAZPrimVtxCut) { //DCA Z cut of photon to primary vertex
1100          if(hPhotonCuts)hPhotonCuts->Fill(cutIndex); //10
1101          return kFALSE;
1102       }
1103    } else {
1104       cutIndex++; //9
1105       cutIndex++; //10
1106    }
1107    cutIndex++; //11
1108
1109    if (photonAOD){
1110                 if (fDoPhotonQualitySelectionCut && photonAOD->GetPhotonQuality() != fPhotonQualityCut){
1111                         if(hPhotonCuts)hPhotonCuts->Fill(cutIndex); //11
1112                         return kFALSE;
1113                 }       
1114    } 
1115    cutIndex++; //12
1116    if(hPhotonCuts)hPhotonCuts->Fill(cutIndex); //11
1117
1118    // Histos after Cuts
1119    if(hInvMassafter)hInvMassafter->Fill(photon->GetMass());
1120    if(hArmenterosafter)hArmenterosafter->Fill(photon->GetArmenterosAlpha(),photon->GetArmenterosQt());
1121    if(hPsiPairDeltaPhiafter)hPsiPairDeltaPhiafter->Fill(deltaPhi,photon->GetPsiPair());
1122    return kTRUE;
1123
1124 }
1125
1126 ///________________________________________________________________________
1127 Bool_t AliConversionCuts::CorrectedTPCClusterCut(AliConversionPhotonBase *photon, AliVEvent * event)
1128 {   //Cut on corrected TPC Cluster Info
1129
1130    AliVTrack * negTrack = GetTrack(event, photon->GetTrackLabelNegative());
1131    AliVTrack * posTrack = GetTrack(event, photon->GetTrackLabelPositive());
1132
1133    if(!negTrack||!posTrack)return kFALSE;
1134
1135    Double_t negclsToF=0;
1136
1137    if (!fUseCorrectedTPCClsInfo ){
1138       if(negTrack->GetTPCNclsF()!=0){
1139          negclsToF = (Double_t)negTrack->GetNcls(1)/(Double_t)negTrack->GetTPCNclsF();}// Ncluster/Nfindablecluster
1140    }
1141    else {
1142       negclsToF = negTrack->GetTPCClusterInfo(2,0,GetFirstTPCRow(photon->GetConversionRadius()));
1143    }
1144
1145    Double_t posclsToF = 0.;
1146    if (!fUseCorrectedTPCClsInfo ){
1147       if(posTrack->GetTPCNclsF()!=0){
1148          posclsToF = (Double_t)posTrack->GetNcls(1)/(Double_t)posTrack->GetTPCNclsF();
1149       }
1150    }else{
1151       posclsToF = posTrack->GetTPCClusterInfo(2,0,GetFirstTPCRow(photon->GetConversionRadius()));
1152    }
1153
1154    if( negclsToF < fMinClsTPCToF || posclsToF < fMinClsTPCToF ){
1155       return kFALSE;
1156    }
1157
1158    return kTRUE;
1159 }
1160
1161 ///________________________________________________________________________
1162 Bool_t AliConversionCuts::PhotonIsSelected(AliConversionPhotonBase *photon, AliVEvent * event)
1163 {
1164    //Selection of Reconstructed Photons
1165
1166    FillPhotonCutIndex(kPhotonIn);
1167
1168    if(event->IsA()==AliESDEvent::Class()) {
1169      if(!SelectV0Finder( ( ((AliESDEvent*)event)->GetV0(photon->GetV0Index()))->GetOnFlyStatus() ) ){
1170          FillPhotonCutIndex(kOnFly);
1171          return kFALSE;
1172       }
1173    }
1174    // else if(event->IsA()==AliAODEvent::Class()) {
1175    //    if(!SelectV0Finder( ( ((AliAODEvent*)event)->GetV0(photon->GetV0Index())) ) ){
1176    //       FillPhotonCutIndex(kOnFly);
1177    //       return kFALSE;
1178    //    }
1179    // }
1180
1181    // Get Tracks
1182    AliVTrack * negTrack = GetTrack(event, photon->GetTrackLabelNegative());
1183    AliVTrack * posTrack = GetTrack(event, photon->GetTrackLabelPositive());
1184
1185    if(!negTrack || !posTrack) {
1186       FillPhotonCutIndex(kNoTracks);
1187       return kFALSE;
1188    }
1189    photon->DeterminePhotonQuality(negTrack,posTrack);
1190    // Track Cuts
1191    if(!TracksAreSelected(negTrack, posTrack)){
1192       FillPhotonCutIndex(kTrackCuts);
1193       return kFALSE;
1194    }
1195    if (hEtaDistV0s)hEtaDistV0s->Fill(photon->GetPhotonEta());
1196    // dEdx Cuts
1197    if(!dEdxCuts(negTrack) || !dEdxCuts(posTrack)) {
1198       FillPhotonCutIndex(kdEdxCuts);
1199       return kFALSE;
1200    }
1201    if (hEtaDistV0sAfterdEdxCuts)hEtaDistV0sAfterdEdxCuts->Fill(photon->GetPhotonEta());
1202    // Photon Cuts
1203    if(!PhotonCuts(photon,event)){
1204       FillPhotonCutIndex(kPhotonCuts);
1205       return kFALSE;
1206    }
1207
1208    // Photon passed cuts
1209    FillPhotonCutIndex(kPhotonOut);
1210    return kTRUE;
1211 }
1212
1213 ///________________________________________________________________________
1214 Bool_t AliConversionCuts::ArmenterosQtCut(AliConversionPhotonBase *photon)
1215 {   // Armenteros Qt Cut
1216    if(fDo2DQt){
1217       if ( !(TMath::Power(photon->GetArmenterosAlpha()/0.95,2)+TMath::Power(photon->GetArmenterosQt()/fQtMax,2) < 1) ){
1218          return kFALSE;
1219       }
1220    } else {
1221       if(photon->GetArmenterosQt()>fQtMax){
1222          return kFALSE;
1223       }
1224    }
1225    return kTRUE;
1226 }
1227
1228
1229 ///________________________________________________________________________
1230 Bool_t AliConversionCuts::AcceptanceCuts(AliConversionPhotonBase *photon) {
1231    // Exclude certain areas for photon reconstruction
1232
1233    Int_t cutIndex=0;
1234    if(hAcceptanceCuts)hAcceptanceCuts->Fill(cutIndex);
1235    cutIndex++;
1236
1237    if(photon->GetConversionRadius()>fMaxR){ // cuts on distance from collision point
1238       if(hAcceptanceCuts)hAcceptanceCuts->Fill(cutIndex);
1239       return kFALSE;
1240    }
1241    cutIndex++;
1242
1243    if(photon->GetConversionRadius()<fMinR){ // cuts on distance from collision point
1244       if(hAcceptanceCuts)hAcceptanceCuts->Fill(cutIndex);
1245       return kFALSE;
1246    }
1247    cutIndex++;
1248
1249    if(photon->GetConversionRadius() <= ((abs(photon->GetConversionZ())*fLineCutZRSlope)-fLineCutZValue)){
1250       if(hAcceptanceCuts)hAcceptanceCuts->Fill(cutIndex);
1251       return kFALSE;
1252    }
1253    else if (fUseEtaMinCut &&  photon->GetConversionRadius() >= ((abs(photon->GetConversionZ())*fLineCutZRSlopeMin)-fLineCutZValueMin )){
1254       if(hAcceptanceCuts)hAcceptanceCuts->Fill(cutIndex);
1255       return kFALSE;
1256    }
1257    cutIndex++;
1258
1259    if(abs(photon->GetConversionZ()) > fMaxZ ){ // cuts out regions where we do not reconstruct
1260       if(hAcceptanceCuts)hAcceptanceCuts->Fill(cutIndex);
1261       return kFALSE;
1262    }
1263    cutIndex++;
1264
1265
1266    if( photon->GetPhotonEta() > (fEtaCut + fEtaShift)    || photon->GetPhotonEta() < (-fEtaCut + fEtaShift) ){
1267       if(hAcceptanceCuts)hAcceptanceCuts->Fill(cutIndex);
1268       return kFALSE;
1269    }
1270    if(fEtaCutMin>-0.1){
1271       if( photon->GetPhotonEta() < (fEtaCutMin + fEtaShift) && photon->GetPhotonEta() > (-fEtaCutMin + fEtaShift) ){
1272          if(hAcceptanceCuts)hAcceptanceCuts->Fill(cutIndex);
1273          return kFALSE;
1274       }
1275    }
1276    cutIndex++;
1277
1278    if(photon->GetPhotonPt()<fPtCut){
1279       if(hAcceptanceCuts)hAcceptanceCuts->Fill(cutIndex);
1280       return kFALSE;
1281    }
1282    cutIndex++;
1283
1284    if(hAcceptanceCuts)hAcceptanceCuts->Fill(cutIndex);
1285
1286    return kTRUE;
1287 }
1288
1289
1290 ///________________________________________________________________________
1291 Bool_t AliConversionCuts::SpecificTrackCuts(AliAODTrack * negTrack, AliAODTrack * posTrack,Int_t &cutIndex) {
1292    // Track Cuts which require AOD/ESD specific implementation
1293
1294    if( !negTrack->IsOn(AliESDtrack::kTPCrefit)  || !posTrack->IsOn(AliESDtrack::kTPCrefit)   )  {
1295       if(hTrackCuts)hTrackCuts->Fill(cutIndex);
1296       return kFALSE;
1297    }
1298    cutIndex++;
1299
1300    AliAODVertex * NegVtxType=negTrack->GetProdVertex();
1301    AliAODVertex * PosVtxType=posTrack->GetProdVertex();
1302    if( (NegVtxType->GetType())==AliAODVertex::kKink || (PosVtxType->GetType())==AliAODVertex::kKink) {
1303       if(hTrackCuts)hTrackCuts->Fill(cutIndex);
1304       return kFALSE;
1305    }
1306    return kTRUE;
1307
1308 }
1309
1310
1311 ///________________________________________________________________________
1312 Bool_t AliConversionCuts::SpecificTrackCuts(AliESDtrack * negTrack, AliESDtrack * posTrack,Int_t &cutIndex) {
1313    // Track Cuts which require AOD/ESD specific implementation
1314
1315    if( !negTrack->IsOn(AliESDtrack::kTPCrefit)  || !posTrack->IsOn(AliESDtrack::kTPCrefit)   )  {
1316       if(hTrackCuts)hTrackCuts->Fill(cutIndex);
1317       return kFALSE;
1318    }
1319    cutIndex++;
1320
1321    if(negTrack->GetKinkIndex(0) > 0  || posTrack->GetKinkIndex(0) > 0 ) {
1322       if(hTrackCuts)hTrackCuts->Fill(cutIndex);
1323       return kFALSE;
1324    }
1325    return kTRUE;
1326 }
1327
1328
1329
1330 ///________________________________________________________________________
1331 Bool_t AliConversionCuts::TracksAreSelected(AliVTrack * negTrack, AliVTrack * posTrack) {
1332    // Track Selection for Photon Reconstruction
1333
1334    Int_t cutIndex=0;
1335    if(hTrackCuts)hTrackCuts->Fill(cutIndex);
1336    cutIndex++;
1337
1338    // avoid like sign
1339    if(negTrack->Charge() == posTrack->Charge()) {
1340       if(hTrackCuts)hTrackCuts->Fill(cutIndex);
1341       return kFALSE;
1342    }
1343    cutIndex++;
1344
1345    // Number of TPC Clusters
1346
1347
1348    if( negTrack->GetNcls(1) < fMinClsTPC || posTrack->GetNcls(1) < fMinClsTPC ) {
1349       if(hTrackCuts)hTrackCuts->Fill(cutIndex);
1350       return kFALSE;
1351    }
1352    cutIndex++;
1353
1354    // Acceptance
1355    if( posTrack->Eta() > (fEtaCut + fEtaShift) || posTrack->Eta() < (-fEtaCut + fEtaShift) ||
1356        negTrack->Eta() > (fEtaCut + fEtaShift) || negTrack->Eta() < (-fEtaCut + fEtaShift) ){
1357       if(hTrackCuts)hTrackCuts->Fill(cutIndex);
1358       return kFALSE;
1359    }
1360    if(fEtaCutMin>-0.1){
1361       if( (posTrack->Eta() < (fEtaCutMin + fEtaShift) && posTrack->Eta() > (-fEtaCutMin + fEtaShift)) ||
1362           (negTrack->Eta() < (fEtaCutMin + fEtaShift) && negTrack->Eta() > (-fEtaCutMin + fEtaShift)) ){
1363          if(hTrackCuts)hTrackCuts->Fill(cutIndex);
1364          return kFALSE;
1365       }
1366    }
1367    cutIndex++;
1368
1369    // Single Pt Cut
1370    if( negTrack->Pt()< fSinglePtCut || posTrack->Pt()< fSinglePtCut){
1371       if(hTrackCuts)hTrackCuts->Fill(cutIndex);
1372       return kFALSE;
1373    }
1374    cutIndex++;
1375
1376    // AOD ESD specific cuts
1377    Bool_t passCuts = kTRUE;
1378
1379    if(negTrack->IsA()==AliAODTrack::Class()) {
1380       passCuts = passCuts * SpecificTrackCuts(static_cast<AliAODTrack*>(negTrack), static_cast<AliAODTrack*>(posTrack),cutIndex);
1381    } else {
1382       passCuts = passCuts * SpecificTrackCuts(static_cast<AliESDtrack*>(negTrack), static_cast<AliESDtrack*>(posTrack),cutIndex);
1383    }
1384
1385    if(!passCuts){
1386       if(hTrackCuts)hTrackCuts->Fill(cutIndex);
1387       return kFALSE;
1388    }
1389    cutIndex++;
1390
1391    if(hTrackCuts)hTrackCuts->Fill(cutIndex);
1392
1393    return kTRUE;
1394
1395 }
1396
1397 ///________________________________________________________________________
1398 Bool_t AliConversionCuts::dEdxCuts(AliVTrack *fCurrentTrack){
1399    // Electron Identification Cuts for Photon reconstruction
1400    if(!fPIDResponse){InitPIDResponse();}// Try to reinitialize PID Response
1401    if(!fPIDResponse){AliError("No PID Response"); return kTRUE;}// if still missing fatal error
1402
1403    Int_t cutIndex=0;
1404    if(hdEdxCuts)hdEdxCuts->Fill(cutIndex);
1405    if(hTPCdEdxSigbefore)hTPCdEdxSigbefore->Fill(fCurrentTrack->P(),fPIDResponse->NumberOfSigmasTPC(fCurrentTrack, AliPID::kElectron));
1406    if(hTPCdEdxbefore)hTPCdEdxbefore->Fill(fCurrentTrack->P(),fCurrentTrack->GetTPCsignal());
1407    cutIndex++;
1408
1409    if(fDodEdxSigmaCut == kTRUE){
1410       // TPC Electron Line
1411       if( fPIDResponse->NumberOfSigmasTPC(fCurrentTrack,AliPID::kElectron)<fPIDnSigmaBelowElectronLine ||
1412           fPIDResponse->NumberOfSigmasTPC(fCurrentTrack,AliPID::kElectron)>fPIDnSigmaAboveElectronLine){
1413
1414          if(hdEdxCuts)hdEdxCuts->Fill(cutIndex);
1415          return kFALSE;
1416       }
1417       cutIndex++;
1418
1419       // TPC Pion Line
1420       if( fCurrentTrack->P()>fPIDMinPnSigmaAbovePionLine && fCurrentTrack->P()<fPIDMaxPnSigmaAbovePionLine ){
1421          if(fPIDResponse->NumberOfSigmasTPC(fCurrentTrack,AliPID::kElectron)>fPIDnSigmaBelowElectronLine &&
1422             fPIDResponse->NumberOfSigmasTPC(fCurrentTrack,AliPID::kElectron)<fPIDnSigmaAboveElectronLine&&
1423             fPIDResponse->NumberOfSigmasTPC(fCurrentTrack,AliPID::kPion)<fPIDnSigmaAbovePionLine){
1424
1425             if(hdEdxCuts)hdEdxCuts->Fill(cutIndex);
1426             return kFALSE;
1427          }
1428       }
1429       cutIndex++;
1430
1431       // High Pt Pion rej
1432       if( fCurrentTrack->P()>fPIDMaxPnSigmaAbovePionLine ){
1433          if(fPIDResponse->NumberOfSigmasTPC(fCurrentTrack,AliPID::kElectron)>fPIDnSigmaBelowElectronLine &&
1434             fPIDResponse->NumberOfSigmasTPC(fCurrentTrack,AliPID::kElectron)<fPIDnSigmaAboveElectronLine &&
1435             fPIDResponse->NumberOfSigmasTPC(fCurrentTrack,AliPID::kPion)<fPIDnSigmaAbovePionLineHighPt){
1436
1437             if(hdEdxCuts)hdEdxCuts->Fill(cutIndex);
1438             return kFALSE;
1439          }
1440       }
1441       cutIndex++;
1442    }
1443    else{cutIndex+=3;}
1444
1445    if(fDoKaonRejectionLowP == kTRUE){
1446       if(fCurrentTrack->P()<fPIDMinPKaonRejectionLowP ){
1447          if( abs(fPIDResponse->NumberOfSigmasTPC(fCurrentTrack,AliPID::kKaon))<fPIDnSigmaAtLowPAroundKaonLine){
1448
1449             if(hdEdxCuts)hdEdxCuts->Fill(cutIndex);
1450             return kFALSE;
1451          }
1452       }
1453    }
1454    cutIndex++;
1455    if(fDoProtonRejectionLowP == kTRUE){
1456       if( fCurrentTrack->P()<fPIDMinPProtonRejectionLowP ){
1457          if( abs(fPIDResponse->NumberOfSigmasTPC(fCurrentTrack,AliPID::kProton))<fPIDnSigmaAtLowPAroundProtonLine){
1458
1459             if(hdEdxCuts)hdEdxCuts->Fill(cutIndex);
1460             return kFALSE;
1461          }
1462       }
1463    }
1464    cutIndex++;
1465
1466    if(fDoPionRejectionLowP == kTRUE){
1467       if( fCurrentTrack->P()<fPIDMinPPionRejectionLowP ){
1468          if( abs(fPIDResponse->NumberOfSigmasTPC(fCurrentTrack,AliPID::kPion))<fPIDnSigmaAtLowPAroundPionLine){
1469
1470             if(hdEdxCuts)hdEdxCuts->Fill(cutIndex);
1471             return kFALSE;
1472          }
1473       }
1474    }
1475    cutIndex++;
1476
1477
1478    // cout<<"Start"<<endl;
1479    // AliPIDResponse::EDetPidStatus status=fPIDResponse->CheckPIDStatus(AliPIDResponse::kTOF,fCurrentTrack);
1480
1481    // if( ( (status & AliVTrack::kTOFout) == AliVTrack::kTOFout ) && ( (status & AliVTrack::kTIME) == AliVTrack::kTIME ))
1482    //    {cout<<"TOF DA"<<endl;}
1483    // if(status == AliPIDResponse::kDetPidOk){
1484    //    Float_t probMis = fPIDResponse->GetTOFMismatchProbability(fCurrentTrack);
1485    //    cout<<"--> "<<probMis<<endl;
1486    //    if(probMis > 0.01){
1487
1488    //    }
1489    // }
1490
1491    if((fCurrentTrack->GetStatus() & AliESDtrack::kTOFpid) && !(fCurrentTrack->GetStatus() & AliESDtrack::kTOFmismatch)){
1492       if(hTOFbefore){
1493          Double_t t0 = fPIDResponse->GetTOFResponse().GetStartTime(fCurrentTrack->P());
1494          Double_t  times[AliPID::kSPECIESC];
1495          fCurrentTrack->GetIntegratedTimes(times);
1496          Double_t TOFsignal = fCurrentTrack->GetTOFsignal();
1497          Double_t dT = TOFsignal - t0 - times[0];
1498          hTOFbefore->Fill(fCurrentTrack->P(),dT);
1499       }
1500       if(hTOFSigbefore) hTOFSigbefore->Fill(fCurrentTrack->P(),fPIDResponse->NumberOfSigmasTOF(fCurrentTrack, AliPID::kElectron));
1501       if(fUseTOFpid){
1502          if(fPIDResponse->NumberOfSigmasTOF(fCurrentTrack, AliPID::kElectron)>fTofPIDnSigmaAboveElectronLine ||
1503             fPIDResponse->NumberOfSigmasTOF(fCurrentTrack, AliPID::kElectron)<fTofPIDnSigmaBelowElectronLine ){
1504             if(hdEdxCuts)hdEdxCuts->Fill(cutIndex);
1505             return kFALSE;
1506          }
1507       }
1508       if(hTOFSigafter)hTOFSigafter->Fill(fCurrentTrack->P(),fPIDResponse->NumberOfSigmasTOF(fCurrentTrack, AliPID::kElectron));
1509    }
1510    cutIndex++;
1511    // Apply TRD PID
1512    if(fDoTRDPID){
1513       if(!fPIDResponse->IdentifiedAsElectronTRD(fCurrentTrack,fPIDTRDEfficiency)){
1514          if(hdEdxCuts)hdEdxCuts->Fill(cutIndex);
1515          return kFALSE;
1516       }
1517    }
1518    cutIndex++;
1519
1520    if(hdEdxCuts)hdEdxCuts->Fill(cutIndex);
1521    if(hTPCdEdxSigafter)hTPCdEdxSigafter->Fill(fCurrentTrack->P(),fPIDResponse->NumberOfSigmasTPC(fCurrentTrack, AliPID::kElectron));
1522    if(hTPCdEdxafter)hTPCdEdxafter->Fill(fCurrentTrack->P(),fCurrentTrack->GetTPCsignal());
1523    
1524    return kTRUE;
1525 }
1526
1527 ///________________________________________________________________________
1528 Bool_t AliConversionCuts::AsymmetryCut(AliConversionPhotonBase * photon,AliVEvent *event) {
1529    // Cut on Energy Assymetry
1530
1531    for(Int_t ii=0;ii<2;ii++){
1532
1533       AliVTrack *track=GetTrack(event,photon->GetTrackLabel(ii));
1534
1535       if( track->P() > fMinPPhotonAsymmetryCut ){
1536          Double_t trackNegAsy=0;
1537          if (photon->GetPhotonP()!=0.){
1538             trackNegAsy= track->P()/photon->GetPhotonP();
1539          }
1540
1541          if( trackNegAsy<fMinPhotonAsymmetry ||trackNegAsy>(1.- fMinPhotonAsymmetry)){
1542             return kFALSE;
1543          }
1544       }
1545
1546    }
1547    return kTRUE;
1548 }
1549
1550 ///________________________________________________________________________
1551 AliVTrack *AliConversionCuts::GetTrack(AliVEvent * event, Int_t label){
1552    //Returns pointer to the track with given ESD label
1553    //(Important for AOD implementation, since Track array in AOD data is different
1554    //from ESD array, but ESD tracklabels are stored in AOD Tracks)
1555
1556    AliESDEvent * esdEvent = dynamic_cast<AliESDEvent*>(event);
1557    if(esdEvent) {
1558       if(label > event->GetNumberOfTracks() ) return NULL;
1559       AliESDtrack * track = esdEvent->GetTrack(label);
1560       return track;
1561
1562    } else {
1563       AliVTrack * track = 0x0;
1564       if(((AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()->GetTask("V0ReaderV1"))->AreAODsRelabeled()){
1565          track = dynamic_cast<AliVTrack*>(event->GetTrack(label));
1566          return track;
1567       }
1568       else{
1569          for(Int_t ii=0; ii<event->GetNumberOfTracks(); ii++) {
1570             track = dynamic_cast<AliVTrack*>(event->GetTrack(ii));
1571             if(track){
1572                if(track->GetID() == label) {
1573                   return track;
1574                }
1575             }
1576          }
1577       }
1578    }
1579    //AliDebug(5,(Form("track not found %d %d",label,event->GetNumberOfTracks()));
1580    return NULL;
1581 }
1582
1583 ///________________________________________________________________________
1584 AliESDtrack *AliConversionCuts::GetESDTrack(AliESDEvent * event, Int_t label){
1585    //Returns pointer to the track with given ESD label
1586    //(Important for AOD implementation, since Track array in AOD data is different
1587    //from ESD array, but ESD tracklabels are stored in AOD Tracks)
1588
1589    if(event) {
1590       if(label > event->GetNumberOfTracks() ) return NULL;
1591       AliESDtrack * track = event->GetTrack(label);
1592       return track;
1593    }
1594    //AliDebug(5,(Form("track not found %d %d",label,event->GetNumberOfTracks()));
1595    return NULL;
1596 }
1597
1598
1599
1600 ///________________________________________________________________________
1601 Bool_t AliConversionCuts::PIDProbabilityCut(AliConversionPhotonBase *photon, AliVEvent * event){
1602    // Cut on Electron Probability for Photon Reconstruction
1603
1604    AliESDEvent * esdEvent = dynamic_cast<AliESDEvent*>(event);
1605
1606    if(esdEvent){
1607
1608       Bool_t iResult=kFALSE;
1609
1610       Double_t *posProbArray = new Double_t[AliPID::kSPECIES];
1611       Double_t *negProbArray = new Double_t[AliPID::kSPECIES];
1612
1613       AliESDtrack* negTrack   = esdEvent->GetTrack(photon->GetTrackLabelNegative());
1614       AliESDtrack* posTrack   = esdEvent->GetTrack(photon->GetTrackLabelPositive());
1615
1616       if(negProbArray && posProbArray){
1617
1618          negTrack->GetTPCpid(negProbArray);
1619          posTrack->GetTPCpid(posProbArray);
1620
1621          if(negProbArray[AliPID::kElectron]>=fPIDProbabilityCutNegativeParticle && posProbArray[AliPID::kElectron]>=fPIDProbabilityCutPositiveParticle){
1622             iResult=kTRUE;
1623          }
1624       }
1625
1626       delete [] posProbArray;
1627       delete [] negProbArray;
1628       return iResult;
1629
1630    } else {
1631       ///Not possible for AODs
1632       return kTRUE;
1633    }
1634
1635
1636
1637
1638 }
1639
1640
1641 ///________________________________________________________________________
1642 Bool_t AliConversionCuts::AcceptanceCut(TParticle *particle, TParticle * ePos,TParticle* eNeg){
1643    // MC Acceptance Cuts
1644    //(Certain areas were excluded for photon reconstruction)
1645
1646    if(particle->R()>fMaxR){
1647       return kFALSE;}
1648
1649    if(ePos->R()>fMaxR){
1650       return kFALSE;
1651    }
1652
1653    if(ePos->R()<fMinR){
1654       return kFALSE;
1655    }
1656
1657    if( ePos->R() <= ((abs(ePos->Vz())*fLineCutZRSlope)-fLineCutZValue)){
1658       return kFALSE;
1659    }
1660    else if (fUseEtaMinCut &&  ePos->R() >= ((abs(ePos->Vz())*fLineCutZRSlopeMin)-fLineCutZValueMin )){
1661       return kFALSE;
1662    }
1663
1664    if(abs(eNeg->Vz()) > fMaxZ ){ // cuts out regions where we do not reconstruct
1665       return kFALSE;
1666    }
1667
1668    if(eNeg->Vz()!=ePos->Vz()||eNeg->R()!=ePos->R()){
1669       return kFALSE;
1670    }
1671
1672    if(abs(ePos->Vz()) > fMaxZ ){ // cuts out regions where we do not reconstruct
1673       return kFALSE;
1674    }
1675
1676
1677    if( particle->Eta() > (fEtaCut + fEtaShift) || particle->Eta() < (-fEtaCut + fEtaShift) ){
1678       return kFALSE;
1679    }
1680    if( ePos->Eta() > (fEtaCut + fEtaShift) || ePos->Eta() < (-fEtaCut + fEtaShift) ){
1681       return kFALSE;
1682    }
1683    if( eNeg->Eta() > (fEtaCut + fEtaShift) || eNeg->Eta() < (-fEtaCut + fEtaShift) ){
1684       return kFALSE;
1685    }
1686    if(fEtaCutMin>-0.1){
1687       if( particle->Eta() < (fEtaCutMin + fEtaShift) && particle->Eta() > (-fEtaCutMin + fEtaShift) ){
1688          return kFALSE;
1689       }
1690       if( ePos->Eta() < (fEtaCutMin + fEtaShift) && ePos->Eta() > (-fEtaCutMin + fEtaShift) ){
1691          return kFALSE;
1692       }
1693       if( eNeg->Eta() < (fEtaCutMin + fEtaShift) && eNeg->Eta() > (-fEtaCutMin + fEtaShift) ){
1694          return kFALSE;
1695       }
1696    }
1697
1698    if( ePos->Pt()< fSinglePtCut ||  eNeg->Pt()< fSinglePtCut){
1699       return kFALSE;
1700    }
1701
1702    if(particle->Pt()<fPtCut){
1703       return kFALSE;
1704    }
1705
1706    return kTRUE;
1707 }
1708 ///________________________________________________________________________
1709 Bool_t AliConversionCuts::UpdateCutString() {
1710    ///Update the cut string (if it has been created yet)
1711
1712    if(fCutString && fCutString->GetString().Length() == kNCuts) {
1713       fCutString->SetString(GetCutNumber());
1714    } else {
1715       return kFALSE;
1716    }
1717    return kTRUE;
1718 }
1719 ///________________________________________________________________________
1720 void AliConversionCuts::LoadReweightingHistosMCFromFile() {
1721
1722   AliInfo("Entering loading of histograms for weighting");
1723   TFile *f = TFile::Open(fPathTrFReweighting.Data());
1724   if(!f){
1725      AliError(Form("file for weighting %s not found",fPathTrFReweighting.Data()));
1726      return;
1727   }
1728   if (fNameHistoReweightingPi0.CompareTo("") != 0 && fDoReweightHistoMCPi0 ){
1729      TH1D *hReweightMCHistPi0temp = (TH1D*)f->Get(fNameHistoReweightingPi0.Data());
1730      hReweightMCHistPi0 = new TH1D(*hReweightMCHistPi0temp);
1731      hReweightMCHistPi0->SetDirectory(0);
1732      if (hReweightMCHistPi0) AliInfo(Form("%s has been loaded from %s", fNameHistoReweightingPi0.Data(),fPathTrFReweighting.Data() ));
1733      else AliWarning(Form("%s not found in %s", fNameHistoReweightingPi0.Data() ,fPathTrFReweighting.Data()));
1734   }
1735   if (fNameFitDataPi0.CompareTo("") != 0 && fDoReweightHistoMCPi0 ){
1736      TF1 *fFitDataPi0temp = (TF1*)f->Get(fNameFitDataPi0.Data());
1737      fFitDataPi0 = new TF1(*fFitDataPi0temp);
1738      if (fFitDataPi0) AliInfo(Form("%s has been loaded from %s", fNameFitDataPi0.Data(),fPathTrFReweighting.Data() ));
1739      else AliWarning(Form("%s not found in %s",fPathTrFReweighting.Data(), fNameFitDataPi0.Data() ));
1740   }
1741
1742   if (fNameHistoReweightingEta.CompareTo("") != 0 && fDoReweightHistoMCEta){
1743      TH1D *hReweightMCHistEtatemp = (TH1D*)f->Get(fNameHistoReweightingEta.Data());
1744      hReweightMCHistEta = new TH1D(*hReweightMCHistEtatemp);
1745      hReweightMCHistEta->SetDirectory(0);
1746      if (hReweightMCHistEta) AliInfo(Form("%s has been loaded from %s", fNameHistoReweightingEta.Data(),fPathTrFReweighting.Data() ));
1747      else AliWarning(Form("%s not found in %s", fNameHistoReweightingEta.Data(),fPathTrFReweighting.Data() ));
1748   }
1749
1750   if (fNameFitDataEta.CompareTo("") != 0 && fDoReweightHistoMCEta){
1751      TF1 *fFitDataEtatemp = (TF1*)f->Get(fNameFitDataEta.Data());
1752      fFitDataEta = new TF1(*fFitDataEtatemp);
1753      if (fFitDataEta) AliInfo(Form("%s has been loaded from %s", fNameFitDataEta.Data(),fPathTrFReweighting.Data() ));
1754      else AliWarning(Form("%s not found in %s", fNameFitDataEta.Data(),fPathTrFReweighting.Data() ));
1755
1756   }
1757   if (fNameHistoReweightingK0s.CompareTo("") != 0 && fDoReweightHistoMCK0s){
1758      TH1D *hReweightMCHistK0stemp = (TH1D*)f->Get(fNameHistoReweightingK0s.Data());
1759      hReweightMCHistK0s = new TH1D(*hReweightMCHistK0stemp);
1760      hReweightMCHistK0s->SetDirectory(0);
1761      if (hReweightMCHistK0s) AliInfo(Form("%s has been loaded from %s", fNameHistoReweightingK0s.Data(),fPathTrFReweighting.Data() ));
1762      else AliWarning(Form("%s not found in %s", fNameHistoReweightingK0s.Data(),fPathTrFReweighting.Data() ));
1763   }
1764
1765   if (fNameFitDataK0s.CompareTo("") != 0 && fDoReweightHistoMCK0s){
1766      TF1 *fFitDataK0stemp = (TF1*)f->Get(fNameFitDataK0s.Data());
1767      fFitDataK0s = new TF1(*fFitDataK0stemp);
1768      if (fFitDataK0s) AliInfo(Form("%s has been loaded from %s", fNameFitDataK0s.Data(),fPathTrFReweighting.Data() ));
1769      else AliWarning(Form("%s not found in %s", fNameFitDataK0s.Data(),fPathTrFReweighting.Data() ));
1770   }
1771   f->Close();
1772   delete f;
1773   
1774 }
1775
1776
1777 ///________________________________________________________________________
1778 Bool_t AliConversionCuts::InitializeCutsFromCutString(const TString analysisCutSelection ) {
1779    // Initialize Cuts from a given Cut string
1780    if(fDoReweightHistoMCPi0 || fDoReweightHistoMCEta || fDoReweightHistoMCK0s) {
1781       AliInfo("Weighting was enabled");
1782       LoadReweightingHistosMCFromFile();
1783    }
1784
1785    AliInfo(Form("Set Photoncut Number: %s",analysisCutSelection.Data()));
1786    if(analysisCutSelection.Length()!=kNCuts) {
1787       AliError(Form("Cut selection has the wrong length! size is %d, number of cuts is %d", analysisCutSelection.Length(), kNCuts));
1788       return kFALSE;
1789    }
1790    if(!analysisCutSelection.IsDigit()){
1791       AliError("Cut selection contains characters");
1792       return kFALSE;
1793    }
1794
1795    const char *cutSelection = analysisCutSelection.Data();
1796 #define ASSIGNARRAY(i)  fCuts[i] = cutSelection[i] - '0'
1797    for(Int_t ii=0;ii<kNCuts;ii++){
1798       ASSIGNARRAY(ii);
1799    }
1800
1801    // Set Individual Cuts
1802    for(Int_t ii=0;ii<kNCuts;ii++){
1803       if(!SetCut(cutIds(ii),fCuts[ii]))return kFALSE;
1804    }
1805
1806    PrintCutsWithValues();
1807
1808    return kTRUE;
1809 }
1810 ///________________________________________________________________________
1811 Bool_t AliConversionCuts::SetCut(cutIds cutID, const Int_t value) {
1812    ///Set individual cut ID
1813
1814    switch (cutID) {
1815
1816    case kv0FinderType:
1817       if( SetV0Finder(value)) {
1818          fCuts[kv0FinderType] = value;
1819          UpdateCutString();
1820          return kTRUE;
1821       } else return kFALSE;
1822
1823    case kededxSigmaCut:
1824       if( SetTPCdEdxCutElectronLine(value)) {
1825          fCuts[kededxSigmaCut] = value;
1826          UpdateCutString();
1827          return kTRUE;
1828       } else return kFALSE;
1829
1830    case kpidedxSigmaCut:
1831       if( SetTPCdEdxCutPionLine(value)) {
1832          fCuts[kpidedxSigmaCut] = value;
1833          UpdateCutString();
1834          return kTRUE;
1835       } else return kFALSE;
1836
1837    case kpiMomdedxSigmaCut:
1838       if( SetMinMomPiondEdxCut(value)) {
1839          fCuts[kpiMomdedxSigmaCut] = value;
1840          UpdateCutString();
1841          return kTRUE;
1842       } else return kFALSE;
1843
1844    case kchi2GammaCut:
1845       if( SetChi2GammaCut(value)) {
1846          fCuts[kchi2GammaCut] = value;
1847          UpdateCutString();
1848          return kTRUE;
1849       } else return kFALSE;
1850
1851    case ksinglePtCut:
1852       if( SetSinglePtCut(value)) {
1853          fCuts[ksinglePtCut] = value;
1854          UpdateCutString();
1855          return kTRUE;
1856       } else return kFALSE;
1857
1858    case kclsTPCCut:
1859       if( SetTPCClusterCut(value)) {
1860          fCuts[kclsTPCCut] = value;
1861          UpdateCutString();
1862          return kTRUE;
1863       } else return kFALSE;
1864
1865    case ketaCut:
1866       if( SetEtaCut(value)) {
1867          fCuts[ketaCut] = value;
1868          UpdateCutString();
1869          return kTRUE;
1870       } else return kFALSE;
1871
1872    case kLowPRejectionSigmaCut:
1873       if( SetLowPRejectionCuts(value)) {
1874          fCuts[kLowPRejectionSigmaCut] = value;
1875          UpdateCutString();
1876          return kTRUE;
1877       } else return kFALSE;
1878
1879    case kQtMaxCut:
1880       if( SetQtMaxCut(value)) {
1881          fCuts[kQtMaxCut] = value;
1882          UpdateCutString();
1883          return kTRUE;
1884       } else return kFALSE;
1885
1886    case kpiMaxMomdedxSigmaCut:
1887       if( SetMaxMomPiondEdxCut(value)) {
1888          fCuts[kpiMaxMomdedxSigmaCut] = value;
1889          UpdateCutString();
1890          return kTRUE;
1891       } else return kFALSE;
1892
1893    case kRCut:
1894       if( SetRCut(value)) {
1895          fCuts[kRCut] = value;
1896          UpdateCutString();
1897          return kTRUE;
1898       } else return kFALSE;
1899
1900    case kremovePileUp:
1901       if( SetRemovePileUp(value)) {
1902          fCuts[kremovePileUp] = value;
1903          UpdateCutString();
1904          return kTRUE;
1905       } else return kFALSE;
1906
1907    case kselectV0AND:
1908       if( SetSelectSpecialTrigger(value)) {
1909          fCuts[kselectV0AND] = value;
1910          UpdateCutString();
1911          return kTRUE;
1912       } else return kFALSE;
1913
1914    case kmultiplicityMethod:
1915       if( SetMultiplicityMethod(value)) {
1916          fCuts[kmultiplicityMethod] = value;
1917          UpdateCutString();
1918          return kTRUE;
1919       } else return kFALSE;
1920
1921    case kisHeavyIon:
1922       if( SetIsHeavyIon(value)) {
1923          fCuts[kisHeavyIon] = value;
1924          UpdateCutString();
1925          return kTRUE;
1926       } else return kFALSE;
1927
1928    case kCentralityMin:
1929       if( SetCentralityMin(value)) {
1930          fCuts[kCentralityMin] = value;
1931          UpdateCutString();
1932          return kTRUE;
1933       } else return kFALSE;
1934
1935    case kCentralityMax:
1936       if( SetCentralityMax(value)) {
1937          fCuts[kCentralityMax] = value;
1938          UpdateCutString();
1939          return kTRUE;
1940       } else return kFALSE;
1941
1942    case kTOFelectronPID:
1943       if( SetTOFElectronPIDCut(value)) {
1944          fCuts[kTOFelectronPID] = value;
1945          UpdateCutString();
1946          return kTRUE;
1947       } else return kFALSE;
1948
1949    case kdoPhotonAsymmetryCut:
1950       if( SetPhotonAsymmetryCut(value)) {
1951          fCuts[kdoPhotonAsymmetryCut] = value;
1952          UpdateCutString();
1953          return kTRUE;
1954       } else return kFALSE;
1955
1956    case kPsiPair:
1957       if( SetPsiPairCut(value)) {
1958          fCuts[kPsiPair] = value;
1959          UpdateCutString();
1960          return kTRUE;
1961       } else return kFALSE;
1962
1963    case kCosPAngle:
1964       if( SetCosPAngleCut(value)) {
1965          fCuts[kCosPAngle] = value;
1966          UpdateCutString();
1967          return kTRUE;
1968       } else return kFALSE;
1969
1970
1971    case kElecShare:
1972       if( SetSharedElectronCut(value)) {
1973          fCuts[kElecShare] = value;
1974          UpdateCutString();
1975          return kTRUE;
1976       } else return kFALSE;
1977
1978    case kToCloseV0s:
1979       if( SetToCloseV0sCut(value)) {
1980          fCuts[kToCloseV0s] = value;
1981          UpdateCutString();
1982          return kTRUE;
1983       } else return kFALSE;
1984
1985    case kExtraSignals:
1986       if( SetRejectExtraSignalsCut(value)) {
1987          fCuts[kExtraSignals] = value;
1988          UpdateCutString();
1989          return kTRUE;
1990       } else return kFALSE;
1991
1992    case kDcaRPrimVtx:
1993       if( SetDCARPhotonPrimVtxCut(value)) {
1994          fCuts[kDcaRPrimVtx] = value;
1995          UpdateCutString();
1996          return kTRUE;
1997       } else return kFALSE;
1998
1999    case kDcaZPrimVtx:
2000       if( SetDCAZPhotonPrimVtxCut(value)) {
2001          fCuts[kDcaZPrimVtx] = value;
2002          UpdateCutString();
2003          return kTRUE;
2004       } else return kFALSE;
2005
2006    case kInPlaneOutOfPlane:
2007    if( SetInPlaneOutOfPlane(value)) {
2008       fCuts[kInPlaneOutOfPlane] = value;
2009       UpdateCutString();
2010       return kTRUE;
2011    } else return kFALSE;
2012
2013
2014       
2015
2016    case kNCuts:
2017       AliError("Cut id out of range");
2018       return kFALSE;
2019    }
2020
2021    AliError("Cut id %d not recognized");
2022    return kFALSE;
2023
2024
2025 }
2026 ///________________________________________________________________________
2027 void AliConversionCuts::PrintCuts() {
2028    // Print out current Cut Selection
2029    for(Int_t ic = 0; ic < kNCuts; ic++) {
2030       printf("%-30s : %d \n", fgkCutNames[ic], fCuts[ic]);
2031    }
2032 }
2033
2034 void AliConversionCuts::PrintCutsWithValues() {
2035    // Print out current Cut Selection with value
2036    if (fIsHeavyIon == 0) {
2037       printf("Running in pp mode \n");
2038       if (fSpecialTrigger == 0){
2039         printf("\t only events triggered by V0OR will be analysed \n");
2040       } else if (fSpecialTrigger == 1){
2041         printf("\t only events triggered by V0AND will be analysed \n");
2042       } else if (fSpecialTrigger == 2){
2043          printf("\t only events where SDD was present will be analysed \n");
2044       } else if (fSpecialTrigger == 3){
2045          printf("\t only events where SDD was present will be analysed and triggered by VOAND\n");
2046       } else if (fSpecialTrigger > 3){   
2047          printf("\t only events triggered by %s \n", fSpecialTriggerName.Data());
2048       }
2049    } else if (fIsHeavyIon == 1){ 
2050       printf("Running in PbPb mode \n");
2051       if (fDetectorCentrality == 0){
2052          printf("\t centrality selection based on V0M \n");
2053       } else if (fDetectorCentrality == 1){
2054          printf("\t centrality selection based on Cl1 \n");
2055       }   
2056       if (fModCentralityClass == 0){
2057         printf("\t %d - %d \n", fCentralityMin*10, fCentralityMax*10);
2058       } else if ( fModCentralityClass == 1){ 
2059         printf("\t %d - %d \n", fCentralityMin*5, fCentralityMax*5);
2060       } else if ( fModCentralityClass == 2){ 
2061         printf("\t %d - %d \n", fCentralityMin*5+45, fCentralityMax*5+45);
2062       } else if (fModCentralityClass == 3){
2063         printf("\t %d - %d, with Track mult in MC as data \n", fCentralityMin*10, fCentralityMax*10);
2064       } else if ( fModCentralityClass == 4){ 
2065         printf("\t %d - %d, with Track mult in MC as data \n", fCentralityMin*5, fCentralityMax*5);
2066       } else if ( fModCentralityClass == 5){ 
2067         printf("\t %d - %d, with Track mult in MC as data \n", fCentralityMin*5+45, fCentralityMax*5+45);
2068       }
2069       if (fSpecialTrigger == 0){
2070         printf("\t only events triggered by kMB, kCentral, kSemiCentral will be analysed \n");
2071       } else if (fSpecialTrigger > 4){   
2072          printf("\t only events triggered by %s \n", fSpecialTriggerName.Data());
2073       }
2074    } else if (fIsHeavyIon == 2){
2075       printf("Running in pPb mode \n");
2076       if (fDetectorCentrality == 0){
2077          printf("\t centrality selection based on V0A \n");
2078       } else if (fDetectorCentrality == 1){
2079          printf("\t centrality selection based on Cl1 \n");
2080       }   
2081       if (fModCentralityClass == 0){
2082         printf("\t %d - %d \n", fCentralityMin*10, fCentralityMax*10);
2083       }
2084       if (fSpecialTrigger == 0){
2085         printf("\t only events triggered by kINT7 will be analysed \n");
2086       } else if (fSpecialTrigger > 4){   
2087          printf("\t only events triggered by %s \n", fSpecialTriggerName.Data());
2088       }
2089    }
2090    printf("Electron cuts: \n");
2091    if (fEtaCutMin > -0.1) printf("\t %3.2f < eta_{e} < %3.2f\n", fEtaCutMin, fEtaCut );
2092      else printf("\t eta_{e} < %3.2f\n", fEtaCut );
2093    printf("\t p_{T,e} > %3.2f\n", fSinglePtCut );
2094    printf("\t %3.2f < n sigma e < %3.2f\n", fPIDnSigmaBelowElectronLine, fPIDnSigmaAboveElectronLine );
2095    
2096    printf("Photon cuts: \n");
2097    printf("\t %3.2f < R_{conv} < %3.2f\n", fMinR, fMaxR );
2098    printf("\t Z_{conv} < %3.2f\n", fMaxZ );
2099    if (fEtaCutMin > -0.1) printf("\t %3.2f < eta_{conv} < %3.2f\n", fEtaCutMin, fEtaCut );
2100      else printf("\t eta_{conv} < %3.2f\n", fEtaCut );
2101    printf("\t p_{T,gamma} > %3.2f\n", fPtCut );  
2102    if (fDo2DQt){
2103           printf("\t 2 dimensional q_{T} cut applied with maximum of %3.2f \n", fQtMax );
2104    } else {
2105           printf("\t 1 dimensional q_{T} cut applied with maximum of %3.2f \n", fQtMax );
2106    }
2107    if (fDo2DPsiPairChi2){
2108           printf("\t 2 dimensional triangle chi^{2} and psi_{pair} cut applied with maximum of chi^{2} = %3.2f and |psi_{pair}| = %3.2f \n", fChi2CutConversion, fPsiPairCut ); 
2109    } else {
2110       printf("\t chi^{2} max cut chi^{2} < %3.2f \n", fChi2CutConversion ); 
2111           printf("\t psi_{pair} max cut |psi_{pair}| < %3.2f \n", fPsiPairCut ); 
2112    }       
2113    printf("\t cos(Theta_{point}) > %3.2f \n", fCosPAngleCut );
2114    printf("\t dca_{R} < %3.2f \n", fDCARPrimVtxCut );
2115    printf("\t dca_{Z} < %3.2f \n", fDCAZPrimVtxCut );
2116    if (fDoPhotonQualitySelectionCut) printf("\t selection based on photon quality with quality %d \n", fPhotonQualityCut );
2117 }
2118
2119 ///________________________________________________________________________
2120 Bool_t AliConversionCuts::SetIsHeavyIon(Int_t isHeavyIon)
2121 {   // Set Cut
2122    switch(isHeavyIon){
2123    case 0:
2124       fIsHeavyIon=0;
2125       break;
2126    case 1:
2127       fIsHeavyIon=1;
2128       fDetectorCentrality=0;
2129       break;
2130    case 2:
2131       fIsHeavyIon=1;
2132       fDetectorCentrality=1;
2133       break;
2134    case 3: //allows to select centrality 0-45% in steps of 5% for V0 Multiplicity
2135       fIsHeavyIon=1;
2136       fDetectorCentrality=0;
2137       fModCentralityClass=1;
2138       break;
2139    case 4: //allows to select centrality 45-90% in steps of 5% for V0 Multiplicity
2140       fIsHeavyIon=1;
2141       fDetectorCentrality=0;
2142       fModCentralityClass=2;
2143       break;
2144    case 5: //strict cut on v0 tracks for MC
2145       fIsHeavyIon=1;
2146       fDetectorCentrality=0;
2147       fModCentralityClass=3;
2148       break;
2149    case 6: //allows to select centrality 0-45% in steps of 5% for track mult
2150       //strict cut on v0 tracks for MC
2151       fIsHeavyIon=1;
2152       fDetectorCentrality=0;
2153       fModCentralityClass=4;
2154       break;
2155    case 7: //allows to select centrality 45-90% in steps of 5% for V0 Multiplicity
2156       //strict cut on v0 tracks for MC
2157       fIsHeavyIon=1;
2158       fDetectorCentrality=0;
2159       fModCentralityClass=5;
2160       break;
2161    case 8:
2162       fIsHeavyIon=2;
2163       fDetectorCentrality=0;
2164       break;
2165    case 9:
2166       fIsHeavyIon=2;
2167       fDetectorCentrality=1;
2168       break;
2169    default:
2170       AliError(Form("SetHeavyIon not defined %d",isHeavyIon));
2171       return kFALSE;
2172    }
2173    return kTRUE;
2174 }
2175
2176 //___________________________________________________________________
2177 Bool_t AliConversionCuts::SetCentralityMin(Int_t minCentrality)
2178 {
2179    // Set Cut
2180    if(minCentrality<0||minCentrality>9){
2181       AliError(Form("minCentrality not defined %d",minCentrality));
2182       return kFALSE;
2183    }
2184
2185    fCentralityMin=minCentrality;
2186    return kTRUE;
2187 }
2188 //___________________________________________________________________
2189 Bool_t AliConversionCuts::SetCentralityMax(Int_t maxCentrality)
2190 {
2191    // Set Cut
2192    if(maxCentrality<0||maxCentrality>9){
2193       AliError(Form("maxCentrality not defined %d",maxCentrality));
2194       return kFALSE;
2195    }
2196    fCentralityMax=maxCentrality;
2197    return kTRUE;
2198 }
2199 ///________________________________________________________________________
2200 Int_t AliConversionCuts::SetSelectSpecialTrigger(Int_t selectSpecialTrigger)
2201 {// Set Cut
2202
2203    switch(selectSpecialTrigger){
2204    case 0:
2205       fSpecialTrigger=0; // dont care
2206       break;
2207    case 1:
2208       fSpecialTrigger=1; // V0AND
2209       break;
2210    case 2:
2211       fSpecialTrigger=2; // with SDD requested
2212       break;
2213    case 3:
2214       fSpecialTrigger=3; // V0AND plus with SDD requested
2215       break;
2216    // allows to run MB & 6 other different trigger classes in parallel with the same photon cut
2217    case 4:
2218       fSpecialTrigger=4; // different trigger class as MB
2219       fTriggerSelectedManually = kTRUE;
2220       break;
2221    case 5:
2222       fSpecialTrigger=4; // different trigger class as MB
2223       fTriggerSelectedManually = kTRUE;
2224       break;
2225    case 6:
2226       fSpecialTrigger=4; // different trigger class as MB
2227       fTriggerSelectedManually = kTRUE;
2228       break;
2229    case 7:
2230       fSpecialTrigger=4; // different trigger class as MB
2231       fTriggerSelectedManually = kTRUE;
2232       break;
2233     case 8:
2234       fSpecialTrigger=4; // different trigger class as MB
2235       fTriggerSelectedManually = kTRUE;
2236       break;
2237     case 9:
2238       fSpecialTrigger=4; // different trigger class as MB
2239       fTriggerSelectedManually = kTRUE;
2240       break;
2241    default:
2242       AliError("Warning: Special Trigger Not known");
2243       return kFALSE;
2244    }
2245    return kTRUE;
2246 }
2247 ///________________________________________________________________________
2248 Bool_t AliConversionCuts::SetMultiplicityMethod(Int_t multiplicityMethod)
2249 {
2250    // Set Cut
2251    fMultiplicityMethod=multiplicityMethod;
2252
2253    // 0 Photon Multiplicity
2254    // 1 TPC Track multiplicity
2255    // 2 V0 Mult
2256    // 3 SPD Mult
2257
2258    return kTRUE;
2259 }
2260 ///________________________________________________________________________
2261 Bool_t AliConversionCuts::SetRemovePileUp(Int_t removePileUp)
2262 {// Set Cut
2263    switch(removePileUp){
2264    case 0:
2265       fRemovePileUp=kFALSE;
2266       break;
2267    case 1:
2268       fRemovePileUp=kTRUE;
2269       break;
2270    default:
2271       AliError("RemovePileUpCut not defined");
2272       return kFALSE;
2273    }
2274    return kTRUE;
2275 }
2276 ///________________________________________________________________________
2277 Bool_t AliConversionCuts::SetRejectExtraSignalsCut(Int_t extraSignal) {
2278
2279    switch(extraSignal){
2280    case 0:
2281       fRejectExtraSignals = 0;
2282       break; // No Rejection
2283    case 1:
2284       fRejectExtraSignals = 1;
2285       break; // MinBias Header
2286    case 2:
2287       fRejectExtraSignals = 2;
2288       break; // User String Array
2289    case 3:
2290       fRejectExtraSignals = 3;
2291       break; // Rejection for Gamma Correction only
2292    default:
2293       AliError(Form("Extra Signal Rejection not defined %d",extraSignal));
2294       return kFALSE;
2295    }
2296    return kTRUE;
2297 }
2298 ///________________________________________________________________________
2299 Bool_t AliConversionCuts::SetV0Finder(Int_t v0FinderType)
2300 {   // Set Cut
2301    switch (v0FinderType){
2302    case 0:  // on fly V0 finder
2303       cout << "have chosen onfly V0" << endl;
2304       fUseOnFlyV0Finder=kTRUE;
2305       break;
2306    case 1:  // offline V0 finder
2307       cout << "have chosen offline V0" << endl;
2308       fUseOnFlyV0Finder=kFALSE;
2309       break;
2310    default:
2311       AliError(Form(" v0FinderType not defined %d",v0FinderType));
2312       return kFALSE;
2313    }
2314    return kTRUE;
2315 }
2316 ///________________________________________________________________________
2317 Bool_t AliConversionCuts::SetEtaCut(Int_t etaCut)
2318 {   // Set Cut
2319
2320    //Set Standard LineCutZValues
2321    fLineCutZValueMin = -2;
2322    fLineCutZValue = 7.;
2323
2324    switch(etaCut){
2325    case 0: // 0.9
2326       fEtaCut     = 0.9;
2327       fLineCutZRSlope = tan(2*atan(exp(-fEtaCut)));
2328       fEtaCutMin     = -0.1;
2329       fLineCutZRSlopeMin = 0.;
2330       break;
2331    case 1:  // 0.6  // changed from 1.2 to 0.6 on 2013.06.10
2332       fEtaCut     = 0.6;
2333       fLineCutZRSlope = tan(2*atan(exp(-fEtaCut)));
2334       fEtaCutMin     = -0.1;
2335       fLineCutZRSlopeMin = 0.;
2336       break;
2337    case 2:  // 1.4
2338       fEtaCut     = 1.4;
2339       fLineCutZRSlope = tan(2*atan(exp(-fEtaCut)));
2340       fEtaCutMin     = -0.1;
2341       fLineCutZRSlopeMin = 0.;
2342       break;
2343    case 3: // 0.65
2344       fEtaCut     = 0.65;
2345       fLineCutZRSlope = tan(2*atan(exp(-fEtaCut)));
2346       fEtaCutMin     = -0.1;
2347       fLineCutZRSlopeMin = 0.;
2348       break;
2349    case 4: // 0.75
2350       fEtaCut     = 0.75;
2351       fLineCutZRSlope = tan(2*atan(exp(-fEtaCut)));
2352       fEtaCutMin     = -0.1;
2353       fLineCutZRSlopeMin = 0.;
2354       break;
2355    case 5: // 0.5
2356       fEtaCut     = 0.5;
2357       fLineCutZRSlope = tan(2*atan(exp(-fEtaCut)));
2358       fEtaCutMin     = -0.1;
2359       fLineCutZRSlopeMin = 0.;
2360       break;
2361    case 6: // 5.
2362       fEtaCut     = 5.;
2363       fLineCutZRSlope = tan(2*atan(exp(-fEtaCut)));
2364       fEtaCutMin     = -0.1;
2365       fLineCutZRSlopeMin = 0.;
2366       break;
2367    case 7:
2368       if (fIsHeavyIon==1){
2369          fEtaCut     = 0.7;
2370          fLineCutZRSlope = tan(2*atan(exp(-fEtaCut)));
2371          fEtaCutMin     = -0.1;
2372          fLineCutZRSlopeMin = 0.;
2373          break;
2374       } else {   
2375          fEtaCut     = 0.3;
2376          fLineCutZRSlope = tan(2*atan(exp(-fEtaCut)));
2377          fEtaCutMin     = -0.1;
2378          fLineCutZRSlopeMin = 0.;
2379          break;
2380       }
2381    // case 8: // 0.1 - 0.8
2382    //    fEtaCut     = 0.9;
2383    //    fLineCutZRSlope = tan(2*atan(exp(-fEtaCut)));
2384    //    fEtaCutMin     = 0.1;
2385    //    fLineCutZRSlopeMin = tan(2*atan(exp(-fEtaCutMin)));
2386    //    break;
2387    case 8: // 0.4
2388       fEtaCut     = 0.4;
2389       fLineCutZRSlope = tan(2*atan(exp(-fEtaCut)));
2390       fEtaCutMin     = -0.1;
2391       fLineCutZRSlopeMin = 0.;
2392       break;
2393    case 9: // 10
2394       fEtaCut     = 10;
2395       fLineCutZRSlope = tan(2*atan(exp(-fEtaCut)));
2396       fEtaCutMin     = -0.1;
2397       fLineCutZRSlopeMin = 0.;
2398       break;
2399    default:
2400       AliError(Form(" EtaCut not defined %d",etaCut));
2401       return kFALSE;
2402    }
2403    return kTRUE;
2404 }
2405 ///________________________________________________________________________
2406 Bool_t AliConversionCuts::SetRCut(Int_t RCut){
2407    // Set Cut
2408    switch(RCut){
2409    case 0:
2410       fMinR=0;
2411       fMaxR = 180.;
2412       break;
2413    case 1:
2414       fMinR=2.8;
2415       fMaxR = 180.;
2416       break;
2417    case 2:
2418       fMinR=5.;
2419       fMaxR = 180.;
2420       break;
2421    case 3:
2422       fMaxR = 70.;
2423       fMinR = 10.;
2424       break;
2425    case 4:
2426       fMaxR = 70.;
2427       fMinR = 5.;
2428       break;
2429       // High purity cuts for PbPb (remove first layers of material)
2430    case 5:
2431       fMaxR = 180.;
2432       fMinR = 10.;
2433       break;
2434    case 6:
2435       fMaxR = 180.;
2436       fMinR = 20.;
2437       break;
2438    case 7:
2439       fMaxR = 180.;
2440       fMinR = 26.;
2441       break;
2442    case 8:
2443       fMaxR = 180.;
2444       fMinR = 12.5;
2445       break;
2446    case 9:
2447       fMaxR = 180.;
2448       fMinR = 7.5;
2449       break;
2450
2451    default:
2452       AliError("RCut not defined");
2453       return kFALSE;
2454    }
2455    return kTRUE;
2456 }
2457 ///________________________________________________________________________
2458 Bool_t AliConversionCuts::SetSinglePtCut(Int_t singlePtCut)
2459 {   // Set Cut
2460    switch(singlePtCut){
2461    case 0: // 0.050 GeV
2462       fSinglePtCut = 0.050;
2463       break;
2464    case 1:  // 0.100 GeV
2465       fSinglePtCut = 0.100;
2466       break;
2467    case 2:  // 0.150 GeV
2468       fSinglePtCut = 0.150;
2469       break;
2470    case 3:  // 0.200 GeV
2471       fSinglePtCut = 0.200;
2472       break;
2473    case 4:  // 0.075 GeV
2474       fSinglePtCut = 0.075;
2475       break;
2476    case 5:  // 0.125 GeV
2477       fSinglePtCut = 0.125;
2478       break;
2479    case 6:  // 0.04 GeV
2480       fSinglePtCut = 0.040;
2481       break;
2482    case 7:  // 0.0 GeV
2483       fSinglePtCut = 0.0;
2484       break;
2485    default:
2486       AliError(Form("singlePtCut not defined %d",singlePtCut));
2487       return kFALSE;
2488    }
2489    return kTRUE;
2490 }
2491 ///________________________________________________________________________
2492 Bool_t AliConversionCuts::SetTPCClusterCut(Int_t clsTPCCut)
2493 {   // Set Cut
2494    switch(clsTPCCut){
2495    case 0: // 0
2496       fMinClsTPC= 0.;
2497       break;
2498    case 1:  // 60
2499       fMinClsTPC= 60.;
2500       break;
2501    case 2:  // 80
2502       fMinClsTPC= 80.;
2503       break;
2504    case 3:  // 100
2505       fMinClsTPC= 100.;
2506       break;
2507    case 4:  // 95% of findable clusters
2508       fMinClsTPCToF= 0.95;
2509       fUseCorrectedTPCClsInfo=1;
2510       break;
2511    case 5:  // 0% of findable clusters
2512       fMinClsTPCToF= 0.0;
2513       fUseCorrectedTPCClsInfo=1;
2514       break;
2515    case 6:  // 70% of findable clusters
2516       fMinClsTPCToF= 0.7;
2517       fUseCorrectedTPCClsInfo=1;
2518       break;
2519    case 7:  // 0% of findable clusters
2520       fMinClsTPCToF= 0.35;
2521       fUseCorrectedTPCClsInfo=0;
2522       break;
2523    case 8:
2524       fMinClsTPCToF= 0.35;
2525       fUseCorrectedTPCClsInfo=1;
2526       break;
2527    case 9:
2528       fMinClsTPCToF= 0.6;
2529       fUseCorrectedTPCClsInfo=1;
2530       break;
2531    default:
2532       AliError(Form("Warning: clsTPCCut not defined %d",clsTPCCut));
2533       return kFALSE;
2534    }
2535    return kTRUE;
2536 }
2537 ///________________________________________________________________________
2538 Bool_t AliConversionCuts::SetTPCdEdxCutElectronLine(Int_t ededxSigmaCut)
2539 {   // Set Cut
2540    switch(ededxSigmaCut){
2541    case 0: // -10,10
2542       fPIDnSigmaBelowElectronLine=-10;
2543       fPIDnSigmaAboveElectronLine=10;
2544       break;
2545    case 1: // -5,5
2546       fPIDnSigmaBelowElectronLine=-5;
2547       fPIDnSigmaAboveElectronLine=5;
2548       break;
2549    case 2: // -3,5
2550       fPIDnSigmaBelowElectronLine=-3;
2551       fPIDnSigmaAboveElectronLine=5;
2552       break;
2553    case 3: // -4,5
2554       fPIDnSigmaBelowElectronLine=-4;
2555       fPIDnSigmaAboveElectronLine=5;
2556       break;
2557    case 4: // -6,7
2558       fPIDnSigmaBelowElectronLine=-6;
2559       fPIDnSigmaAboveElectronLine=7;
2560       break;
2561    case 5: // -4,4
2562       fPIDnSigmaBelowElectronLine=-4;
2563       fPIDnSigmaAboveElectronLine=4;
2564       break;
2565    case 6: // -2.5,4
2566       fPIDnSigmaBelowElectronLine=-2.5;
2567       fPIDnSigmaAboveElectronLine=4;
2568       break;
2569    case 7: // -2,3.5
2570       fPIDnSigmaBelowElectronLine=-2;
2571       fPIDnSigmaAboveElectronLine=3.5;
2572       break;
2573    default:
2574       AliError("TPCdEdxCutElectronLine not defined");
2575       return kFALSE;
2576
2577    }
2578    return kTRUE;
2579 }
2580 ///________________________________________________________________________
2581 Bool_t AliConversionCuts::SetTPCdEdxCutPionLine(Int_t pidedxSigmaCut)
2582 {   // Set Cut
2583
2584    switch(pidedxSigmaCut){
2585    case 0:  // -10
2586       fPIDnSigmaAbovePionLine=-10;
2587       fPIDnSigmaAbovePionLineHighPt=-10;
2588       break;
2589    case 1:   // 0
2590       fPIDnSigmaAbovePionLine=0;
2591       fPIDnSigmaAbovePionLineHighPt=-10;
2592       break;
2593    case 2:  // 1
2594       fPIDnSigmaAbovePionLine=1;
2595       fPIDnSigmaAbovePionLineHighPt=-10;
2596       break;
2597    case 3:  // 1
2598       fPIDnSigmaAbovePionLine=2.5;
2599       fPIDnSigmaAbovePionLineHighPt=-10;
2600       break;
2601    case 4:  // 1
2602       fPIDnSigmaAbovePionLine=0.5;
2603       fPIDnSigmaAbovePionLineHighPt=-10;
2604       break;
2605    case 5:  // 1
2606       fPIDnSigmaAbovePionLine=2.;
2607       fPIDnSigmaAbovePionLineHighPt=-10;
2608       break;
2609    case 6:  // 1
2610       fPIDnSigmaAbovePionLine=2.;
2611       fPIDnSigmaAbovePionLineHighPt=0.5;
2612       break;
2613    case 7:  // 1
2614       fPIDnSigmaAbovePionLine=3.5;
2615       fPIDnSigmaAbovePionLineHighPt=-10;
2616       break;
2617    case 8:  // 1
2618       fPIDnSigmaAbovePionLine=2.;
2619       fPIDnSigmaAbovePionLineHighPt=1.;
2620       break;
2621    case 9:
2622       fPIDnSigmaAbovePionLine=3.0; // We need a bit less tight cut on dE/dx
2623       fPIDnSigmaAbovePionLineHighPt=-10;
2624       break;
2625    default:
2626       AliError(Form("Warning: pidedxSigmaCut not defined %d",pidedxSigmaCut));
2627       return kFALSE;
2628    }
2629    return kTRUE;
2630 }
2631 ///________________________________________________________________________
2632 Bool_t AliConversionCuts::SetMinMomPiondEdxCut(Int_t piMomdedxSigmaCut)
2633 {   // Set Cut
2634    switch(piMomdedxSigmaCut){
2635    case 0:  // 0.5 GeV
2636       fPIDMinPnSigmaAbovePionLine=0.5;
2637       break;
2638    case 1:  // 1. GeV
2639       fPIDMinPnSigmaAbovePionLine=1.;
2640       break;
2641    case 2:  // 1.5 GeV
2642       fPIDMinPnSigmaAbovePionLine=1.5;
2643       break;
2644    case 3:  // 20.0 GeV
2645       fPIDMinPnSigmaAbovePionLine=20.;
2646       break;
2647    case 4:  // 50.0 GeV
2648       fPIDMinPnSigmaAbovePionLine=50.;
2649       break;
2650    case 5:  // 0.3 GeV
2651       fPIDMinPnSigmaAbovePionLine=0.3;
2652       break;
2653    case 6:  // 0.25 GeV
2654       fPIDMinPnSigmaAbovePionLine=0.25;
2655       break;
2656    case 7:  // 0.4 GeV
2657       fPIDMinPnSigmaAbovePionLine=0.4;
2658       break;
2659    case 8:  // 0.2 GeV
2660       fPIDMinPnSigmaAbovePionLine=0.2;
2661       break;
2662    default:
2663       AliError(Form("piMomdedxSigmaCut not defined %d",piMomdedxSigmaCut));
2664       return kFALSE;
2665    }
2666    return kTRUE;
2667 }
2668 ///________________________________________________________________________
2669 Bool_t AliConversionCuts::SetMaxMomPiondEdxCut(Int_t piMaxMomdedxSigmaCut)
2670 {   // Set Cut
2671    switch(piMaxMomdedxSigmaCut){
2672    case 0:  // 100. GeV
2673       fPIDMaxPnSigmaAbovePionLine=100.;
2674       break;
2675    case 1:  // 5. GeV
2676       fPIDMaxPnSigmaAbovePionLine=5.;
2677       break;
2678    case 2:  // 4. GeV
2679       fPIDMaxPnSigmaAbovePionLine=4.;
2680       break;
2681    case 3:  // 3.5 GeV
2682       fPIDMaxPnSigmaAbovePionLine=3.5;
2683       break;
2684    case 4:  // 3. GeV
2685       fPIDMaxPnSigmaAbovePionLine=3.;
2686       break;
2687    case 5:  // 7. GeV
2688       fPIDMaxPnSigmaAbovePionLine=7.;
2689       break;
2690    default:
2691       AliError(Form("piMaxMomdedxSigmaCut not defined %d",piMaxMomdedxSigmaCut));
2692       return kFALSE;
2693    }
2694    return kTRUE;
2695 }
2696 ///________________________________________________________________________
2697 Bool_t AliConversionCuts::SetLowPRejectionCuts(Int_t LowPRejectionSigmaCut)
2698 {   // Set Cut
2699    switch(LowPRejectionSigmaCut){
2700    case 0:  //
2701       fPIDnSigmaAtLowPAroundKaonLine=0;
2702       fPIDnSigmaAtLowPAroundProtonLine=0;
2703       fPIDnSigmaAtLowPAroundPionLine=0;
2704       fDoKaonRejectionLowP = kFALSE;
2705       fDoProtonRejectionLowP = kFALSE;
2706       fDoPionRejectionLowP = kFALSE;
2707       fPIDMinPPionRejectionLowP = fPIDMinPnSigmaAbovePionLine;
2708       break;
2709    case 1:  //
2710       fPIDnSigmaAtLowPAroundKaonLine=0.5;
2711       fPIDnSigmaAtLowPAroundProtonLine=0.5;
2712       fPIDnSigmaAtLowPAroundPionLine=0.5;
2713       fDoKaonRejectionLowP = kTRUE;
2714       fDoProtonRejectionLowP = kTRUE;
2715       fDoPionRejectionLowP = kTRUE;
2716       fPIDMinPPionRejectionLowP = fPIDMinPnSigmaAbovePionLine;
2717       break;
2718    case 2:  //
2719       fPIDnSigmaAtLowPAroundKaonLine=1;
2720       fPIDnSigmaAtLowPAroundProtonLine=1;
2721       fPIDnSigmaAtLowPAroundPionLine=1;
2722       fDoKaonRejectionLowP = kTRUE;
2723       fDoProtonRejectionLowP = kTRUE;
2724       fDoPionRejectionLowP = kTRUE;
2725       fPIDMinPPionRejectionLowP = fPIDMinPnSigmaAbovePionLine;
2726       break;
2727    case 3:  //
2728       fPIDnSigmaAtLowPAroundKaonLine=2.;
2729       fPIDnSigmaAtLowPAroundProtonLine=2.;
2730       fPIDnSigmaAtLowPAroundPionLine=2.;
2731       fDoKaonRejectionLowP = kTRUE;
2732       fDoProtonRejectionLowP = kTRUE;
2733       fDoPionRejectionLowP = kTRUE;
2734       fPIDMinPPionRejectionLowP = fPIDMinPnSigmaAbovePionLine;
2735       break;
2736    case 4:  //
2737       fPIDnSigmaAtLowPAroundKaonLine=0.;
2738       fPIDnSigmaAtLowPAroundProtonLine=0.;
2739       fPIDnSigmaAtLowPAroundPionLine=1;
2740       fDoKaonRejectionLowP = kFALSE;
2741       fDoProtonRejectionLowP = kFALSE;
2742       fDoPionRejectionLowP = kTRUE;
2743       fPIDMinPPionRejectionLowP = fPIDMinPnSigmaAbovePionLine;
2744       break;
2745    case 5:  //
2746       fPIDnSigmaAtLowPAroundKaonLine=0.;
2747       fPIDnSigmaAtLowPAroundProtonLine=0.;
2748       fPIDnSigmaAtLowPAroundPionLine=1.5;
2749       fDoKaonRejectionLowP = kFALSE;
2750       fDoProtonRejectionLowP = kFALSE;
2751       fDoPionRejectionLowP = kTRUE;
2752       fPIDMinPPionRejectionLowP = fPIDMinPnSigmaAbovePionLine;
2753       break;
2754    case 6:  //
2755       fPIDnSigmaAtLowPAroundKaonLine=0.;
2756       fPIDnSigmaAtLowPAroundProtonLine=0.;
2757       fPIDnSigmaAtLowPAroundPionLine=2.;
2758       fDoKaonRejectionLowP = kFALSE;
2759       fDoProtonRejectionLowP = kFALSE;
2760       fDoPionRejectionLowP = kTRUE;
2761       fPIDMinPPionRejectionLowP = fPIDMinPnSigmaAbovePionLine;
2762       break;
2763    case 7:  //
2764       fPIDnSigmaAtLowPAroundKaonLine=0.;
2765       fPIDnSigmaAtLowPAroundProtonLine=0.;
2766       fPIDnSigmaAtLowPAroundPionLine=0.5;
2767       fDoKaonRejectionLowP = kFALSE;
2768       fDoProtonRejectionLowP = kFALSE;
2769       fDoPionRejectionLowP = kTRUE;
2770       fPIDMinPPionRejectionLowP = fPIDMinPnSigmaAbovePionLine;
2771       break;
2772    default:
2773       AliError(Form("LowPRejectionSigmaCut not defined %d",LowPRejectionSigmaCut));
2774       return kFALSE;
2775    }
2776    return kTRUE;
2777 }
2778 ///________________________________________________________________________
2779 Bool_t AliConversionCuts::SetTOFElectronPIDCut(Int_t TOFelectronPID){
2780    // Set Cut
2781    switch(TOFelectronPID){
2782    case 0: // no cut
2783       fUseTOFpid = kFALSE;
2784       fTofPIDnSigmaBelowElectronLine=-100;
2785       fTofPIDnSigmaAboveElectronLine=100;
2786       break;
2787    case 1: // -7,7
2788       fUseTOFpid = kTRUE;
2789       fTofPIDnSigmaBelowElectronLine=-7;
2790       fTofPIDnSigmaAboveElectronLine=7;
2791       break;
2792    case 2: // -5,5
2793       fUseTOFpid = kTRUE;
2794       fTofPIDnSigmaBelowElectronLine=-5;
2795       fTofPIDnSigmaAboveElectronLine=5;
2796       break;
2797    case 3: // -3,5
2798       fUseTOFpid = kTRUE;
2799       fTofPIDnSigmaBelowElectronLine=-3;
2800       fTofPIDnSigmaAboveElectronLine=5;
2801       break;
2802    case 4: // -2,3
2803       fUseTOFpid = kTRUE;
2804       fTofPIDnSigmaBelowElectronLine=-2;
2805       fTofPIDnSigmaAboveElectronLine=3;
2806       break;
2807    case 5: // -3,3
2808       fUseTOFpid = kTRUE;
2809       fTofPIDnSigmaBelowElectronLine=-3;
2810       fTofPIDnSigmaAboveElectronLine=3;
2811       break;
2812    default:
2813       AliError(Form("TOFElectronCut not defined %d",TOFelectronPID));
2814       return kFALSE;
2815    }
2816    return kTRUE;
2817 }
2818 ///________________________________________________________________________
2819 Bool_t AliConversionCuts::SetQtMaxCut(Int_t QtMaxCut)
2820 {   // Set Cut
2821    switch(QtMaxCut){
2822    case 0: //
2823       fQtMax=1.;
2824       fDoQtGammaSelection=kFALSE;
2825       fDo2DQt=kFALSE;
2826       break;
2827    case 1:
2828       fQtMax=0.1;
2829       fDo2DQt=kFALSE;
2830       break;
2831    case 2:
2832       fQtMax=0.07;
2833       fDo2DQt=kFALSE;
2834       break;
2835    case 3:
2836       fQtMax=0.05;
2837       fDo2DQt=kFALSE;
2838       break;
2839    case 4:
2840       fQtMax=0.03;
2841       fDo2DQt=kFALSE;
2842       break;
2843    case 5:
2844       fQtMax=0.02;
2845       fDo2DQt=kFALSE;
2846       break;
2847    case 6:
2848       fQtMax=0.02;
2849       fDo2DQt=kTRUE;
2850       break;
2851    case 7:
2852       fQtMax=0.15;
2853       fDo2DQt=kFALSE;
2854       break;
2855    case 8:
2856       fQtMax=0.05;
2857       fDo2DQt=kTRUE;
2858       break;   
2859    case 9:
2860       fQtMax=0.03;
2861       fDo2DQt=kTRUE;
2862       break;      
2863    default:
2864       AliError(Form("Warning: QtMaxCut not defined %d",QtMaxCut));
2865       return kFALSE;
2866    }
2867    return kTRUE;
2868 }
2869 ///________________________________________________________________________
2870 Bool_t AliConversionCuts::SetChi2GammaCut(Int_t chi2GammaCut)
2871 {   // Set Cut
2872
2873    switch(chi2GammaCut){
2874    case 0: // 100
2875       fChi2CutConversion = 100.;
2876       break;
2877    case 1:  // 50
2878       fChi2CutConversion = 50.;
2879       break;
2880    case 2:  // 30
2881       fChi2CutConversion = 30.;
2882       break;
2883    case 3:
2884       fChi2CutConversion = 200.;
2885       break;
2886    case 4:
2887       fChi2CutConversion = 500.;
2888       break;
2889    case 5:
2890       fChi2CutConversion = 100000.;
2891       break;
2892    case 6:
2893       fChi2CutConversion = 5.;
2894       break;
2895    case 7:
2896       fChi2CutConversion = 10.;
2897       break;
2898    case 8:
2899       fChi2CutConversion = 20.;
2900       break;
2901    case 9:
2902       fChi2CutConversion = 15.;
2903       break;
2904    default:
2905       AliError(Form("Warning: Chi2GammaCut not defined %d",chi2GammaCut));
2906       return kFALSE;
2907    }
2908    return kTRUE;
2909 }
2910 ///________________________________________________________________________
2911 Bool_t AliConversionCuts::SetPsiPairCut(Int_t psiCut) {
2912
2913    switch(psiCut) {
2914    case 0:
2915       fPsiPairCut = 10000; //
2916       break;
2917    case 1:
2918       fPsiPairCut = 0.1; //
2919       break;
2920    case 2:
2921       fPsiPairCut = 0.05; // Standard
2922       break;
2923    case 3:
2924       fPsiPairCut = 0.035; //
2925       break;
2926    case 4:
2927       fPsiPairCut = 0.2; //
2928       break;   
2929    case 5:
2930       fPsiPairCut = 0.1; //
2931       fDo2DPsiPairChi2 = kTRUE;
2932       break;
2933    case 6:
2934       fPsiPairCut = 0.05; //
2935       fDo2DPsiPairChi2 = kTRUE;
2936       break;
2937    case 7:
2938       fPsiPairCut = 0.035; //
2939       fDo2DPsiPairChi2 = kTRUE;
2940       break;
2941    case 8:
2942       fPsiPairCut = 0.2; //
2943       fDo2DPsiPairChi2 = kTRUE; //
2944       break;
2945    case 9:
2946       fPsiPairCut = 0.5; //
2947       break;
2948    default:
2949       AliError(Form("PsiPairCut not defined %d",psiCut));
2950       return kFALSE;
2951    }
2952
2953    return kTRUE;
2954 }
2955 ///________________________________________________________________________
2956 Bool_t AliConversionCuts::SetPhotonAsymmetryCut(Int_t doPhotonAsymmetryCut){
2957    // Set Cut
2958    switch(doPhotonAsymmetryCut){
2959    case 0:
2960       fDoPhotonAsymmetryCut=0;
2961       fMinPPhotonAsymmetryCut=100.;
2962       fMinPhotonAsymmetry=0.;
2963       break;
2964    case 1:
2965       fDoPhotonAsymmetryCut=1;
2966       fMinPPhotonAsymmetryCut=3.5;
2967       fMinPhotonAsymmetry=0.04;
2968       break;
2969    case 2:
2970       fDoPhotonAsymmetryCut=1;
2971       fMinPPhotonAsymmetryCut=3.5;
2972       fMinPhotonAsymmetry=0.06;
2973       break;
2974    case 3:
2975       fDoPhotonAsymmetryCut=1;
2976       fMinPPhotonAsymmetryCut=0.0;
2977       fMinPhotonAsymmetry=0.05;
2978       break; 
2979    default:
2980       AliError(Form("PhotonAsymmetryCut not defined %d",doPhotonAsymmetryCut));
2981       return kFALSE;
2982    }
2983    fCuts[kdoPhotonAsymmetryCut]=doPhotonAsymmetryCut;
2984    return kTRUE;
2985 }
2986 ///________________________________________________________________________
2987 Bool_t AliConversionCuts::SetCosPAngleCut(Int_t cosCut) {
2988
2989    switch(cosCut){
2990    case 0:
2991       fCosPAngleCut = -1; 
2992       break;
2993    case 1:
2994       fCosPAngleCut = 0; 
2995       break;
2996    case 2:
2997       fCosPAngleCut = 0.5; 
2998       break;
2999    case 3:
3000       fCosPAngleCut = 0.75; 
3001       break;
3002    case 4:
3003       fCosPAngleCut = 0.85; 
3004       break;
3005    case 5:
3006       fCosPAngleCut = 0.88; 
3007       break;
3008    case 6:
3009       fCosPAngleCut = 0.9;
3010       break;
3011    case 7:
3012       fCosPAngleCut = 0.95;
3013       break;
3014    default:
3015       AliError(Form("Cosine Pointing Angle cut not defined %d",cosCut));
3016       return kFALSE;
3017    }
3018
3019    return kTRUE;
3020 }
3021 ///________________________________________________________________________
3022 Bool_t AliConversionCuts::SetSharedElectronCut(Int_t sharedElec) {
3023
3024         switch(sharedElec){
3025         case 0:
3026                 fDoSharedElecCut = kFALSE;
3027                 fDoPhotonQualitySelectionCut = kFALSE;
3028                 fPhotonQualityCut = 0;
3029                 break;
3030         case 1:
3031                 fDoSharedElecCut = kTRUE;
3032                 fDoPhotonQualitySelectionCut = kFALSE;
3033                 fPhotonQualityCut = 0;
3034                 break;
3035         case 2:
3036                 fDoSharedElecCut = kFALSE;
3037                 fDoPhotonQualitySelectionCut = kTRUE;
3038                 fPhotonQualityCut = 1;
3039                 break;
3040         case 3:
3041                 fDoSharedElecCut = kFALSE;
3042                 fDoPhotonQualitySelectionCut = kTRUE;     
3043                 fPhotonQualityCut = 2;
3044                 break;
3045         case 4:
3046                 fDoSharedElecCut = kFALSE;
3047                 fDoPhotonQualitySelectionCut = kTRUE;     
3048                 fPhotonQualityCut = 3;
3049                 break;
3050         default:
3051                 AliError(Form("Shared Electron Cut not defined %d",sharedElec));        
3052                 return kFALSE;
3053    }
3054
3055    return kTRUE;
3056 }
3057 ///________________________________________________________________________
3058 Bool_t AliConversionCuts::SetToCloseV0sCut(Int_t toClose) {
3059
3060    switch(toClose){
3061    case 0:
3062       fDoToCloseV0sCut = kFALSE;
3063       fminV0Dist = 250;
3064       break;
3065    case 1:
3066       fDoToCloseV0sCut = kTRUE;
3067       fminV0Dist = 1;
3068       break;
3069    case 2:
3070       fDoToCloseV0sCut = kTRUE;
3071       fminV0Dist = 2;
3072       break;
3073    case 3:
3074       fDoToCloseV0sCut = kTRUE;
3075       fminV0Dist = 3;
3076       break;
3077    default:
3078       AliError(Form("Shared Electron Cut not defined %d",toClose));
3079       return kFALSE;
3080    }
3081    return kTRUE;
3082 }
3083 ///________________________________________________________________________
3084 Bool_t AliConversionCuts::SetTRDElectronCut(Int_t TRDElectronCut)
3085 {   // Set Cut
3086    switch(TRDElectronCut){
3087    case 0:
3088       fDoTRDPID=kFALSE;
3089       break;
3090    case 1:
3091       fDoTRDPID=kTRUE;
3092       fPIDTRDEfficiency=0.1;
3093       break;
3094    case 8:
3095       fDoTRDPID=kTRUE;
3096       fPIDTRDEfficiency=0.8;
3097       break;
3098    case 9:
3099       fDoTRDPID=kTRUE;
3100       fPIDTRDEfficiency=0.9;
3101       break;
3102    default:
3103       AliError(Form("TRDElectronCut not defined %d",TRDElectronCut));
3104       return kFALSE;
3105    }
3106
3107    return kTRUE;
3108 }
3109
3110 ///________________________________________________________________________
3111 Bool_t AliConversionCuts::SetDCAZPhotonPrimVtxCut(Int_t DCAZPhotonPrimVtx){
3112    // Set Cut
3113    switch(DCAZPhotonPrimVtx){
3114    case 0:  //
3115       fDCAZPrimVtxCut   = 1000;
3116       break;
3117    case 1:  //
3118       fDCAZPrimVtxCut   = 10;
3119       break;
3120    case 2:  //
3121       fDCAZPrimVtxCut   = 5;
3122       break;
3123    case 3:  //
3124       fDCAZPrimVtxCut   = 4;
3125       break;
3126    case 4:  //
3127       fDCAZPrimVtxCut   = 3;
3128       break;
3129    case 5:  //
3130       fDCAZPrimVtxCut   = 2.5;
3131       break;
3132    case 6:  //
3133       fDCAZPrimVtxCut   = 2;
3134       break;
3135    case 7:  //
3136       fDCAZPrimVtxCut   = 1.5;
3137       break;
3138    case 8:  //
3139       fDCAZPrimVtxCut   = 1;
3140       break;
3141    case 9:  //
3142       fDCAZPrimVtxCut   = 0.5;
3143       break;
3144    default:
3145       cout<<"Warning: DCAZPhotonPrimVtx not defined "<<DCAZPhotonPrimVtx<<endl;
3146       return kFALSE;
3147    }
3148    return kTRUE;
3149 }
3150
3151 ///________________________________________________________________________
3152 Bool_t AliConversionCuts::SetDCARPhotonPrimVtxCut(Int_t DCARPhotonPrimVtx){
3153    // Set Cut
3154    switch(DCARPhotonPrimVtx){
3155    case 0:  //
3156       fDCARPrimVtxCut   = 1000;
3157       break;
3158    case 1:  //
3159       fDCARPrimVtxCut   = 10;
3160       break;
3161    case 2:  //
3162       fDCARPrimVtxCut   = 5;
3163       break;
3164    case 3:  //
3165       fDCARPrimVtxCut   = 4;
3166       break;
3167    case 4:  //
3168       fDCARPrimVtxCut   = 3;
3169       break;
3170    case 5:  //
3171       fDCARPrimVtxCut   = 2.5;
3172       break;
3173    case 6:  //
3174       fDCARPrimVtxCut   = 2;
3175       break;
3176    case 7:  //
3177       fDCARPrimVtxCut   = 1.5;
3178       break;
3179    case 8:  //
3180       fDCARPrimVtxCut   = 1;
3181       break;
3182    case 9:  //
3183       fDCARPrimVtxCut   = 0.5;
3184       break;
3185    default:
3186       cout<<"Warning: DCARPhotonPrimVtx not defined "<<DCARPhotonPrimVtx<<endl;
3187       return kFALSE;
3188    }
3189    return kTRUE;
3190 }
3191
3192 ///________________________________________________________________________
3193 Bool_t AliConversionCuts::SetInPlaneOutOfPlane(Int_t inOutPlane){
3194    // Set Cut
3195    switch(inOutPlane){
3196    case 0:  //
3197       fInPlaneOutOfPlane = 0; // No Event Plane
3198       break;
3199    case 1:  //
3200       fInPlaneOutOfPlane = 1; // In-Plane
3201       break;
3202    case 2:  //
3203       fInPlaneOutOfPlane = 2; // Out-Of-Plane
3204       break;
3205    default:
3206       cout<<"Warning: In-Plane or Out-Of-Plane not defined "<<inOutPlane<<endl;
3207       return kFALSE;
3208    }
3209    return kTRUE;
3210 }
3211
3212
3213 //-------------------------------------------------------------
3214 Double_t AliConversionCuts::GetCentrality(AliVEvent *event)
3215 {   // Get Event Centrality
3216
3217    AliESDEvent *esdEvent=dynamic_cast<AliESDEvent*>(event);
3218    if(esdEvent){
3219       AliCentrality *fESDCentrality=(AliCentrality*)esdEvent->GetCentrality();
3220
3221       if(fDetectorCentrality==0){
3222          if (fIsHeavyIon==2){
3223             return fESDCentrality->GetCentralityPercentile("V0A"); // default for pPb
3224          } else{
3225             return fESDCentrality->GetCentralityPercentile("V0M"); // default
3226          }
3227       }
3228       if(fDetectorCentrality==1){
3229          return fESDCentrality->GetCentralityPercentile("CL1");
3230       }
3231    }
3232
3233    AliAODEvent *aodEvent=dynamic_cast<AliAODEvent*>(event);
3234    if(aodEvent){
3235       if(aodEvent->GetHeader()){return aodEvent->GetHeader()->GetCentrality();}
3236    }
3237
3238    return -1;
3239 }
3240 //-------------------------------------------------------------
3241 Bool_t AliConversionCuts::IsCentralitySelected(AliVEvent *event, AliVEvent *fMCEvent)
3242 {   // Centrality Selection
3243    if(!fIsHeavyIon)return kTRUE;
3244
3245    if(fCentralityMin == fCentralityMax ) return kTRUE;//0-100%
3246    else if(fCentralityMax==0) fCentralityMax=10; //CentralityRange = fCentralityMin-100%
3247
3248    Double_t centrality=GetCentrality(event);
3249    if(centrality<0)return kFALSE;
3250
3251    Int_t centralityC=0;
3252    if (fModCentralityClass == 0){
3253       centralityC= Int_t(centrality/10);
3254       if(centralityC >= fCentralityMin && centralityC < fCentralityMax)
3255          return kTRUE;
3256       else return kFALSE;
3257    }
3258    else if (fModCentralityClass ==1){
3259       centralityC= Int_t(centrality);
3260       if(centralityC >= fCentralityMin*5 && centralityC < fCentralityMax*5){
3261          return kTRUE;
3262       } else return kFALSE;
3263    }
3264    else if (fModCentralityClass ==2){
3265       centralityC= Int_t(centrality);
3266       if(centralityC >= ((fCentralityMin*5)+45) && centralityC < ((fCentralityMax*5)+45))
3267          return kTRUE;
3268       else return kFALSE;
3269    }
3270
3271    Int_t nprimaryTracks = ((AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()->GetTask("V0ReaderV1"))->GetNumberOfPrimaryTracks();
3272    Int_t PrimaryTracks10[10][2] =
3273       {
3274          {9999,9999}, //  0
3275          {1210, 928}, // 10
3276          { 817, 658}, // 20
3277          { 536, 435}, // 30
3278          { 337, 276}, // 40
3279          { 197, 162}, // 50
3280          { 106, 100}, // 60
3281          {  51,  44}, // 70
3282          {  21,  18}, // 80
3283          {   0,   0}  // 90
3284       };
3285    Int_t PrimaryTracks5a[10][2] =
3286       {
3287          {9999,9999}, // 0
3288          {1485,1168}, // 5
3289          {1210, 928}, // 10
3290          { 995, 795}, // 15
3291          { 817, 658}, // 20
3292          { 666, 538}, // 25
3293          { 536, 435}, // 30
3294          { 428, 350}, // 35
3295          { 337, 276}, // 40
3296          { 260, 214}  // 45
3297       };
3298    Int_t PrimaryTracks5b[10][2] =
3299       {
3300          { 260, 214}, // 45
3301          { 197, 162}, // 50
3302          { 147, 125}, // 55
3303          { 106, 100}, // 60
3304          {  75,  63}, // 65
3305          {  51,  44}, // 70
3306          {  34,  29}, // 75
3307          {  21,  18}, // 80
3308          {  13,  11}, // 85
3309          {   0,   0}  // 90
3310       };
3311
3312    Int_t column = -1;
3313    if(event->IsA()==AliESDEvent::Class()) column = 0;
3314    if(event->IsA()==AliAODEvent::Class()) column = 1;
3315
3316    if (fModCentralityClass == 3){
3317       if(fMCEvent){
3318          if(nprimaryTracks > PrimaryTracks10[fCentralityMax][column] && nprimaryTracks <= PrimaryTracks10[fCentralityMin][column])
3319             return kTRUE;
3320          else return kFALSE;
3321       }
3322       else{
3323          centralityC= Int_t(centrality/10);
3324          if(centralityC >= fCentralityMin && centralityC < fCentralityMax)
3325             return kTRUE;
3326          else return kFALSE;
3327       }
3328    }
3329    else if (fModCentralityClass ==4){
3330       if(fMCEvent){
3331          if(nprimaryTracks > PrimaryTracks5a[fCentralityMax][column] && nprimaryTracks <= PrimaryTracks5a[fCentralityMin][column])
3332             return kTRUE;
3333          else return kFALSE;
3334       }
3335       else{
3336          centralityC= Int_t(centrality);
3337          if(centralityC >= fCentralityMin*5 && centralityC < fCentralityMax*5){
3338             return kTRUE;
3339          } else return kFALSE;
3340       }
3341    }
3342    else if (fModCentralityClass ==5){
3343       if(fMCEvent){
3344          if(nprimaryTracks > PrimaryTracks5b[fCentralityMax][column] && nprimaryTracks <= PrimaryTracks5b[fCentralityMin][column])
3345             return kTRUE;
3346          else return kFALSE;
3347       }
3348       else{
3349          centralityC= Int_t(centrality);
3350          if(centralityC >= ((fCentralityMin*5)+45) && centralityC < ((fCentralityMax*5)+45))
3351             return kTRUE;
3352          else return kFALSE;
3353       }
3354    }
3355
3356    return kFALSE;
3357 }
3358 ///________________________________________________________________________
3359 Bool_t AliConversionCuts::VertexZCut(AliVEvent *event){
3360    // Cut on z position of primary vertex
3361    Double_t fVertexZ=event->GetPrimaryVertex()->GetZ();
3362    Double_t fVertexZSPD = 0;
3363    AliESDEvent *fESDEvent=dynamic_cast<AliESDEvent*>(event);
3364    if(fESDEvent){
3365       fVertexZSPD = fESDEvent->GetPrimaryVertexSPD()->GetZ();
3366    } 
3367    AliAODEvent *fAODEvent=dynamic_cast<AliAODEvent*>(event);
3368    if(fAODEvent){
3369       fVertexZSPD = fAODEvent->GetPrimaryVertexSPD()->GetZ();
3370    }
3371   
3372    if(abs(fVertexZ)>fMaxVertexZ)return kFALSE;
3373
3374    TString periodName = ((AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()
3375                                                 ->GetTask("V0ReaderV1"))->GetPeriodName();
3376    if (periodName.CompareTo("LHC11h")==0){
3377      if (abs(fVertexZ-fVertexZSPD) > 0.1) return kFALSE;
3378    }                                            
3379    if (fIsHeavyIon == 2){
3380      if(fUtils->IsFirstEventInChunk(event)) return kFALSE;
3381      if(!fUtils->IsVertexSelected2013pA(event)) return kFALSE;
3382      if(fUtils->IsPileUpEvent(event)) return kFALSE;
3383    }
3384
3385    return kTRUE;
3386 }
3387 ///________________________________________________________________________
3388
3389 Int_t AliConversionCuts::GetNumberOfContributorsVtx(AliVEvent *event){
3390    // returns number of contributors to the vertex
3391
3392    AliESDEvent *fESDEvent=dynamic_cast<AliESDEvent*>(event);
3393    if(fESDEvent){
3394       if (fESDEvent->GetPrimaryVertex() != NULL){
3395          if(fESDEvent->GetPrimaryVertex()->GetNContributors()>0) {
3396 //     cout << "accepted global" << fESDEvent->GetEventNumberInFile() << " with NCont: " << fESDEvent->GetPrimaryVertex()->GetNContributors() << endl;
3397             return fESDEvent->GetPrimaryVertex()->GetNContributors();
3398          }
3399       }
3400
3401       if(fESDEvent->GetPrimaryVertexSPD() !=NULL){
3402          if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors()>0) {
3403 //     cout << "accepted SPD" << fESDEvent->GetEventNumberInFile() << " with NCont: " << fESDEvent->GetPrimaryVertexSPD()->GetNContributors() << endl;
3404             return fESDEvent->GetPrimaryVertexSPD()->GetNContributors();
3405          }  else {
3406             AliWarning(Form("Number of contributors from bad vertex type:: %s",fESDEvent->GetPrimaryVertex()->GetName()));
3407 //            cout << "rejected " << fESDEvent->GetEventNumberInFile() << endl;
3408             return 0;
3409          }
3410       }
3411    }
3412
3413    AliAODEvent *fAODEvent=dynamic_cast<AliAODEvent*>(event);
3414    if(fAODEvent){
3415       if (fAODEvent->GetPrimaryVertex() != NULL){
3416          if(fAODEvent->GetPrimaryVertex()->GetNContributors()>0) {
3417             return fAODEvent->GetPrimaryVertex()->GetNContributors();
3418          }
3419       }
3420       if(fAODEvent->GetPrimaryVertexSPD() !=NULL){
3421          if(fAODEvent->GetPrimaryVertexSPD()->GetNContributors()>0) {
3422             return fAODEvent->GetPrimaryVertexSPD()->GetNContributors();
3423          } else {
3424             AliWarning(Form("Number of contributors from bad vertex type:: %s",fAODEvent->GetPrimaryVertex()->GetName()));
3425             return 0;
3426          }
3427       }
3428    }
3429   // cout << "rejected " << fESDEvent->GetEventNumberInFile() << endl;
3430    return 0;
3431 }
3432
3433 ///________________________________________________________________________
3434
3435 Bool_t AliConversionCuts::IsTriggerSelected(AliVEvent *fInputEvent)
3436 {
3437
3438    AliInputEventHandler *fInputHandler=(AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
3439
3440    UInt_t isSelected = AliVEvent::kAny;
3441    TString periodName = ((AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()->GetTask("V0ReaderV1"))->GetPeriodName();
3442 //    cout <<   periodName.Data() << endl;
3443    
3444    if (fInputHandler==NULL) return kFALSE;
3445    if( fInputHandler->GetEventSelection() || fInputEvent->IsA()==AliAODEvent::Class()) {
3446       if (!fTriggerSelectedManually){
3447          if (fPreSelCut) fOfflineTriggerMask = AliVEvent::kAny;
3448          else {
3449             if (fIsHeavyIon == 1) fOfflineTriggerMask = AliVEvent::kMB | AliVEvent::kCentral | AliVEvent::kSemiCentral;
3450             else if (fIsHeavyIon == 2) fOfflineTriggerMask = AliVEvent::kINT7;
3451             else if (periodName.CompareTo("LHC11c") == 0 || periodName.CompareTo("LHC11d") == 0 || periodName.CompareTo("LHC11e") == 0 || periodName.CompareTo("LHC11f") == 0 || periodName.CompareTo("LHC11g") == 0  || periodName.CompareTo("LHC12a") == 0 || periodName.CompareTo("LHC12b") == 0 || periodName.CompareTo("LHC12c") == 0 || periodName.CompareTo("LHC12d") == 0 || periodName.CompareTo("LHC11f") == 0  || periodName.CompareTo("LHC13g") == 0 ) fOfflineTriggerMask = AliVEvent::kINT7;      
3452             else fOfflineTriggerMask = AliVEvent::kMB;
3453          }
3454       }
3455       // Get the actual offline trigger mask for the event and AND it with the
3456       // requested mask. If no mask requested select by default the event.
3457 //       if (fPreSelCut) cout << "Trigger selected from outside: "<< fTriggerSelectedManually <<"\t Offline Trigger mask for Precut: " << fOfflineTriggerMask << endl;
3458 //       else cout << "Trigger selected from outside: "<< fTriggerSelectedManually <<"\t Offline Trigger mask: " << fOfflineTriggerMask << endl;
3459
3460       if (fOfflineTriggerMask)
3461          isSelected = fOfflineTriggerMask & fInputHandler->IsEventSelected();
3462    }
3463    fIsSDDFired = !(fInputHandler->IsEventSelected() & AliVEvent::kFastOnly);
3464
3465    // Fill Histogram
3466    if(hTriggerClass){
3467       if (fIsSDDFired) hTriggerClass->Fill(33);
3468       if (fInputHandler->IsEventSelected() & AliVEvent::kMB)hTriggerClass->Fill(0);
3469       if (fInputHandler->IsEventSelected() & AliVEvent::kINT7)hTriggerClass->Fill(1);
3470       if (fInputHandler->IsEventSelected() & AliVEvent::kMUON)hTriggerClass->Fill(2);
3471       if (fInputHandler->IsEventSelected() & AliVEvent::kHighMult)hTriggerClass->Fill(3);
3472       if (fInputHandler->IsEventSelected() & AliVEvent::kEMC1)hTriggerClass->Fill(4);
3473       if (fInputHandler->IsEventSelected() & AliVEvent::kCINT5)hTriggerClass->Fill(5);
3474       if (fInputHandler->IsEventSelected() & AliVEvent::kCMUS5)hTriggerClass->Fill(6);
3475       if (fInputHandler->IsEventSelected() & AliVEvent::kMUSPB)hTriggerClass->Fill(6);
3476       if (fInputHandler->IsEventSelected() & AliVEvent::kMUSH7)hTriggerClass->Fill(7);
3477       if (fInputHandler->IsEventSelected() & AliVEvent::kMUSHPB)hTriggerClass->Fill(7);
3478       if (fInputHandler->IsEventSelected() & AliVEvent::kMUL7)hTriggerClass->Fill(8);
3479       if (fInputHandler->IsEventSelected() & AliVEvent::kMuonLikePB)hTriggerClass->Fill(8);
3480       if (fInputHandler->IsEventSelected() & AliVEvent::kMUU7)hTriggerClass->Fill(9);
3481       if (fInputHandler->IsEventSelected() & AliVEvent::kMuonUnlikePB)hTriggerClass->Fill(9);
3482       if (fInputHandler->IsEventSelected() & AliVEvent::kEMC7)hTriggerClass->Fill(10);
3483       if (fInputHandler->IsEventSelected() & AliVEvent::kEMC8)hTriggerClass->Fill(10);
3484       if (fInputHandler->IsEventSelected() & AliVEvent::kMUS7)hTriggerClass->Fill(11);
3485       if (fInputHandler->IsEventSelected() & AliVEvent::kPHI1)hTriggerClass->Fill(12);
3486       if (fInputHandler->IsEventSelected() & AliVEvent::kPHI7)hTriggerClass->Fill(13);
3487       if (fInputHandler->IsEventSelected() & AliVEvent::kPHI8)hTriggerClass->Fill(13);
3488       if (fInputHandler->IsEventSelected() & AliVEvent::kPHOSPb)hTriggerClass->Fill(13);
3489       if (fInputHandler->IsEventSelected() & AliVEvent::kEMCEJE)hTriggerClass->Fill(14);
3490       if (fInputHandler->IsEventSelected() & AliVEvent::kEMCEGA)hTriggerClass->Fill(15);
3491       if (fInputHandler->IsEventSelected() & AliVEvent::kCentral)hTriggerClass->Fill(16);
3492       if (fInputHandler->IsEventSelected() & AliVEvent::kSemiCentral)hTriggerClass->Fill(17);
3493       if (fInputHandler->IsEventSelected() & AliVEvent::kDG5)hTriggerClass->Fill(18);
3494       if (fInputHandler->IsEventSelected() & AliVEvent::kZED)hTriggerClass->Fill(19);
3495       if (fInputHandler->IsEventSelected() & AliVEvent::kSPI7)hTriggerClass->Fill(20);
3496       if (fInputHandler->IsEventSelected() & AliVEvent::kSPI)hTriggerClass->Fill(20);
3497       if (fInputHandler->IsEventSelected() & AliVEvent::kINT8)hTriggerClass->Fill(21);
3498       if (fInputHandler->IsEventSelected() & AliVEvent::kMuonSingleLowPt8)hTriggerClass->Fill(22);
3499       if (fInputHandler->IsEventSelected() & AliVEvent::kMuonSingleHighPt8)hTriggerClass->Fill(23);
3500       if (fInputHandler->IsEventSelected() & AliVEvent::kMuonLikeLowPt8)hTriggerClass->Fill(24);
3501       if (fInputHandler->IsEventSelected() & AliVEvent::kMuonUnlikeLowPt8)hTriggerClass->Fill(25);
3502       if (fInputHandler->IsEventSelected() & AliVEvent::kMuonUnlikeLowPt0)hTriggerClass->Fill(26);
3503       if (fInputHandler->IsEventSelected() & AliVEvent::kUserDefined)hTriggerClass->Fill(27);
3504       if (fInputHandler->IsEventSelected() & AliVEvent::kTRD)hTriggerClass->Fill(28);
3505       if (fInputHandler->IsEventSelected() & AliVEvent::kFastOnly)hTriggerClass->Fill(29);
3506       if (fInputHandler->IsEventSelected() & AliVEvent::kAnyINT)hTriggerClass->Fill(30);
3507       if (fInputHandler->IsEventSelected() & AliVEvent::kAny)hTriggerClass->Fill(31);
3508       if (!fInputHandler->IsEventSelected()) hTriggerClass->Fill(34);
3509    }
3510
3511    if(hTriggerClassSelected && isSelected){
3512       if (!fIsSDDFired) hTriggerClassSelected->Fill(33);
3513       if (fInputHandler->IsEventSelected() & AliVEvent::kMB)hTriggerClassSelected->Fill(0);
3514       if (fInputHandler->IsEventSelected() & AliVEvent::kINT7)hTriggerClassSelected->Fill(1);
3515       if (fInputHandler->IsEventSelected() & AliVEvent::kMUON)hTriggerClassSelected->Fill(2);
3516       if (fInputHandler->IsEventSelected() & AliVEvent::kHighMult)hTriggerClassSelected->Fill(3);
3517       if (fInputHandler->IsEventSelected() & AliVEvent::kEMC1)hTriggerClassSelected->Fill(4);
3518       if (fInputHandler->IsEventSelected() & AliVEvent::kCINT5)hTriggerClassSelected->Fill(5);
3519       if (fInputHandler->IsEventSelected() & AliVEvent::kCMUS5)hTriggerClassSelected->Fill(6);
3520       if (fInputHandler->IsEventSelected() & AliVEvent::kMUSPB)hTriggerClassSelected->Fill(6);
3521       if (fInputHandler->IsEventSelected() & AliVEvent::kMUSH7)hTriggerClassSelected->Fill(7);
3522       if (fInputHandler->IsEventSelected() & AliVEvent::kMUSHPB)hTriggerClassSelected->Fill(7);
3523       if (fInputHandler->IsEventSelected() & AliVEvent::kMUL7)hTriggerClassSelected->Fill(8);
3524       if (fInputHandler->IsEventSelected() & AliVEvent::kMuonLikePB)hTriggerClassSelected->Fill(8);
3525       if (fInputHandler->IsEventSelected() & AliVEvent::kMUU7)hTriggerClassSelected->Fill(9);
3526       if (fInputHandler->IsEventSelected() & AliVEvent::kMuonUnlikePB)hTriggerClassSelected->Fill(9);
3527       if (fInputHandler->IsEventSelected() & AliVEvent::kEMC7)hTriggerClassSelected->Fill(10);
3528       if (fInputHandler->IsEventSelected() & AliVEvent::kEMC8)hTriggerClassSelected->Fill(10);
3529       if (fInputHandler->IsEventSelected() & AliVEvent::kMUS7)hTriggerClassSelected->Fill(11);
3530       if (fInputHandler->IsEventSelected() & AliVEvent::kPHI1)hTriggerClassSelected->Fill(12);
3531       if (fInputHandler->IsEventSelected() & AliVEvent::kPHI7)hTriggerClassSelected->Fill(13);
3532       if (fInputHandler->IsEventSelected() & AliVEvent::kPHI8)hTriggerClassSelected->Fill(13);
3533       if (fInputHandler->IsEventSelected() & AliVEvent::kPHOSPb)hTriggerClassSelected->Fill(13);
3534       if (fInputHandler->IsEventSelected() & AliVEvent::kEMCEJE)hTriggerClassSelected->Fill(14);
3535       if (fInputHandler->IsEventSelected() & AliVEvent::kEMCEGA)hTriggerClassSelected->Fill(15);
3536       if (fInputHandler->IsEventSelected() & AliVEvent::kCentral)hTriggerClassSelected->Fill(16);
3537       if (fInputHandler->IsEventSelected() & AliVEvent::kSemiCentral)hTriggerClassSelected->Fill(17);
3538       if (fInputHandler->IsEventSelected() & AliVEvent::kDG5)hTriggerClassSelected->Fill(18);
3539       if (fInputHandler->IsEventSelected() & AliVEvent::kZED)hTriggerClassSelected->Fill(19);
3540       if (fInputHandler->IsEventSelected() & AliVEvent::kSPI7)hTriggerClassSelected->Fill(20);
3541       if (fInputHandler->IsEventSelected() & AliVEvent::kSPI)hTriggerClassSelected->Fill(20);
3542       if (fInputHandler->IsEventSelected() & AliVEvent::kINT8)hTriggerClassSelected->Fill(21);
3543       if (fInputHandler->IsEventSelected() & AliVEvent::kMuonSingleLowPt8)hTriggerClassSelected->Fill(22);
3544       if (fInputHandler->IsEventSelected() & AliVEvent::kMuonSingleHighPt8)hTriggerClassSelected->Fill(23);
3545       if (fInputHandler->IsEventSelected() & AliVEvent::kMuonLikeLowPt8)hTriggerClassSelected->Fill(24);
3546       if (fInputHandler->IsEventSelected() & AliVEvent::kMuonUnlikeLowPt8)hTriggerClassSelected->Fill(25);
3547       if (fInputHandler->IsEventSelected() & AliVEvent::kMuonUnlikeLowPt0)hTriggerClassSelected->Fill(26);
3548       if (fInputHandler->IsEventSelected() & AliVEvent::kUserDefined)hTriggerClassSelected->Fill(27);
3549       if (fInputHandler->IsEventSelected() & AliVEvent::kTRD)hTriggerClassSelected->Fill(28);
3550       if (fInputHandler->IsEventSelected() & AliVEvent::kFastOnly)hTriggerClassSelected->Fill(29);
3551       if (fInputHandler->IsEventSelected() & AliVEvent::kAnyINT)hTriggerClassSelected->Fill(30);
3552       if (fInputHandler->IsEventSelected() & AliVEvent::kAny)hTriggerClassSelected->Fill(31);
3553    }
3554
3555    if(!isSelected)return kFALSE;
3556
3557    return kTRUE;
3558
3559 }
3560
3561 ///________________________________________________________________________
3562 Int_t AliConversionCuts::GetFirstTPCRow(Double_t radius){
3563    // Get first TPC row
3564    Int_t firstTPCRow = 0;
3565    Double_t radiusI = 84.8;
3566    Double_t radiusO = 134.6;
3567    Double_t radiusOB = 198.;
3568    Double_t rSizeI = 0.75;
3569    Double_t rSizeO = 1.;
3570    Double_t rSizeOB = 1.5;
3571    Int_t nClsI = 63;
3572    Int_t nClsIO = 127;
3573
3574    if(radius <= radiusI){
3575       return firstTPCRow;
3576    }
3577    if(radius>radiusI && radius<=radiusO){
3578       firstTPCRow = (Int_t)((radius-radiusI)/rSizeI);
3579    }
3580    if(radius>radiusO && radius<=radiusOB){
3581       firstTPCRow = (Int_t)(nClsI+(radius-radiusO)/rSizeO);
3582    }
3583
3584    if(radius>radiusOB){
3585       firstTPCRow =(Int_t)(nClsIO+(radius-radiusOB)/rSizeOB);
3586    }
3587
3588    return firstTPCRow;
3589 }
3590
3591 Bool_t AliConversionCuts::CosinePAngleCut(const AliConversionPhotonBase * photon, AliVEvent * event) const {
3592    ///Check if passes cosine of pointing angle cut
3593    if(GetCosineOfPointingAngle(photon, event) < fCosPAngleCut){
3594       return kFALSE;
3595    }
3596    return kTRUE;
3597 }
3598
3599 Double_t AliConversionCuts::GetCosineOfPointingAngle( const AliConversionPhotonBase * photon, AliVEvent * event) const{
3600    // calculates the pointing angle of the recalculated V0
3601
3602    Double_t momV0[3] = {0,0,0};
3603    if(event->IsA()==AliESDEvent::Class()){
3604       AliESDEvent *esdEvent = dynamic_cast<AliESDEvent*>(event);
3605       if(!esdEvent) return -999;
3606       AliESDv0 *v0 = esdEvent->GetV0(photon->GetV0Index());
3607       if(!v0) return -999;
3608       v0->GetPxPyPz(momV0[0],momV0[1],momV0[2]);
3609    }
3610    if(event->IsA()==AliAODEvent::Class()){
3611       momV0[0] = photon->GetPx();
3612       momV0[1] = photon->GetPy();
3613       momV0[2] = photon->GetPz();
3614    }
3615
3616    //Double_t momV0[3] = { photon->GetPx(), photon->GetPy(), photon->GetPz() }; //momentum of the V0
3617    Double_t PosV0[3] = { photon->GetConversionX() - event->GetPrimaryVertex()->GetX(),
3618                          photon->GetConversionY() - event->GetPrimaryVertex()->GetY(),
3619                          photon->GetConversionZ() - event->GetPrimaryVertex()->GetZ() }; //Recalculated V0 Position vector
3620
3621    Double_t momV02 = momV0[0]*momV0[0] + momV0[1]*momV0[1] + momV0[2]*momV0[2];
3622    Double_t PosV02 = PosV0[0]*PosV0[0] + PosV0[1]*PosV0[1] + PosV0[2]*PosV0[2];
3623
3624
3625    Double_t cosinePointingAngle = -999;
3626    if(momV02*PosV02 > 0.0)
3627       cosinePointingAngle = (PosV0[0]*momV0[0] +  PosV0[1]*momV0[1] + PosV0[2]*momV0[2] ) / TMath::Sqrt(momV02 * PosV02);
3628
3629    return cosinePointingAngle;
3630 }
3631
3632 ///________________________________________________________________________
3633 Bool_t AliConversionCuts::PsiPairCut(const AliConversionPhotonBase * photon) const {
3634
3635    if (fDo2DPsiPairChi2){
3636       if (abs(photon->GetPsiPair()) < -fPsiPairCut/fChi2CutConversion*photon->GetChi2perNDF() + fPsiPairCut ){  
3637          return kTRUE;
3638       } else {
3639          return kFALSE;
3640       }    
3641    } else {
3642       if(abs(photon->GetPsiPair()) > fPsiPairCut){
3643          return kFALSE;}
3644       else{return kTRUE;}
3645    } 
3646 }
3647
3648 ///________________________________________________________________________
3649 TString AliConversionCuts::GetCutNumber(){
3650    // returns TString with current cut number
3651    TString a(kNCuts);
3652    for(Int_t ii=0;ii<kNCuts;ii++){
3653       a.Append(Form("%d",fCuts[ii]));
3654    }
3655    return a;
3656 }
3657
3658 ///________________________________________________________________________
3659 void AliConversionCuts::FillElectonLabelArray(AliAODConversionPhoton* photon, Int_t nV0){
3660
3661    Int_t posLabel = photon->GetTrackLabelPositive();
3662    Int_t negLabel = photon->GetTrackLabelNegative();
3663
3664    fElectronLabelArray[nV0*2] = posLabel;
3665    fElectronLabelArray[(nV0*2)+1] = negLabel;
3666 }
3667 ///________________________________________________________________________
3668 Bool_t AliConversionCuts::RejectSharedElectronV0s(AliAODConversionPhoton* photon, Int_t nV0, Int_t nV0s){
3669
3670    Int_t posLabel = photon->GetTrackLabelPositive();
3671    Int_t negLabel = photon->GetTrackLabelNegative();
3672
3673    for(Int_t i = 0; i<nV0s*2;i++){
3674       if(i==nV0*2)     continue;
3675       if(i==(nV0*2)+1) continue;
3676       if(fElectronLabelArray[i] == posLabel){
3677          return kFALSE;}
3678       if(fElectronLabelArray[i] == negLabel){
3679          return kFALSE;}
3680    }
3681
3682    return kTRUE;
3683 }
3684 ///________________________________________________________________________
3685 Bool_t AliConversionCuts::RejectToCloseV0s(AliAODConversionPhoton* photon, TList *photons, Int_t nV0){
3686
3687
3688    Double_t posX = photon->GetConversionX();
3689    Double_t posY = photon->GetConversionY();
3690    Double_t posZ = photon->GetConversionZ();
3691
3692    for(Int_t i = 0;i<photons->GetEntries();i++){
3693       if(nV0 == i) continue;
3694       AliAODConversionPhoton *photonComp = (AliAODConversionPhoton*) photons->At(i);
3695       Double_t posCompX = photonComp->GetConversionX();
3696       Double_t posCompY = photonComp->GetConversionY();
3697       Double_t posCompZ = photonComp->GetConversionZ();
3698
3699       Double_t dist = pow((posX - posCompX),2)+pow((posY - posCompY),2)+pow((posZ - posCompZ),2);
3700
3701       if(dist < fminV0Dist*fminV0Dist){
3702          if(photon->GetChi2perNDF() < photonComp->GetChi2perNDF()) return kTRUE;
3703          else {
3704             return kFALSE;}
3705       }
3706
3707    }
3708    return kTRUE;
3709 }
3710 ///________________________________________________________________________
3711 void AliConversionCuts::GetNotRejectedParticles(Int_t rejection, TList *HeaderList, AliVEvent *MCEvent){
3712
3713         TString periodName = ((AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()->GetTask("V0ReaderV1"))->GetPeriodName();
3714
3715
3716         if(fNotRejectedStart){
3717                 delete[] fNotRejectedStart;
3718                 fNotRejectedStart = NULL;
3719         }
3720         if(fNotRejectedEnd){
3721                 delete[] fNotRejectedEnd;
3722                 fNotRejectedEnd = NULL;
3723         }
3724         if(fGeneratorNames){
3725                 delete[] fGeneratorNames;
3726                 fGeneratorNames = NULL;
3727         }
3728
3729         if(rejection == 0) return; // No Rejection
3730
3731         AliGenCocktailEventHeader *cHeader = 0x0;
3732         AliAODMCHeader *cHeaderAOD = 0x0;
3733         Bool_t headerFound = kFALSE;
3734         AliStack *fMCStack = 0x0;
3735         TClonesArray *fMCStackAOD = 0x0;
3736         if(MCEvent->IsA()==AliMCEvent::Class()){
3737                 cHeader = dynamic_cast<AliGenCocktailEventHeader*>(dynamic_cast<AliMCEvent*>(MCEvent)->GenEventHeader());
3738                 if(cHeader) headerFound = kTRUE;
3739                 fMCStack = dynamic_cast<AliStack*>(dynamic_cast<AliMCEvent*>(MCEvent)->Stack());
3740         }
3741         if(MCEvent->IsA()==AliAODEvent::Class()){ // MCEvent is a AODEvent in case of AOD
3742                 cHeaderAOD = dynamic_cast<AliAODMCHeader*>(MCEvent->FindListObject(AliAODMCHeader::StdBranchName()));
3743                 fMCStackAOD = dynamic_cast<TClonesArray*>(MCEvent->FindListObject(AliAODMCParticle::StdBranchName()));
3744                 
3745                 
3746                 if(cHeaderAOD) headerFound = kTRUE;
3747         }
3748
3749         if(headerFound){
3750                 TList *genHeaders = 0x0;
3751                 if(cHeader) genHeaders = cHeader->GetHeaders();
3752                 if(cHeaderAOD){
3753                         genHeaders = cHeaderAOD->GetCocktailHeaders();
3754                         if(genHeaders->GetEntries()==1){
3755                                 SetRejectExtraSignalsCut(0);
3756                                 return;
3757                         }
3758                 }
3759                 AliGenEventHeader* gh = 0;
3760                 fnHeaders = 0;
3761                 Int_t firstindexA = 0;
3762                 Int_t lastindexA =  -1;
3763                 if(rejection == 1 || rejection == 3) fnHeaders = 1; // MinBiasHeader
3764                 if(rejection == 2){ // TList of Headers Names
3765                         for(Int_t i = 0; i<genHeaders->GetEntries();i++){
3766                                 gh = (AliGenEventHeader*)genHeaders->At(i);
3767                                 TString GeneratorName = gh->GetName();
3768                                 lastindexA = lastindexA + gh->NProduced();
3769 //                              cout << i << "\t" << GeneratorName.Data() << endl;
3770                                 for(Int_t j = 0; j<HeaderList->GetEntries();j++){
3771                                 TString GeneratorInList = ((TObjString*)HeaderList->At(j))->GetString();
3772                                 if(GeneratorName.CompareTo(GeneratorInList) == 0){
3773                                         if (GeneratorInList.CompareTo("PARAM") == 0 || GeneratorInList.CompareTo("BOX") == 0 ){
3774                                                 if(fMCStack){
3775                                                         if (fMCStack->Particle(firstindexA)->GetPdgCode() == 111 || fMCStack->Particle(firstindexA)->GetPdgCode() == 221 ) {
3776                                                                 if (periodName.CompareTo("LHC14a1b")==0 || periodName.CompareTo("LHC14a1b")==0 ){
3777                                                                         if (gh->NProduced() > 10 && (fMCStack->Particle(firstindexA+10)->GetPdgCode() == 111 || fMCStack->Particle(firstindexA+10)->GetPdgCode() == 221 )){
3778                                                                                 fnHeaders++;
3779                                                                                 continue;
3780                                                                         }       
3781                                                                 } else {        
3782                                                                         fnHeaders++;
3783                                                                         continue;
3784                                                                 }       
3785                                                         }
3786                                                 }   
3787                                                 if ( fMCStackAOD){
3788                                                         AliAODMCParticle *aodMCParticle = static_cast<AliAODMCParticle*>(fMCStackAOD->At(firstindexA));
3789                                                         if (  aodMCParticle->GetPdgCode() == 111 || aodMCParticle->GetPdgCode() == 221 ){
3790                                                                 if (periodName.CompareTo("LHC14a1b")==0 || periodName.CompareTo("LHC14a1b")==0 ){
3791                                                                         if (gh->NProduced() > 10 && (fMCStack->Particle(firstindexA+10)->GetPdgCode() == 111 || fMCStack->Particle(firstindexA+10)->GetPdgCode() == 221 )){
3792                                                                                 fnHeaders++;
3793                                                                                 continue;
3794                                                                         }       
3795                                                                 } else {
3796                                                                         fnHeaders++;
3797                                                                         continue;
3798                                                                 }       
3799                                                         }   
3800                                                 }
3801                                         }
3802                                         fnHeaders++;
3803                                         continue;
3804                                 }
3805                                 }
3806                                 firstindexA = firstindexA + gh->NProduced();
3807                         }
3808                 }
3809
3810                 fNotRejectedStart = new Int_t[fnHeaders];
3811                 fNotRejectedEnd = new Int_t[fnHeaders];
3812                 fGeneratorNames = new TString[fnHeaders];
3813
3814                 if(rejection == 1 || rejection == 3){
3815                         fNotRejectedStart[0] = 0;
3816                         fNotRejectedEnd[0] = ((AliGenEventHeader*)genHeaders->At(0))->NProduced()-1;
3817                         fGeneratorNames[0] = ((AliGenEventHeader*)genHeaders->At(0))->GetName();
3818                         return;
3819                 }
3820
3821                 Int_t firstindex = 0;
3822                 Int_t lastindex =  -1;
3823                 Int_t number = 0;
3824                 for(Int_t i = 0; i<genHeaders->GetEntries();i++){
3825                         gh = (AliGenEventHeader*)genHeaders->At(i);
3826                         TString GeneratorName = gh->GetName();
3827                         lastindex = lastindex + gh->NProduced();
3828                         for(Int_t j = 0; j<HeaderList->GetEntries();j++){
3829                                 TString GeneratorInList = ((TObjString*)HeaderList->At(j))->GetString();
3830                                 if(GeneratorName.CompareTo(GeneratorInList) == 0){
3831                                         if (GeneratorInList.CompareTo("PARAM") == 0 || GeneratorInList.CompareTo("BOX") == 0 ){
3832                                                 if(fMCStack){
3833                                                         if (fMCStack->Particle(firstindex)->GetPdgCode() == 111 || fMCStack->Particle(firstindex)->GetPdgCode() == 221 ) {
3834                                                                 if (periodName.CompareTo("LHC14a1b")==0 || periodName.CompareTo("LHC14a1b")==0 ){
3835                                                                         if (gh->NProduced() > 10 && (fMCStack->Particle(firstindexA+10)->GetPdgCode() == 111 || fMCStack->Particle(firstindexA+10)->GetPdgCode() == 221 )){
3836                                                                                 fNotRejectedStart[number] = firstindex;
3837                                                                                 fNotRejectedEnd[number] = lastindex;
3838                                                                                 fGeneratorNames[number] = GeneratorName;
3839                                                                                 number++;
3840                                                                                 continue;
3841                                                                         }       
3842                                                                 } else {
3843                                                                         fNotRejectedStart[number] = firstindex;
3844                                                                         fNotRejectedEnd[number] = lastindex;
3845                                                                         fGeneratorNames[number] = GeneratorName;
3846                                                                         number++;
3847                                                                         continue;
3848                                                                 }       
3849                                                         }
3850                                                 }   
3851                                                 if ( fMCStackAOD){
3852                                                         AliAODMCParticle *aodMCParticle = static_cast<AliAODMCParticle*>(fMCStackAOD->At(firstindex));
3853                                                         if (  aodMCParticle->GetPdgCode() == 111 || aodMCParticle->GetPdgCode() == 221 ){
3854                                                                 if (periodName.CompareTo("LHC14a1b")==0 || periodName.CompareTo("LHC14a1b")==0 ){
3855                                                                         if (gh->NProduced() > 10 && (fMCStack->Particle(firstindexA+10)->GetPdgCode() == 111 || fMCStack->Particle(firstindexA+10)->GetPdgCode() == 221 )){
3856                                                                                 fNotRejectedStart[number] = firstindex;
3857                                                                                 fNotRejectedEnd[number] = lastindex;
3858                                                                                 fGeneratorNames[number] = GeneratorName;
3859                                                                                 number++;
3860                                                                                 continue;
3861                                                                         } 
3862                                                                 } else {
3863                                                                         fNotRejectedStart[number] = firstindex;
3864                                                                         fNotRejectedEnd[number] = lastindex;
3865                                                                         fGeneratorNames[number] = GeneratorName;
3866                                                                         number++;
3867                                                                         continue;       
3868                                                                 }       
3869                                                         }   
3870                                                 }
3871                                                         
3872                                         } else {
3873                                                 fNotRejectedStart[number] = firstindex;
3874                                                 fNotRejectedEnd[number] = lastindex;
3875                                                 fGeneratorNames[number] = GeneratorName;
3876                 //                cout << "Number of particles produced for: " << i << "\t" << GeneratorName.Data() << "\t" << lastindex-firstindex+1 << endl;
3877                                                 number++;
3878                                                 continue;
3879                                         }
3880                                 }
3881                         }
3882                         firstindex = firstindex + gh->NProduced();
3883                 }
3884         } else { // No Cocktail Header Found
3885                 fNotRejectedStart = new Int_t[1];
3886                 fNotRejectedEnd = new Int_t[1];
3887
3888                 fnHeaders = 1;
3889                 fNotRejectedStart[0] = 0;
3890                 fNotRejectedEnd[0] = static_cast<AliMCEvent*>(MCEvent)->Stack()->GetNprimary()-1;
3891                 fGeneratorNames = new TString[1];
3892                 fGeneratorNames[0] = "NoCocktailGeneratorFound";
3893
3894                 AliGenPythiaEventHeader *mcHeaderPythia = dynamic_cast<AliGenPythiaEventHeader*>(static_cast<AliMCEvent*>(MCEvent)->GenEventHeader());
3895                 if (mcHeaderPythia) fGeneratorNames[0] = "NoCocktailGeneratorFound_Pythia";
3896                 AliGenDPMjetEventHeader *mcHeaderPhojet = dynamic_cast<AliGenDPMjetEventHeader*>(static_cast<AliMCEvent*>(MCEvent)->GenEventHeader());
3897                 if (mcHeaderPhojet) fGeneratorNames[0] = "NoCocktailGeneratorFound_Phojet";
3898                 AliGenHijingEventHeader *mcHeaderHijing = dynamic_cast<AliGenHijingEventHeader*>(static_cast<AliMCEvent*>(MCEvent)->GenEventHeader());
3899                 if (mcHeaderHijing) fGeneratorNames[0] = "NoCocktailGeneratorFound_Hijing";
3900
3901                 SetRejectExtraSignalsCut(0);
3902         }
3903 }
3904
3905 //_________________________________________________________________________
3906 Int_t AliConversionCuts::IsParticleFromBGEvent(Int_t index, AliStack *MCStack, AliVEvent *InputEvent){
3907
3908    // Not Accepted == kFALSE == 0
3909    //     Accepted ==  kTRUE == 1
3910    //  FirstHeader ==  kTRUE == 3
3911    if(index < 0) return 0; // No Particle
3912
3913    Int_t accepted = 0;
3914    if(!InputEvent || InputEvent->IsA()==AliESDEvent::Class()){
3915       if( index >= MCStack->GetNprimary()){ // Secondary Particle
3916          if( ((TParticle*)MCStack->Particle(index))->GetMother(0) < 0) return 1; // Secondary Particle without Mother??
3917          return IsParticleFromBGEvent(((TParticle*)MCStack->Particle(index))->GetMother(0),MCStack,InputEvent);
3918       }
3919       for(Int_t i = 0;i<fnHeaders;i++){
3920          if(index >= fNotRejectedStart[i] && index <= fNotRejectedEnd[i]){
3921             accepted = 1;
3922             if(i == 0) accepted = 2; // MB Header
3923          }
3924       }
3925    }
3926    else if(InputEvent->IsA()==AliAODEvent::Class()){
3927       TClonesArray *AODMCTrackArray = dynamic_cast<TClonesArray*>(InputEvent->FindListObject(AliAODMCParticle::StdBranchName()));
3928       AliAODMCParticle *aodMCParticle = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(index));
3929       if(!aodMCParticle) return 1; // Photon Without a Mother ? --> Accepted
3930       if(!aodMCParticle->IsPrimary()){
3931          if( aodMCParticle->GetMother() < 0) return 1;// Secondary Particle without Mother??
3932          return IsParticleFromBGEvent(aodMCParticle->GetMother(),MCStack,InputEvent);
3933       }
3934       index = abs(static_cast<AliAODMCParticle*>(AODMCTrackArray->At(index))->GetLabel());
3935       for(Int_t i = 0;i<fnHeaders;i++){
3936          if(index >= fNotRejectedStart[i] && index <= fNotRejectedEnd[i]){
3937             accepted = 1;
3938             if(i == 0) accepted = 2; // MB Header
3939          }
3940       }
3941    }
3942
3943    return accepted;
3944 }
3945
3946 //_________________________________________________________________________
3947 Int_t AliConversionCuts::IsEventAcceptedByConversionCut(AliConversionCuts *ReaderCuts, AliVEvent *InputEvent, AliMCEvent *MCEvent, Int_t isHeavyIon){
3948
3949    if ( !IsTriggerSelected(InputEvent) )
3950       return 3;
3951
3952    if(isHeavyIon != 0 && !(IsCentralitySelected(InputEvent,MCEvent)))
3953       return 1; // Check Centrality --> Not Accepted => eventQuality = 1
3954       
3955       
3956    if(isHeavyIon == 0 && GetIsFromPileup()){
3957       if(InputEvent->IsPileupFromSPD(3,0.8,3.,2.,5.)){
3958
3959          return 6; // Check Pileup --> Not Accepted => eventQuality = 6
3960       }
3961    }
3962
3963    Bool_t hasV0And = ReaderCuts->HasV0AND();
3964    Bool_t isSDDFired = ReaderCuts->IsSDDFired();
3965    if( (IsSpecialTrigger() == 2 || IsSpecialTrigger() == 3) && !isSDDFired && !MCEvent)
3966       return 7; // With SDD requested but no fired
3967
3968    if( (IsSpecialTrigger() == 1 || IsSpecialTrigger() == 3) && !hasV0And)
3969       return 8; // V0AND requested but no fired
3970
3971    if(hCentrality)hCentrality->Fill(GetCentrality(InputEvent));
3972    if(hCentralityVsNumberOfPrimaryTracks)
3973       hCentralityVsNumberOfPrimaryTracks->Fill(GetCentrality(InputEvent),
3974                                              ((AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()
3975                                                 ->GetTask("V0ReaderV1"))->GetNumberOfPrimaryTracks());     
3976
3977    return 0;
3978 }
3979
3980 //_________________________________________________________________________
3981 Float_t AliConversionCuts::GetWeightForMeson(TString period, Int_t index, AliStack *MCStack, AliVEvent *InputEvent){
3982    if (!(period.CompareTo("LHC12f1a") == 0 || period.CompareTo("LHC12f1b") == 0  || period.CompareTo("LHC12i3") == 0 || period.CompareTo("LHC11a10a") == 0 || period.CompareTo("LHC11a10b") == 0 || period.CompareTo("LHC11a10b_bis") == 0 || period.CompareTo("LHC11a10a_bis") == 0 || period.CompareTo("LHC11a10b_plus") == 0 || period.Contains("LHC13d2")|| period.Contains("LHC14a1") || 
3983    period.CompareTo("LHC13e7") == 0 || period.Contains("LHC13b2_efix"))) return 1.;
3984
3985    Int_t kCaseGen = 0;
3986    for (Int_t i = 0; i < fnHeaders; i++){
3987       if (index >= fNotRejectedStart[i] && index < fNotRejectedEnd[i]+1){
3988          if (fGeneratorNames[i].CompareTo("Pythia") == 0){
3989             kCaseGen = 1;
3990          } else if (fGeneratorNames[i].CompareTo("DPMJET") == 0){
3991             kCaseGen = 2;
3992          } else if (fGeneratorNames[i].CompareTo("HIJING") == 0 ||
3993                     fGeneratorNames[i].CompareTo("Hijing") == 0 ||
3994                     fGeneratorNames[i].Contains("hijing")){
3995             kCaseGen = 3;
3996          } else if (fGeneratorNames[i].CompareTo("BOX") == 0){
3997              kCaseGen = 4;
3998          } else if (fGeneratorNames[i].CompareTo("PARAM") == 0){
3999             kCaseGen = 5;
4000          } else if (fGeneratorNames[i].CompareTo("NoCocktailGeneratorFound") == 0){
4001             kCaseGen = 6;
4002          } else if (fGeneratorNames[i].CompareTo("NoCocktailGeneratorFound_Pythia") == 0){
4003             kCaseGen = 1;
4004          } else if (fGeneratorNames[i].CompareTo("NoCocktailGeneratorFound_Phojet") == 0){
4005             kCaseGen = 2;
4006          } else if (fGeneratorNames[i].CompareTo("NoCocktailGeneratorFound_Hijing") == 0){
4007             kCaseGen = 3;
4008          }
4009          if (period.Contains("LHC13d2") || period.CompareTo("LHC13e7") == 0 || period.Contains("LHC13b2_efix")  || period.Contains("LHC14a1")  ){
4010             kCaseGen = 3;
4011          }
4012       }
4013    }
4014    if (kCaseGen == 0) return 1;
4015
4016
4017    Double_t mesonPt = 0;
4018    Double_t mesonMass = 0;
4019    Int_t PDGCode = 0;
4020    if(!InputEvent || InputEvent->IsA()==AliESDEvent::Class()){
4021       mesonPt = ((TParticle*)MCStack->Particle(index))->Pt();
4022       mesonMass = ((TParticle*)MCStack->Particle(index))->GetCalcMass();
4023       PDGCode = ((TParticle*)MCStack->Particle(index))->GetPdgCode();
4024    }
4025    else if(InputEvent->IsA()==AliAODEvent::Class()){
4026       TClonesArray *AODMCTrackArray = dynamic_cast<TClonesArray*>(InputEvent->FindListObject(AliAODMCParticle::StdBranchName()));
4027       AliAODMCParticle *aodMCParticle = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(index));
4028       mesonPt = aodMCParticle->Pt();
4029       mesonMass = aodMCParticle->GetCalcMass();
4030       PDGCode = aodMCParticle->GetPdgCode();
4031    }
4032
4033    Float_t functionResultMC = 1.;
4034    if (kCaseGen == 1){ // Pythia 6
4035       Float_t dNdyMC = 2.1462;
4036       Float_t nMC = 7.06055;
4037       Float_t tMC = 0.12533;
4038       if ( PDGCode ==  111){
4039          dNdyMC = 2.1462;
4040          nMC = 7.06055;
4041          tMC = 0.12533;
4042       } else if ( PDGCode ==  221){
4043          dNdyMC = 0.2357;
4044          nMC = 5.9105;
4045          tMC = 0.1525;
4046       }
4047       functionResultMC = dNdyMC / ( 2 * TMath::Pi())*(nMC-1.)*(nMC-2.) / (nMC*tMC*(nMC*tMC+mesonMass*(nMC-2.)))  * TMath::Power(1.+(TMath::Sqrt(mesonPt*mesonPt+mesonMass*mesonMass)-mesonMass)/(nMC*tMC), -nMC);
4048    } else if (kCaseGen == 2){ // Phojet
4049       Float_t dNdyMC = 2.35978;
4050       Float_t nMC = 6.81795;
4051       Float_t tMC = 0.11492;
4052       if ( PDGCode ==  111){
4053          dNdyMC = 2.35978;
4054          nMC = 6.81795;
4055          tMC = 0.11492;
4056       } else if ( PDGCode ==  221){
4057          dNdyMC = 0.3690;
4058          nMC = 5.55809;
4059          tMC = 0.13387;
4060       }
4061       functionResultMC = dNdyMC / ( 2 * TMath::Pi())*(nMC-1.)*(nMC-2.) / (nMC*tMC*(nMC*tMC+mesonMass*(nMC-2.)))  * TMath::Power(1.+(TMath::Sqrt(mesonPt*mesonPt+mesonMass*mesonMass)-mesonMass)/(nMC*tMC), -nMC);
4062    } else if (kCaseGen == 4){ // BOX generators pp
4063 //       functionResultMC = 1./sqrt(1.-mesonMass*mesonMass/((mesonMass*mesonMass+mesonPt*mesonPt)*cosh(mesonY)*cosh(mesonY)));
4064       Float_t a = 0.23437;
4065       Float_t b = 5.6661;
4066       Float_t c = -1430.5863;
4067       Float_t d = -0.6966624;
4068       Float_t e = 252.3742;
4069       if ( PDGCode ==  111){
4070          a = 0.23437;
4071          b = 5.6661;
4072          c = -1430.5863;
4073          d = -0.6966624;
4074          e = 252.3742;
4075       } else if ( PDGCode ==  221){
4076          a = 0.10399;
4077          b = 4.35311;
4078          c = -12.17723;
4079          d = -0.01172;
4080          e =1.85140;
4081       }
4082       functionResultMC = a*TMath::Power(mesonPt,-1.*(b+c/(TMath::Power(mesonPt,d)+e)))*1./mesonPt *1./1.6 *1./(2.* TMath::Pi());
4083 //       cout << functionResultMC << endl;
4084    } else if (kCaseGen == 3 ){ // HIJING
4085       if ( PDGCode ==  111 && fDoReweightHistoMCPi0 && hReweightMCHistPi0!= 0x0){
4086          functionResultMC = hReweightMCHistPi0->Interpolate(mesonPt);
4087       }
4088       if ( PDGCode ==  221 && fDoReweightHistoMCEta && hReweightMCHistEta!= 0x0){
4089          functionResultMC = hReweightMCHistEta->Interpolate(mesonPt);
4090       }
4091       if ( PDGCode ==  310 && fDoReweightHistoMCK0s && hReweightMCHistK0s!= 0x0){
4092          functionResultMC = hReweightMCHistK0s->Interpolate(mesonPt);
4093       }
4094    }
4095
4096    Float_t functionResultData = 1;
4097    if (kCaseGen == 1 || kCaseGen == 2 || kCaseGen == 4 ){
4098       Float_t dNdyData = 2.2328;
4099       Float_t nData = 7.1473;
4100       Float_t tData = 0.1346;
4101       if ( PDGCode ==  111){
4102          dNdyData = 2.2328;
4103          nData = 7.1473;
4104          tData = 0.1346;
4105       } else if ( PDGCode ==  221){
4106          dNdyData = 0.38992; //be careful this fit is not optimal, eta in data still has problems
4107          nData = 5.72778;
4108          tData = 0.13835;
4109       }
4110       functionResultData = dNdyData / ( 2 * TMath::Pi())*(nData-1.)*(nData-2.) / (nData*tData*(nData*tData+mesonMass*(nData-2.)))  * TMath::Power(1.+(TMath::Sqrt(mesonPt*mesonPt+mesonMass*mesonMass)-mesonMass)/(nData*tData), -nData);
4111 //       cout << functionResultData << endl;
4112    } else {
4113       if ( PDGCode ==  111 && fDoReweightHistoMCPi0 && fFitDataPi0!= 0x0){
4114          functionResultData = fFitDataPi0->Eval(mesonPt);
4115       }
4116       if ( PDGCode ==  221 && fDoReweightHistoMCEta && fFitDataEta!= 0x0){
4117          functionResultData = fFitDataEta->Eval(mesonPt);
4118       }
4119       if ( PDGCode ==  310 && fDoReweightHistoMCK0s && fFitDataK0s!= 0x0){
4120          functionResultData = fFitDataK0s->Eval(mesonPt);
4121       }
4122
4123    }
4124
4125    Double_t weight = 1;
4126    if (PDGCode ==  111 || PDGCode ==  221){
4127       if (functionResultData != 0. && functionResultMC != 0. && isfinite(functionResultData) && isfinite(functionResultMC)){
4128          weight = functionResultData/functionResultMC;
4129          if ( kCaseGen == 3){
4130             if (PDGCode ==  111){ 
4131                if (!(fDoReweightHistoMCPi0 && hReweightMCHistPi0!= 0x0 && PDGCode ==  111)){
4132                   weight = 1.;
4133                }
4134             } 
4135             if (PDGCode ==  221){ 
4136                if (!(fDoReweightHistoMCEta && hReweightMCHistEta!= 0x0 && PDGCode ==  221)){
4137                   weight = 1.;
4138                }
4139             }
4140          }
4141          if (!isfinite(functionResultData)) weight = 1.;
4142          if (!isfinite(weight)) weight = 1.;
4143       }
4144    } else if (PDGCode ==  310 && functionResultMC != 0 && isfinite(functionResultMC)){
4145         weight = functionResultMC;
4146    }
4147
4148 //    if (fModCentralityClass == 0 && fCentralityMin == 4 && fCentralityMax == 6 && PDGCode ==  111){
4149 //        cout << period.Data() << "\t" << kCaseGen << "\t" <<fModCentralityClass<< "\t" <<fCentralityMin<< "\t" <<fCentralityMax << "\t" << mesonPt << "\t" <<mesonMass<< "\t"<<functionResultData << "\t"<< functionResultMC << "\t" << weight <<endl;
4150 //    }
4151    return weight;
4152 }
4153 ///________________________________________________________________________
4154 AliConversionCuts* AliConversionCuts::GetStandardCuts2010PbPb(){
4155     //Create and return standard 2010 PbPb cuts
4156     AliConversionCuts *cuts=new AliConversionCuts("StandardCuts2010PbPb","StandardCuts2010PbPb");
4157     if(!cuts->InitializeCutsFromCutString("100000204209297002322000000")){
4158    cout<<"Warning: Initialization of Standardcuts2010PbPb failed"<<endl;}
4159     return cuts;
4160 }
4161
4162 ///________________________________________________________________________
4163 AliConversionCuts* AliConversionCuts::GetStandardCuts2010pp(){
4164     //Create and return standard 2010 PbPb cuts
4165     AliConversionCuts *cuts=new AliConversionCuts("StandardCuts2010pp","StandardCuts2010pp");
4166     if(!cuts->InitializeCutsFromCutString("000001100209366300380000000")){
4167    cout<<"Warning: Initialization of Standardcuts2010pp failed"<<endl;}
4168     return cuts;
4169 }
4170 ///________________________________________________________________________
4171 void AliConversionCuts::GetCorrectEtaShiftFromPeriod(TString periodName){
4172
4173    if(periodName.CompareTo("LHC12g") == 0 || //pilot run 2012
4174       periodName.CompareTo("LHC13b") == 0 || //mainly minimum bias
4175       periodName.CompareTo("LHC13c") == 0 || //mainly minimum bias
4176       periodName.CompareTo("LHC13d") == 0 || //mainly triggered
4177       periodName.CompareTo("LHC13e") == 0 || //mainly triggered
4178       periodName.CompareTo("LHC13c3") == 0 || //MC Starlight, anchor LHC13d+e
4179       periodName.CompareTo("LHC13c2") == 0 || //MC Starlight, coherent J/Psi, UPC muon anchor LHC13d+e
4180       periodName.CompareTo("LHC13b4") == 0 || //MC Pythia 6 (Jet-Jet), anchor LHC13b
4181       periodName.CompareTo("LHC13b2_fix_1") == 0 || //MC DPMJET, anchr LHC13b+c
4182       periodName.CompareTo("LHC13b2_efix_p1") == 0 || //MC DPMJET, anchr LHC13b+c
4183       periodName.CompareTo("LHC13b2_efix_p2") == 0 || //MC DPMJET, anchr LHC13b+c
4184       periodName.CompareTo("LHC13b2_efix_p3") == 0 || //MC DPMJET, anchr LHC13b+c
4185       periodName.CompareTo("LHC13b2_efix_p4") == 0 || //MC DPMJET, anchr LHC13b+c
4186       periodName.CompareTo("LHC13e7") == 0 || //MC DPMJET, anchr LHC13b+c
4187       periodName.CompareTo("LHC13b3") == 0 || //MC HIJING, weighted to number of events per run, anchor LHC13b
4188       periodName.CompareTo("LHC13b2") == 0 ||  // MC DPMJET, wrong energy, anchor LHC13b
4189       periodName.CompareTo("LHC13b2_plus") == 0 || // MC DPMJET, weighted to number event per run, anchor LHC13b
4190       periodName.CompareTo("LHC13c1_bis") == 0 || // MC AMPT fast generation, pT hardbin, anchor ?
4191       periodName.CompareTo("LHC13c1") == 0 || // MC AMPT fast generation, anchor ?
4192       periodName.CompareTo("LHC13b1") == 0 || // MC DPMJET, fragments, with fixed label 0, anchor LHC12g
4193       periodName.CompareTo("LHC12g4b_fix") == 0 || // MC DPMJET, with fixed label 0, anchor LHC12g
4194       periodName.CompareTo("LHC12g1_fix") == 0 || // MC ?, with fixed label 0, anchor LHC12g
4195       periodName.CompareTo("LHC12g4c") == 0 || // MC DPMJET, shifted vertex runs, anchor LHC12g
4196       periodName.CompareTo("LHC12h6") == 0 || // MC muon cocktail, anchor LHC12g
4197       periodName.CompareTo("LHC12g4b") == 0 || // MC DPMJET 3rd iteration, anchor LHC12g
4198       periodName.CompareTo("LHC12g4a") == 0 || // MC DPMJET improved, anchor LHC12g
4199       periodName.CompareTo("LHC12g4") == 0 || // MC DPMJET, anchor LHC12g
4200       periodName.CompareTo("LHC12g5") == 0 || // MC PHOJET, anchor LHC12g
4201       periodName.CompareTo("LHC12g2") == 0 || // MC Starlight background, anchor LHC12g
4202       periodName.CompareTo("LHC12g1") == 0 ) // MC ?, anchor LHC12g
4203       {
4204          printf(" Gamma Conversion Cuts %s :: pPb Run doing Eta Shift of %f \n\n",(GetCutNumber()).Data(),-0.465);
4205          SetEtaShift(-0.465);
4206       }
4207    else if(periodName.CompareTo("LHC13f") == 0 ||
4208            periodName.CompareTo("LHC13c6b") == 0 ||// MC Jpsi -> mumu, anchor LHC13f
4209            periodName.CompareTo("LHC13c5") == 0 || //MC Starlight, gamma gamma UPC muon, anchor LHC13f
4210            periodName.CompareTo("LHC13c4") == 0 )//MC Starlight, coherent JPsi, UPC muon, anchor LHC13f
4211       {
4212          printf(" Gamma Conversion Cuts %s :: Pbp Run doing Eta Shift of %f \n\n",(GetCutNumber()).Data(),0.465);
4213          SetEtaShift(+0.465);
4214       }
4215    else printf(" Gamma Conversion Cuts %s :: Automatic Eta Shift requested but Period is not known -> No Shift \n\n",(GetCutNumber()).Data());
4216 }
4217
4218 ///________________________________________________________________________
4219 Bool_t AliConversionCuts::InPlaneOutOfPlaneCut(Double_t photonPhi, Double_t eventPlaneAngle, Bool_t fill){
4220    
4221    //GetPhotonPhi() 0-2 Pi  //eventPlaneAngle -1pi-1pi
4222    eventPlaneAngle=eventPlaneAngle+TMath::Pi();
4223    Double_t gammaToEPAngle = eventPlaneAngle-photonPhi;
4224    if(gammaToEPAngle < 0) gammaToEPAngle=gammaToEPAngle+2*TMath::Pi();
4225    gammaToEPAngle = gammaToEPAngle-TMath::Pi(); // angle from -pi +pi
4226
4227    if(!fInPlaneOutOfPlane){
4228       if(fill&&hEventPlanePhi)hEventPlanePhi->Fill(gammaToEPAngle);
4229       return kTRUE;
4230    }
4231    else if(fInPlaneOutOfPlane == 1){
4232       if(abs(gammaToEPAngle)<=0.25*TMath::Pi() || abs(gammaToEPAngle)>=0.75*TMath::Pi()){
4233          if(fill&&hEventPlanePhi)hEventPlanePhi->Fill(gammaToEPAngle);
4234          return kTRUE;
4235       }
4236       else return kFALSE;
4237    }
4238    else if(fInPlaneOutOfPlane == 2){
4239       if(abs(gammaToEPAngle)>0.25*TMath::Pi() && abs(gammaToEPAngle)<0.75*TMath::Pi()){
4240          if(fill&&hEventPlanePhi)hEventPlanePhi->Fill(gammaToEPAngle);
4241          return kTRUE;
4242       }
4243       else return kFALSE;
4244    }
4245    return kFALSE;
4246
4247 }
4248