Replace GetMass for M(), smearing for virtual photons added, histogramas for eta...
[u/mrichter/AliRoot.git] / PWGGA / GammaConv / AliConversionMesonCuts.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 "AliConversionMesonCuts.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 "AliStack.h"
36 #include "AliAODConversionMother.h"
37 #include "TObjString.h"
38 #include "AliAODEvent.h"
39 #include "AliESDEvent.h"
40 #include "AliCentrality.h"
41 #include "TList.h"
42 #include "TPDGCode.h"
43 #include "TDatabasePDG.h"
44 #include "AliAODMCParticle.h"
45
46 class iostream;
47
48 using namespace std;
49
50 ClassImp(AliConversionMesonCuts)
51
52
53 const char* AliConversionMesonCuts::fgkCutNames[AliConversionMesonCuts::kNCuts] = {
54    "MesonKind", //0
55    "BackgroundScheme", //1
56    "NumberOfBGEvents", //2
57    "DegreesForRotationMethod", //3
58    "RapidityMesonCut", //4
59    "RCut", //5
60    "AlphaMesonCut", //6
61    "Chi2MesonCut", //7
62    "SharedElectronCuts", //8
63    "RejectToCloseV0s", //9
64    "UseMCPSmearing", //10
65    "DcaGammaGamma", //11
66    "DcaRPrimVtx", //12
67    "DcaZPrimVtx" //13
68 };
69
70
71 //________________________________________________________________________
72 AliConversionMesonCuts::AliConversionMesonCuts(const char *name,const char *title) :
73    AliAnalysisCuts(name,title),
74    fHistograms(NULL),
75    fMesonKind(0),
76    fMaxR(200),
77    fChi2CutMeson(1000),
78    fAlphaMinCutMeson(0),
79    fAlphaCutMeson(1),
80    fRapidityCutMeson(1),
81    fUseRotationMethodInBG(kFALSE),
82    fDoBG(kTRUE),
83    fdoBGProbability(kFALSE),
84    fUseTrackMultiplicityForBG(kFALSE),
85    fnDegreeRotationPMForBG(0),
86    fNumberOfBGEvents(0),
87    fOpeningAngle(0.005),
88    fDoToCloseV0sCut(kFALSE),
89    fminV0Dist(200.),
90    fDoSharedElecCut(kFALSE),
91    fUseMCPSmearing(kFALSE),
92    fPBremSmearing(0),
93    fPSigSmearing(0),
94    fPSigSmearingCte(0),
95    fBrem(NULL),
96    fRandom(0),
97    fFAlphaCut(0),
98    fAlphaPtDepCut(kFALSE),
99    fElectronLabelArraySize(500),
100    fElectronLabelArray(NULL),
101    fDCAGammaGammaCut(1000),
102    fDCAZMesonPrimVtxCut(1000),
103    fDCARMesonPrimVtxCut(1000),
104    fBackgroundHandler(0),
105    fCutString(NULL),
106    hMesonCuts(NULL),
107    hMesonBGCuts(NULL),
108    hDCAGammaGammaMesonBefore(NULL),
109    hDCAZMesonPrimVtxBefore(NULL),
110    hDCARMesonPrimVtxBefore(NULL),
111    hDCAGammaGammaMesonAfter(NULL),
112    hDCAZMesonPrimVtxAfter(NULL),
113    hDCARMesonPrimVtxAfter(NULL)
114
115 {
116    for(Int_t jj=0;jj<kNCuts;jj++){fCuts[jj]=0;}
117    fCutString=new TObjString((GetCutNumber()).Data());
118    fElectronLabelArray = new Int_t[fElectronLabelArraySize];
119    if (fBrem == NULL){
120       fBrem = new TF1("fBrem","pow(-log(x),[0]/log(2.0)-1.0)/TMath::Gamma([0]/log(2.0))",0.00001,0.999999999);
121       // tests done with 1.0e-14
122       fBrem->SetParameter(0,fPBremSmearing);
123       fBrem->SetNpx(100000);
124    }
125
126 }
127
128 //________________________________________________________________________
129 AliConversionMesonCuts::AliConversionMesonCuts(const AliConversionMesonCuts &ref) :
130    AliAnalysisCuts(ref),
131    fHistograms(NULL),
132    fMesonKind(ref.fMesonKind),
133    fMaxR(ref.fMaxR),
134    fChi2CutMeson(ref.fChi2CutMeson),
135    fAlphaMinCutMeson(ref.fAlphaMinCutMeson),
136    fAlphaCutMeson(ref.fAlphaCutMeson),
137    fRapidityCutMeson(ref.fRapidityCutMeson),
138    fUseRotationMethodInBG(ref.fUseRotationMethodInBG),
139    fDoBG(ref.fDoBG),
140    fdoBGProbability(ref.fdoBGProbability),
141    fUseTrackMultiplicityForBG(ref.fUseTrackMultiplicityForBG),
142    fnDegreeRotationPMForBG(ref.fnDegreeRotationPMForBG),
143    fNumberOfBGEvents(ref. fNumberOfBGEvents),
144    fOpeningAngle(ref.fOpeningAngle),
145    fDoToCloseV0sCut(ref.fDoToCloseV0sCut),
146    fminV0Dist(ref.fminV0Dist),
147    fDoSharedElecCut(ref.fDoSharedElecCut),
148    fUseMCPSmearing(ref.fUseMCPSmearing),
149    fPBremSmearing(ref.fPBremSmearing),
150    fPSigSmearing(ref.fPSigSmearing),
151    fPSigSmearingCte(ref.fPSigSmearingCte),
152    fBrem(NULL),
153    fRandom(ref.fRandom),
154    fFAlphaCut(NULL),
155    fAlphaPtDepCut(ref.fAlphaPtDepCut),
156    fElectronLabelArraySize(ref.fElectronLabelArraySize),
157    fElectronLabelArray(NULL),
158    fDCAGammaGammaCut(ref.fDCAGammaGammaCut),
159    fDCAZMesonPrimVtxCut(ref.fDCAZMesonPrimVtxCut),
160    fDCARMesonPrimVtxCut(ref.fDCARMesonPrimVtxCut),
161    fBackgroundHandler(ref.fBackgroundHandler),
162    fCutString(NULL),
163    hMesonCuts(NULL),
164    hMesonBGCuts(NULL),
165    hDCAGammaGammaMesonBefore(NULL),
166    hDCAZMesonPrimVtxBefore(NULL),
167    hDCARMesonPrimVtxBefore(NULL),
168    hDCAGammaGammaMesonAfter(NULL),
169    hDCAZMesonPrimVtxAfter(NULL),
170    hDCARMesonPrimVtxAfter(NULL)
171 {
172    // Copy Constructor
173    for(Int_t jj=0;jj<kNCuts;jj++){fCuts[jj]=ref.fCuts[jj];}
174    fCutString=new TObjString((GetCutNumber()).Data());
175    fElectronLabelArray = new Int_t[fElectronLabelArraySize];
176    if (fBrem == NULL)fBrem = (TF1*)ref.fBrem->Clone("fBrem");
177    // Histograms are not copied, if you need them, call InitCutHistograms
178 }
179
180
181 //________________________________________________________________________
182 AliConversionMesonCuts::~AliConversionMesonCuts() {
183    // Destructor
184    //Deleting fHistograms leads to seg fault it it's added to output collection of a task
185    // if(fHistograms)
186    //   delete fHistograms;
187    // fHistograms = NULL;
188    if(fCutString != NULL){
189       delete fCutString;
190       fCutString = NULL;
191    }
192    if(fElectronLabelArray){
193       delete fElectronLabelArray;
194       fElectronLabelArray = NULL;
195    }
196
197 }
198
199 //________________________________________________________________________
200 void AliConversionMesonCuts::InitCutHistograms(TString name, Bool_t additionalHists){
201
202    // Initialize Cut Histograms for QA (only initialized and filled if function is called)
203    TH1::AddDirectory(kFALSE);
204
205    if(fHistograms != NULL){
206       delete fHistograms;
207       fHistograms=NULL;
208    }
209    if(fHistograms==NULL){
210       fHistograms=new TList();
211       fHistograms->SetOwner(kTRUE);
212       if(name=="")fHistograms->SetName(Form("ConvMesonCuts_%s",GetCutNumber().Data()));
213       else fHistograms->SetName(Form("%s_%s",name.Data(),GetCutNumber().Data()));
214    }
215
216    // Meson Cuts
217    hMesonCuts=new TH1F(Form("MesonCuts %s",GetCutNumber().Data()),"MesonCuts",13,-0.5,12.5);
218    hMesonCuts->GetXaxis()->SetBinLabel(1,"in");
219    hMesonCuts->GetXaxis()->SetBinLabel(2,"undef rapidity");
220    hMesonCuts->GetXaxis()->SetBinLabel(3,"rapidity cut");
221    hMesonCuts->GetXaxis()->SetBinLabel(4,"opening angle");
222    hMesonCuts->GetXaxis()->SetBinLabel(5,"alpha max");
223    hMesonCuts->GetXaxis()->SetBinLabel(6,"alpha min");
224    hMesonCuts->GetXaxis()->SetBinLabel(7,"dca gamma gamma");
225    hMesonCuts->GetXaxis()->SetBinLabel(8,"dca R prim Vtx");
226    hMesonCuts->GetXaxis()->SetBinLabel(9,"dca Z prim Vtx");
227    hMesonCuts->GetXaxis()->SetBinLabel(10,"out");
228    fHistograms->Add(hMesonCuts);
229
230    hMesonBGCuts=new TH1F(Form("MesonBGCuts %s",GetCutNumber().Data()),"MesonBGCuts",13,-0.5,12.5);
231    hMesonBGCuts->GetXaxis()->SetBinLabel(1,"in");
232    hMesonBGCuts->GetXaxis()->SetBinLabel(2,"undef rapidity");
233    hMesonBGCuts->GetXaxis()->SetBinLabel(3,"rapidity cut");
234    hMesonBGCuts->GetXaxis()->SetBinLabel(4,"opening angle");
235    hMesonBGCuts->GetXaxis()->SetBinLabel(5,"alpha max");
236    hMesonBGCuts->GetXaxis()->SetBinLabel(6,"alpha min");
237    hMesonBGCuts->GetXaxis()->SetBinLabel(7,"dca gamma gamma");
238    hMesonBGCuts->GetXaxis()->SetBinLabel(8,"dca R prim Vtx");
239    hMesonBGCuts->GetXaxis()->SetBinLabel(9,"dca Z prim Vtx");
240    hMesonBGCuts->GetXaxis()->SetBinLabel(10,"out");
241
242    fHistograms->Add(hMesonBGCuts);
243
244    if (additionalHists){
245       hDCAGammaGammaMesonBefore=new TH1F(Form("DCAGammaGammaMeson Before %s",GetCutNumber().Data()),"DCAGammaGammaMeson Before",200,0,10);
246       fHistograms->Add(hDCAGammaGammaMesonBefore);
247
248       hDCARMesonPrimVtxBefore=new TH1F(Form("DCARMesonPrimVtx Before %s",GetCutNumber().Data()),"DCARMesonPrimVtx Before",200,0,10);
249       fHistograms->Add(hDCARMesonPrimVtxBefore);
250
251       hDCAZMesonPrimVtxBefore=new TH1F(Form("DCAZMesonPrimVtx Before %s",GetCutNumber().Data()),"DCAZMesonPrimVtx Before",401,-10,10);
252       fHistograms->Add(hDCAZMesonPrimVtxBefore);
253
254    }
255
256    hDCAGammaGammaMesonAfter=new TH1F(Form("DCAGammaGammaMeson After %s",GetCutNumber().Data()),"DCAGammaGammaMeson After",200,0,10);
257    fHistograms->Add(hDCAGammaGammaMesonAfter);
258
259    hDCAZMesonPrimVtxAfter=new TH2F(Form("InvMassDCAZMesonPrimVtx After %s",GetCutNumber().Data()),"InvMassDCAZMesonPrimVtx After",800,0,0.8,401,-10,10);
260    fHistograms->Add(hDCAZMesonPrimVtxAfter);
261
262    hDCARMesonPrimVtxAfter=new TH1F(Form("DCARMesonPrimVtx After %s",GetCutNumber().Data()),"DCARMesonPrimVtx After",200,0,10);
263    fHistograms->Add(hDCARMesonPrimVtxAfter);
264
265    TH1::AddDirectory(kTRUE);
266 }
267
268 //________________________________________________________________________
269 Bool_t AliConversionMesonCuts::MesonIsSelectedMC(TParticle *fMCMother,AliStack *fMCStack, Double_t fRapidityShift){
270    // Returns true for all pions within acceptance cuts for decay into 2 photons
271    // If bMCDaughtersInAcceptance is selected, it requires in addition that both daughter photons are within acceptance cuts
272
273    if(!fMCStack)return kFALSE;
274
275    if(fMCMother->GetPdgCode()==111 || fMCMother->GetPdgCode()==221){
276       if(fMCMother->R()>fMaxR)  return kFALSE; // cuts on distance from collision point
277
278       Double_t rapidity = 10.;
279       if(fMCMother->Energy() - fMCMother->Pz() == 0 || fMCMother->Energy() + fMCMother->Pz() == 0){
280          rapidity=8.-fRapidityShift;
281       } else{
282          rapidity = 0.5*(TMath::Log((fMCMother->Energy()+fMCMother->Pz()) / (fMCMother->Energy()-fMCMother->Pz())))-fRapidityShift;
283       }
284
285       // Rapidity Cut
286       if(abs(rapidity)>fRapidityCutMeson)return kFALSE;
287
288       // Select only -> 2y decay channel
289       if(fMCMother->GetNDaughters()!=2)return kFALSE;
290
291       for(Int_t i=0;i<2;i++){
292          TParticle *MDaughter=fMCStack->Particle(fMCMother->GetDaughter(i));
293          // Is Daughter a Photon?
294          if(MDaughter->GetPdgCode()!=22)return kFALSE;
295          // Is Photon in Acceptance?
296          //   if(bMCDaughtersInAcceptance){
297          //     if(!PhotonIsSelectedMC(MDaughter,fMCStack)){return kFALSE;}
298          //   }
299       }
300       return kTRUE;
301    }
302    return kFALSE;
303 }
304 //________________________________________________________________________
305 Bool_t AliConversionMesonCuts::MesonIsSelectedAODMC(AliAODMCParticle *MCMother,TClonesArray *AODMCArray, Double_t fRapidityShift){
306    // Returns true for all pions within acceptance cuts for decay into 2 photons
307    // If bMCDaughtersInAcceptance is selected, it requires in addition that both daughter photons are within acceptance cuts
308
309    if(!AODMCArray)return kFALSE;
310
311    if(MCMother->GetPdgCode()==111 || MCMother->GetPdgCode()==221){
312       Double_t rMeson = sqrt( (MCMother->Xv()*MCMother->Xv()) + (MCMother->Yv()*MCMother->Yv()) ) ;
313       if(rMeson>fMaxR)  return kFALSE; // cuts on distance from collision point
314
315       Double_t rapidity = 10.;
316       if(MCMother->E() - MCMother->Pz() == 0 || MCMother->E() + MCMother->Pz() == 0){
317          rapidity=8.-fRapidityShift;
318       } else{
319          rapidity = 0.5*(TMath::Log((MCMother->E()+MCMother->Pz()) / (MCMother->E()-MCMother->Pz())))-fRapidityShift;
320       }
321
322       // Rapidity Cut
323       if(abs(rapidity)>fRapidityCutMeson)return kFALSE;
324
325       // Select only -> 2y decay channel
326       if(MCMother->GetNDaughters()!=2)return kFALSE;
327
328       for(Int_t i=0;i<2;i++){
329          AliAODMCParticle *MDaughter=static_cast<AliAODMCParticle*>(AODMCArray->At(MCMother->GetDaughter(i)));
330          // Is Daughter a Photon?
331          if(MDaughter->GetPdgCode()!=22)return kFALSE;
332          // Is Photon in Acceptance?
333          //   if(bMCDaughtersInAcceptance){
334          //     if(!PhotonIsSelectedMC(MDaughter,fMCStack)){return kFALSE;}
335          //   }
336       }
337       return kTRUE;
338    }
339    return kFALSE;
340 }
341
342 //________________________________________________________________________
343 Bool_t AliConversionMesonCuts::MesonIsSelectedMCDalitz(TParticle *fMCMother,AliStack *fMCStack, Int_t &labelelectron, Int_t &labelpositron, Int_t &labelgamma, Double_t fRapidityShift){
344
345    // Returns true for all pions within acceptance cuts for decay into 2 photons
346    // If bMCDaughtersInAcceptance is selected, it requires in addition that both daughter photons are within acceptance cuts
347
348    if( !fMCStack )return kFALSE;
349
350    if(  fMCMother->GetPdgCode() != 111 && fMCMother->GetPdgCode() != 221 ) return kFALSE;
351
352    if(  fMCMother->R()>fMaxR ) return kFALSE; // cuts on distance from collision point
353
354    Double_t rapidity = 10.;
355
356    if( fMCMother->Energy() - fMCMother->Pz() == 0 || fMCMother->Energy() + fMCMother->Pz() == 0 ){
357       rapidity=8.-fRapidityShift;
358    }
359    else{
360       rapidity = 0.5*(TMath::Log((fMCMother->Energy()+fMCMother->Pz()) / (fMCMother->Energy()-fMCMother->Pz())))-fRapidityShift;
361    }
362
363    // Rapidity Cut
364    if( abs(rapidity) > fRapidityCutMeson )return kFALSE;
365
366    // Select only -> Dalitz decay channel
367    if( fMCMother->GetNDaughters() != 3 )return kFALSE;
368
369    TParticle *positron = 0x0;
370    TParticle *electron = 0x0;
371    TParticle    *gamma = 0x0;
372
373    for(Int_t index= fMCMother->GetFirstDaughter();index<= fMCMother->GetLastDaughter();index++){
374
375       TParticle* temp = (TParticle*)fMCStack->Particle( index );
376
377       switch( temp->GetPdgCode() ) {
378       case ::kPositron:
379          positron      =  temp;
380          labelpositron = index;
381          break;
382       case ::kElectron:
383          electron      =  temp;
384          labelelectron = index;
385          break;
386       case ::kGamma:
387          gamma         =  temp;
388          labelgamma    = index;
389          break;
390       }
391    }
392
393    if( positron && electron && gamma) return kTRUE;
394    return kFALSE;
395 }
396
397 //________________________________________________________________________
398 Bool_t AliConversionMesonCuts::MesonIsSelectedMCEtaPiPlPiMiGamma(TParticle *fMCMother,AliStack *fMCStack, Int_t &labelNegPion, Int_t &labelPosPion, Int_t &labelGamma, Double_t fRapidityShift){
399
400         // Returns true for all pions within acceptance cuts for decay into 2 photons
401         // If bMCDaughtersInAcceptance is selected, it requires in addition that both daughter photons are within acceptance cuts
402
403         if( !fMCStack )return kFALSE;
404
405         if( fMCMother->GetPdgCode() != 221 ) return kFALSE;
406
407         if( fMCMother->R()>fMaxR ) return kFALSE; // cuts on distance from collision point
408
409         Double_t rapidity = 10.;
410
411         if( fMCMother->Energy() - fMCMother->Pz() == 0 || fMCMother->Energy() + fMCMother->Pz() == 0 ){
412                 rapidity=8.-fRapidityShift;
413         }
414         else{
415                 rapidity = 0.5*(TMath::Log((fMCMother->Energy()+fMCMother->Pz()) / (fMCMother->Energy()-fMCMother->Pz())))-fRapidityShift;
416         }
417
418         // Rapidity Cut
419         if( abs(rapidity) > fRapidityCutMeson )return kFALSE;
420
421         // Select only -> Dalitz decay channel
422         if( fMCMother->GetNDaughters() != 3 )return kFALSE;
423
424         TParticle *posPion = 0x0;
425         TParticle *negPion = 0x0;
426         TParticle    *gamma = 0x0;
427
428         for(Int_t index= fMCMother->GetFirstDaughter();index<= fMCMother->GetLastDaughter();index++){
429
430                 TParticle* temp = (TParticle*)fMCStack->Particle( index );
431
432                 switch( temp->GetPdgCode() ) {
433                 case 211:
434                         posPion      =  temp;
435                         labelPosPion = index;
436                         break;
437                 case -211:
438                         negPion      =  temp;
439                         labelNegPion = index;
440                         break;
441                 case ::kGamma:
442                         gamma         =  temp;
443                         labelGamma    = index;
444                         break;
445                 }
446         }
447
448         if( posPion && negPion && gamma) return kTRUE;
449         return kFALSE;
450 }
451
452 //________________________________________________________________________
453 Bool_t AliConversionMesonCuts::MesonIsSelectedMCChiC(TParticle *fMCMother,AliStack *fMCStack,Int_t & labelelectronChiC, Int_t & labelpositronChiC, Int_t & labelgammaChiC, Double_t fRapidityShift){
454    // Returns true for all ChiC within acceptance cuts for decay into JPsi + gamma -> e+ + e- + gamma
455    // If bMCDaughtersInAcceptance is selected, it requires in addition that both daughter photons are within acceptance cuts
456
457    if(!fMCStack)return kFALSE;
458    // if(fMCMother->GetPdgCode()==20443 ){
459    //    return kFALSE;
460    // }
461    if(fMCMother->GetPdgCode()==10441 || fMCMother->GetPdgCode()==10443 || fMCMother->GetPdgCode()==445 ){
462       if(fMCMother->R()>fMaxR)  return kFALSE; // cuts on distance from collision point
463
464       Double_t rapidity = 10.;
465       if(fMCMother->Energy() - fMCMother->Pz() == 0 || fMCMother->Energy() + fMCMother->Pz() == 0){
466          rapidity=8.-fRapidityShift;
467       }
468       else{
469          rapidity = 0.5*(TMath::Log((fMCMother->Energy()+fMCMother->Pz()) / (fMCMother->Energy()-fMCMother->Pz())))-fRapidityShift;
470       }
471
472       // Rapidity Cut
473       if(abs(rapidity)>fRapidityCutMeson)return kFALSE;
474
475       // Select only -> ChiC radiative (JPsi+gamma) decay channel
476       if(fMCMother->GetNDaughters()!=2)return kFALSE;
477
478       TParticle *jpsi   = 0x0;
479       TParticle *gamma  = 0x0;
480       TParticle *positron = 0x0;
481       TParticle *electron = 0x0;
482
483       Int_t labeljpsiChiC = -1;
484
485       for(Int_t index= fMCMother->GetFirstDaughter();index<= fMCMother->GetLastDaughter();index++){
486
487          TParticle* temp = (TParticle*)fMCStack->Particle( index );
488
489          switch( temp->GetPdgCode() ) {
490          case 443:
491             jpsi =  temp;
492             labeljpsiChiC = index;
493             break;
494          case 22:
495             gamma    =  temp;
496             labelgammaChiC = index;
497             break;
498          }
499       }
500
501       if ( !jpsi || ! gamma) return kFALSE;
502       if(jpsi->GetNDaughters()!=2)return kFALSE;
503
504
505       for(Int_t index= jpsi->GetFirstDaughter();index<= jpsi->GetLastDaughter();index++){
506          TParticle* temp = (TParticle*)fMCStack->Particle( index );
507          switch( temp->GetPdgCode() ) {
508          case -11:
509             electron =  temp;
510             labelelectronChiC = index;
511             break;
512          case 11:
513             positron =  temp;
514             labelpositronChiC = index;
515             break;
516          }
517       }
518       if( !electron || !positron) return kFALSE;
519       if( positron && electron && gamma) return kTRUE;
520    }
521    return kFALSE;
522 }
523
524
525
526 ///________________________________________________________________________
527 Bool_t AliConversionMesonCuts::MesonIsSelected(AliAODConversionMother *pi0,Bool_t IsSignal, Double_t fRapidityShift)
528 {
529
530    // Selection of reconstructed Meson candidates
531    // Use flag IsSignal in order to fill Fill different
532    // histograms for Signal and Background
533    TH1 *hist=0x0;
534
535    if(IsSignal){hist=hMesonCuts;}
536    else{hist=hMesonBGCuts;}
537
538    Int_t cutIndex=0;
539
540    if(hist)hist->Fill(cutIndex);
541    cutIndex++;
542
543    // Undefined Rapidity -> Floating Point exception
544    if((pi0->E()+pi0->Pz())/(pi0->E()-pi0->Pz())<=0){
545       if(hist)hist->Fill(cutIndex);
546       cutIndex++;
547 //        cout << "undefined rapidity" << endl;
548       return kFALSE;
549    }
550    else{
551       // PseudoRapidity Cut --> But we cut on Rapidity !!!
552       cutIndex++;
553       if(abs(pi0->Rapidity()-fRapidityShift)>fRapidityCutMeson){
554          if(hist)hist->Fill(cutIndex);
555 //               cout << abs(pi0->Rapidity()-fRapidityShift) << ">" << fRapidityCutMeson << endl;
556          return kFALSE;
557       }
558    }
559    cutIndex++;
560
561    // Opening Angle Cut
562    //fOpeningAngle=2*TMath::ATan(0.134/pi0->P());// physical minimum opening angle
563    if( pi0->GetOpeningAngle() < fOpeningAngle){
564 //        cout << pi0->GetOpeningAngle() << "<" << fOpeningAngle << endl; 
565       if(hist)hist->Fill(cutIndex);
566       return kFALSE;
567    }
568    cutIndex++;
569
570    if ( fAlphaPtDepCut == kTRUE ) {
571  
572         fAlphaCutMeson = fFAlphaCut->Eval( pi0->Pt() );
573    }
574    
575    
576    // Alpha Max Cut
577    if(pi0->GetAlpha()>fAlphaCutMeson){
578 //         cout << pi0->GetAlpha() << ">" << fAlphaCutMeson << endl; 
579       if(hist)hist->Fill(cutIndex);
580       return kFALSE;
581    }
582    cutIndex++;
583
584    // Alpha Min Cut
585    if(pi0->GetAlpha()<fAlphaMinCutMeson){
586 //        cout << pi0->GetAlpha() << "<" << fAlphaMinCutMeson << endl; 
587       if(hist)hist->Fill(cutIndex);
588       return kFALSE;
589    }
590    cutIndex++;
591
592    if (hDCAGammaGammaMesonBefore)hDCAGammaGammaMesonBefore->Fill(pi0->GetDCABetweenPhotons());
593    if (hDCARMesonPrimVtxBefore)hDCARMesonPrimVtxBefore->Fill(pi0->GetDCARMotherPrimVtx());
594
595    if (pi0->GetDCABetweenPhotons() > fDCAGammaGammaCut){
596 //        cout << pi0->GetDCABetweenPhotons() << ">" << fDCAGammaGammaCut << endl; 
597       if(hist)hist->Fill(cutIndex);
598       return kFALSE;
599    }
600    cutIndex++;
601
602    if (pi0->GetDCARMotherPrimVtx() > fDCARMesonPrimVtxCut){
603 //         cout << pi0->GetDCARMotherPrimVtx() << ">" << fDCARMesonPrimVtxCut << endl; 
604       if(hist)hist->Fill(cutIndex);
605       return kFALSE;
606    }
607    cutIndex++;
608
609
610    if (hDCAZMesonPrimVtxBefore)hDCAZMesonPrimVtxBefore->Fill(pi0->GetDCAZMotherPrimVtx());
611
612    if (abs(pi0->GetDCAZMotherPrimVtx()) > fDCAZMesonPrimVtxCut){
613 //        cout << pi0->GetDCAZMotherPrimVtx() << ">" << fDCAZMesonPrimVtxCut << endl; 
614       if(hist)hist->Fill(cutIndex);
615       return kFALSE;
616    }
617    cutIndex++;
618
619
620    if (hDCAGammaGammaMesonAfter)hDCAGammaGammaMesonAfter->Fill(pi0->GetDCABetweenPhotons());
621    if (hDCARMesonPrimVtxAfter)hDCARMesonPrimVtxAfter->Fill(pi0->GetDCARMotherPrimVtx());
622    if (hDCAZMesonPrimVtxAfter)hDCAZMesonPrimVtxAfter->Fill(pi0->M(),pi0->GetDCAZMotherPrimVtx());
623
624    if(hist)hist->Fill(cutIndex);
625    return kTRUE;
626 }
627
628
629
630 ///________________________________________________________________________
631 ///________________________________________________________________________
632 Bool_t AliConversionMesonCuts::UpdateCutString() {
633    ///Update the cut string (if it has been created yet)
634
635    if(fCutString && fCutString->GetString().Length() == kNCuts) {
636       fCutString->SetString(GetCutNumber());
637    } else {
638       return kFALSE;
639    }
640    return kTRUE;
641 }
642
643 ///________________________________________________________________________
644 Bool_t AliConversionMesonCuts::InitializeCutsFromCutString(const TString analysisCutSelection ) {
645    // Initialize Cuts from a given Cut string
646    AliInfo(Form("Set Meson Cutnumber: %s",analysisCutSelection.Data()));
647    if(analysisCutSelection.Length()!=kNCuts) {
648       AliError(Form("Cut selection has the wrong length! size is %d, number of cuts is %d", analysisCutSelection.Length(), kNCuts));
649       return kFALSE;
650    }
651    if(!analysisCutSelection.IsDigit()){
652       AliError("Cut selection contains characters");
653       return kFALSE;
654    }
655
656    const char *cutSelection = analysisCutSelection.Data();
657 #define ASSIGNARRAY(i)  fCuts[i] = cutSelection[i] - '0'
658    for(Int_t ii=0;ii<kNCuts;ii++){
659       ASSIGNARRAY(ii);
660    }
661
662    // Set Individual Cuts
663    for(Int_t ii=0;ii<kNCuts;ii++){
664       if(!SetCut(cutIds(ii),fCuts[ii]))return kFALSE;
665    }
666
667    //PrintCuts();
668    return kTRUE;
669 }
670 ///________________________________________________________________________
671 Bool_t AliConversionMesonCuts::SetCut(cutIds cutID, const Int_t value) {
672    ///Set individual cut ID
673
674    //cout << "Updating cut  " << fgkCutNames[cutID] << " (" << cutID << ") to " << value << endl;
675    switch (cutID) {
676    case kMesonKind:
677       if( SetMesonKind(value)) {
678          fCuts[kMesonKind] = value;
679          UpdateCutString();
680          return kTRUE;
681       } else return kFALSE;
682    case kchi2MesonCut:
683       if( SetChi2MesonCut(value)) {
684          fCuts[kchi2MesonCut] = value;
685          UpdateCutString();
686          return kTRUE;
687       } else return kFALSE;
688    case kalphaMesonCut:
689       if( SetAlphaMesonCut(value)) {
690          fCuts[kalphaMesonCut] = value;
691          UpdateCutString();
692          return kTRUE;
693       } else return kFALSE;
694    case kRCut:
695       if( SetRCut(value)) {
696          fCuts[kRCut] = value;
697          UpdateCutString();
698          return kTRUE;
699       } else return kFALSE;
700
701    case kRapidityMesonCut:
702       if( SetRapidityMesonCut(value)) {
703          fCuts[kRapidityMesonCut] = value;
704          UpdateCutString();
705          return kTRUE;
706       } else return kFALSE;
707
708    case kBackgroundScheme:
709       if( SetBackgroundScheme(value)) {
710          fCuts[kBackgroundScheme] = value;
711          UpdateCutString();
712          return kTRUE;
713       } else return kFALSE;
714
715    case kDegreesForRotationMethod:
716       if( SetNDegreesForRotationMethod(value)) {
717          fCuts[kDegreesForRotationMethod] = value;
718          UpdateCutString();
719          return kTRUE;
720       } else return kFALSE;
721
722    case kNumberOfBGEvents:
723       if( SetNumberOfBGEvents(value)) {
724          fCuts[kNumberOfBGEvents] = value;
725          UpdateCutString();
726          return kTRUE;
727       } else return kFALSE;
728
729    case kuseMCPSmearing:
730       if( SetMCPSmearing(value)) {
731          fCuts[kuseMCPSmearing] = value;
732          UpdateCutString();
733          return kTRUE;
734       } else return kFALSE;
735    case kElecShare:
736       if( SetSharedElectronCut(value)) {
737          fCuts[kElecShare] = value;
738          UpdateCutString();
739          return kTRUE;
740       } else return kFALSE;
741    case kToCloseV0s:
742       if( SetToCloseV0sCut(value)) {
743          fCuts[kToCloseV0s] = value;
744          UpdateCutString();
745          return kTRUE;
746       } else return kFALSE;
747    case kDcaGammaGamma:
748       if( SetDCAGammaGammaCut(value)) {
749          fCuts[kDcaGammaGamma] = value;
750          UpdateCutString();
751          return kTRUE;
752       } else return kFALSE;
753    case kDcaZPrimVtx:
754       if( SetDCAZMesonPrimVtxCut(value)) {
755          fCuts[kDcaZPrimVtx] = value;
756          UpdateCutString();
757          return kTRUE;
758       } else return kFALSE;
759    case kDcaRPrimVtx:
760       if( SetDCARMesonPrimVtxCut(value)) {
761          fCuts[kDcaRPrimVtx] = value;
762          UpdateCutString();
763          return kTRUE;
764       } else return kFALSE;
765
766    case kNCuts:
767       cout << "Error:: Cut id out of range"<< endl;
768       return kFALSE;
769    }
770
771    cout << "Error:: Cut id " << cutID << " not recognized "<< endl;
772    return kFALSE;
773
774 }
775
776
777 ///________________________________________________________________________
778 void AliConversionMesonCuts::PrintCuts() {
779    // Print out current Cut Selection
780    for(Int_t ic = 0; ic < kNCuts; ic++) {
781       printf("%-30s : %d \n", fgkCutNames[ic], fCuts[ic]);
782    }
783 }
784
785 ///________________________________________________________________________
786 Bool_t AliConversionMesonCuts::SetMesonKind(Int_t mesonKind){
787    // Set Cut
788    switch(mesonKind){
789    case 0:
790       fMesonKind=0;
791       break;
792    case 1:
793       fMesonKind=1;;
794       break;
795    default:
796       cout<<"Warning: Meson kind not defined"<<mesonKind<<endl;
797       return kFALSE;
798    }
799    return kTRUE;
800 }
801
802 ///________________________________________________________________________
803 Bool_t AliConversionMesonCuts::SetRCut(Int_t rCut){
804    // Set Cut
805    switch(rCut){
806    case 0:
807       fMaxR = 180.;
808       break;
809    case 1:
810       fMaxR = 180.;
811       break;
812    case 2:
813       fMaxR = 180.;
814       break;
815    case 3:
816       fMaxR = 70.;
817       break;
818    case 4:
819       fMaxR = 70.;
820       break;
821    case 5:
822       fMaxR = 180.;
823       break;
824       // High purity cuts for PbPb
825    case 9:
826       fMaxR = 180.;
827       break;
828    default:
829       cout<<"Warning: rCut not defined"<<rCut<<endl;
830       return kFALSE;
831    }
832    return kTRUE;
833 }
834
835 ///________________________________________________________________________
836 Bool_t AliConversionMesonCuts::SetChi2MesonCut(Int_t chi2MesonCut){
837    // Set Cut
838    switch(chi2MesonCut){
839    case 0:  // 100.
840       fChi2CutMeson = 100.;
841       break;
842    case 1:  // 50.
843       fChi2CutMeson = 50.;
844       break;
845    case 2:  // 30.
846       fChi2CutMeson = 30.;
847       break;
848    case 3:
849       fChi2CutMeson = 200.;
850       break;
851    case 4:
852       fChi2CutMeson = 500.;
853       break;
854    case 5:
855       fChi2CutMeson = 1000.;
856       break;
857    default:
858       cout<<"Warning: Chi2MesonCut not defined "<<chi2MesonCut<<endl;
859       return kFALSE;
860    }
861    return kTRUE;
862 }
863
864
865 ///________________________________________________________________________
866 Bool_t AliConversionMesonCuts::SetAlphaMesonCut(Int_t alphaMesonCut)
867 {   // Set Cut
868    switch(alphaMesonCut){
869    case 0:      // 0- 0.7
870       fAlphaMinCutMeson  = 0.0;
871       fAlphaCutMeson     = 0.7;
872       fAlphaPtDepCut = kFALSE;
873       break;
874    case 1:      // Updated 31 October 2013 before 0.0 - 0.5
875       if( fFAlphaCut ) delete fFAlphaCut;
876       fFAlphaCut= new TF1("fFAlphaCut","[0]*tanh([1]*x)",0.,100.);
877       fFAlphaCut->SetParameter(0,0.7);
878       fFAlphaCut->SetParameter(1,1.2);
879       fAlphaMinCutMeson  =  0.0;
880       fAlphaCutMeson     = -1.0;
881       fAlphaPtDepCut = kTRUE;
882       break;
883    case 2:      // Updated 31 October 2013 before 0.5-1  
884       if( fFAlphaCut ) delete fFAlphaCut;
885       fFAlphaCut= new TF1("fFAlphaCut","[0]*tanh([1]*x)",0.,100.);
886       fFAlphaCut->SetParameter(0,0.8);
887       fFAlphaCut->SetParameter(1,1.2);
888       fAlphaMinCutMeson  =  0.0;
889       fAlphaCutMeson     = -1.0;
890       fAlphaPtDepCut = kTRUE;
891       break;
892    case 3:      // 0.0-1
893       fAlphaMinCutMeson  = 0.0;
894       fAlphaCutMeson     = 1.;
895       fAlphaPtDepCut = kFALSE;
896       break;
897    case 4:      // 0-0.65
898       fAlphaMinCutMeson  = 0.0;
899       fAlphaCutMeson     = 0.65;
900       fAlphaPtDepCut = kFALSE;
901       break;
902    case 5:      // 0-0.75
903       fAlphaMinCutMeson  = 0.0;
904       fAlphaCutMeson     = 0.75;
905       fAlphaPtDepCut = kFALSE;
906       break;
907    case 6:      // 0-0.8
908       fAlphaMinCutMeson  = 0.0;
909       fAlphaCutMeson     = 0.8;
910       fAlphaPtDepCut = kFALSE;
911       break;
912    case 7:      // 0.0-0.85
913       fAlphaMinCutMeson  = 0.0;
914       fAlphaCutMeson     = 0.85;
915       fAlphaPtDepCut = kFALSE;
916       break;
917    case 8:      // 0.0-0.6
918       fAlphaMinCutMeson  = 0.0;
919       fAlphaCutMeson     = 0.6;
920       fAlphaPtDepCut = kFALSE;
921       break;
922    case 9: // Updated 11 November 2013 before 0.0 - 0.3
923       if( fFAlphaCut ) delete fFAlphaCut;
924       fFAlphaCut= new TF1("fFAlphaCut","[0]*tanh([1]*x)",0.,100.);
925       fFAlphaCut->SetParameter(0,0.65);
926       fFAlphaCut->SetParameter(1,1.2);
927       fAlphaMinCutMeson  =  0.0;
928       fAlphaCutMeson     = -1.0;
929       fAlphaPtDepCut = kTRUE;
930       break;
931    default:
932       cout<<"Warning: AlphaMesonCut not defined "<<alphaMesonCut<<endl;
933       return kFALSE;
934    }
935    return kTRUE;
936 }
937
938 ///________________________________________________________________________
939 Bool_t AliConversionMesonCuts::SetRapidityMesonCut(Int_t RapidityMesonCut){
940    // Set Cut
941    switch(RapidityMesonCut){
942    case 0:  // changed from 0.9 to 1.35
943       fRapidityCutMeson   = 1.35;
944       break;
945    case 1:  //
946       fRapidityCutMeson   = 0.8;
947       break;
948    case 2:  //
949       fRapidityCutMeson   = 0.7;
950       break;
951    case 3:  //
952       fRapidityCutMeson   = 0.6;
953       break;
954    case 4:  //
955       fRapidityCutMeson   = 0.5;
956       break;
957    case 5:  //
958       fRapidityCutMeson   = 0.85;
959       break;
960    case 6:  //
961       fRapidityCutMeson   = 0.75;
962       break;
963    case 7:  //
964       fRapidityCutMeson   = 0.3;
965       break;
966    case 8:  //changed, before 0.35
967       fRapidityCutMeson   = 0.25;
968       break;
969    case 9:  //
970       fRapidityCutMeson   = 0.4;
971       break;
972    default:
973       cout<<"Warning: RapidityMesonCut not defined "<<RapidityMesonCut<<endl;
974       return kFALSE;
975    }
976    return kTRUE;
977 }
978
979
980 ///________________________________________________________________________
981 Bool_t AliConversionMesonCuts::SetBackgroundScheme(Int_t BackgroundScheme){
982    // Set Cut
983    switch(BackgroundScheme){
984    case 0: //Rotation
985       fUseRotationMethodInBG=kTRUE;
986       fdoBGProbability=kFALSE;
987       break;
988    case 1: // mixed event with V0 multiplicity
989       fUseRotationMethodInBG=kFALSE;
990       fUseTrackMultiplicityForBG=kFALSE;
991       fdoBGProbability=kFALSE;
992       break;
993    case 2: // mixed event with track multiplicity
994       fUseRotationMethodInBG=kFALSE;
995       fUseTrackMultiplicityForBG=kTRUE;
996       fdoBGProbability=kFALSE;
997       break;
998    case 3: //Rotation
999       fUseRotationMethodInBG=kTRUE;
1000       fdoBGProbability=kTRUE;
1001       break;
1002    case 4: //No BG calculation
1003       cout << "no BG calculation should be done" << endl;
1004       fUseRotationMethodInBG=kFALSE;
1005       fdoBGProbability=kFALSE;
1006       fDoBG=kFALSE;
1007       break;
1008    case 5: //Rotation
1009       fUseRotationMethodInBG=kTRUE;
1010       fdoBGProbability=kFALSE;
1011       fBackgroundHandler = 1;
1012       break;
1013    case 6: // mixed event with V0 multiplicity
1014       fUseRotationMethodInBG=kFALSE;
1015       fUseTrackMultiplicityForBG=kFALSE;
1016       fdoBGProbability=kFALSE;
1017       fBackgroundHandler = 1;
1018       break;
1019    case 7: // mixed event with track multiplicity
1020       fUseRotationMethodInBG=kFALSE;
1021       fUseTrackMultiplicityForBG=kTRUE;
1022       fdoBGProbability=kFALSE;
1023       fBackgroundHandler = 1;
1024       break;
1025    case 8: //Rotation
1026       fUseRotationMethodInBG=kTRUE;
1027       fdoBGProbability=kTRUE;
1028       fBackgroundHandler = 1;
1029       break;
1030    default:
1031       cout<<"Warning: BackgroundScheme not defined "<<BackgroundScheme<<endl;
1032       return kFALSE;
1033    }
1034    return kTRUE;
1035 }
1036
1037
1038 ///________________________________________________________________________
1039 Bool_t AliConversionMesonCuts::SetNDegreesForRotationMethod(Int_t DegreesForRotationMethod){
1040    // Set Cut
1041    switch(DegreesForRotationMethod){
1042    case 0:
1043       fnDegreeRotationPMForBG = 5;
1044       break;
1045    case 1:
1046       fnDegreeRotationPMForBG = 10;
1047       break;
1048    case 2:
1049       fnDegreeRotationPMForBG = 15;
1050       break;
1051    case 3:
1052       fnDegreeRotationPMForBG = 20;
1053       break;
1054    default:
1055       cout<<"Warning: DegreesForRotationMethod not defined "<<DegreesForRotationMethod<<endl;
1056       return kFALSE;
1057    }
1058    fCuts[kDegreesForRotationMethod]=DegreesForRotationMethod;
1059    return kTRUE;
1060 }
1061
1062 ///________________________________________________________________________
1063 Bool_t AliConversionMesonCuts::SetNumberOfBGEvents(Int_t NumberOfBGEvents){
1064    // Set Cut
1065    switch(NumberOfBGEvents){
1066    case 0:
1067       fNumberOfBGEvents = 5;
1068       break;
1069    case 1:
1070       fNumberOfBGEvents = 10;
1071       break;
1072    case 2:
1073       fNumberOfBGEvents = 15;
1074       break;
1075    case 3:
1076       fNumberOfBGEvents = 20;
1077       break;
1078    case 4:
1079       fNumberOfBGEvents = 2;
1080       break;
1081    case 5:
1082       fNumberOfBGEvents = 50;
1083       break;
1084    case 6:
1085       fNumberOfBGEvents = 80;
1086       break;
1087    case 7:
1088       fNumberOfBGEvents = 100;
1089       break;
1090    default:
1091       cout<<"Warning: NumberOfBGEvents not defined "<<NumberOfBGEvents<<endl;
1092       return kFALSE;
1093    }
1094    return kTRUE;
1095 }
1096 ///________________________________________________________________________
1097 Bool_t AliConversionMesonCuts::SetSharedElectronCut(Int_t sharedElec) {
1098
1099    switch(sharedElec){
1100    case 0:
1101       fDoSharedElecCut = kFALSE;
1102       break;
1103    case 1:
1104       fDoSharedElecCut = kTRUE;
1105       break;
1106    default:
1107       cout<<"Warning: Shared Electron Cut not defined "<<sharedElec<<endl;
1108       return kFALSE;
1109    }
1110
1111    return kTRUE;
1112 }
1113
1114 ///________________________________________________________________________
1115 Bool_t AliConversionMesonCuts::SetToCloseV0sCut(Int_t toClose) {
1116
1117    switch(toClose){
1118    case 0:
1119       fDoToCloseV0sCut = kFALSE;
1120       fminV0Dist = 250;
1121       break;
1122    case 1:
1123       fDoToCloseV0sCut = kTRUE;
1124       fminV0Dist = 1;
1125       break;
1126    case 2:
1127       fDoToCloseV0sCut = kTRUE;
1128       fminV0Dist = 2;
1129       break;
1130    case 3:
1131       fDoToCloseV0sCut = kTRUE;
1132       fminV0Dist = 3;
1133       break;
1134    default:
1135       cout<<"Warning: Shared Electron Cut not defined "<<toClose<<endl;
1136       return kFALSE;
1137    }
1138    return kTRUE;
1139 }
1140
1141 ///________________________________________________________________________
1142 Bool_t AliConversionMesonCuts::SetMCPSmearing(Int_t useMCPSmearing)
1143 {// Set Cut
1144    switch(useMCPSmearing){
1145    case 0:
1146       fUseMCPSmearing=0;
1147       fPBremSmearing=1.;
1148       fPSigSmearing=0.;
1149       fPSigSmearingCte=0.;
1150       break;
1151    case 1:
1152       fUseMCPSmearing=1;
1153       fPBremSmearing=1.0e-14;
1154       fPSigSmearing=0.;
1155       fPSigSmearingCte=0.;
1156       break;
1157    case 2:
1158       fUseMCPSmearing=1;
1159       fPBremSmearing=1.0e-15;
1160       fPSigSmearing=0.0;
1161       fPSigSmearingCte=0.;
1162       break;
1163    case 3:
1164       fUseMCPSmearing=1;
1165       fPBremSmearing=1.;
1166       fPSigSmearing=0.003;
1167       fPSigSmearingCte=0.002;
1168       break;
1169    case 4:
1170       fUseMCPSmearing=1;
1171       fPBremSmearing=1.;
1172       fPSigSmearing=0.003;
1173       fPSigSmearingCte=0.007;
1174       break;
1175    case 5:
1176       fUseMCPSmearing=1;
1177       fPBremSmearing=1.;
1178       fPSigSmearing=0.003;
1179       fPSigSmearingCte=0.016;
1180       break;
1181    case 6:
1182       fUseMCPSmearing=1;
1183       fPBremSmearing=1.;
1184       fPSigSmearing=0.007;
1185       fPSigSmearingCte=0.016;
1186       break;
1187    case 7:
1188       fUseMCPSmearing=1;
1189       fPBremSmearing=1.0e-16;
1190       fPSigSmearing=0.0;
1191       fPSigSmearingCte=0.;
1192       break;
1193    case 8:
1194       fUseMCPSmearing=1;
1195       fPBremSmearing=1.;
1196       fPSigSmearing=0.007;
1197       fPSigSmearingCte=0.014;
1198       break;
1199    case 9:
1200       fUseMCPSmearing=1;
1201       fPBremSmearing=1.;
1202       fPSigSmearing=0.007;
1203       fPSigSmearingCte=0.011;
1204       break;
1205
1206    default:
1207       AliError("Warning: UseMCPSmearing not defined");
1208       return kFALSE;
1209    }
1210    return kTRUE;
1211 }
1212
1213
1214 ///________________________________________________________________________
1215 Bool_t AliConversionMesonCuts::SetDCAGammaGammaCut(Int_t DCAGammaGamma){
1216    // Set Cut
1217    switch(DCAGammaGamma){
1218    case 0:  //
1219       fDCAGammaGammaCut   = 1000;
1220       break;
1221    case 1:  //
1222       fDCAGammaGammaCut   = 10;
1223       break;
1224    case 2:  //
1225       fDCAGammaGammaCut   = 5;
1226       break;
1227    case 3:  //
1228       fDCAGammaGammaCut   = 4;
1229       break;
1230    case 4:  //
1231       fDCAGammaGammaCut   = 3;
1232       break;
1233    case 5:  //
1234       fDCAGammaGammaCut   = 2.5;
1235       break;
1236    case 6:  //
1237       fDCAGammaGammaCut   = 2;
1238       break;
1239    case 7:  //
1240       fDCAGammaGammaCut   = 1.5;
1241       break;
1242    case 8:  //
1243       fDCAGammaGammaCut   = 1;
1244       break;
1245    case 9:  //
1246       fDCAGammaGammaCut   = 0.5;
1247       break;
1248    default:
1249       cout<<"Warning: DCAGammaGamma not defined "<<DCAGammaGamma<<endl;
1250       return kFALSE;
1251    }
1252    return kTRUE;
1253 }
1254
1255 ///________________________________________________________________________
1256 Bool_t AliConversionMesonCuts::SetDCAZMesonPrimVtxCut(Int_t DCAZMesonPrimVtx){
1257    // Set Cut
1258    switch(DCAZMesonPrimVtx){
1259    case 0:  //
1260       fDCAZMesonPrimVtxCut   = 1000;
1261       break;
1262    case 1:  //
1263       fDCAZMesonPrimVtxCut   = 10;
1264       break;
1265    case 2:  //
1266       fDCAZMesonPrimVtxCut   = 5;
1267       break;
1268    case 3:  //
1269       fDCAZMesonPrimVtxCut   = 4;
1270       break;
1271    case 4:  //
1272       fDCAZMesonPrimVtxCut   = 3;
1273       break;
1274    case 5:  //
1275       fDCAZMesonPrimVtxCut   = 2.5;
1276       break;
1277    case 6:  //
1278       fDCAZMesonPrimVtxCut   = 2;
1279       break;
1280    case 7:  //
1281       fDCAZMesonPrimVtxCut   = 1.5;
1282       break;
1283    case 8:  //
1284       fDCAZMesonPrimVtxCut   = 1;
1285       break;
1286    case 9:  //
1287       fDCAZMesonPrimVtxCut   = 0.5;
1288       break;
1289    default:
1290       cout<<"Warning: DCAZMesonPrimVtx not defined "<<DCAZMesonPrimVtx<<endl;
1291       return kFALSE;
1292    }
1293    return kTRUE;
1294 }
1295
1296 ///________________________________________________________________________
1297 Bool_t AliConversionMesonCuts::SetDCARMesonPrimVtxCut(Int_t DCARMesonPrimVtx){
1298    // Set Cut
1299    switch(DCARMesonPrimVtx){
1300    case 0:  //
1301       fDCARMesonPrimVtxCut   = 1000;
1302       break;
1303    case 1:  //
1304       fDCARMesonPrimVtxCut   = 10;
1305       break;
1306    case 2:  //
1307       fDCARMesonPrimVtxCut   = 5;
1308       break;
1309    case 3:  //
1310       fDCARMesonPrimVtxCut   = 4;
1311       break;
1312    case 4:  //
1313       fDCARMesonPrimVtxCut   = 3;
1314       break;
1315    case 5:  //
1316       fDCARMesonPrimVtxCut   = 2.5;
1317       break;
1318    case 6:  //
1319       fDCARMesonPrimVtxCut   = 2;
1320       break;
1321    case 7:  //
1322       fDCARMesonPrimVtxCut   = 1.5;
1323       break;
1324    case 8:  //
1325       fDCARMesonPrimVtxCut   = 1;
1326       break;
1327    case 9:  //
1328       fDCARMesonPrimVtxCut   = 0.5;
1329       break;
1330    default:
1331       cout<<"Warning: DCARMesonPrimVtx not defined "<<DCARMesonPrimVtx<<endl;
1332       return kFALSE;
1333    }
1334    return kTRUE;
1335 }
1336
1337
1338 ///________________________________________________________________________
1339 TString AliConversionMesonCuts::GetCutNumber(){
1340    // returns TString with current cut number
1341    TString a(kNCuts);
1342    for(Int_t ii=0;ii<kNCuts;ii++){
1343       a.Append(Form("%d",fCuts[ii]));
1344    }
1345    return a;
1346 }
1347
1348 ///________________________________________________________________________
1349 void AliConversionMesonCuts::FillElectonLabelArray(AliAODConversionPhoton* photon, Int_t nV0){
1350
1351    Int_t posLabel = photon->GetTrackLabelPositive();
1352    Int_t negLabel = photon->GetTrackLabelNegative();
1353
1354    fElectronLabelArray[nV0*2] = posLabel;
1355    fElectronLabelArray[(nV0*2)+1] = negLabel;
1356 }
1357
1358 ///________________________________________________________________________
1359 Bool_t AliConversionMesonCuts::RejectSharedElectronV0s(AliAODConversionPhoton* photon, Int_t nV0, Int_t nV0s){
1360
1361    Int_t posLabel = photon->GetTrackLabelPositive();
1362    Int_t negLabel = photon->GetTrackLabelNegative();
1363
1364    for(Int_t i = 0; i<nV0s*2;i++){
1365       if(i==nV0*2)     continue;
1366       if(i==(nV0*2)+1) continue;
1367       if(fElectronLabelArray[i] == posLabel){
1368          return kFALSE;}
1369       if(fElectronLabelArray[i] == negLabel){
1370          return kFALSE;}
1371    }
1372
1373    return kTRUE;
1374 }
1375 ///________________________________________________________________________
1376 Bool_t AliConversionMesonCuts::RejectToCloseV0s(AliAODConversionPhoton* photon, TList *photons, Int_t nV0){
1377    Double_t posX = photon->GetConversionX();
1378    Double_t posY = photon->GetConversionY();
1379    Double_t posZ = photon->GetConversionZ();
1380
1381    for(Int_t i = 0;i<photons->GetEntries();i++){
1382       if(nV0 == i) continue;
1383       AliAODConversionPhoton *photonComp = (AliAODConversionPhoton*) photons->At(i);
1384       Double_t posCompX = photonComp->GetConversionX();
1385       Double_t posCompY = photonComp->GetConversionY();
1386       Double_t posCompZ = photonComp->GetConversionZ();
1387
1388       Double_t dist = pow((posX - posCompX),2)+pow((posY - posCompY),2)+pow((posZ - posCompZ),2);
1389
1390       if(dist < fminV0Dist*fminV0Dist){
1391          if(photon->GetChi2perNDF() < photonComp->GetChi2perNDF()) return kTRUE;
1392          else {
1393             return kFALSE;}
1394       }
1395
1396    }
1397    return kTRUE;
1398 }
1399
1400 ///________________________________________________________________________
1401 void AliConversionMesonCuts::SmearParticle(AliAODConversionPhoton* photon)
1402 {
1403
1404    if (photon==NULL) return;
1405    Double_t facPBrem = 1.;
1406    Double_t facPSig = 0.;
1407
1408    Double_t phi=0.;
1409    Double_t theta=0.;
1410    Double_t P=0.;
1411
1412
1413    P=photon->P();
1414    phi=photon->Phi();
1415    if( photon->P()!=0){
1416       theta=acos( photon->Pz()/ photon->P());
1417    }
1418
1419    if( fPSigSmearing != 0. || fPSigSmearingCte!=0. ){
1420       facPSig = TMath::Sqrt(fPSigSmearingCte*fPSigSmearingCte+fPSigSmearing*fPSigSmearing*P*P)*fRandom.Gaus(0.,1.);
1421    }
1422
1423    if( fPBremSmearing != 1.){
1424       if(fBrem!=NULL){
1425          facPBrem = fBrem->GetRandom();
1426       }
1427    }
1428
1429    photon->SetPx(facPBrem* (1+facPSig)* P*sin(theta)*cos(phi)) ;
1430    photon->SetPy(facPBrem* (1+facPSig)* P*sin(theta)*sin(phi)) ;
1431    photon->SetPz(facPBrem* (1+facPSig)* P*cos(theta)) ;
1432    photon->SetE(photon->P());
1433 }
1434 ///________________________________________________________________________
1435 void AliConversionMesonCuts::SmearVirtualPhoton(AliAODConversionPhoton* photon)
1436 {
1437
1438    if (photon==NULL) return;
1439    Double_t facPBrem = 1.;
1440    Double_t facPSig = 0.;
1441
1442    Double_t phi=0.;
1443    Double_t theta=0.;
1444    Double_t P=0.;
1445
1446
1447    P=photon->P();
1448    phi=photon->Phi();
1449    if( photon->P()!=0){
1450       theta=acos( photon->Pz()/ photon->P());
1451    }
1452
1453    if( fPSigSmearing != 0. || fPSigSmearingCte!=0. ){
1454       facPSig = TMath::Sqrt(fPSigSmearingCte*fPSigSmearingCte+fPSigSmearing*fPSigSmearing*P*P)*fRandom.Gaus(0.,1.);
1455    }
1456
1457    if( fPBremSmearing != 1.){
1458       if(fBrem!=NULL){
1459          facPBrem = fBrem->GetRandom();
1460       }
1461    }
1462
1463    photon->SetPx(facPBrem* (1+facPSig)* P*sin(theta)*cos(phi)) ;
1464    photon->SetPy(facPBrem* (1+facPSig)* P*sin(theta)*sin(phi)) ;
1465    photon->SetPz(facPBrem* (1+facPSig)* P*cos(theta)) ;
1466    
1467    
1468 }
1469