updated task to make it possible to run on the grid
[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 class iostream;
43
44 using namespace std;
45
46 ClassImp(AliConversionMesonCuts)
47
48
49 const char* AliConversionMesonCuts::fgkCutNames[AliConversionMesonCuts::kNCuts] = {
50    "MesonKind",
51    "BackgroundScheme",
52    "NumberOfBGEvents",
53    "DegreesForRotationMethod",
54    "RapidityMesonCut",
55    "RCut",
56    "AlphaMesonCut",
57    "Chi2MesonCut",
58    "SharedElectronCuts",
59    "RejectToCloseV0s",
60    "UseMCPSmearing",
61 };
62
63
64 //________________________________________________________________________
65 AliConversionMesonCuts::AliConversionMesonCuts(const char *name,const char *title) :
66    AliAnalysisCuts(name,title),
67    fHistograms(NULL),
68    fMesonKind(0),
69    fMaxR(200),
70    fChi2CutMeson(1000),
71    fAlphaMinCutMeson(0),
72    fAlphaCutMeson(1),
73    fRapidityCutMeson(1),
74    fUseRotationMethodInBG(kFALSE),
75    fDoBG(kTRUE),
76    fdoBGProbability(kFALSE),
77    fUseTrackMultiplicityForBG(kFALSE),
78    fnDegreeRotationPMForBG(0),
79    fNumberOfBGEvents(0),
80    fOpeningAngle(0.005),
81    fDoToCloseV0sCut(kFALSE),
82    fminV0Dist(200.),
83    fDoSharedElecCut(kFALSE),
84    fUseMCPSmearing(kFALSE),
85    fPBremSmearing(0),
86    fPSigSmearing(0),
87    fPSigSmearingCte(0),
88    fBrem(NULL),
89    fRandom(0),
90    fElectronLabelArraySize(500),
91    fElectronLabelArray(NULL),
92    fBackgroundHandler(0),
93    fCutString(NULL),
94    hMesonCuts(NULL),
95    hMesonBGCuts(NULL)
96    
97 {
98    for(Int_t jj=0;jj<kNCuts;jj++){fCuts[jj]=0;}
99    fCutString=new TObjString((GetCutNumber()).Data());
100    fElectronLabelArray = new Int_t[fElectronLabelArraySize];
101    if (fBrem == NULL){
102       fBrem = new TF1("fBrem","pow(-log(x),[0]/log(2.0)-1.0)/TMath::Gamma([0]/log(2.0))",0.00001,0.999999999);
103       // tests done with 1.0e-14
104       fBrem->SetParameter(0,fPBremSmearing);
105       fBrem->SetNpx(100000);
106    }
107
108 }
109
110 //________________________________________________________________________
111 AliConversionMesonCuts::AliConversionMesonCuts(const AliConversionMesonCuts &ref) :
112    AliAnalysisCuts(ref),
113    fHistograms(NULL),
114    fMesonKind(ref.fMesonKind),
115    fMaxR(ref.fMaxR),
116    fChi2CutMeson(ref.fChi2CutMeson),
117    fAlphaMinCutMeson(ref.fAlphaMinCutMeson),
118    fAlphaCutMeson(ref.fAlphaCutMeson),
119    fRapidityCutMeson(ref.fRapidityCutMeson),
120    fUseRotationMethodInBG(ref.fUseRotationMethodInBG),
121    fDoBG(ref.fDoBG),
122    fdoBGProbability(ref.fdoBGProbability),
123    fUseTrackMultiplicityForBG(ref.fUseTrackMultiplicityForBG),
124    fnDegreeRotationPMForBG(ref.fnDegreeRotationPMForBG),
125    fNumberOfBGEvents(ref. fNumberOfBGEvents),
126    fOpeningAngle(ref.fOpeningAngle),
127    fDoToCloseV0sCut(ref.fDoToCloseV0sCut),
128    fminV0Dist(ref.fminV0Dist),
129    fDoSharedElecCut(ref.fDoSharedElecCut),
130    fUseMCPSmearing(ref.fUseMCPSmearing),
131    fPBremSmearing(ref.fPBremSmearing),
132    fPSigSmearing(ref.fPSigSmearing),
133    fPSigSmearingCte(ref.fPSigSmearingCte),
134    fBrem(NULL),
135    fRandom(ref.fRandom),
136    fElectronLabelArraySize(ref.fElectronLabelArraySize),
137    fElectronLabelArray(NULL),
138    fBackgroundHandler(ref.fBackgroundHandler),
139    fCutString(NULL),
140    hMesonCuts(NULL),
141    hMesonBGCuts(NULL)
142    
143 {
144    // Copy Constructor
145    for(Int_t jj=0;jj<kNCuts;jj++){fCuts[jj]=ref.fCuts[jj];}
146    fCutString=new TObjString((GetCutNumber()).Data());
147    fElectronLabelArray = new Int_t[fElectronLabelArraySize];
148    if (fBrem == NULL)fBrem = (TF1*)ref.fBrem->Clone("fBrem");
149    // Histograms are not copied, if you need them, call InitCutHistograms
150 }
151
152
153 //________________________________________________________________________
154 AliConversionMesonCuts::~AliConversionMesonCuts() {
155    // Destructor
156    //Deleting fHistograms leads to seg fault it it's added to output collection of a task
157    // if(fHistograms)
158    //   delete fHistograms;
159    // fHistograms = NULL;
160    if(fCutString != NULL){
161       delete fCutString;
162       fCutString = NULL;
163    }
164    if(fElectronLabelArray){
165       delete fElectronLabelArray;
166       fElectronLabelArray = NULL;
167    }
168
169 }
170
171 //________________________________________________________________________
172 void AliConversionMesonCuts::InitCutHistograms(TString name){
173
174    // Initialize Cut Histograms for QA (only initialized and filled if function is called)
175
176    if(fHistograms != NULL){
177       delete fHistograms;
178       fHistograms=NULL;
179    }
180    if(fHistograms==NULL){
181       fHistograms=new TList();
182       if(name=="")fHistograms->SetName(Form("ConvMesonCuts_%s",GetCutNumber().Data()));
183       else fHistograms->SetName(Form("%s_%s",name.Data(),GetCutNumber().Data()));
184    }
185
186    // Meson Cuts
187    hMesonCuts=new TH1F(Form("MesonCuts %s",GetCutNumber().Data()),"MesonCuts",10,-0.5,9.5);
188    hMesonCuts->GetXaxis()->SetBinLabel(1,"in");
189    hMesonCuts->GetXaxis()->SetBinLabel(2,"undef rapidity");
190    hMesonCuts->GetXaxis()->SetBinLabel(3,"rapidity cut");
191    hMesonCuts->GetXaxis()->SetBinLabel(4,"opening angle");
192    hMesonCuts->GetXaxis()->SetBinLabel(5,"alpha max");
193    hMesonCuts->GetXaxis()->SetBinLabel(6,"alpha min");
194    hMesonCuts->GetXaxis()->SetBinLabel(7,"out");
195    fHistograms->Add(hMesonCuts);
196
197    hMesonBGCuts=new TH1F(Form("MesonBGCuts %s",GetCutNumber().Data()),"MesonBGCuts",10,-0.5,9.5);
198    hMesonBGCuts->GetXaxis()->SetBinLabel(1,"in");
199    hMesonBGCuts->GetXaxis()->SetBinLabel(2,"undef rapidity");
200    hMesonBGCuts->GetXaxis()->SetBinLabel(3,"rapidity cut");
201    hMesonBGCuts->GetXaxis()->SetBinLabel(4,"opening angle");
202    hMesonBGCuts->GetXaxis()->SetBinLabel(5,"alpha max");
203    hMesonBGCuts->GetXaxis()->SetBinLabel(6,"alpha min");
204    hMesonBGCuts->GetXaxis()->SetBinLabel(7,"out");
205    fHistograms->Add(hMesonBGCuts);
206 }
207
208
209 //________________________________________________________________________
210 Bool_t AliConversionMesonCuts::MesonIsSelectedMC(TParticle *fMCMother,AliStack *fMCStack){
211    // Returns true for all pions within acceptance cuts for decay into 2 photons
212    // If bMCDaughtersInAcceptance is selected, it requires in addition that both daughter photons are within acceptance cuts
213    
214    if(!fMCStack)return kFALSE;
215    
216    if(fMCMother->GetPdgCode()==111 || fMCMother->GetPdgCode()==221){                    
217       if(fMCMother->R()>fMaxR)  return kFALSE; // cuts on distance from collision point
218       
219       Double_t rapidity = 10.;
220       if(fMCMother->Energy() - fMCMother->Pz() == 0 || fMCMother->Energy() + fMCMother->Pz() == 0){
221          rapidity=8.;
222       } else{
223          rapidity = 0.5*(TMath::Log((fMCMother->Energy()+fMCMother->Pz()) / (fMCMother->Energy()-fMCMother->Pz())));
224       } 
225       
226       // Rapidity Cut
227       if(abs(rapidity)>fRapidityCutMeson)return kFALSE;
228       
229       // Select only -> 2y decay channel
230       if(fMCMother->GetNDaughters()!=2)return kFALSE;
231       
232       for(Int_t i=0;i<2;i++){
233          TParticle *MDaughter=fMCStack->Particle(fMCMother->GetDaughter(i));
234          // Is Daughter a Photon?
235          if(MDaughter->GetPdgCode()!=22)return kFALSE;
236          // Is Photon in Acceptance?
237          //   if(bMCDaughtersInAcceptance){
238          //     if(!PhotonIsSelectedMC(MDaughter,fMCStack)){return kFALSE;}
239          //   }
240       }
241       return kTRUE;
242    }
243    return kFALSE;
244 }
245
246 //________________________________________________________________________
247 Bool_t AliConversionMesonCuts::MesonIsSelectedMCDalitz(TParticle *fMCMother,AliStack *fMCStack){
248    // Returns true for all pions within acceptance cuts for decay into 2 photons
249    // If bMCDaughtersInAcceptance is selected, it requires in addition that both daughter photons are within acceptance cuts
250
251    if(!fMCStack)return kFALSE;
252         
253    if(fMCMother->GetPdgCode()==111 || fMCMother->GetPdgCode()==221){
254                 
255       if(fMCMother->R()>fMaxR)  return kFALSE; // cuts on distance from collision point
256
257       Double_t rapidity = 10.;
258       if(fMCMother->Energy() - fMCMother->Pz() == 0 || fMCMother->Energy() + fMCMother->Pz() == 0){
259          rapidity=8.;
260       }
261       else{
262          rapidity = 0.5*(TMath::Log((fMCMother->Energy()+fMCMother->Pz()) / (fMCMother->Energy()-fMCMother->Pz())));
263       } 
264                 
265       // Rapidity Cut
266       if(abs(rapidity)>fRapidityCutMeson)return kFALSE;
267
268       // Select only -> Dalitz decay channel
269       if(fMCMother->GetNDaughters()!=3)return kFALSE;
270
271       Int_t daughterPDGs[3] = {0,0,0};
272       Int_t index = 0;
273
274       //                iParticle->GetFirstDaughter(); idxPi0 <= iParticle->GetLastDaughter()
275
276       for(Int_t i=fMCMother->GetFirstDaughter(); i<= fMCMother->GetLastDaughter();i++){
277          TParticle *MDaughter=fMCStack->Particle(i);
278          // Is Daughter a Photon or an electron?
279          daughterPDGs[index]=MDaughter->GetPdgCode();
280          index++;
281       }
282       for (Int_t j=0;j<2;j++){
283
284          for (Int_t i=0;i<2;i++){
285             if (daughterPDGs[i] > daughterPDGs[i+1]){
286                Int_t interMed = daughterPDGs[i] ; 
287                daughterPDGs[i] = daughterPDGs[i+1];
288                daughterPDGs[i+1] = interMed;
289             }
290          }
291       }
292       if (daughterPDGs[0] == -11 && daughterPDGs[1] == 11 && daughterPDGs[2] == 22) return kTRUE;
293    }
294    return kFALSE;
295 }
296
297
298 ///________________________________________________________________________
299 Bool_t AliConversionMesonCuts::MesonIsSelected(AliAODConversionMother *pi0,Bool_t IsSignal)
300 {
301    // Selection of reconstructed Meson candidates
302    // Use flag IsSignal in order to fill Fill different
303    // histograms for Signal and Background
304    TH1 *hist=0x0;
305
306    if(IsSignal){hist=hMesonCuts;}
307    else{hist=hMesonBGCuts;}
308
309    Int_t cutIndex=0;
310
311    if(hist)hist->Fill(cutIndex);
312    cutIndex++;
313
314    // Undefined Rapidity -> Floating Point exception
315    if((pi0->E()+pi0->Pz())/(pi0->E()-pi0->Pz())<=0){
316       if(hist)hist->Fill(cutIndex);
317       cutIndex++;
318       return kFALSE;
319    }
320    else{
321       // PseudoRapidity Cut --> But we cut on Rapidity !!!
322       cutIndex++;
323       if(abs(pi0->Rapidity())>fRapidityCutMeson){
324          if(hist)hist->Fill(cutIndex);
325          return kFALSE;
326       }
327    }
328    cutIndex++;
329
330    // Opening Angle Cut
331    //fOpeningAngle=2*TMath::ATan(0.134/pi0->P());// physical minimum opening angle
332    if( pi0->GetOpeningAngle() < fOpeningAngle){
333       if(hist)hist->Fill(cutIndex);
334       return kFALSE;
335    }
336    cutIndex++;
337
338    // Alpha Max Cut
339    if(pi0->GetAlpha()>fAlphaCutMeson){
340       if(hist)hist->Fill(cutIndex);
341       return kFALSE;
342    }
343    cutIndex++;
344
345    // Alpha Min Cut
346    if(pi0->GetAlpha()<fAlphaMinCutMeson){
347       if(hist)hist->Fill(cutIndex);
348       return kFALSE;
349    }
350    cutIndex++;
351
352    if(hist)hist->Fill(cutIndex);
353    return kTRUE;
354 }
355
356
357
358 ///________________________________________________________________________
359 ///________________________________________________________________________
360 Bool_t AliConversionMesonCuts::UpdateCutString() {
361    ///Update the cut string (if it has been created yet)
362
363    if(fCutString && fCutString->GetString().Length() == kNCuts) {
364       //         cout << "Updating cut id in spot number " << cutID << " to " << value << endl;
365       fCutString->SetString(GetCutNumber());
366    } else {
367       //         cout << "fCutString not yet initialized, will not be updated" << endl;
368       return kFALSE;
369    }
370    //   cout << fCutString->GetString().Data() << endl;
371    return kTRUE;
372 }
373
374 ///________________________________________________________________________
375 Bool_t AliConversionMesonCuts::InitializeCutsFromCutString(const TString analysisCutSelection ) {
376    // Initialize Cuts from a given Cut string
377    AliInfo(Form("Set Meson Cutnumber: %s",analysisCutSelection.Data()));
378    if(analysisCutSelection.Length()!=kNCuts) {
379       AliError(Form("Cut selection has the wrong length! size is %d, number of cuts is %d", analysisCutSelection.Length(), kNCuts));
380       return kFALSE;
381    }
382    if(!analysisCutSelection.IsDigit()){
383       AliError("Cut selection contains characters");
384       return kFALSE;
385    }
386         
387    const char *cutSelection = analysisCutSelection.Data();
388 #define ASSIGNARRAY(i)  fCuts[i] = cutSelection[i] - '0'
389    for(Int_t ii=0;ii<kNCuts;ii++){
390       ASSIGNARRAY(ii);
391    }
392
393    // Set Individual Cuts
394    for(Int_t ii=0;ii<kNCuts;ii++){
395       if(!SetCut(cutIds(ii),fCuts[ii]))return kFALSE;
396    }
397
398    //PrintCuts();
399    return kTRUE;
400 }
401 ///________________________________________________________________________
402 Bool_t AliConversionMesonCuts::SetCut(cutIds cutID, const Int_t value) {
403    ///Set individual cut ID
404
405    //cout << "Updating cut  " << fgkCutNames[cutID] << " (" << cutID << ") to " << value << endl;
406    switch (cutID) {
407    case kMesonKind:
408       if( SetMesonKind(value)) {
409          fCuts[kMesonKind] = value;
410          UpdateCutString();
411          return kTRUE;
412       } else return kFALSE;
413    case kchi2MesonCut:
414       if( SetChi2MesonCut(value)) {
415          fCuts[kchi2MesonCut] = value;
416          UpdateCutString();
417          return kTRUE;
418       } else return kFALSE;
419    case kalphaMesonCut:
420       if( SetAlphaMesonCut(value)) {
421          fCuts[kalphaMesonCut] = value;
422          UpdateCutString();
423          return kTRUE;
424       } else return kFALSE;
425    case kRCut:
426       if( SetRCut(value)) {
427          fCuts[kRCut] = value;
428          UpdateCutString();
429          return kTRUE;
430       } else return kFALSE;
431
432    case kRapidityMesonCut:
433       if( SetRapidityMesonCut(value)) {
434          fCuts[kRapidityMesonCut] = value;
435          UpdateCutString();
436          return kTRUE;
437       } else return kFALSE;
438
439    case kBackgroundScheme:
440       if( SetBackgroundScheme(value)) {
441          fCuts[kBackgroundScheme] = value;
442          UpdateCutString();
443          return kTRUE;
444       } else return kFALSE;
445
446    case kDegreesForRotationMethod:
447       if( SetNDegreesForRotationMethod(value)) {
448          fCuts[kDegreesForRotationMethod] = value;
449          UpdateCutString();
450          return kTRUE;
451       } else return kFALSE;
452
453    case kNumberOfBGEvents:
454       if( SetNumberOfBGEvents(value)) {
455          fCuts[kNumberOfBGEvents] = value;
456          UpdateCutString();
457          return kTRUE;
458       } else return kFALSE;
459
460    case kuseMCPSmearing:
461       if( SetMCPSmearing(value)) {
462          fCuts[kuseMCPSmearing] = value;
463          UpdateCutString();
464          return kTRUE;
465       } else return kFALSE;
466    case kElecShare:
467       if( SetSharedElectronCut(value)) {
468          fCuts[kElecShare] = value;
469          UpdateCutString();
470          return kTRUE;
471       } else return kFALSE;
472    case kToCloseV0s:
473       if( SetToCloseV0sCut(value)) {
474          fCuts[kToCloseV0s] = value;
475          UpdateCutString();
476          return kTRUE;
477       } else return kFALSE;
478
479    case kNCuts:
480       cout << "Error:: Cut id out of range"<< endl;
481       return kFALSE;
482    }
483   
484    cout << "Error:: Cut id " << cutID << " not recognized "<< endl;
485    return kFALSE;
486   
487 }
488
489
490 ///________________________________________________________________________
491 void AliConversionMesonCuts::PrintCuts() {
492    // Print out current Cut Selection
493    for(Int_t ic = 0; ic < kNCuts; ic++) {
494       printf("%-30s : %d \n", fgkCutNames[ic], fCuts[ic]);
495    }
496 }
497
498 ///________________________________________________________________________
499 Bool_t AliConversionMesonCuts::SetMesonKind(Int_t mesonKind){
500    // Set Cut
501    switch(mesonKind){
502    case 0:
503       fMesonKind=0;
504       break;
505    case 1:
506       fMesonKind=1;;
507       break;
508    default:
509       cout<<"Warning: Meson kind not defined"<<mesonKind<<endl;
510       return kFALSE;
511    }
512    return kTRUE;
513 }
514
515 ///________________________________________________________________________
516 Bool_t AliConversionMesonCuts::SetRCut(Int_t rCut){
517    // Set Cut
518    switch(rCut){
519    case 0:
520       fMaxR = 180.;
521       break;
522    case 1:
523       fMaxR = 180.;
524       break;
525    case 2:
526       fMaxR = 180.;
527       break;
528    case 3:
529       fMaxR = 70.;                      
530       break;
531    case 4:
532       fMaxR = 70.;
533       break;
534    case 5:
535       fMaxR = 180.;
536       break;
537       // High purity cuts for PbPb
538    case 9:
539       fMaxR = 180.;
540       break;
541    default:
542       cout<<"Warning: rCut not defined"<<rCut<<endl;
543       return kFALSE;
544    }
545    return kTRUE;
546 }
547
548 ///________________________________________________________________________
549 Bool_t AliConversionMesonCuts::SetChi2MesonCut(Int_t chi2MesonCut){
550    // Set Cut
551    switch(chi2MesonCut){
552    case 0:  // 100.
553       fChi2CutMeson = 100.;
554       break;
555    case 1:  // 50.
556       fChi2CutMeson = 50.;
557       break;
558    case 2:  // 30.
559       fChi2CutMeson = 30.;
560       break;
561    case 3:
562       fChi2CutMeson = 200.;
563       break;
564    case 4:
565       fChi2CutMeson = 500.;
566       break;
567    case 5:
568       fChi2CutMeson = 1000.;
569       break;
570    default:
571       cout<<"Warning: Chi2MesonCut not defined "<<chi2MesonCut<<endl;
572       return kFALSE;
573    }
574    return kTRUE;
575 }
576
577
578 ///________________________________________________________________________
579 Bool_t AliConversionMesonCuts::SetAlphaMesonCut(Int_t alphaMesonCut)
580 {   // Set Cut
581    switch(alphaMesonCut){
582    case 0:      // 0- 0.7
583       fAlphaMinCutMeson  = 0.0;
584       fAlphaCutMeson     = 0.7;
585       break;
586    case 1:      // 0-0.5
587       fAlphaMinCutMeson  = 0.0;
588       fAlphaCutMeson     = 0.5;
589       break;
590    case 2:      // 0.5-1
591       fAlphaMinCutMeson  = 0.5;
592       fAlphaCutMeson     = 1.;
593       break;
594    case 3:      // 0.0-1
595       fAlphaMinCutMeson  = 0.0;
596       fAlphaCutMeson     = 1.;
597       break;
598    case 4:      // 0-0.65
599       fAlphaMinCutMeson  = 0.0;
600       fAlphaCutMeson     = 0.65;
601       break;
602    case 5:      // 0-0.75
603       fAlphaMinCutMeson  = 0.0;
604       fAlphaCutMeson     = 0.75;
605       break;
606    case 6:      // 0-0.8
607       fAlphaMinCutMeson  = 0.0;
608       fAlphaCutMeson     = 0.8;
609       break;
610    case 7:      // 0.0-0.85
611       fAlphaMinCutMeson  = 0.0;
612       fAlphaCutMeson     = 0.85;
613       break;
614    case 8:      // 0.0-0.6
615       fAlphaMinCutMeson  = 0.0;
616       fAlphaCutMeson     = 0.6;
617       break;
618    case 9:      // 0.0-0.3
619       fAlphaMinCutMeson  = 0.0;
620       fAlphaCutMeson     = 0.3;
621       break;
622    default:
623       cout<<"Warning: AlphaMesonCut not defined "<<alphaMesonCut<<endl;
624       return kFALSE;
625    }
626    return kTRUE;
627 }
628
629 ///________________________________________________________________________
630 Bool_t AliConversionMesonCuts::SetRapidityMesonCut(Int_t RapidityMesonCut){ 
631    // Set Cut
632    switch(RapidityMesonCut){
633    case 0:  //
634       fRapidityCutMeson   = 0.9;
635       break;
636    case 1:  //
637       fRapidityCutMeson   = 0.8;
638       break;
639    case 2:  //
640       fRapidityCutMeson   = 0.7;
641       break;
642    case 3:  //
643       fRapidityCutMeson   = 0.6;
644       break;
645    case 4:  //
646       fRapidityCutMeson   = 0.5;
647       break;
648    case 5:  //
649       fRapidityCutMeson   = 0.85;
650       break;
651    case 6:  //
652       fRapidityCutMeson   = 0.75;
653       break;
654
655    default:
656       cout<<"Warning: RapidityMesonCut not defined "<<RapidityMesonCut<<endl;
657       return kFALSE;
658    }
659    return kTRUE;
660 }
661
662
663 ///________________________________________________________________________
664 Bool_t AliConversionMesonCuts::SetBackgroundScheme(Int_t BackgroundScheme){
665    // Set Cut
666    switch(BackgroundScheme){
667    case 0: //Rotation
668       fUseRotationMethodInBG=kTRUE;
669       fdoBGProbability=kFALSE;
670       break;
671    case 1: // mixed event with V0 multiplicity
672       fUseRotationMethodInBG=kFALSE;
673       fUseTrackMultiplicityForBG=kFALSE;
674       fdoBGProbability=kFALSE;
675       break;
676    case 2: // mixed event with track multiplicity
677       fUseRotationMethodInBG=kFALSE;
678       fUseTrackMultiplicityForBG=kTRUE;
679       fdoBGProbability=kFALSE;
680       break;
681    case 3: //Rotation
682       fUseRotationMethodInBG=kTRUE;
683       fdoBGProbability=kTRUE;
684       break;
685    case 4: //No BG calculation
686       cout << "no BG calculation should be done" << endl;
687       fUseRotationMethodInBG=kFALSE;
688       fdoBGProbability=kFALSE;
689       fDoBG=kFALSE;
690       break;
691    case 5: //Rotation
692       fUseRotationMethodInBG=kTRUE;
693       fdoBGProbability=kFALSE;
694       fBackgroundHandler = 1;
695       break;
696    case 6: // mixed event with V0 multiplicity
697       fUseRotationMethodInBG=kFALSE;
698       fUseTrackMultiplicityForBG=kFALSE;
699       fdoBGProbability=kFALSE;
700       fBackgroundHandler = 1;
701       break;
702    case 7: // mixed event with track multiplicity
703       fUseRotationMethodInBG=kFALSE;
704       fUseTrackMultiplicityForBG=kTRUE;
705       fdoBGProbability=kFALSE;
706       fBackgroundHandler = 1;
707       break;
708    case 8: //Rotation
709       fUseRotationMethodInBG=kTRUE;
710       fdoBGProbability=kTRUE;
711       fBackgroundHandler = 1;
712       break;
713    default:
714       cout<<"Warning: BackgroundScheme not defined "<<BackgroundScheme<<endl;
715       return kFALSE;
716    }
717    return kTRUE;
718 }
719
720
721 ///________________________________________________________________________
722 Bool_t AliConversionMesonCuts::SetNDegreesForRotationMethod(Int_t DegreesForRotationMethod){
723    // Set Cut
724    switch(DegreesForRotationMethod){
725    case 0:
726       fnDegreeRotationPMForBG = 5;
727       break;
728    case 1:
729       fnDegreeRotationPMForBG = 10;
730       break;
731    case 2:
732       fnDegreeRotationPMForBG = 15;
733       break;
734    case 3:
735       fnDegreeRotationPMForBG = 20;
736       break;
737    default:
738       cout<<"Warning: DegreesForRotationMethod not defined "<<DegreesForRotationMethod<<endl;
739       return kFALSE;
740    }
741    fCuts[kDegreesForRotationMethod]=DegreesForRotationMethod;
742    return kTRUE;
743 }
744
745 ///________________________________________________________________________
746 Bool_t AliConversionMesonCuts::SetNumberOfBGEvents(Int_t NumberOfBGEvents){
747    // Set Cut
748    switch(NumberOfBGEvents){
749    case 0:
750       fNumberOfBGEvents = 5;
751       break;
752    case 1:
753       fNumberOfBGEvents = 10;
754       break;
755    case 2:
756       fNumberOfBGEvents = 15;
757       break;
758    case 3:
759       fNumberOfBGEvents = 20;
760       break;
761    case 4:
762       fNumberOfBGEvents = 2;
763       break;
764    case 5:
765       fNumberOfBGEvents = 50;
766       break;
767    case 6:
768       fNumberOfBGEvents = 80;
769       break;
770    case 7:
771       fNumberOfBGEvents = 100;
772       break;
773    default:
774       cout<<"Warning: NumberOfBGEvents not defined "<<NumberOfBGEvents<<endl;
775       return kFALSE;
776    }
777    return kTRUE;
778 }
779 ///________________________________________________________________________
780 Bool_t AliConversionMesonCuts::SetSharedElectronCut(Int_t sharedElec) {
781
782    switch(sharedElec){
783    case 0:
784       fDoSharedElecCut = kFALSE;
785       break;
786    case 1:
787       fDoSharedElecCut = kTRUE;
788       break;
789    default:
790       cout<<"Warning: Shared Electron Cut not defined "<<sharedElec<<endl;
791       return kFALSE;
792    }
793         
794    return kTRUE;
795 }
796
797 ///________________________________________________________________________
798 Bool_t AliConversionMesonCuts::SetToCloseV0sCut(Int_t toClose) {
799
800    switch(toClose){
801    case 0:
802       fDoToCloseV0sCut = kFALSE;
803       fminV0Dist = 250;
804       break;
805    case 1:
806       fDoToCloseV0sCut = kTRUE;
807       fminV0Dist = 1;
808       break;
809    case 2:
810       fDoToCloseV0sCut = kTRUE;
811       fminV0Dist = 2;
812       break;
813    case 3:
814       fDoToCloseV0sCut = kTRUE;
815       fminV0Dist = 3;
816       break;
817    default:
818       cout<<"Warning: Shared Electron Cut not defined "<<toClose<<endl;
819       return kFALSE;
820    }
821    return kTRUE;
822 }
823
824 ///________________________________________________________________________
825 Bool_t AliConversionMesonCuts::SetMCPSmearing(Int_t useMCPSmearing)
826 {// Set Cut
827    switch(useMCPSmearing){
828    case 0:
829       fUseMCPSmearing=0;
830       fPBremSmearing=1.;
831       fPSigSmearing=0.;
832       fPSigSmearingCte=0.;
833       break;
834    case 1:
835       fUseMCPSmearing=1;
836       fPBremSmearing=1.0e-14;
837       fPSigSmearing=0.;
838       fPSigSmearingCte=0.;
839       break;
840    case 2:
841       fUseMCPSmearing=1;
842       fPBremSmearing=1.0e-15;
843       fPSigSmearing=0.0;
844       fPSigSmearingCte=0.;
845       break;
846    case 3:
847       fUseMCPSmearing=1;
848       fPBremSmearing=1.;
849       fPSigSmearing=0.003;
850       fPSigSmearingCte=0.002;
851       break;
852    case 4:
853       fUseMCPSmearing=1;
854       fPBremSmearing=1.;
855       fPSigSmearing=0.003;
856       fPSigSmearingCte=0.007;
857       break;
858    case 5:
859       fUseMCPSmearing=1;
860       fPBremSmearing=1.;
861       fPSigSmearing=0.003;
862       fPSigSmearingCte=0.016;
863       break;
864    case 6:
865       fUseMCPSmearing=1;
866       fPBremSmearing=1.;
867       fPSigSmearing=0.007;
868       fPSigSmearingCte=0.016;
869       break;
870    case 7:
871       fUseMCPSmearing=1;
872       fPBremSmearing=1.0e-16;
873       fPSigSmearing=0.0;
874       fPSigSmearingCte=0.;
875       break;
876    case 8:
877       fUseMCPSmearing=1;
878       fPBremSmearing=1.;
879       fPSigSmearing=0.007;
880       fPSigSmearingCte=0.014;
881       break;
882    case 9:
883       fUseMCPSmearing=1;
884       fPBremSmearing=1.;
885       fPSigSmearing=0.007;
886       fPSigSmearingCte=0.011;
887       break;
888
889    default:
890       AliError("Warning: UseMCPSmearing not defined");
891       return kFALSE;
892    }
893    return kTRUE;
894 }
895 ///________________________________________________________________________
896 TString AliConversionMesonCuts::GetCutNumber(){
897    // returns TString with current cut number
898    TString a(kNCuts);
899    for(Int_t ii=0;ii<kNCuts;ii++){
900       a.Append(Form("%d",fCuts[ii]));
901    }
902    return a;
903 }
904
905 ///________________________________________________________________________
906 void AliConversionMesonCuts::FillElectonLabelArray(AliAODConversionPhoton* photon, Int_t nV0){
907         
908    Int_t posLabel = photon->GetTrackLabelPositive();
909    Int_t negLabel = photon->GetTrackLabelNegative();
910         
911    fElectronLabelArray[nV0*2] = posLabel;
912    fElectronLabelArray[(nV0*2)+1] = negLabel;
913 }
914
915 ///________________________________________________________________________
916 Bool_t AliConversionMesonCuts::RejectSharedElectronV0s(AliAODConversionPhoton* photon, Int_t nV0, Int_t nV0s){
917
918    Int_t posLabel = photon->GetTrackLabelPositive();
919    Int_t negLabel = photon->GetTrackLabelNegative();
920
921    for(Int_t i = 0; i<nV0s*2;i++){
922       if(i==nV0*2)     continue;
923       if(i==(nV0*2)+1) continue;
924       if(fElectronLabelArray[i] == posLabel){
925          return kFALSE;}
926       if(fElectronLabelArray[i] == negLabel){
927          return kFALSE;}
928    }
929
930    return kTRUE;
931 }
932 ///________________________________________________________________________
933 Bool_t AliConversionMesonCuts::RejectToCloseV0s(AliAODConversionPhoton* photon, TList *photons, Int_t nV0){
934    Double_t posX = photon->GetConversionX();
935    Double_t posY = photon->GetConversionY();
936    Double_t posZ = photon->GetConversionZ();
937
938    for(Int_t i = 0;i<photons->GetEntries();i++){
939       if(nV0 == i) continue;
940       AliAODConversionPhoton *photonComp = (AliAODConversionPhoton*) photons->At(i);
941       Double_t posCompX = photonComp->GetConversionX();
942       Double_t posCompY = photonComp->GetConversionY();
943       Double_t posCompZ = photonComp->GetConversionZ();
944                 
945       Double_t dist = pow((posX - posCompX),2)+pow((posY - posCompY),2)+pow((posZ - posCompZ),2);
946
947       if(dist < fminV0Dist*fminV0Dist){
948          if(photon->GetChi2perNDF() < photonComp->GetChi2perNDF()) return kTRUE;
949          else {
950             return kFALSE;}
951       }
952                 
953    }
954    return kTRUE;
955 }
956
957 ///________________________________________________________________________
958 void AliConversionMesonCuts::SmearParticle(AliAODConversionPhoton* photon)
959 {
960    Double_t facPBrem = 1.;
961    Double_t facPSig = 0.;
962
963    Double_t phi=0.;
964    Double_t theta=0.;
965    Double_t P=0.;
966
967    
968    P=photon->P();
969    phi=photon->Phi();
970    if( photon->P()!=0){
971       theta=acos( photon->Pz()/ photon->P());
972    }
973
974    if( fPSigSmearing != 0. || fPSigSmearingCte!=0. ){ 
975       facPSig = TMath::Sqrt(fPSigSmearingCte*fPSigSmearingCte+fPSigSmearing*fPSigSmearing*P*P)*fRandom.Gaus(0.,1.);
976    }
977         
978    if( fPBremSmearing != 1.){
979       if(fBrem!=NULL){
980          facPBrem = fBrem->GetRandom();
981       }
982    }
983
984    photon->SetPx(facPBrem* (1+facPSig)* P*sin(theta)*cos(phi)) ;
985    photon->SetPy(facPBrem* (1+facPSig)* P*sin(theta)*sin(phi)) ;
986    photon->SetPz(facPBrem* (1+facPSig)* P*cos(theta)) ;
987    photon->SetE(photon->P());
988 }