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