[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),
  fEventPlane(0.),
  fMethodEP(0)
{
  //empty constructor
}
//________________________________________________
AliHFAfterBurner::AliHFAfterBurner(Int_t decChannel):
  fSigv2(0.1),
  fBkgv2(0.2),
  fUseNewton(kTRUE),
  fPrecisionNewton(0.0005),
  fDecChannel(decChannel),
  fSignal(0),
  fEventPlane(0.),
  fMethodEP(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),
  fEventPlane(source.fEventPlane),
  fMethodEP(source.fMethodEP)
{
  //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;
  fEventPlane=source.fEventPlane;
  fMethodEP=source.fMethodEP;
  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=0x0;
  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;
  }
  if(!pdgdaughters) return 0.;

  lab = d->MatchToMC(pdgmother,mcArray,nProngs,pdgdaughters);
  delete [] 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),
4609bdb6	45	fSignal(0),
	46	fEventPlane(0.),
	47	fMethodEP(0)
ed352b10	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),
4609bdb6	58	fSignal(0),
	59	fEventPlane(0.),
	60	fMethodEP(0)
ed352b10	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),
4609bdb6	72	fSignal(source.fSignal),
	73	fEventPlane(source.fEventPlane),
	74	fMethodEP(source.fMethodEP)
ed352b10	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;
4609bdb6	91	fEventPlane=source.fEventPlane;
4609bdb6	92	fMethodEP=source.fMethodEP;
ed352b10	93	return *this;
	94	}
	95	//______________________________________________________________________________
	96	AliHFAfterBurner::~AliHFAfterBurner(){
262b44fa	97	// destructor
ed352b10	98	}
	99	//______________________________________________________________________________
	100	Double_t AliHFAfterBurner::GetNewAngle(AliAODRecoDecayHF d,TClonesArray mcArray){
262b44fa	101	// modify the phi angle of teh tracks
ed352b10	102	Int_t lab=-1;
4609bdb6	103	Int_t *pdgdaughters=0x0;
ed352b10	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	}
4609bdb6	134	if(!pdgdaughters) return 0.;
4609bdb6	135
ed352b10	136	lab = d->MatchToMC(pdgmother,mcArray,nProngs,pdgdaughters);
374eb67b	137	delete [] pdgdaughters;
374eb67b	138
ed352b10	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){
262b44fa	175	// modifies the phi angle to after-burn the given v2
ed352b10	176	Double_t evplane = fEventPlane;
	177	if(fUseNewton){
	178	return NewtonMethodv2(phi,v2);
	179	}else{
	180	return phi-v2TMath::Sin(2(phi-evplane));
	181	}
	182	}
	183	//______________________________________________________________________________
	184	Double_t AliHFAfterBurner::NewtonMethodv2(Double_t phi,Double_t v2,Double_t phi0){
262b44fa	185	// modifies the phi angle to after-burn the given v2 using newton method
ed352b10	186	Double_t eventplane = fEventPlane;
	187	Double_t phi1 = phi-(phi+v2TMath::Sin(2.(phi-eventplane))-phi0)/(1.+2.v2TMath::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){
262b44fa	196	// SetMCv2
ed352b10	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){
262b44fa	227	// get phi angle in the range 0 - 2*pi
ed352b10	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){
262b44fa	239	// set the decay channel
ed352b10	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	}