[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / charmFlow / AliHFAfterBurner.cxx

/**************************************************************************
 * Copyright(c) 1998-2010, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/////////////////////////////////////////////////////////////
//
// AliHFAfterBurner introduces a v2 modulation in MC for 
// the D mesons v2 analysis with event plane method
// Authors: Giacomo Ortona, ortona@to.infn.it
// 
/////////////////////////////////////////////////////////////

/* $Id$ */

#include <TDatabasePDG.h>
#include <TVector3.h>
#include <TRandom.h>
#include <TRandom3.h>
#include <AliAODEvent.h>
#include <AliAODMCParticle.h>
#include "AliAODRecoDecay.h"
#include "AliAODRecoDecayHF.h"
#include "AliHFAfterBurner.h"

ClassImp(AliHFAfterBurner)

//________________________________________________
AliHFAfterBurner::AliHFAfterBurner():
  fSigv2(0),
  fBkgv2(0),
  fUseNewton(kTRUE),
  fPrecisionNewton(0),
  fDecChannel(0),
  fSignal(0)
{
  //empty constructor
}
//________________________________________________
AliHFAfterBurner::AliHFAfterBurner(Int_t decChannel):
  fSigv2(0.1),
  fBkgv2(0.2),
  fUseNewton(kTRUE),
  fPrecisionNewton(0.0005),
  fDecChannel(decChannel),
  fSignal(0)
{
  //default constructor
}
//______________________________________________________________________________
AliHFAfterBurner::AliHFAfterBurner(const AliHFAfterBurner &source):
  TObject(source),
  fSigv2(source.fSigv2),
  fBkgv2(source.fBkgv2),
  fUseNewton(source.fUseNewton),
  fPrecisionNewton(source.fPrecisionNewton),
  fDecChannel(source.fDecChannel),
  fSignal(source.fSignal)
{
  //copy constructor
}
//______________________________________________________________________________
AliHFAfterBurner &AliHFAfterBurner::operator=(const AliHFAfterBurner &source)
{
  //assignment operator
  if(&source == this) return *this;

  TObject::operator=(source);
  fSigv2 = source.fSigv2;
  fBkgv2=source.fBkgv2;
  fUseNewton=source.fUseNewton;
  fPrecisionNewton=source.fPrecisionNewton;
  fDecChannel=source.fDecChannel;
  fSignal=source.fSignal;
  return *this;
}
//______________________________________________________________________________
AliHFAfterBurner::~AliHFAfterBurner(){
  // destructor
}
//______________________________________________________________________________
Double_t AliHFAfterBurner::GetNewAngle(AliAODRecoDecayHF *d,TClonesArray *mcArray){
  // modify the phi angle of teh tracks
  Int_t lab=-1;
  Int_t *pdgdaughters;
  Int_t pdgmother=0;
  Int_t nProngs=0;
  switch(fDecChannel){
  case 0:
    //D+
    nProngs=3;
    pdgmother=411;
    pdgdaughters=new Int_t[3];
    pdgdaughters[0]=211;//pi
    pdgdaughters[1]=321;//K
    pdgdaughters[2]=211;//pi
    break;
  case 1:
    //D0
    nProngs=2;
    pdgmother=421;
    pdgdaughters=new Int_t[2];
    pdgdaughters[0]=211;//pi 
    pdgdaughters[1]=321;//K
    break;
  case 2:
    //D*
    nProngs=3;
    pdgmother=413;
    pdgdaughters=new Int_t[3];
    pdgdaughters[1]=211;//pi
    pdgdaughters[0]=321;//K
    pdgdaughters[2]=211;//pi (soft?)
    break;
  }
  lab = d->MatchToMC(pdgmother,mcArray,nProngs,pdgdaughters);
  Double_t phi=-999.;
  if(lab>=0){
    fSignal=kTRUE;
    phi=GetPhi(d->Phi(),fSigv2);
  }
  else {//phi=NewtonMethodv2(phi,fBkgv2,eventplane);
    //background
    fSignal=kFALSE;
    //const Int_t nProngs=d->GetNProngs();
    Float_t phidau[nProngs];
    for(Int_t ipr=0;ipr<nProngs;ipr++){
      phidau[ipr]=(Float_t)d->PhiProng(ipr);
      //AliAODTrack *trk = (AliAODTrack*)d->GetDaughter(ipr);
      Int_t labdau=(Int_t)d->GetProngID(ipr);
      if(labdau<0)continue;
      AliAODMCParticle *mcpart= (AliAODMCParticle*)mcArray->At(labdau);
      if(!mcpart)continue;
      Int_t laborig=TMath::Abs(CheckOrigin(mcpart,mcArray));
      if(laborig>=0){//charm
	mcpart= (AliAODMCParticle*)mcArray->At(laborig);
	if(mcpart)phidau[ipr]=GetPhi(mcpart->Phi(),fSigv2);
      }else{//not charm
	phidau[ipr]=GetPhi(phidau[ipr],fBkgv2);
      }
    }
    Float_t py=0,px=0;
    for(Int_t ipr=0;ipr<nProngs;ipr++){
      py+=d->PtProng(ipr)*TMath::Sin(phidau[ipr]);
      px+=d->PtProng(ipr)*TMath::Cos(phidau[ipr]);
    }
    phi=TMath::Pi()+TMath::ATan2(-py,-px);
  }
  return GetPhi02Pi(phi);
}
//______________________________________________________________________________
Double_t AliHFAfterBurner::GetPhi(Double_t phi,Float_t v2){
  // modifies the phi angle to after-burn the given v2
  Double_t evplane = fEventPlane;
  if(fUseNewton){
    return NewtonMethodv2(phi,v2);
  }else{
    return phi-v2*TMath::Sin(2*(phi-evplane));
  }
}
//______________________________________________________________________________
Double_t AliHFAfterBurner::NewtonMethodv2(Double_t phi,Double_t v2,Double_t phi0){
  // modifies the phi angle to after-burn the given v2 using newton method
  Double_t eventplane = fEventPlane;
  Double_t phi1 = phi-(phi+v2*TMath::Sin(2.*(phi-eventplane))-phi0)/(1.+2.*v2*TMath::Cos(2.*(phi-eventplane)));
  if(TMath::Abs(phi/phi1-1.)<fPrecisionNewton){
    return phi1;
  }else {
    return NewtonMethodv2(phi1,v2,phi0);
  }
}
//______________________________________________________________________________
void AliHFAfterBurner::SetMCv2(Float_t v2sig,Float_t v2bkg){
  // SetMCv2
  if(v2sig>=0)fSigv2=v2sig;
  if(v2bkg>=0)fBkgv2=v2bkg;
}
//______________________________________________________________________________
Int_t AliHFAfterBurner::CheckOrigin(const AliAODMCParticle* mcPart,TClonesArray *arrayMC)const{
  
  //
  // checking whether the mother of the particle is a charm
  //
  Float_t charmmother=kFALSE;
  Int_t pdgGranma = 0;
  Int_t mother = -1;
  mother = mcPart->GetMother();
  Int_t istep = 0;
  while (mother >=0 ){
    istep++;
    AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
    if(!mcGranma) break;
    pdgGranma = mcGranma->GetPdgCode();
    Int_t abspdgGranma = TMath::Abs(pdgGranma);
    if ((abspdgGranma > 400 && abspdgGranma < 500) || (abspdgGranma > 4000 && abspdgGranma < 5000)) {
      charmmother=kTRUE;
    }
    mother = mcGranma->GetMother();
  }
  if(!charmmother)mother=-1;
  return mother;
}
//________________________________________________________________________
Float_t AliHFAfterBurner::GetPhi02Pi(Float_t phi){
  // get phi angle in the range 0 - 2*pi
  Float_t result=phi;
  while(result<0){
    result=result+2.*TMath::Pi();
  }
  while(result>TMath::Pi()*2.){
    result=result-2.*TMath::Pi();
  }
  return result;
}
//________________________________________________________________________
void AliHFAfterBurner::SetDecChannel(Int_t decch){
  // set the decay channel
  if(decch>2){
    AliWarning("Invalid decay channel");
    return;
  }
  fDecChannel=decch;
}
//________________________________________________________________________
void AliHFAfterBurner::SetEventPlaneMethod(Int_t method){
  //Only Random EP supported by now, feature to generate EP from histos to be added
  fMethodEP=method;
  return;
}
//________________________________________________________________________
void AliHFAfterBurner::SetEventPlane(){
  //Only Random EP supported by now, feature to generate EP from histos to be added
  TRandom3 *g = new TRandom3(0);
  fEventPlane=g->Rndm()*TMath::Pi();
  delete g;g=0x0;
  return;
}
Commit	Line	Data
ed352b10	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(){
262b44fa	89	// destructor
ed352b10	90	}
	91	//______________________________________________________________________________
	92	Double_t AliHFAfterBurner::GetNewAngle(AliAODRecoDecayHF d,TClonesArray mcArray){
262b44fa	93	// modify the phi angle of teh tracks
ed352b10	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){
262b44fa	163	// modifies the phi angle to after-burn the given v2
ed352b10	164	Double_t evplane = fEventPlane;
	165	if(fUseNewton){
	166	return NewtonMethodv2(phi,v2);
	167	}else{
	168	return phi-v2TMath::Sin(2(phi-evplane));
	169	}
	170	}
	171	//______________________________________________________________________________
	172	Double_t AliHFAfterBurner::NewtonMethodv2(Double_t phi,Double_t v2,Double_t phi0){
262b44fa	173	// modifies the phi angle to after-burn the given v2 using newton method
ed352b10	174	Double_t eventplane = fEventPlane;
	175	Double_t phi1 = phi-(phi+v2TMath::Sin(2.(phi-eventplane))-phi0)/(1.+2.v2TMath::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){
262b44fa	184	// SetMCv2
ed352b10	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){
262b44fa	215	// get phi angle in the range 0 - 2*pi
ed352b10	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){
262b44fa	227	// set the decay channel
ed352b10	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	}