PWGHF/vertexingHF/charmFlow/AliHFAfterBurner.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2010, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   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 // AliHFAfterBurner introduces a v2 modulation in MC for
  19 // the D mesons v2 analysis with event plane method
  20 // Authors: Giacomo Ortona, ortona@to.infn.it
  21 //
  22 /////////////////////////////////////////////////////////////
  23
  24 /* $Id$ */
  25
  26 #include <TDatabasePDG.h>
  27 #include <TVector3.h>
  28 #include <TRandom.h>
  29 #include <TRandom3.h>
  30 #include <AliAODEvent.h>
  31 #include <AliAODMCParticle.h>
  32 #include "AliAODRecoDecay.h"
  33 #include "AliAODRecoDecayHF.h"
  34 #include "AliHFAfterBurner.h"
  35
  36 ClassImp(AliHFAfterBurner)
  37
  38 //________________________________________________
  39 AliHFAfterBurner::AliHFAfterBurner():
  40   fSigv2(0),
  41   fBkgv2(0),
  42   fUseNewton(kTRUE),
  43   fPrecisionNewton(0),
  44   fDecChannel(0),
  45   fSignal(0),
  46   fEventPlane(0.),
  47   fMethodEP(0)
  48 {
  49   //empty constructor
  50 }
  51 //________________________________________________
  52 AliHFAfterBurner::AliHFAfterBurner(Int_t decChannel):
  53   fSigv2(0.1),
  54   fBkgv2(0.2),
  55   fUseNewton(kTRUE),
  56   fPrecisionNewton(0.0005),
  57   fDecChannel(decChannel),
  58   fSignal(0),
  59   fEventPlane(0.),
  60   fMethodEP(0)
  61 {
  62   //default constructor
  63 }
  64 //______________________________________________________________________________
  65 AliHFAfterBurner::AliHFAfterBurner(const AliHFAfterBurner &source):
  66   TObject(source),
  67   fSigv2(source.fSigv2),
  68   fBkgv2(source.fBkgv2),
  69   fUseNewton(source.fUseNewton),
  70   fPrecisionNewton(source.fPrecisionNewton),
  71   fDecChannel(source.fDecChannel),
  72   fSignal(source.fSignal),
  73   fEventPlane(source.fEventPlane),
  74   fMethodEP(source.fMethodEP)
  75 {
  76   //copy constructor
  77 }
  78 //______________________________________________________________________________
  79 AliHFAfterBurner &AliHFAfterBurner::operator=(const AliHFAfterBurner &source)
  80 {
  81   //assignment operator
  82   if(&source == this) return *this;
  83
  84   TObject::operator=(source);
  85   fSigv2 = source.fSigv2;
  86   fBkgv2=source.fBkgv2;
  87   fUseNewton=source.fUseNewton;
  88   fPrecisionNewton=source.fPrecisionNewton;
  89   fDecChannel=source.fDecChannel;
  90   fSignal=source.fSignal;
  91   fEventPlane=source.fEventPlane;
  92   fMethodEP=source.fMethodEP;
  93   return *this;
  94 }
  95 //______________________________________________________________________________
  96 AliHFAfterBurner::~AliHFAfterBurner(){
  97   // destructor
  98 }
  99 //______________________________________________________________________________
 100 Double_t AliHFAfterBurner::GetNewAngle(AliAODRecoDecayHF *d,TClonesArray *mcArray){
 101   // modify the phi angle of teh tracks
 102   Int_t lab=-1;
 103   Int_t *pdgdaughters=0x0;
 104   Int_t pdgmother=0;
 105   Int_t nProngs=0;
 106   switch(fDecChannel){
 107   case 0:
 108     //D+
 109     nProngs=3;
 110     pdgmother=411;
 111     pdgdaughters=new Int_t[3];
 112     pdgdaughters[0]=211;//pi
 113     pdgdaughters[1]=321;//K
 114     pdgdaughters[2]=211;//pi
 115     break;
 116   case 1:
 117     //D0
 118     nProngs=2;
 119     pdgmother=421;
 120     pdgdaughters=new Int_t[2];
 121     pdgdaughters[0]=211;//pi
 122     pdgdaughters[1]=321;//K
 123     break;
 124   case 2:
 125     //D*
 126     nProngs=3;
 127     pdgmother=413;
 128     pdgdaughters=new Int_t[3];
 129     pdgdaughters[1]=211;//pi
 130     pdgdaughters[0]=321;//K
 131     pdgdaughters[2]=211;//pi (soft?)
 132     break;
 133   }
 134   if(!pdgdaughters) return 0.;
 135
 136   lab = d->MatchToMC(pdgmother,mcArray,nProngs,pdgdaughters);
 137   delete [] pdgdaughters;
 138
 139   Double_t phi=-999.;
 140   if(lab>=0){
 141     fSignal=kTRUE;
 142     phi=GetPhi(d->Phi(),fSigv2);
 143   }
 144   else {//phi=NewtonMethodv2(phi,fBkgv2,eventplane);
 145     //background
 146     fSignal=kFALSE;
 147     //const Int_t nProngs=d->GetNProngs();
 148     Float_t phidau[nProngs];
 149     for(Int_t ipr=0;ipr<nProngs;ipr++){
 150       phidau[ipr]=(Float_t)d->PhiProng(ipr);
 151       //AliAODTrack *trk = (AliAODTrack*)d->GetDaughter(ipr);
 152       Int_t labdau=(Int_t)d->GetProngID(ipr);
 153       if(labdau<0)continue;
 154       AliAODMCParticle *mcpart= (AliAODMCParticle*)mcArray->At(labdau);
 155       if(!mcpart)continue;
 156       Int_t laborig=TMath::Abs(CheckOrigin(mcpart,mcArray));
 157       if(laborig>=0){//charm
 158         mcpart= (AliAODMCParticle*)mcArray->At(laborig);
 159         if(mcpart)phidau[ipr]=GetPhi(mcpart->Phi(),fSigv2);
 160       }else{//not charm
 161         phidau[ipr]=GetPhi(phidau[ipr],fBkgv2);
 162       }
 163     }
 164     Float_t py=0,px=0;
 165     for(Int_t ipr=0;ipr<nProngs;ipr++){
 166       py+=d->PtProng(ipr)*TMath::Sin(phidau[ipr]);
 167       px+=d->PtProng(ipr)*TMath::Cos(phidau[ipr]);
 168     }
 169     phi=TMath::Pi()+TMath::ATan2(-py,-px);
 170   }
 171   return GetPhi02Pi(phi);
 172 }
 173 //______________________________________________________________________________
 174 Double_t AliHFAfterBurner::GetPhi(Double_t phi,Float_t v2){
 175   // modifies the phi angle to after-burn the given v2
 176   Double_t evplane = fEventPlane;
 177   if(fUseNewton){
 178     return NewtonMethodv2(phi,v2);
 179   }else{
 180     return phi-v2*TMath::Sin(2*(phi-evplane));
 181   }
 182 }
 183 //______________________________________________________________________________
 184 Double_t AliHFAfterBurner::NewtonMethodv2(Double_t phi,Double_t v2,Double_t phi0){
 185   // modifies the phi angle to after-burn the given v2 using newton method
 186   Double_t eventplane = fEventPlane;
 187   Double_t phi1 = phi-(phi+v2*TMath::Sin(2.*(phi-eventplane))-phi0)/(1.+2.*v2*TMath::Cos(2.*(phi-eventplane)));
 188   if(TMath::Abs(phi/phi1-1.)<fPrecisionNewton){
 189     return phi1;
 190   }else {
 191     return NewtonMethodv2(phi1,v2,phi0);
 192   }
 193 }
 194 //______________________________________________________________________________
 195 void AliHFAfterBurner::SetMCv2(Float_t v2sig,Float_t v2bkg){
 196   // SetMCv2
 197   if(v2sig>=0)fSigv2=v2sig;
 198   if(v2bkg>=0)fBkgv2=v2bkg;
 199 }
 200 //______________________________________________________________________________
 201 Int_t AliHFAfterBurner::CheckOrigin(const AliAODMCParticle* mcPart,TClonesArray *arrayMC)const{
 202
 203   //
 204   // checking whether the mother of the particle is a charm
 205   //
 206   Float_t charmmother=kFALSE;
 207   Int_t pdgGranma = 0;
 208   Int_t mother = -1;
 209   mother = mcPart->GetMother();
 210   Int_t istep = 0;
 211   while (mother >=0 ){
 212     istep++;
 213     AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
 214     if(!mcGranma) break;
 215     pdgGranma = mcGranma->GetPdgCode();
 216     Int_t abspdgGranma = TMath::Abs(pdgGranma);
 217     if ((abspdgGranma > 400 && abspdgGranma < 500) || (abspdgGranma > 4000 && abspdgGranma < 5000)) {
 218       charmmother=kTRUE;
 219     }
 220     mother = mcGranma->GetMother();
 221   }
 222   if(!charmmother)mother=-1;
 223   return mother;
 224 }
 225 //________________________________________________________________________
 226 Float_t AliHFAfterBurner::GetPhi02Pi(Float_t phi){
 227   // get phi angle in the range 0 - 2*pi
 228   Float_t result=phi;
 229   while(result<0){
 230     result=result+2.*TMath::Pi();
 231   }
 232   while(result>TMath::Pi()*2.){
 233     result=result-2.*TMath::Pi();
 234   }
 235   return result;
 236 }
 237 //________________________________________________________________________
 238 void AliHFAfterBurner::SetDecChannel(Int_t decch){
 239   // set the decay channel
 240   if(decch>2){
 241     AliWarning("Invalid decay channel");
 242     return;
 243   }
 244   fDecChannel=decch;
 245 }
 246 //________________________________________________________________________
 247 void AliHFAfterBurner::SetEventPlaneMethod(Int_t method){
 248   //Only Random EP supported by now, feature to generate EP from histos to be added
 249   fMethodEP=method;
 250   return;
 251 }
 252 //________________________________________________________________________
 253 void AliHFAfterBurner::SetEventPlane(){
 254   //Only Random EP supported by now, feature to generate EP from histos to be added
 255   TRandom3 *g = new TRandom3(0);
 256   fEventPlane=g->Rndm()*TMath::Pi();
 257   delete g;g=0x0;
 258   return;
 259 }