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