#include "AliAODRecoDecayHF2Prong.h"
#include "AliAODRecoDecayHF4Prong.h"
#include "AliAODRecoCascadeHF.h"
-#include "AliAODHFUtil.h"
#include "AliAnalysisTaskSE.h"
#include "AliRDHFCutsDplustoKpipi.h"
fCentUpLimit(100),
fNMassBins(200),
fReadMC(kFALSE),
+ fUseAfterBurner(kFALSE),
fDecChannel(0),
+ fAfterBurner(0),
fUseV0EP(kFALSE),
fV0EPorder(2),
fMinCentr(10),
fCentUpLimit(100),
fNMassBins(200),
fReadMC(kFALSE),
+ fUseAfterBurner(kFALSE),
fDecChannel(decaychannel),
fUseV0EP(kFALSE),
fV0EPorder(2),
pdg=4122;
break;
}
+ fAfterBurner = new AliHFAfterBurner(fDecChannel);
for(Int_t i=0;i<6;i++)fHistvzero[i]=(TH2D*)histPar[i]->Clone();
for(Int_t i=0;i<6;i++)if(!fHistvzero[i])printf("No VZERO histograms!\n");
if(pdg==413) SetMassLimits((Float_t)0.135,(Float_t)0.165);
fHistvzero[i]=0x0;
}
}
-
+ if(fAfterBurner){
+ delete fAfterBurner;
+ fAfterBurner=0;
+ }
}
//_________________________________________________________________
if(fDebug > 1) printf("AnalysisTaskSEHFv2::UserCreateOutputObjects() \n");
- fhEventsInfo=new TH1F(GetOutputSlot(1)->GetContainer()->GetName(), "Number of AODs scanned",7,-0.5,6.5);
+ fhEventsInfo=new TH1F(GetOutputSlot(1)->GetContainer()->GetName(), "Number of AODs scanned",8,-0.5,7.5);
fhEventsInfo->GetXaxis()->SetBinLabel(1,"nEventsAnal");
fhEventsInfo->GetXaxis()->SetBinLabel(2,"nEvSelected");
fhEventsInfo->GetXaxis()->SetBinLabel(3,"nCandidatesSelected");
fhEventsInfo->GetXaxis()->SetBinLabel(5,"Pile-up Rej");
fhEventsInfo->GetXaxis()->SetBinLabel(6,"N. of 0SMH");
fhEventsInfo->GetXaxis()->SetBinLabel(7,Form("Ev Sel in Centr %.0f-%.0f%s",fRDCuts->GetMinCentrality(),fRDCuts->GetMaxCentrality(),"%"));
+ fhEventsInfo->GetXaxis()->SetBinLabel(8,"mismatch lab");
fhEventsInfo->GetXaxis()->SetNdivisions(1,kFALSE);
}
}
TString centrbinname=Form("centr%d_%d",icentr-5,icentr);
-
- if(fUseV0EP){
- rpangleevent=GetEventPlaneFromV0(aod);
+ if(fReadMC){
+ fUseV0EP=kTRUE;
+ TRandom3 *g = new TRandom3(0);
+ rpangleevent=g->Rndm()*TMath::Pi();
+ delete g;g=0x0;
eventplane=rpangleevent;
+ if(fUseAfterBurner)fAfterBurner->SetEventPlane((Double_t)rpangleevent);
}else{
- // event plane and resolution
- //--------------------------------------------------------------------------
- // extracting Q vectors and calculating v2 + resolution
- pl = aod->GetHeader()->GetEventplaneP();
- if(!pl){
- AliError("AliAnalysisTaskSEHFv2::UserExec:no eventplane! v2 analysis without eventplane not possible!\n");
- return;
+ if(fUseV0EP){
+ rpangleevent=GetEventPlaneFromV0(aod);
+ eventplane=rpangleevent;
+ }else{
+ // event plane and resolution
+ //--------------------------------------------------------------------------
+ // extracting Q vectors and calculating v2 + resolution
+ pl = aod->GetHeader()->GetEventplaneP();
+ if(!pl){
+ AliError("AliAnalysisTaskSEHFv2::UserExec:no eventplane! v2 analysis without eventplane not possible!\n");
+ return;
+ }
+ q = pl->GetQVector();
+ rpangleevent = pl->GetEventplane("Q"); // reaction plane angle without autocorrelations removal
+ Double_t deltaPsi = pl->GetQsubRes();
+ if(TMath::Abs(deltaPsi)>TMath::Pi()/2.){
+ if(deltaPsi>0.) deltaPsi-=TMath::Pi();
+ else deltaPsi +=TMath::Pi();
+ } // difference of subevents reaction plane angle cannot be bigger than phi/2
+ Double_t planereso = TMath::Cos(2.*deltaPsi); // reaction plane resolution
+ //--------------------------------------------------------------------------
+ ((TH1F*)fOutput->FindObject(Form("hEvPlaneReso%s",centrbinname.Data())))->Fill(planereso);
}
- q = pl->GetQVector();
- rpangleevent = pl->GetEventplane("Q"); // reaction plane angle without autocorrelations removal
- Double_t deltaPsi = pl->GetQsubRes();
- if(TMath::Abs(deltaPsi)>TMath::Pi()/2.){
- if(deltaPsi>0.) deltaPsi-=TMath::Pi();
- else deltaPsi +=TMath::Pi();
- } // difference of subevents reaction plane angle cannot be bigger than phi/2
- Double_t planereso = TMath::Cos(2.*deltaPsi); // reaction plane resolution
- //--------------------------------------------------------------------------
- ((TH1F*)fOutput->FindObject(Form("hEvPlaneReso%s",centrbinname.Data())))->Fill(planereso);
}
((TH1F*)fOutput->FindObject(Form("hEvPlane%s",centrbinname.Data())))->Fill(rpangleevent);
Float_t phi=d->Phi();
((TH1F*)fOutput->FindObject(Form("hPhi_pt%d%s",ptbin,centrbinname.Data())))->Fill(phi);
+ if(fReadMC&&fUseAfterBurner)phi=fAfterBurner->GetNewAngle(d,arrayMC);
Float_t deltaphi=GetPhi0Pi(phi-eventplane);
//fill the histograms with the appropriate method
if(fReadMC){
Int_t lab=-1;
- lab = d->MatchToMC(411,arrayMC,3,pdgdaughters);
+ if(fUseAfterBurner){
+ Bool_t isSignal=fAfterBurner->GetIsSignal();
+ if(isSignal)lab=10;
+ }else {
+ lab = d->MatchToMC(411,arrayMC,3,pdgdaughters);
+ }
if(lab>=0){ //signal
((TH1F*)fOutput->FindObject(Form("hSgn_pt%dphi%dcentr%d_%d",ptbin,phibin,icentr-5,icentr)))->Fill(masses[0]);
((TH2F*)fOutput->FindObject(Form("hMc2phiS_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
if(inte>0)flowTrack->SetWeight(flowTrack->Weight()/inte);
}
if(fDebug>15)printf("EPfromV0 flow tracks weights done\n");
-
+
AliFlowVector qvec=flowEvent.GetQ(fV0EPorder);
Double_t angleEP=(1./(Double_t)fV0EPorder)*qvec.Phi();
if(fDebug>15)printf("EPfromV0 phi %f\n",angleEP);
#include "AliAnalysisTaskSE.h"
#include "AliAnalysisVertexingHF.h"
+#include "AliHFAfterBurner.h"
class TH1F;
class TH2D;
void SetV0EventPlaneOrder(Int_t n){fV0EPorder=n;}
void SetMinCentrality(Int_t mincentr){fMinCentr=mincentr;}
void SetMaxCentrality(Int_t maxcentr){fMaxCentr=maxcentr;}
+ void SetUseAfterBurner(Bool_t ab){fUseAfterBurner=ab;}
+ void SetAfterBurner(AliHFAfterBurner *ab){fAfterBurner=ab;}
Float_t GetUpperMassLimit()const {return fUpmasslimit;}
Float_t GetLowerMassLimit()const {return fLowmasslimit;}
Float_t GetPhi0Pi(Float_t phi);
Float_t GetLowerCentLimit()const {return fCentLowLimit;}
Float_t GetUpperCentLimit()const {return fCentUpLimit;}
+ AliHFAfterBurner *GetAfterBurner()const {return fAfterBurner;}
// Implementation of interface methods
virtual void UserCreateOutputObjects();
virtual void LocalInit();// {Init();}
Float_t fCentUpLimit; //upper centrality limit
Int_t fNMassBins; //number of bins in the mass histograms
Bool_t fReadMC; //flag for access to MC
+ Bool_t fUseAfterBurner; //enable afterburning
Int_t fDecChannel; //decay channel identifier
+ AliHFAfterBurner *fAfterBurner;//Afterburner options
Bool_t fUseV0EP; //flag to select EP method
Int_t fV0EPorder; //harmonic for VZERO event plane
Int_t fMinCentr; //minimum centrality
--- /dev/null
+/**************************************************************************
+ * 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(){
+
+}
+//______________________________________________________________________________
+Double_t AliHFAfterBurner::GetNewAngle(AliAODRecoDecayHF *d,TClonesArray *mcArray){
+ 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){
+ 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){
+ 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){
+ 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){
+ 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){
+ 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;
+}
--- /dev/null
+#ifndef ALIANALYSISTASKAFTERBURNER_H
+#define ALIANALYSISTASKAFTERBURNER_H
+
+/* Copyright(c) 1998-2010, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+//*************************************************************************
+// AliAnalysisTaskSEHFv2 gives the needed tools for the D
+// mesons v2 analysis
+// Authors: Giacomo Ortona, ortona@to.infn.it,
+//
+//*************************************************************************
+
+/* $Id$ */
+
+class TH1F;
+class TH2D;
+class AliRDHFCuts;
+class TVector2;
+
+class AliHFAfterBurner : public TObject
+{
+ public:
+
+ AliHFAfterBurner();
+ AliHFAfterBurner(Int_t decChannel);
+
+ virtual ~AliHFAfterBurner();
+
+ Double_t GetNewAngle(AliAODRecoDecayHF *d,TClonesArray *arrayMC);
+ Double_t NewtonMethodv2(Double_t phi,Double_t v2,Double_t phi0);
+ Double_t NewtonMethodv2(Double_t phi,Double_t v2){
+ return NewtonMethodv2(phi,v2,phi);}
+ void SetMCv2(Float_t v2sig,Float_t v2bkg);
+ Int_t CheckOrigin(const AliAODMCParticle* mcPart,TClonesArray *arrayMC) const;
+ Double_t GetPhi(Double_t phi,Float_t v2);
+ Float_t GetPhi02Pi(Float_t phi);
+ void SetPrecision(Bool_t prec){fPrecisionNewton=prec;}
+ void SetUseNewton(Bool_t newt){fUseNewton=newt;}
+ void SetDecChannel(Int_t decch);
+ void SetEventPlane(Double_t ep){fEventPlane=ep;}
+ void SetEventPlane();
+ void SetEventPlaneMethod(Int_t method);
+
+ Bool_t GetIsSignal(){return fSignal;}
+ Double_t GetEventPlane(){return fEventPlane;}
+ private:
+ AliHFAfterBurner(const AliHFAfterBurner &source);
+ AliHFAfterBurner& operator=(const AliHFAfterBurner& source);
+
+ Float_t fSigv2; //v2 for signal
+ Float_t fBkgv2; //v2 for background
+ Bool_t fUseNewton; //Switch between Newton and fast methods
+ Double_t fPrecisionNewton;//precision to stop Newton method
+ Int_t fDecChannel; //D2H decay channel
+ Bool_t fSignal; //kTRUE if candidate is signal
+ Double_t fEventPlane; //event plane
+ Int_t fMethodEP; //switch between different EP method calculation
+
+ ClassDef(AliHFAfterBurner,1); // AliFlowAB class
+};
+#endif
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff