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::MesonIsSelectedMCPiPlPiMiPiZero(TParticle *fMCMother,AliStack *fMCStack, Int_t &labelNegPion, Int_t &labelPosPion, Int_t &labelNeutPion, Double_t fRapidityShift){
 454
 455         // Returns true for all pions within acceptance cuts for decay into 2 photons
 456         // If bMCDaughtersInAcceptance is selected, it requires in addition that both daughter photons are within acceptance cuts
 457
 458         if( !fMCStack )return kFALSE;
 459
 460         if( !(fMCMother->GetPdgCode() == 221 || fMCMother->GetPdgCode() == 223) ) return kFALSE;
 461
 462         if( fMCMother->R()>fMaxR ) return kFALSE; // cuts on distance from collision point
 463
 464         Double_t rapidity = 10.;
 465
 466         if( fMCMother->Energy() - fMCMother->Pz() == 0 || fMCMother->Energy() + fMCMother->Pz() == 0 ){
 467                 rapidity=8.-fRapidityShift;
 468         }
 469         else{
 470                 rapidity = 0.5*(TMath::Log((fMCMother->Energy()+fMCMother->Pz()) / (fMCMother->Energy()-fMCMother->Pz())))-fRapidityShift;
 471         }
 472
 473         // Rapidity Cut
 474         if( abs(rapidity) > fRapidityCutMeson )return kFALSE;
 475
 476         // Select only -> pi+ pi- pi0
 477         if( fMCMother->GetNDaughters() != 3 )return kFALSE;
 478
 479         TParticle *posPion = 0x0;
 480         TParticle *negPion = 0x0;
 481         TParticle *neutPion = 0x0;
 482
 483         for(Int_t index= fMCMother->GetFirstDaughter();index<= fMCMother->GetLastDaughter();index++){
 484
 485                 TParticle* temp = (TParticle*)fMCStack->Particle( index );
 486                 switch( temp->GetPdgCode() ) {
 487                 case 211:
 488                         posPion      =  temp;
 489                         labelPosPion = index;
 490                         break;
 491                 case -211:
 492                         negPion      =  temp;
 493                         labelNegPion = index;
 494                         break;
 495                 case 111:
 496                         neutPion         =  temp;
 497                         labelNeutPion    = index;
 498                         break;
 499                 }
 500         }
 501
 502         if( posPion && negPion && neutPion ) return kTRUE;
 503         return kFALSE;
 504 }
 505
 506
 507 //________________________________________________________________________
 508 Bool_t AliConversionMesonCuts::MesonIsSelectedMCChiC(TParticle *fMCMother,AliStack *fMCStack,Int_t & labelelectronChiC, Int_t & labelpositronChiC, Int_t & labelgammaChiC, Double_t fRapidityShift){
 509    // Returns true for all ChiC within acceptance cuts for decay into JPsi + gamma -> e+ + e- + gamma
 510    // If bMCDaughtersInAcceptance is selected, it requires in addition that both daughter photons are within acceptance cuts
 511
 512    if(!fMCStack)return kFALSE;
 513    // if(fMCMother->GetPdgCode()==20443 ){
 514    //    return kFALSE;
 515    // }
 516    if(fMCMother->GetPdgCode()==10441 || fMCMother->GetPdgCode()==10443 || fMCMother->GetPdgCode()==445 ){
 517       if(fMCMother->R()>fMaxR)  return kFALSE; // cuts on distance from collision point
 518
 519       Double_t rapidity = 10.;
 520       if(fMCMother->Energy() - fMCMother->Pz() == 0 || fMCMother->Energy() + fMCMother->Pz() == 0){
 521          rapidity=8.-fRapidityShift;
 522       }
 523       else{
 524          rapidity = 0.5*(TMath::Log((fMCMother->Energy()+fMCMother->Pz()) / (fMCMother->Energy()-fMCMother->Pz())))-fRapidityShift;
 525       }
 526
 527       // Rapidity Cut
 528       if(abs(rapidity)>fRapidityCutMeson)return kFALSE;
 529
 530       // Select only -> ChiC radiative (JPsi+gamma) decay channel
 531       if(fMCMother->GetNDaughters()!=2)return kFALSE;
 532
 533       TParticle *jpsi   = 0x0;
 534       TParticle *gamma  = 0x0;
 535       TParticle *positron = 0x0;
 536       TParticle *electron = 0x0;
 537
 538       Int_t labeljpsiChiC = -1;
 539
 540       for(Int_t index= fMCMother->GetFirstDaughter();index<= fMCMother->GetLastDaughter();index++){
 541
 542          TParticle* temp = (TParticle*)fMCStack->Particle( index );
 543
 544          switch( temp->GetPdgCode() ) {
 545          case 443:
 546             jpsi =  temp;
 547             labeljpsiChiC = index;
 548             break;
 549          case 22:
 550             gamma    =  temp;
 551             labelgammaChiC = index;
 552             break;
 553          }
 554       }
 555
 556       if ( !jpsi || ! gamma) return kFALSE;
 557       if(jpsi->GetNDaughters()!=2)return kFALSE;
 558
 559
 560       for(Int_t index= jpsi->GetFirstDaughter();index<= jpsi->GetLastDaughter();index++){
 561          TParticle* temp = (TParticle*)fMCStack->Particle( index );
 562          switch( temp->GetPdgCode() ) {
 563          case -11:
 564             electron =  temp;
 565             labelelectronChiC = index;
 566             break;
 567          case 11:
 568             positron =  temp;
 569             labelpositronChiC = index;
 570             break;
 571          }
 572       }
 573       if( !electron || !positron) return kFALSE;
 574       if( positron && electron && gamma) return kTRUE;
 575    }
 576    return kFALSE;
 577 }
 578
 579 ///________________________________________________________________________
 580 Bool_t AliConversionMesonCuts::MesonIsSelected(AliAODConversionMother *pi0,Bool_t IsSignal, Double_t fRapidityShift)
 581 {
 582
 583    // Selection of reconstructed Meson candidates
 584    // Use flag IsSignal in order to fill Fill different
 585    // histograms for Signal and Background
 586    TH1 *hist=0x0;
 587
 588    if(IsSignal){hist=hMesonCuts;}
 589    else{hist=hMesonBGCuts;}
 590
 591    Int_t cutIndex=0;
 592
 593    if(hist)hist->Fill(cutIndex);
 594    cutIndex++;
 595
 596    // Undefined Rapidity -> Floating Point exception
 597    if((pi0->E()+pi0->Pz())/(pi0->E()-pi0->Pz())<=0){
 598       if(hist)hist->Fill(cutIndex);
 599       cutIndex++;
 600 //        cout << "undefined rapidity" << endl;
 601       return kFALSE;
 602    }
 603    else{
 604       // PseudoRapidity Cut --> But we cut on Rapidity !!!
 605       cutIndex++;
 606       if(abs(pi0->Rapidity()-fRapidityShift)>fRapidityCutMeson){
 607          if(hist)hist->Fill(cutIndex);
 608 //               cout << abs(pi0->Rapidity()-fRapidityShift) << ">" << fRapidityCutMeson << endl;
 609          return kFALSE;
 610       }
 611    }
 612    cutIndex++;
 613
 614    // Opening Angle Cut
 615    //fOpeningAngle=2*TMath::ATan(0.134/pi0->P());// physical minimum opening angle
 616    if( pi0->GetOpeningAngle() < fOpeningAngle){
 617 //        cout << pi0->GetOpeningAngle() << "<" << fOpeningAngle << endl;
 618       if(hist)hist->Fill(cutIndex);
 619       return kFALSE;
 620    }
 621    cutIndex++;
 622
 623    if ( fAlphaPtDepCut == kTRUE ) {
 624
 625         fAlphaCutMeson = fFAlphaCut->Eval( pi0->Pt() );
 626    }
 627
 628
 629    // Alpha Max Cut
 630    if(pi0->GetAlpha()>fAlphaCutMeson){
 631 //         cout << pi0->GetAlpha() << ">" << fAlphaCutMeson << endl;
 632       if(hist)hist->Fill(cutIndex);
 633       return kFALSE;
 634    }
 635    cutIndex++;
 636
 637    // Alpha Min Cut
 638    if(pi0->GetAlpha()<fAlphaMinCutMeson){
 639 //        cout << pi0->GetAlpha() << "<" << fAlphaMinCutMeson << endl;
 640       if(hist)hist->Fill(cutIndex);
 641       return kFALSE;
 642    }
 643    cutIndex++;
 644
 645    if (hDCAGammaGammaMesonBefore)hDCAGammaGammaMesonBefore->Fill(pi0->GetDCABetweenPhotons());
 646    if (hDCARMesonPrimVtxBefore)hDCARMesonPrimVtxBefore->Fill(pi0->GetDCARMotherPrimVtx());
 647
 648    if (pi0->GetDCABetweenPhotons() > fDCAGammaGammaCut){
 649 //        cout << pi0->GetDCABetweenPhotons() << ">" << fDCAGammaGammaCut << endl;
 650       if(hist)hist->Fill(cutIndex);
 651       return kFALSE;
 652    }
 653    cutIndex++;
 654
 655    if (pi0->GetDCARMotherPrimVtx() > fDCARMesonPrimVtxCut){
 656 //         cout << pi0->GetDCARMotherPrimVtx() << ">" << fDCARMesonPrimVtxCut << endl;
 657       if(hist)hist->Fill(cutIndex);
 658       return kFALSE;
 659    }
 660    cutIndex++;
 661
 662
 663    if (hDCAZMesonPrimVtxBefore)hDCAZMesonPrimVtxBefore->Fill(pi0->GetDCAZMotherPrimVtx());
 664
 665    if (abs(pi0->GetDCAZMotherPrimVtx()) > fDCAZMesonPrimVtxCut){
 666 //        cout << pi0->GetDCAZMotherPrimVtx() << ">" << fDCAZMesonPrimVtxCut << endl;
 667       if(hist)hist->Fill(cutIndex);
 668       return kFALSE;
 669    }
 670    cutIndex++;
 671
 672
 673    if (hDCAGammaGammaMesonAfter)hDCAGammaGammaMesonAfter->Fill(pi0->GetDCABetweenPhotons());
 674    if (hDCARMesonPrimVtxAfter)hDCARMesonPrimVtxAfter->Fill(pi0->GetDCARMotherPrimVtx());
 675    if (hDCAZMesonPrimVtxAfter)hDCAZMesonPrimVtxAfter->Fill(pi0->M(),pi0->GetDCAZMotherPrimVtx());
 676
 677    if(hist)hist->Fill(cutIndex);
 678    return kTRUE;
 679 }
 680
 681
 682
 683 ///________________________________________________________________________
 684 ///________________________________________________________________________
 685 Bool_t AliConversionMesonCuts::UpdateCutString() {
 686    ///Update the cut string (if it has been created yet)
 687
 688    if(fCutString && fCutString->GetString().Length() == kNCuts) {
 689       fCutString->SetString(GetCutNumber());
 690    } else {
 691       return kFALSE;
 692    }
 693    return kTRUE;
 694 }
 695
 696 ///________________________________________________________________________
 697 Bool_t AliConversionMesonCuts::InitializeCutsFromCutString(const TString analysisCutSelection ) {
 698    // Initialize Cuts from a given Cut string
 699    AliInfo(Form("Set Meson Cutnumber: %s",analysisCutSelection.Data()));
 700    if(analysisCutSelection.Length()!=kNCuts) {
 701       AliError(Form("Cut selection has the wrong length! size is %d, number of cuts is %d", analysisCutSelection.Length(), kNCuts));
 702       return kFALSE;
 703    }
 704    if(!analysisCutSelection.IsDigit()){
 705       AliError("Cut selection contains characters");
 706       return kFALSE;
 707    }
 708
 709    const char *cutSelection = analysisCutSelection.Data();
 710 #define ASSIGNARRAY(i)  fCuts[i] = cutSelection[i] - '0'
 711    for(Int_t ii=0;ii<kNCuts;ii++){
 712       ASSIGNARRAY(ii);
 713    }
 714
 715    // Set Individual Cuts
 716    for(Int_t ii=0;ii<kNCuts;ii++){
 717       if(!SetCut(cutIds(ii),fCuts[ii]))return kFALSE;
 718    }
 719
 720    PrintCutsWithValues();
 721    return kTRUE;
 722 }
 723 ///________________________________________________________________________
 724 Bool_t AliConversionMesonCuts::SetCut(cutIds cutID, const Int_t value) {
 725    ///Set individual cut ID
 726
 727    //cout << "Updating cut  " << fgkCutNames[cutID] << " (" << cutID << ") to " << value << endl;
 728    switch (cutID) {
 729    case kMesonKind:
 730       if( SetMesonKind(value)) {
 731          fCuts[kMesonKind] = value;
 732          UpdateCutString();
 733          return kTRUE;
 734       } else return kFALSE;
 735    case kchi2MesonCut:
 736       if( SetChi2MesonCut(value)) {
 737          fCuts[kchi2MesonCut] = value;
 738          UpdateCutString();
 739          return kTRUE;
 740       } else return kFALSE;
 741    case kalphaMesonCut:
 742       if( SetAlphaMesonCut(value)) {
 743          fCuts[kalphaMesonCut] = value;
 744          UpdateCutString();
 745          return kTRUE;
 746       } else return kFALSE;
 747    case kRCut:
 748       if( SetRCut(value)) {
 749          fCuts[kRCut] = value;
 750          UpdateCutString();
 751          return kTRUE;
 752       } else return kFALSE;
 753
 754    case kRapidityMesonCut:
 755       if( SetRapidityMesonCut(value)) {
 756          fCuts[kRapidityMesonCut] = value;
 757          UpdateCutString();
 758          return kTRUE;
 759       } else return kFALSE;
 760
 761    case kBackgroundScheme:
 762       if( SetBackgroundScheme(value)) {
 763          fCuts[kBackgroundScheme] = value;
 764          UpdateCutString();
 765          return kTRUE;
 766       } else return kFALSE;
 767
 768    case kDegreesForRotationMethod:
 769       if( SetNDegreesForRotationMethod(value)) {
 770          fCuts[kDegreesForRotationMethod] = value;
 771          UpdateCutString();
 772          return kTRUE;
 773       } else return kFALSE;
 774
 775    case kNumberOfBGEvents:
 776       if( SetNumberOfBGEvents(value)) {
 777          fCuts[kNumberOfBGEvents] = value;
 778          UpdateCutString();
 779          return kTRUE;
 780       } else return kFALSE;
 781
 782    case kuseMCPSmearing:
 783       if( SetMCPSmearing(value)) {
 784          fCuts[kuseMCPSmearing] = value;
 785          UpdateCutString();
 786          return kTRUE;
 787       } else return kFALSE;
 788    case kElecShare:
 789       if( SetSharedElectronCut(value)) {
 790          fCuts[kElecShare] = value;
 791          UpdateCutString();
 792          return kTRUE;
 793       } else return kFALSE;
 794    case kToCloseV0s:
 795       if( SetToCloseV0sCut(value)) {
 796          fCuts[kToCloseV0s] = value;
 797          UpdateCutString();
 798          return kTRUE;
 799       } else return kFALSE;
 800    case kDcaGammaGamma:
 801       if( SetDCAGammaGammaCut(value)) {
 802          fCuts[kDcaGammaGamma] = value;
 803          UpdateCutString();
 804          return kTRUE;
 805       } else return kFALSE;
 806    case kDcaZPrimVtx:
 807       if( SetDCAZMesonPrimVtxCut(value)) {
 808          fCuts[kDcaZPrimVtx] = value;
 809          UpdateCutString();
 810          return kTRUE;
 811       } else return kFALSE;
 812    case kDcaRPrimVtx:
 813       if( SetDCARMesonPrimVtxCut(value)) {
 814          fCuts[kDcaRPrimVtx] = value;
 815          UpdateCutString();
 816          return kTRUE;
 817       } else return kFALSE;
 818
 819    case kNCuts:
 820       cout << "Error:: Cut id out of range"<< endl;
 821       return kFALSE;
 822    }
 823
 824    cout << "Error:: Cut id " << cutID << " not recognized "<< endl;
 825    return kFALSE;
 826
 827 }
 828
 829
 830 ///________________________________________________________________________
 831 void AliConversionMesonCuts::PrintCuts() {
 832    // Print out current Cut Selection
 833    for(Int_t ic = 0; ic < kNCuts; ic++) {
 834       printf("%-30s : %d \n", fgkCutNames[ic], fCuts[ic]);
 835    }
 836 }
 837
 838 ///________________________________________________________________________
 839 void AliConversionMesonCuts::PrintCutsWithValues() {
 840    // Print out current Cut Selection with values
 841         printf("\nMeson cutnumber \n");
 842         for(Int_t ic = 0; ic < kNCuts; ic++) {
 843                 printf("%d",fCuts[ic]);
 844         }
 845         printf("\n\n");
 846
 847         printf("Meson cuts \n");
 848         printf("\t |y| < %3.2f \n", fRapidityCutMeson);
 849         printf("\t theta_{open} < %3.2f\n", fOpeningAngle);
 850         if (!fAlphaPtDepCut) printf("\t %3.2f < alpha < %3.2f\n", fAlphaMinCutMeson, fAlphaCutMeson);
 851         printf("\t dca_{gamma,gamma} > %3.2f\n", fDCAGammaGammaCut);
 852         printf("\t dca_{R, prim Vtx} > %3.2f\n", fDCARMesonPrimVtxCut);
 853         printf("\t dca_{Z, prim Vtx} > %3.2f\n\n", fDCAZMesonPrimVtxCut);
 854
 855          printf("Meson BG settings \n");
 856          if (!fDoBG){
 857                 if (!fUseRotationMethodInBG  & !fUseTrackMultiplicityForBG) printf("\t BG scheme: mixing V0 mult \n");
 858                 if (!fUseRotationMethodInBG  & fUseTrackMultiplicityForBG) printf("\t BG scheme: mixing track mult \n");
 859                 if (fUseRotationMethodInBG )printf("\t BG scheme: rotation \n");
 860                 if (fdoBGProbability) printf("\t -> use BG probability \n");
 861                 if (fBackgroundHandler) printf("\t -> use new BG handler \n");
 862                 printf("\t depth of pool: %d\n", fNumberOfBGEvents);
 863                 if (fUseRotationMethodInBG )printf("\t degree's for BG rotation: %d\n", fnDegreeRotationPMForBG);
 864          }
 865
 866 }
 867
 868
 869 ///________________________________________________________________________
 870 Bool_t AliConversionMesonCuts::SetMesonKind(Int_t mesonKind){
 871    // Set Cut
 872    switch(mesonKind){
 873    case 0:
 874       fMesonKind=0;
 875       break;
 876    case 1:
 877       fMesonKind=1;;
 878       break;
 879    default:
 880       cout<<"Warning: Meson kind not defined"<<mesonKind<<endl;
 881       return kFALSE;
 882    }
 883    return kTRUE;
 884 }
 885
 886 ///________________________________________________________________________
 887 Bool_t AliConversionMesonCuts::SetRCut(Int_t rCut){
 888    // Set Cut
 889    switch(rCut){
 890    case 0:
 891       fMaxR = 180.;
 892       break;
 893    case 1:
 894       fMaxR = 180.;
 895       break;
 896    case 2:
 897       fMaxR = 180.;
 898       break;
 899    case 3:
 900       fMaxR = 70.;
 901       break;
 902    case 4:
 903       fMaxR = 70.;
 904       break;
 905    case 5:
 906       fMaxR = 180.;
 907       break;
 908       // High purity cuts for PbPb
 909    case 9:
 910       fMaxR = 180.;
 911       break;
 912    default:
 913       cout<<"Warning: rCut not defined"<<rCut<<endl;
 914       return kFALSE;
 915    }
 916    return kTRUE;
 917 }
 918
 919 ///________________________________________________________________________
 920 Bool_t AliConversionMesonCuts::SetChi2MesonCut(Int_t chi2MesonCut){
 921    // Set Cut
 922    switch(chi2MesonCut){
 923    case 0:  // 100.
 924       fChi2CutMeson = 100.;
 925       break;
 926    case 1:  // 50.
 927       fChi2CutMeson = 50.;
 928       break;
 929    case 2:  // 30.
 930       fChi2CutMeson = 30.;
 931       break;
 932    case 3:
 933       fChi2CutMeson = 200.;
 934       break;
 935    case 4:
 936       fChi2CutMeson = 500.;
 937       break;
 938    case 5:
 939       fChi2CutMeson = 1000.;
 940       break;
 941    default:
 942       cout<<"Warning: Chi2MesonCut not defined "<<chi2MesonCut<<endl;
 943       return kFALSE;
 944    }
 945    return kTRUE;
 946 }
 947
 948
 949 ///________________________________________________________________________
 950 Bool_t AliConversionMesonCuts::SetAlphaMesonCut(Int_t alphaMesonCut)
 951 {   // Set Cut
 952    switch(alphaMesonCut){
 953    case 0:      // 0- 0.7
 954       fAlphaMinCutMeson  = 0.0;
 955       fAlphaCutMeson     = 0.7;
 956       fAlphaPtDepCut = kFALSE;
 957       break;
 958    case 1:      // Updated 31 October 2013 before 0.0 - 0.5
 959       if( fFAlphaCut ) delete fFAlphaCut;
 960       fFAlphaCut= new TF1("fFAlphaCut","[0]*tanh([1]*x)",0.,100.);
 961       fFAlphaCut->SetParameter(0,0.7);
 962       fFAlphaCut->SetParameter(1,1.2);
 963       fAlphaMinCutMeson  =  0.0;
 964       fAlphaCutMeson     = -1.0;
 965       fAlphaPtDepCut = kTRUE;
 966       break;
 967    case 2:      // Updated 31 October 2013 before 0.5-1
 968       if( fFAlphaCut ) delete fFAlphaCut;
 969       fFAlphaCut= new TF1("fFAlphaCut","[0]*tanh([1]*x)",0.,100.);
 970       fFAlphaCut->SetParameter(0,0.8);
 971       fFAlphaCut->SetParameter(1,1.2);
 972       fAlphaMinCutMeson  =  0.0;
 973       fAlphaCutMeson     = -1.0;
 974       fAlphaPtDepCut = kTRUE;
 975       break;
 976    case 3:      // 0.0-1
 977       fAlphaMinCutMeson  = 0.0;
 978       fAlphaCutMeson     = 1.;
 979       fAlphaPtDepCut = kFALSE;
 980       break;
 981    case 4:      // 0-0.65
 982       fAlphaMinCutMeson  = 0.0;
 983       fAlphaCutMeson     = 0.65;
 984       fAlphaPtDepCut = kFALSE;
 985       break;
 986    case 5:      // 0-0.75
 987       fAlphaMinCutMeson  = 0.0;
 988       fAlphaCutMeson     = 0.75;
 989       fAlphaPtDepCut = kFALSE;
 990       break;
 991    case 6:      // 0-0.8
 992       fAlphaMinCutMeson  = 0.0;
 993       fAlphaCutMeson     = 0.8;
 994       fAlphaPtDepCut = kFALSE;
 995       break;
 996    case 7:      // 0.0-0.85
 997       fAlphaMinCutMeson  = 0.0;
 998       fAlphaCutMeson     = 0.85;
 999       fAlphaPtDepCut = kFALSE;
1000       break;
1001    case 8:      // 0.0-0.6
1002       fAlphaMinCutMeson  = 0.0;
1003       fAlphaCutMeson     = 0.6;
1004       fAlphaPtDepCut = kFALSE;
1005       break;
1006    case 9: // Updated 11 November 2013 before 0.0 - 0.3
1007       if( fFAlphaCut ) delete fFAlphaCut;
1008       fFAlphaCut= new TF1("fFAlphaCut","[0]*tanh([1]*x)",0.,100.);
1009       fFAlphaCut->SetParameter(0,0.65);
1010       fFAlphaCut->SetParameter(1,1.2);
1011       fAlphaMinCutMeson  =  0.0;
1012       fAlphaCutMeson     = -1.0;
1013       fAlphaPtDepCut = kTRUE;
1014       break;
1015    default:
1016       cout<<"Warning: AlphaMesonCut not defined "<<alphaMesonCut<<endl;
1017       return kFALSE;
1018    }
1019    return kTRUE;
1020 }
1021
1022 ///________________________________________________________________________
1023 Bool_t AliConversionMesonCuts::SetRapidityMesonCut(Int_t RapidityMesonCut){
1024    // Set Cut
1025    switch(RapidityMesonCut){
1026    case 0:  // changed from 0.9 to 1.35
1027       fRapidityCutMeson   = 1.35;
1028       break;
1029    case 1:  //
1030       fRapidityCutMeson   = 0.8;
1031       break;
1032    case 2:  //
1033       fRapidityCutMeson   = 0.7;
1034       break;
1035    case 3:  //
1036       fRapidityCutMeson   = 0.6;
1037       break;
1038    case 4:  //
1039       fRapidityCutMeson   = 0.5;
1040       break;
1041    case 5:  //
1042       fRapidityCutMeson   = 0.85;
1043       break;
1044    case 6:  //
1045       fRapidityCutMeson   = 0.75;
1046       break;
1047    case 7:  //
1048       fRapidityCutMeson   = 0.3;
1049       break;
1050    case 8:  //changed, before 0.35
1051       fRapidityCutMeson   = 0.25;
1052       break;
1053    case 9:  //
1054       fRapidityCutMeson   = 0.4;
1055       break;
1056    default:
1057       cout<<"Warning: RapidityMesonCut not defined "<<RapidityMesonCut<<endl;
1058       return kFALSE;
1059    }
1060    return kTRUE;
1061 }
1062
1063
1064 ///________________________________________________________________________
1065 Bool_t AliConversionMesonCuts::SetBackgroundScheme(Int_t BackgroundScheme){
1066    // Set Cut
1067    switch(BackgroundScheme){
1068    case 0: //Rotation
1069       fUseRotationMethodInBG=kTRUE;
1070       fdoBGProbability=kFALSE;
1071       break;
1072    case 1: // mixed event with V0 multiplicity
1073       fUseRotationMethodInBG=kFALSE;
1074       fUseTrackMultiplicityForBG=kFALSE;
1075       fdoBGProbability=kFALSE;
1076       break;
1077    case 2: // mixed event with track multiplicity
1078       fUseRotationMethodInBG=kFALSE;
1079       fUseTrackMultiplicityForBG=kTRUE;
1080       fdoBGProbability=kFALSE;
1081       break;
1082    case 3: //Rotation
1083       fUseRotationMethodInBG=kTRUE;
1084       fdoBGProbability=kTRUE;
1085       break;
1086    case 4: //No BG calculation
1087       cout << "no BG calculation should be done" << endl;
1088       fUseRotationMethodInBG=kFALSE;
1089       fdoBGProbability=kFALSE;
1090       fDoBG=kFALSE;
1091       break;
1092    case 5: //Rotation
1093       fUseRotationMethodInBG=kTRUE;
1094       fdoBGProbability=kFALSE;
1095       fBackgroundHandler = 1;
1096       break;
1097    case 6: // mixed event with V0 multiplicity
1098       fUseRotationMethodInBG=kFALSE;
1099       fUseTrackMultiplicityForBG=kFALSE;
1100       fdoBGProbability=kFALSE;
1101       fBackgroundHandler = 1;
1102       break;
1103    case 7: // mixed event with track multiplicity
1104       fUseRotationMethodInBG=kFALSE;
1105       fUseTrackMultiplicityForBG=kTRUE;
1106       fdoBGProbability=kFALSE;
1107       fBackgroundHandler = 1;
1108       break;
1109    case 8: //Rotation
1110       fUseRotationMethodInBG=kTRUE;
1111       fdoBGProbability=kTRUE;
1112       fBackgroundHandler = 1;
1113       break;
1114    default:
1115       cout<<"Warning: BackgroundScheme not defined "<<BackgroundScheme<<endl;
1116       return kFALSE;
1117    }
1118    return kTRUE;
1119 }
1120
1121
1122 ///________________________________________________________________________
1123 Bool_t AliConversionMesonCuts::SetNDegreesForRotationMethod(Int_t DegreesForRotationMethod){
1124    // Set Cut
1125    switch(DegreesForRotationMethod){
1126    case 0:
1127       fnDegreeRotationPMForBG = 5;
1128       break;
1129    case 1:
1130       fnDegreeRotationPMForBG = 10;
1131       break;
1132    case 2:
1133       fnDegreeRotationPMForBG = 15;
1134       break;
1135    case 3:
1136       fnDegreeRotationPMForBG = 20;
1137       break;
1138    default:
1139       cout<<"Warning: DegreesForRotationMethod not defined "<<DegreesForRotationMethod<<endl;
1140       return kFALSE;
1141    }
1142    fCuts[kDegreesForRotationMethod]=DegreesForRotationMethod;
1143    return kTRUE;
1144 }
1145
1146 ///________________________________________________________________________
1147 Bool_t AliConversionMesonCuts::SetNumberOfBGEvents(Int_t NumberOfBGEvents){
1148    // Set Cut
1149    switch(NumberOfBGEvents){
1150    case 0:
1151       fNumberOfBGEvents = 5;
1152       break;
1153    case 1:
1154       fNumberOfBGEvents = 10;
1155       break;
1156    case 2:
1157       fNumberOfBGEvents = 15;
1158       break;
1159    case 3:
1160       fNumberOfBGEvents = 20;
1161       break;
1162    case 4:
1163       fNumberOfBGEvents = 2;
1164       break;
1165    case 5:
1166       fNumberOfBGEvents = 50;
1167       break;
1168    case 6:
1169       fNumberOfBGEvents = 80;
1170       break;
1171    case 7:
1172       fNumberOfBGEvents = 100;
1173       break;
1174    default:
1175       cout<<"Warning: NumberOfBGEvents not defined "<<NumberOfBGEvents<<endl;
1176       return kFALSE;
1177    }
1178    return kTRUE;
1179 }
1180 ///________________________________________________________________________
1181 Bool_t AliConversionMesonCuts::SetSharedElectronCut(Int_t sharedElec) {
1182
1183    switch(sharedElec){
1184    case 0:
1185       fDoSharedElecCut = kFALSE;
1186       break;
1187    case 1:
1188       fDoSharedElecCut = kTRUE;
1189       break;
1190    default:
1191       cout<<"Warning: Shared Electron Cut not defined "<<sharedElec<<endl;
1192       return kFALSE;
1193    }
1194
1195    return kTRUE;
1196 }
1197
1198 ///________________________________________________________________________
1199 Bool_t AliConversionMesonCuts::SetToCloseV0sCut(Int_t toClose) {
1200
1201    switch(toClose){
1202    case 0:
1203       fDoToCloseV0sCut = kFALSE;
1204       fminV0Dist = 250;
1205       break;
1206    case 1:
1207       fDoToCloseV0sCut = kTRUE;
1208       fminV0Dist = 1;
1209       break;
1210    case 2:
1211       fDoToCloseV0sCut = kTRUE;
1212       fminV0Dist = 2;
1213       break;
1214    case 3:
1215       fDoToCloseV0sCut = kTRUE;
1216       fminV0Dist = 3;
1217       break;
1218    default:
1219       cout<<"Warning: Shared Electron Cut not defined "<<toClose<<endl;
1220       return kFALSE;
1221    }
1222    return kTRUE;
1223 }
1224
1225 ///________________________________________________________________________
1226 Bool_t AliConversionMesonCuts::SetMCPSmearing(Int_t useMCPSmearing)
1227 {// Set Cut
1228    switch(useMCPSmearing){
1229    case 0:
1230       fUseMCPSmearing=0;
1231       fPBremSmearing=1.;
1232       fPSigSmearing=0.;
1233       fPSigSmearingCte=0.;
1234       break;
1235    case 1:
1236       fUseMCPSmearing=1;
1237       fPBremSmearing=1.0e-14;
1238       fPSigSmearing=0.;
1239       fPSigSmearingCte=0.;
1240       break;
1241    case 2:
1242       fUseMCPSmearing=1;
1243       fPBremSmearing=1.0e-15;
1244       fPSigSmearing=0.0;
1245       fPSigSmearingCte=0.;
1246       break;
1247    case 3:
1248       fUseMCPSmearing=1;
1249       fPBremSmearing=1.;
1250       fPSigSmearing=0.003;
1251       fPSigSmearingCte=0.002;
1252       break;
1253    case 4:
1254       fUseMCPSmearing=1;
1255       fPBremSmearing=1.;
1256       fPSigSmearing=0.003;
1257       fPSigSmearingCte=0.007;
1258       break;
1259    case 5:
1260       fUseMCPSmearing=1;
1261       fPBremSmearing=1.;
1262       fPSigSmearing=0.003;
1263       fPSigSmearingCte=0.016;
1264       break;
1265    case 6:
1266       fUseMCPSmearing=1;
1267       fPBremSmearing=1.;
1268       fPSigSmearing=0.007;
1269       fPSigSmearingCte=0.016;
1270       break;
1271    case 7:
1272       fUseMCPSmearing=1;
1273       fPBremSmearing=1.0e-16;
1274       fPSigSmearing=0.0;
1275       fPSigSmearingCte=0.;
1276       break;
1277    case 8:
1278       fUseMCPSmearing=1;
1279       fPBremSmearing=1.;
1280       fPSigSmearing=0.007;
1281       fPSigSmearingCte=0.014;
1282       break;
1283    case 9:
1284       fUseMCPSmearing=1;
1285       fPBremSmearing=1.;
1286       fPSigSmearing=0.007;
1287       fPSigSmearingCte=0.011;
1288       break;
1289
1290    default:
1291       AliError("Warning: UseMCPSmearing not defined");
1292       return kFALSE;
1293    }
1294    return kTRUE;
1295 }
1296
1297
1298 ///________________________________________________________________________
1299 Bool_t AliConversionMesonCuts::SetDCAGammaGammaCut(Int_t DCAGammaGamma){
1300    // Set Cut
1301    switch(DCAGammaGamma){
1302    case 0:  //
1303       fDCAGammaGammaCut   = 1000;
1304       break;
1305    case 1:  //
1306       fDCAGammaGammaCut   = 10;
1307       break;
1308    case 2:  //
1309       fDCAGammaGammaCut   = 5;
1310       break;
1311    case 3:  //
1312       fDCAGammaGammaCut   = 4;
1313       break;
1314    case 4:  //
1315       fDCAGammaGammaCut   = 3;
1316       break;
1317    case 5:  //
1318       fDCAGammaGammaCut   = 2.5;
1319       break;
1320    case 6:  //
1321       fDCAGammaGammaCut   = 2;
1322       break;
1323    case 7:  //
1324       fDCAGammaGammaCut   = 1.5;
1325       break;
1326    case 8:  //
1327       fDCAGammaGammaCut   = 1;
1328       break;
1329    case 9:  //
1330       fDCAGammaGammaCut   = 0.5;
1331       break;
1332    default:
1333       cout<<"Warning: DCAGammaGamma not defined "<<DCAGammaGamma<<endl;
1334       return kFALSE;
1335    }
1336    return kTRUE;
1337 }
1338
1339 ///________________________________________________________________________
1340 Bool_t AliConversionMesonCuts::SetDCAZMesonPrimVtxCut(Int_t DCAZMesonPrimVtx){
1341    // Set Cut
1342    switch(DCAZMesonPrimVtx){
1343    case 0:  //
1344       fDCAZMesonPrimVtxCut   = 1000;
1345       break;
1346    case 1:  //
1347       fDCAZMesonPrimVtxCut   = 10;
1348       break;
1349    case 2:  //
1350       fDCAZMesonPrimVtxCut   = 5;
1351       break;
1352    case 3:  //
1353       fDCAZMesonPrimVtxCut   = 4;
1354       break;
1355    case 4:  //
1356       fDCAZMesonPrimVtxCut   = 3;
1357       break;
1358    case 5:  //
1359       fDCAZMesonPrimVtxCut   = 2.5;
1360       break;
1361    case 6:  //
1362       fDCAZMesonPrimVtxCut   = 2;
1363       break;
1364    case 7:  //
1365       fDCAZMesonPrimVtxCut   = 1.5;
1366       break;
1367    case 8:  //
1368       fDCAZMesonPrimVtxCut   = 1;
1369       break;
1370    case 9:  //
1371       fDCAZMesonPrimVtxCut   = 0.5;
1372       break;
1373    default:
1374       cout<<"Warning: DCAZMesonPrimVtx not defined "<<DCAZMesonPrimVtx<<endl;
1375       return kFALSE;
1376    }
1377    return kTRUE;
1378 }
1379
1380 ///________________________________________________________________________
1381 Bool_t AliConversionMesonCuts::SetDCARMesonPrimVtxCut(Int_t DCARMesonPrimVtx){
1382    // Set Cut
1383    switch(DCARMesonPrimVtx){
1384    case 0:  //
1385       fDCARMesonPrimVtxCut   = 1000;
1386       break;
1387    case 1:  //
1388       fDCARMesonPrimVtxCut   = 10;
1389       break;
1390    case 2:  //
1391       fDCARMesonPrimVtxCut   = 5;
1392       break;
1393    case 3:  //
1394       fDCARMesonPrimVtxCut   = 4;
1395       break;
1396    case 4:  //
1397       fDCARMesonPrimVtxCut   = 3;
1398       break;
1399    case 5:  //
1400       fDCARMesonPrimVtxCut   = 2.5;
1401       break;
1402    case 6:  //
1403       fDCARMesonPrimVtxCut   = 2;
1404       break;
1405    case 7:  //
1406       fDCARMesonPrimVtxCut   = 1.5;
1407       break;
1408    case 8:  //
1409       fDCARMesonPrimVtxCut   = 1;
1410       break;
1411    case 9:  //
1412       fDCARMesonPrimVtxCut   = 0.5;
1413       break;
1414    default:
1415       cout<<"Warning: DCARMesonPrimVtx not defined "<<DCARMesonPrimVtx<<endl;
1416       return kFALSE;
1417    }
1418    return kTRUE;
1419 }
1420
1421
1422 ///________________________________________________________________________
1423 TString AliConversionMesonCuts::GetCutNumber(){
1424    // returns TString with current cut number
1425    TString a(kNCuts);
1426    for(Int_t ii=0;ii<kNCuts;ii++){
1427       a.Append(Form("%d",fCuts[ii]));
1428    }
1429    return a;
1430 }
1431
1432 ///________________________________________________________________________
1433 void AliConversionMesonCuts::FillElectonLabelArray(AliAODConversionPhoton* photon, Int_t nV0){
1434
1435    Int_t posLabel = photon->GetTrackLabelPositive();
1436    Int_t negLabel = photon->GetTrackLabelNegative();
1437
1438    fElectronLabelArray[nV0*2] = posLabel;
1439    fElectronLabelArray[(nV0*2)+1] = negLabel;
1440 }
1441
1442 ///________________________________________________________________________
1443 Bool_t AliConversionMesonCuts::RejectSharedElectronV0s(AliAODConversionPhoton* photon, Int_t nV0, Int_t nV0s){
1444
1445    Int_t posLabel = photon->GetTrackLabelPositive();
1446    Int_t negLabel = photon->GetTrackLabelNegative();
1447
1448    for(Int_t i = 0; i<nV0s*2;i++){
1449       if(i==nV0*2)     continue;
1450       if(i==(nV0*2)+1) continue;
1451       if(fElectronLabelArray[i] == posLabel){
1452          return kFALSE;}
1453       if(fElectronLabelArray[i] == negLabel){
1454          return kFALSE;}
1455    }
1456
1457    return kTRUE;
1458 }
1459 ///________________________________________________________________________
1460 Bool_t AliConversionMesonCuts::RejectToCloseV0s(AliAODConversionPhoton* photon, TList *photons, Int_t nV0){
1461    Double_t posX = photon->GetConversionX();
1462    Double_t posY = photon->GetConversionY();
1463    Double_t posZ = photon->GetConversionZ();
1464
1465    for(Int_t i = 0;i<photons->GetEntries();i++){
1466       if(nV0 == i) continue;
1467       AliAODConversionPhoton *photonComp = (AliAODConversionPhoton*) photons->At(i);
1468       Double_t posCompX = photonComp->GetConversionX();
1469       Double_t posCompY = photonComp->GetConversionY();
1470       Double_t posCompZ = photonComp->GetConversionZ();
1471
1472       Double_t dist = pow((posX - posCompX),2)+pow((posY - posCompY),2)+pow((posZ - posCompZ),2);
1473
1474       if(dist < fminV0Dist*fminV0Dist){
1475          if(photon->GetChi2perNDF() < photonComp->GetChi2perNDF()) return kTRUE;
1476          else {
1477             return kFALSE;}
1478       }
1479
1480    }
1481    return kTRUE;
1482 }
1483
1484 ///________________________________________________________________________
1485 void AliConversionMesonCuts::SmearParticle(AliAODConversionPhoton* photon)
1486 {
1487
1488    if (photon==NULL) return;
1489    Double_t facPBrem = 1.;
1490    Double_t facPSig = 0.;
1491
1492    Double_t phi=0.;
1493    Double_t theta=0.;
1494    Double_t P=0.;
1495
1496
1497    P=photon->P();
1498    phi=photon->Phi();
1499    if( photon->P()!=0){
1500       theta=acos( photon->Pz()/ photon->P());
1501    }
1502
1503    if( fPSigSmearing != 0. || fPSigSmearingCte!=0. ){
1504       facPSig = TMath::Sqrt(fPSigSmearingCte*fPSigSmearingCte+fPSigSmearing*fPSigSmearing*P*P)*fRandom.Gaus(0.,1.);
1505    }
1506
1507    if( fPBremSmearing != 1.){
1508       if(fBrem!=NULL){
1509          facPBrem = fBrem->GetRandom();
1510       }
1511    }
1512
1513    photon->SetPx(facPBrem* (1+facPSig)* P*sin(theta)*cos(phi)) ;
1514    photon->SetPy(facPBrem* (1+facPSig)* P*sin(theta)*sin(phi)) ;
1515    photon->SetPz(facPBrem* (1+facPSig)* P*cos(theta)) ;
1516    photon->SetE(photon->P());
1517 }
1518 ///________________________________________________________________________
1519 void AliConversionMesonCuts::SmearVirtualPhoton(AliAODConversionPhoton* photon)
1520 {
1521
1522    if (photon==NULL) return;
1523    Double_t facPBrem = 1.;
1524    Double_t facPSig = 0.;
1525
1526    Double_t phi=0.;
1527    Double_t theta=0.;
1528    Double_t P=0.;
1529
1530
1531    P=photon->P();
1532    phi=photon->Phi();
1533    if( photon->P()!=0){
1534       theta=acos( photon->Pz()/ photon->P());
1535    }
1536
1537    if( fPSigSmearing != 0. || fPSigSmearingCte!=0. ){
1538       facPSig = TMath::Sqrt(fPSigSmearingCte*fPSigSmearingCte+fPSigSmearing*fPSigSmearing*P*P)*fRandom.Gaus(0.,1.);
1539    }
1540
1541    if( fPBremSmearing != 1.){
1542       if(fBrem!=NULL){
1543          facPBrem = fBrem->GetRandom();
1544       }
1545    }
1546
1547    photon->SetPx(facPBrem* (1+facPSig)* P*sin(theta)*cos(phi)) ;
1548    photon->SetPy(facPBrem* (1+facPSig)* P*sin(theta)*sin(phi)) ;
1549    photon->SetPz(facPBrem* (1+facPSig)* P*cos(theta)) ;
1550
1551 }
1552 ///________________________________________________________________________
1553 TLorentzVector AliConversionMesonCuts::SmearElectron(TLorentzVector particle)
1554 {
1555
1556    //if (particle==0) return;
1557    Double_t facPBrem = 1.;
1558    Double_t facPSig = 0.;
1559
1560    Double_t phi=0.;
1561    Double_t theta=0.;
1562    Double_t P=0.;
1563
1564
1565    P=particle.P();
1566
1567
1568    phi=particle.Phi();
1569    if (phi < 0.) phi += 2. * TMath::Pi();
1570
1571    if( particle.P()!=0){
1572       theta=acos( particle.Pz()/ particle.P());
1573    }
1574
1575
1576    Double_t fPSigSmearingHalf    =  fPSigSmearing  / 2.0;  //The parameter was set for gammas with 2 particles and here we have just one electron
1577    Double_t sqrtfPSigSmearingCteHalf =  fPSigSmearingCte / 2.0 ;  //The parameter was set for gammas with 2 particles and here we have just one electron
1578
1579
1580
1581    if( fPSigSmearingHalf != 0. || sqrtfPSigSmearingCteHalf!=0. ){
1582       facPSig = TMath::Sqrt(sqrtfPSigSmearingCteHalf*sqrtfPSigSmearingCteHalf+fPSigSmearingHalf*fPSigSmearingHalf*P*P)*fRandom.Gaus(0.,1.);
1583    }
1584
1585    if( fPBremSmearing != 1.){
1586       if(fBrem!=NULL){
1587          facPBrem = fBrem->GetRandom();
1588       }
1589    }
1590
1591    TLorentzVector SmearedParticle;
1592
1593    SmearedParticle.SetXYZM( facPBrem* (1+facPSig)* P*sin(theta)*cos(phi) , facPBrem* (1+facPSig)* P*sin(theta)*sin(phi)  ,
1594                           facPBrem* (1+facPSig)* P*cos(theta) , TDatabasePDG::Instance()->GetParticle(  ::kElectron   )->Mass()) ;
1595
1596    //particle.SetPx(facPBrem* (1+facPSig)* P*sin(theta)*cos(phi)) ;
1597    //particle.SetPy(facPBrem* (1+facPSig)* P*sin(theta)*sin(phi)) ;
1598    //particle.SetPz(facPBrem* (1+facPSig)* P*cos(theta)) ;
1599
1600    return SmearedParticle;
1601
1602 }
1603