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