* about the suitability of this software for any purpose. It is *
* provided "as is" without express or implied warranty. *
**************************************************************************/
-
-/* $Id$ */
-
-//
-//
-// Base class for DStar in Jets Analysis
-//
-// The D* (+ and -) is reconstructed inside jets. Two different cuts are
-// implemented:
-//
-// 1) C.Ivan D* cuts adapted for correlation analysis
-// 2) Topological cut enforcing the correlation D0-softPion pt + relaxed
-// CosThetaStar. This second should be better for correlations.
-//
-// USAGE:
//
-// The analysis is performed separately for D*+ and D*-. A flag in the .C is
-// taking care to decide which analysis.
//
-// The cut number 2 can be activaded with a flag in the .C (not active in this version)
-// Cuts 1) are the default.
-//
-// The task requires reconstructed jets in the AODs
+// Base class for DStar in Jets Analysis
//
//-----------------------------------------------------------------------
// Author A.Grelli
#include <TDatabasePDG.h>
#include <TParticle.h>
#include <TVector3.h>
+#include "TROOT.h"
#include "AliAnalysisTaskSEDStarJets.h"
+#include "AliRDHFCutsDStartoKpipi.h"
#include "AliStack.h"
#include "AliMCEvent.h"
#include "AliAODMCHeader.h"
-#include "AliAnalysisManager.h"
#include "AliAODHandler.h"
+#include "AliAnalysisManager.h"
#include "AliLog.h"
#include "AliAODVertex.h"
#include "AliAODJet.h"
#include "AliAODRecoDecay.h"
-#include "AliAODRecoDecayHF.h"
+#include "AliAODRecoCascadeHF.h"
#include "AliAODRecoDecayHF2Prong.h"
#include "AliESDtrack.h"
#include "AliAODMCParticle.h"
//__________________________________________________________________________
AliAnalysisTaskSEDStarJets::AliAnalysisTaskSEDStarJets() :
AliAnalysisTaskSE(),
- fCountReco(0),
- fCountRecoAcc(0),
- fCountRecoITSClusters(0),
- fCountRecoPPR(0),
- fCountDStar(0),
fEvents(0),
- fMinITSClusters(0),
- fComputeD0(kTRUE),
- fUseMCInfo(kTRUE),
- ftopologicalCut(kFALSE),
+ fchargeFrCorr(0),
+ fUseMCInfo(kTRUE),
fRequireNormalization(kTRUE),
- fLorentzTrack1(0,0,0,0),
- fLorentzTrack2(0,0,0,0),
- fLorentzTrack3(0,0,0,0),
- fLorentzTrack4(0,0,0,0),
fOutput(0),
- fD0ptvsSoftPtSignal(0),
- fD0ptvsSoftPt(0),
+ fCuts(0),
ftrigger(0),
fPtPion(0),
fInvMass(0),
- fRECOPtDStar(0),
+ fRECOPtDStar(0),
+ fRECOPtBkg(0),
fDStar(0),
fDiff(0),
fDiffSideBand(0),
fPhiBkg(0),
fTrueDiff(0),
fResZ(0),
- fResZBkg(0),
- fcharmpt(0),
- fdstarE(0),
+ fResZBkg(0),
fEjet(0),
fPhijet(0),
fEtaJet(0),
- fdstarpt(0)
+ theMCFF(0),
+ fDphiD0Dstar(0),
+ fPtJet(0)
{
//
// Default ctor
//
}
//___________________________________________________________________________
-AliAnalysisTaskSEDStarJets::AliAnalysisTaskSEDStarJets(const Char_t* name) :
+AliAnalysisTaskSEDStarJets::AliAnalysisTaskSEDStarJets(const Char_t* name, AliRDHFCutsDStartoKpipi* cuts) :
AliAnalysisTaskSE(name),
- fCountReco(0),
- fCountRecoAcc(0),
- fCountRecoITSClusters(0),
- fCountRecoPPR(0),
- fCountDStar(0),
fEvents(0),
- fMinITSClusters(0),
- fComputeD0(kTRUE),
+ fchargeFrCorr(0),
fUseMCInfo(kTRUE),
- ftopologicalCut(kFALSE),
fRequireNormalization(kTRUE),
- fLorentzTrack1(0,0,0,0),
- fLorentzTrack2(0,0,0,0),
- fLorentzTrack3(0,0,0,0),
- fLorentzTrack4(0,0,0,0),
fOutput(0),
- fD0ptvsSoftPtSignal(0),
- fD0ptvsSoftPt(0),
+ fCuts(0),
ftrigger(0),
fPtPion(0),
fInvMass(0),
- fRECOPtDStar(0),
+ fRECOPtDStar(0),
+ fRECOPtBkg(0),
fDStar(0),
fDiff(0),
fDiffSideBand(0),
fTrueDiff(0),
fResZ(0),
fResZBkg(0),
- fcharmpt(0),
- fdstarE(0),
fEjet(0),
fPhijet(0),
fEtaJet(0),
- fdstarpt(0)
+ theMCFF(0),
+ fDphiD0Dstar(0),
+ fPtJet(0)
{
//
// Constructor. Initialization of Inputs and Outputs
//
+ fCuts=cuts;
Info("AliAnalysisTaskSEDStarJets","Calling Constructor");
- DefineOutput(1,TList::Class());
-}
-
-//___________________________________________________________________________
-AliAnalysisTaskSEDStarJets& AliAnalysisTaskSEDStarJets::operator=(const AliAnalysisTaskSEDStarJets& c)
-{
- //
- // Assignment operator
- //
- if (this!=&c) {
- AliAnalysisTaskSE::operator=(c) ;
- }
-
- return *this;
+ DefineOutput(1,TList::Class()); // histos
+ DefineOutput(2,AliRDHFCutsDStartoKpipi::Class()); // my cuts
}
-
//___________________________________________________________________________
-AliAnalysisTaskSEDStarJets::AliAnalysisTaskSEDStarJets(const AliAnalysisTaskSEDStarJets& c) :
- AliAnalysisTaskSE(c),
- fCountReco(c.fCountReco),
- fCountRecoAcc(c.fCountRecoAcc),
- fCountRecoITSClusters(c.fCountRecoITSClusters),
- fCountRecoPPR(c.fCountRecoPPR),
- fCountDStar(c.fCountDStar),
- fEvents(c.fEvents),
- fMinITSClusters(c.fMinITSClusters),
- fComputeD0(c.fComputeD0),
- fUseMCInfo(c.fUseMCInfo),
- ftopologicalCut(c.ftopologicalCut),
- fRequireNormalization(c.fRequireNormalization),
- fLorentzTrack1(c.fLorentzTrack1),
- fLorentzTrack2(c.fLorentzTrack2),
- fLorentzTrack3(c.fLorentzTrack3),
- fLorentzTrack4(c.fLorentzTrack4),
- fOutput(c.fOutput),
- fD0ptvsSoftPtSignal(c.fD0ptvsSoftPtSignal),
- fD0ptvsSoftPt(c.fD0ptvsSoftPt),
- ftrigger(c.ftrigger),
- fPtPion(c.fPtPion),
- fInvMass(c.fInvMass),
- fRECOPtDStar(c.fRECOPtDStar),
- fDStar(c.fDStar),
- fDiff(c.fDiff),
- fDiffSideBand(c.fDiffSideBand),
- fDStarMass(c.fDStarMass),
- fPhi(c.fPhi),
- fPhiBkg(c.fPhiBkg),
- fTrueDiff(c.fTrueDiff),
- fResZ(c.fResZ),
- fResZBkg(c.fResZBkg),
- fcharmpt(c.fcharmpt),
- fdstarE(c.fdstarE),
- fEjet(c.fEjet),
- fPhijet(c.fPhijet),
- fEtaJet(c.fEtaJet),
- fdstarpt(c.fdstarpt)
-
-{
+AliAnalysisTaskSEDStarJets::~AliAnalysisTaskSEDStarJets() {
//
- // Copy Constructor
+ // destructor
//
-}
-//___________________________________________________________________________
-AliAnalysisTaskSEDStarJets::~AliAnalysisTaskSEDStarJets() {
- // destructor
-
Info("~AliAnalysisTaskSEDStarJets","Calling Destructor");
if (fOutput) {
delete fOutput;
fOutput = 0;
- }
+ }
+
+ if (fCuts) {
+ delete fCuts;
+ fCuts = 0;
+ }
+
+ if (ftrigger) { delete ftrigger; ftrigger = 0;}
+ if (fPtPion) { delete fPtPion; fPtPion = 0;}
+ if (fInvMass) { delete fInvMass; fInvMass = 0;}
+ if (fRECOPtDStar) { delete fRECOPtDStar; fRECOPtDStar = 0;}
+ if (fRECOPtBkg) { delete fRECOPtBkg; fRECOPtBkg = 0;}
+ if (fDStar) { delete fDStar; fDStar = 0;}
+ if (fDiff) { delete fDiff; fDiff = 0;}
+ if (fDiffSideBand) { delete fDiffSideBand; fDiffSideBand = 0;}
+ if (fDStarMass) { delete fDStarMass; fDStarMass = 0;}
+ if (fPhi) { delete fPhi; fPhi = 0;}
+ if (fPhiBkg) { delete fPhiBkg; fPhiBkg = 0;}
+ if (fTrueDiff){ delete fTrueDiff; fTrueDiff = 0;}
+ if (fResZ) { delete fResZ; fResZ = 0;}
+ if (fResZBkg) { delete fResZBkg; fResZBkg = 0;}
+ if (fEjet) { delete fEjet; fEjet = 0;}
+ if (fPhijet) { delete fPhijet; fPhijet = 0;}
+ if (fEtaJet) { delete fEtaJet; fEtaJet = 0;}
+ if (theMCFF) { delete theMCFF; theMCFF = 0;}
+ if (fDphiD0Dstar) { delete fDphiD0Dstar; fDphiD0Dstar = 0;}
+ if (fPtJet) { delete fPtJet; fPtJet = 0;}
+}
+
+//___________________________________________________________
+void AliAnalysisTaskSEDStarJets::Init(){
+ //
+ // Initialization
+ //
+ if(fDebug > 1) printf("AnalysisTaskSEDStarJets::Init() \n");
+ AliRDHFCutsDStartoKpipi* copyfCuts=new AliRDHFCutsDStartoKpipi(*fCuts);
+ // Post the cuts
+ PostData(2,copyfCuts);
+
+ return;
}
//_________________________________________________
Error("UserExec","NO EVENT FOUND!");
return;
}
-
+
+ fEvents++;
+ AliInfo(Form("Event %d",fEvents));
+ if (fEvents%10000 ==0) AliInfo(Form("Event %d",fEvents));
+
// Load the event
AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
-
- TClonesArray *arrayVerticesHF=0;
-
+
+ TClonesArray *arrayDStartoD0pi=0;
+
if(!aodEvent && AODEvent() && IsStandardAOD()) {
// In case there is an AOD handler writing a standard AOD, use the AOD
// event in memory rather than the input (ESD) event.
if(aodHandler->GetExtensions()) {
AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
AliAODEvent *aodFromExt = ext->GetAOD();
- arrayVerticesHF=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
+ arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
}
} else {
- arrayVerticesHF=(TClonesArray*)aodEvent->GetList()->FindObject("D0toKpi");
+ arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar");
}
- if (!arrayVerticesHF){
+ if (!arrayDStartoD0pi){
AliInfo("Could not find array of HF vertices, skipping the event");
return;
- }else AliDebug(2, Form("Found %d vertices",arrayVerticesHF->GetEntriesFast()));
+ }else AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast()));
// fix for temporary bug in ESDfilter
// the AODs with null vertex pointer didn't pass the PhysSel
if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;
- fEvents++;
- AliDebug(2,Form("Event %d",fEvents));
- if (fEvents%10000 ==0) AliDebug(2,Form("Event %d",fEvents));
-
-
// Simulate a jet triggered sample
TClonesArray *arrayofJets = (TClonesArray*)aodEvent->GetJets();
if(aodEvent->GetNJets()<=0) return;
- AliInfo("found a jet: processing D* in jet analysis");
-
- // AOD primary vertex
- AliAODVertex *prVtx = (AliAODVertex*)aodEvent->GetPrimaryVertex();
- Double_t primaryPos[3];
- prVtx->GetXYZ(primaryPos);
- Bool_t vtxFlag = kTRUE;
- TString title=prVtx->GetTitle();
- if(!title.Contains("VertexerTracks")) vtxFlag=kFALSE;
// counters for efficiencies
Int_t icountReco = 0;
- Int_t icountRecoAcc = 0;
- Int_t icountRecoITSClusters = 0;
- Int_t icountRecoPPR = 0;
- Int_t fiDstar = 0;
- Int_t fDStarD0 = 0;
-
- // TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
-
- TClonesArray *mcArray = 0;
-
- if(fUseMCInfo){
- mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
- //loop on the MC event - some basic MC info on D*, D0 and soft pion
- if(!mcArray) {
- printf("AliAnalysisTaskSEDStarSpectra::UserExec: MC particles branch not found!\n");
- return;
- }
-
- for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) {
- AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
- if (!mcPart) {
- AliWarning("Particle not found in tree, skipping");
- continue;
- }
-
- // charm pt
- if(TMath::Abs(mcPart->GetPdgCode())==4){
- fcharmpt->Fill(mcPart->Pt());
- }
-
- // fill energy and pt for D* in acceptance with correct prongs
- Bool_t isOk = DstarInMC(mcPart,mcArray);
-
- if (isOk){ //D*
- AliDebug(2, "Found a DStar in MC with correct prongs and in acceptance");
- fdstarE ->Fill(mcPart->E());
- fdstarpt->Fill(mcPart->Pt());
-
- // check the MC-Acceptance level cuts
- // since standard CF functions are not applicable, using Kine Cuts on daughters
-
- Int_t daughter0 = mcPart->GetDaughter(0);
- Int_t daughter1 = mcPart->GetDaughter(1);
-
- AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughter0,daughter1));
- AliAODMCParticle* mcPartDaughter0 = 0;
- AliAODMCParticle* mcPartDaughter1 = 0;
-
- mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter0));
- if(!mcPartDaughter0) {
- printf("AliAnalysisTaskSEDStarSpectra::UserExec: MC particle daugter 0 not found!\n");
- return;
- }
- mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter1));
- if(!mcPartDaughter1) {
- printf("AliAnalysisTaskSEDStarSpectra::UserExec: MC particle daugter 1 not found!\n");
- return;
- }
-
- Double_t eta0 = mcPartDaughter0->Eta();
- Double_t eta1 = mcPartDaughter1->Eta();
- Double_t y0 = mcPartDaughter0->Y();
- Double_t y1 = mcPartDaughter1->Y();
- Double_t pt0 = mcPartDaughter0->Pt();
- Double_t pt1 = mcPartDaughter1->Pt();
-
- AliDebug(2, Form("D* Daughter 0: eta = %f, y = %f, pt = %f", eta0, y0, pt0));
- AliDebug(2, Form("D* Daughter 1: eta = %f, y = %f, pt = %f", eta1, y1, pt1));
-
- Int_t daughD00 = 0;
- Int_t daughD01 = 0;
-
- // D0 daughters - do not need to check D0-kpi, already done
-
- if(TMath::Abs(mcPartDaughter0->GetPdgCode())==421){
- daughD00 = mcPartDaughter0->GetDaughter(0);
- daughD01 = mcPartDaughter0->GetDaughter(1);
- }else{
- daughD00 = mcPartDaughter1->GetDaughter(0);
- daughD01 = mcPartDaughter1->GetDaughter(1);
- }
-
- AliAODMCParticle* mcD0PartDaughter0 = 0;
- AliAODMCParticle* mcD0PartDaughter1 = 0;
- mcD0PartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughD00));
- mcD0PartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughD01));
-
- if (!mcD0PartDaughter0 || !mcD0PartDaughter1) {
- AliWarning("At least one Daughter Particle not found in tree, but it should be, this check was already done...");
- return;
- }
-
- // D0 daughters - needed for acceptance
- Double_t pD0pt0 = mcD0PartDaughter0->Pt();
- Double_t pD0pt1 = mcD0PartDaughter1->Pt();
- Double_t pD0eta0 = mcD0PartDaughter0->Eta();
- Double_t pD0eta1 = mcD0PartDaughter1->Eta();
-
- // ACCEPTANCE REQUESTS ---------
-
- // soft pion
- Bool_t daught1inAcceptance = (TMath::Abs(eta1) <= 0.9 && pt1 > 0.05);
- // Do daughters
- Bool_t theD0daught0inAcceptance = (TMath::Abs(pD0eta0) <= 0.9 && pD0pt0 >= 0.1);
- Bool_t theD0daught1inAcceptance = (TMath::Abs(pD0eta1) <= 0.9 && pD0pt1 >= 0.1);
-
- if (daught1inAcceptance && theD0daught0inAcceptance && theD0daught1inAcceptance) {
-
- AliDebug(2, "Daughter particles in acceptance");
-
- // check on the vertex
- if (vtxFlag){
- printf("Vertex cut passed 2\n");
- fDStarD0++;
- Bool_t refitFlag = kTRUE;
- for (Int_t iaod =0; iaod<aodEvent->GetNumberOfTracks(); iaod++){
- AliAODTrack *track = (AliAODTrack*)aodEvent->GetTrack(iaod);
-
- // refit only for D0 daughters
- if ((track->GetLabel() == daughD00) || (track->GetLabel() == daughD01)) {
- if(!(track->GetStatus()&AliESDtrack::kTPCrefit) || !(track->GetStatus()&AliESDtrack::kITSrefit)) {
- refitFlag = kFALSE;
- }
- }
- }
- if (refitFlag){
- printf("Refit cut passed\n");
- }
- else{
- AliDebug(3,"Refit cut not passed\n");
- }
- }
- else{
- AliDebug(3,"Vertex cut not passed\n");
- }
- }
- else{
- AliDebug(3,"Acceptance cut not passed\n");
- }
- }
- }
-
- AliDebug(2, Form("Found %i MC particles that are D* in Kpipi and satisfy Acc cuts!!",fDStarD0));
- fCountDStar += fDStarD0;
-
- } // End of MC
-
- // Now perform the D* in jet reconstruction
-
// Normalization factor
if(fRequireNormalization){
ftrigger->Fill(1);
}
+
+ //D* and D0 prongs needed to MatchToMC method
+ Int_t pdgDgDStartoD0pi[2]={421,211};
+ Int_t pdgDgD0toKpi[2]={321,211};
+
+ Double_t max =0;
+ Double_t ejet = 0;
+ Double_t phiJet = 0;
+ Double_t etaJet = 0;
+ Double_t ptjet = 0;
- Int_t nJets = 0; // for reco D0 countings
-
- for (Int_t iJets = 0; iJets<arrayofJets->GetEntriesFast(); iJets++) { // main loop on jets
-
+ //loop over jets and consider only the leading jey in the event
+ for (Int_t iJets = 0; iJets<arrayofJets->GetEntriesFast(); iJets++) {
AliAODJet* jet = (AliAODJet*)arrayofJets->At(iJets);
-
+
//jets variables
- Double_t ejet = jet->E();
- Double_t phiJet = jet->Phi();
- Double_t etaJet = jet->Eta();
+ ejet = jet->E();
+ if(ejet>max){
+ max = jet->E();
+ phiJet = jet->Phi();
+ etaJet = jet->Eta();
+ ptjet = jet->Pt();
+ }
+
// fill energy, eta and phi of the jet
fEjet ->Fill(ejet);
fPhijet ->Fill(phiJet);
fEtaJet ->Fill(etaJet);
+ fPtJet->Fill(ptjet);
+ }
+
+ //loop over D* candidates
+ for (Int_t iDStartoD0pi = 0; iDStartoD0pi<arrayDStartoD0pi->GetEntriesFast(); iDStartoD0pi++) {
- nJets++;
+ // D* candidates
+ AliAODRecoCascadeHF* dstarD0pi = (AliAODRecoCascadeHF*)arrayDStartoD0pi->At(iDStartoD0pi);
+ AliAODRecoDecayHF2Prong* theD0particle = (AliAODRecoDecayHF2Prong*)dstarD0pi->Get2Prong();
+
+ Double_t finvM =-999;
+ Double_t finvMDStar =-999;
+ Double_t dPhi =-999;
+ Bool_t isDStar =0;
+ Int_t mcLabel = -9999;
+
+ // find associated MC particle for D* ->D0toKpi
+ if(fUseMCInfo){
+ TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
+ if(!mcArray) {
+ printf("AliAnalysisTaskSEDStarSpectra::UserExec: MC particles not found!\n");
+ return;
+ }
+ mcLabel = dstarD0pi->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,mcArray);
+
+ if(mcLabel>=0) isDStar = 1;
+ if(mcLabel>0){
+ Double_t zMC =-999;
+ AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcLabel));
+ //fragmentation function in mc
+ zMC= FillMCFF(mcPart,mcArray,mcLabel);
+ if(zMC>0) theMCFF->Fill(zMC);
+ }
+ }
- // loop over the tracks to search for candidates soft pions
- for (Int_t i=0; i<aodEvent->GetNTracks(); i++){
-
- AliAODTrack* aodTrack = aodEvent->GetTrack(i);
- Double_t vD0[4] = {0.,0.,0.,0.};
- Double_t vD0bar[4] ={0.,0.,0.,0.};
+ // soft pion
+ AliAODTrack *track2 = (AliAODTrack*)dstarD0pi->GetBachelor();
+ Double_t pt = dstarD0pi->Pt();
- Int_t pCharge = aodTrack->Charge();
+ // track quality cuts
+ Int_t isTkSelected = fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kTracks); // quality cuts on tracks
+ if(!isTkSelected) continue;
- // few selections on soft pion
- Bool_t tPrimCand = aodTrack->IsPrimaryCandidate(); // is it primary?
+ // region of interest + topological cuts + PID
+ if(!fCuts->IsInFiducialAcceptance(dstarD0pi->Pt(),dstarD0pi->YDstar())) continue;
+ Int_t isSelected=fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kCandidate); //selected
+ if (!isSelected) continue;
+
+ // fill histos
+ finvM = dstarD0pi->InvMassD0();
+ fInvMass->Fill(finvM);
- if(aodTrack->Pt()>= 5 || aodTrack->Pt()< 0.08 ) continue; //cut on soft pion pt VERY large
- //~ D*s of pt >80GeV with a soft pion of 5GeV!
- if(TMath::Abs(aodTrack->Eta())>0.9) continue;
+ //DStar invariant mass
+ finvMDStar = dstarD0pi->InvMassDstarKpipi();
+
+ Double_t EGjet = 0;
+ Double_t dStarMom = dstarD0pi->P();
+ Double_t phiDStar = dstarD0pi->Phi();
+ Double_t phiD0 = theD0particle->Phi();
+ //check suggested by Federico
+ Double_t dPhiD0Dstar = phiD0 - phiDStar;
+
+ dPhi = phiJet - phiDStar;
- // if D*+ analysis tha D0 and pi+ otherwise pi-
- if(fComputeD0 && pCharge!= 1 ) continue;
- if(!fComputeD0 && pCharge!= -1 ) continue;
+ //plot right dphi
+ if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi());
+ if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi());
+
+ Double_t corrFactorCharge = (ejet/100)*20;
+ EGjet = ejet + corrFactorCharge;
+
+ // fill D* candidates
+ Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
+ if(finvM >= (mPDGD0-0.05) && finvM <=(mPDGD0+0.05)){ // ~3 sigma (sigma=17MeV, conservative)
+
+ if(isDStar == 1) {
+ fDphiD0Dstar->Fill(dPhiD0Dstar);
+ fDStarMass->Fill(finvMDStar);
+ fTrueDiff->Fill(finvMDStar-finvM);
+ }
+ if(isDStar == 0) fDphiD0Dstar->Fill(dPhiD0Dstar); // angle between D0 and D*
+
+ fDStar->Fill(finvMDStar);
+ fDiff ->Fill(finvMDStar-finvM);
+
+ Double_t mPDGDstar=TDatabasePDG::Instance()->GetParticle(413)->Mass();
+ Double_t invmassDelta = dstarD0pi->DeltaInvMass();
- if(tPrimCand && arrayVerticesHF->GetEntriesFast()>0){ // isPion and is Primary, no PID for now
+ // now the dphi signal and the fragmentation function
+ if (TMath::Abs(invmassDelta-(mPDGDstar-mPDGD0))<0.0019){
+ //fill candidates D* and soft pion reco pt
- // label to the candidate soft pion
- Int_t pLabel = -1;
- if(fUseMCInfo) pLabel = aodTrack->GetLabel();
-
- // prepare the TLorentz vector for the pion
- Float_t pionMass = TDatabasePDG::Instance()->GetParticle(211)->Mass();
- fLorentzTrack3.SetPxPyPzE(aodTrack->Px(),aodTrack->Py(),aodTrack->Pz(),aodTrack->E(pionMass));
+ fRECOPtDStar->Fill(pt);
+ fPtPion->Fill(track2->Pt());
- // search for the D0
- for (Int_t iVertex = 0; iVertex<arrayVerticesHF->GetEntriesFast(); iVertex++) {
-
- Double_t invM = 0;
- Double_t invMDStar = 0;
- Double_t dPhi = 0;
-
- AliAODRecoDecayHF2Prong* vtx = (AliAODRecoDecayHF2Prong*)arrayVerticesHF->At(iVertex);
-
- Double_t pt = vtx->Pt();
-
- // acceptance variables
- Bool_t acceptanceProng0 = (TMath::Abs(vtx->EtaProng(0))<= 0.9 && vtx->PtProng(0) >= 0.1);
- Bool_t acceptanceProng1 = (TMath::Abs(vtx->EtaProng(1))<= 0.9 && vtx->PtProng(1) >= 0.1);
- Bool_t acceptanceSoftPi = (TMath::Abs(aodTrack->Eta())<= 0.9 && aodTrack->Pt() >= 0.05);
-
- Int_t pdgDgD0toKpi[2]={321,211};
-
- Int_t mcLabel =-1;
-
- Bool_t isDStar = kFALSE; // just to count
-
- //switch of MC for real data
-
- if(fUseMCInfo){
- mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
- if(!mcArray) {
- printf("AliAnalysisTaskSEDStarSpectra::UserExec: MC particles not found!\n");
- return;
- }
- mcLabel = vtx->MatchToMC(421, mcArray,2,pdgDgD0toKpi) ; //MC D0
- // matching to MC D*
- if(mcLabel !=-1 && pLabel!=-1 && nJets ==1) { // count only once in case of multijets
-
- // search for a D0 and a pi with mother and check it is a D*
- AliAODMCParticle* theMCpion = (AliAODMCParticle*)mcArray->At(pLabel);
- Int_t motherMCPion = theMCpion->GetMother();
- AliAODMCParticle* theMCD0 = (AliAODMCParticle*)mcArray->At(mcLabel);
- Int_t motherMCD0 = theMCD0->GetMother();
-
- if(motherMCPion!=-1 && (motherMCD0 == motherMCPion)){
- AliAODMCParticle* mcMother = (AliAODMCParticle*)mcArray->At(motherMCPion);
- if(TMath::Abs(mcMother->GetPdgCode()) == 413 && TMath::Abs(theMCpion->GetPdgCode())==211){
- isDStar = kTRUE;
- fiDstar++;
- }
- }
- }
- }
-
-
- if (acceptanceProng0 && acceptanceProng1 && acceptanceSoftPi) {
-
- AliDebug(2,"D0 reco daughters in acceptance");
- if(isDStar && nJets==1) icountRecoAcc++;
-
- // D0 daughters
- AliAODTrack *track0 = (AliAODTrack*)vtx->GetDaughter(0);
- AliAODTrack *track1 = (AliAODTrack*)vtx->GetDaughter(1);
-
- // check for ITS refit (already required at the ESD filter level )
- Bool_t kRefitITS = kTRUE;
-
- if((!(track0->GetStatus()&AliESDtrack::kITSrefit)|| (!(track1->GetStatus()&AliESDtrack::kITSrefit)))) {
- kRefitITS = kFALSE;
- }
-
- Int_t ncls0=0,ncls1=0,ncls2=0;
-
- for(Int_t l=0;l<6;l++) {
- if(TESTBIT(track0->GetITSClusterMap(),l)) ncls0++;
- if(TESTBIT(track1->GetITSClusterMap(),l)) ncls1++;
- if(TESTBIT(aodTrack->GetITSClusterMap(),l)) ncls2++;
- }
-
- // clusters in ITS for D0 daugthers and soft pion
- if (ncls0 >= fMinITSClusters && ncls1 >= fMinITSClusters && ncls2>=4 && kRefitITS) {
-
- if(isDStar && nJets==1) icountRecoITSClusters++;
-
- // D0 cutting varibles
- Double_t cosThetaStar = 9999.;
- Double_t pTpi = 0.;
- Double_t pTK = 0.;
- Double_t d01 = 0;
- Double_t d00 = 0;
-
- // D0, D0bar
- if (fComputeD0){
- cosThetaStar = vtx->CosThetaStarD0();
- pTpi = vtx->PtProng(0);
- pTK = vtx->PtProng(1);
- d01 = vtx->Getd0Prong(0)*1E4;
- d00 = vtx->Getd0Prong(1)*1E4;
- }else{
- cosThetaStar = vtx->CosThetaStarD0bar();
- pTpi = vtx->PtProng(1);
- pTK = vtx->PtProng(0);
- d01 = vtx->Getd0Prong(1)*1E4;
- d00 = vtx->Getd0Prong(0)*1E4;
- }
-
- AliDebug(2,"D0 reco daughters in acceptance");
-
- Double_t dca = vtx->GetDCA()*1E4;
- Double_t d0d0 = vtx->Prodd0d0()*1E8;
-
- TVector3 mom(vtx->Px(),vtx->Py(),vtx->Pz());
- TVector3 flight((vtx->Xv())- primaryPos[0],(vtx->Yv())- primaryPos[1],(vtx->Zv())- primaryPos[2]);
- Double_t pta = mom.Angle(flight);
- Double_t cosPointingAngle = TMath::Cos(pta);
-
- // Crstian Ivan D* cuts. Multidimensional optimization
- Double_t cuts[7] = {9999999., 1.1, 0., 9999999.,999999., 9999999., 0.};
-
- if (pt <= 1){ // first bin not optimized
- cuts[0] = 400;
- cuts[1] = 0.8;
- cuts[2] = 0.21;
- cuts[3] = 500;
- cuts[4] = 500;
- cuts[5] = -20000;
- cuts[6] = 0.6;
- }
- else if (pt > 1 && pt <= 2){
- cuts[0] = 200;
- cuts[1] = 0.7;
- cuts[2] = 0.8;
- cuts[3] = 210;
- cuts[4] = 210;
- cuts[5] = -20000;
- cuts[6] = 0.9;
- }
- else if (pt > 2 && pt <= 3){
- cuts[0] = 400;
- cuts[1] = 0.8;
- cuts[2] = 0.8;
- cuts[3] = 420;
- cuts[4] = 350;
- cuts[5] = -8500;
- cuts[6] = 0.9;
- }
- else if (pt > 3 && pt <= 5){
- cuts[0] = 160;
- cuts[1] = 1.0;
- cuts[2] = 1.2;
- cuts[3] = 420;
- cuts[4] = 560;
- cuts[5] = -8500;
- cuts[6] = 0.9;
- }
- else if (pt > 5){
- cuts[0] = 800;
- cuts[1] = 1.0;
- cuts[2] = 1.2;
- cuts[3] = 700;
- cuts[4] = 700;
- cuts[5] = 10000;
- cuts[6] = 0.9;
- }
- // apply D0 cuts
-
- if (dca < cuts[0]
- && TMath::Abs(cosThetaStar) < cuts[1]
- && pTpi > cuts[2]
- && pTK > cuts[2]
- && TMath::Abs(d01) < cuts[3]
- && TMath::Abs(d00) < cuts[4]
- && d0d0 < cuts[5]
- && cosPointingAngle > cuts[6]
- ){
-
- if(fComputeD0){ // D0 from default
-
- if(vtx->ChargeProng(0)==-1){
- fLorentzTrack1.SetPxPyPzE( vtx->PxProng(0),vtx->PyProng(0), vtx->PzProng(0),vtx->EProng(0,321) );
- fLorentzTrack2.SetPxPyPzE( vtx->PxProng(1),vtx->PyProng(1), vtx->PzProng(1),vtx->EProng(1,211) );
- }else{
- fLorentzTrack1.SetPxPyPzE( vtx->PxProng(0),vtx->PyProng(0), vtx->PzProng(0),vtx->EProng(0,211) );
- fLorentzTrack2.SetPxPyPzE( vtx->PxProng(1),vtx->PyProng(1), vtx->PzProng(1),vtx->EProng(1,321) );
- }
-
- vD0[0] = (fLorentzTrack1+fLorentzTrack2).Px();
- vD0[1] = (fLorentzTrack1+fLorentzTrack2).Py();
- vD0[2] = (fLorentzTrack1+fLorentzTrack2).Pz();
- vD0[3] = (fLorentzTrack1+fLorentzTrack2).E();
-
- fLorentzTrack4.SetPxPyPzE(vD0[0],vD0[1],vD0[2],vD0[3]);
-
- }else{ // D0bar analysis
-
- if(vtx->ChargeProng(0)==-1){
- fLorentzTrack1.SetPxPyPzE( vtx->PxProng(0),vtx->PyProng(0), vtx->PzProng(0),vtx->EProng(0,211) );
- fLorentzTrack2.SetPxPyPzE( vtx->PxProng(1),vtx->PyProng(1), vtx->PzProng(1),vtx->EProng(1,321) );
- }else{
- fLorentzTrack1.SetPxPyPzE( vtx->PxProng(0),vtx->PyProng(0), vtx->PzProng(0),vtx->EProng(0,321) );
- fLorentzTrack2.SetPxPyPzE( vtx->PxProng(1),vtx->PyProng(1), vtx->PzProng(1),vtx->EProng(1,211) );
- }
-
- vD0bar[0] = (fLorentzTrack1+fLorentzTrack2).Px();
- vD0bar[1] = (fLorentzTrack1+fLorentzTrack2).Py();
- vD0bar[2] = (fLorentzTrack1+fLorentzTrack2).Pz();
- vD0bar[3] = (fLorentzTrack1+fLorentzTrack2).E();
-
- fLorentzTrack4.SetPxPyPzE(vD0bar[0],vD0bar[1],vD0bar[2],vD0bar[3]);
- }
-
- // D0-D0bar related variables
-
- invM = GetInvariantMass(fLorentzTrack1,fLorentzTrack2);
- if(nJets==1) fInvMass->Fill(invM);
-
- if(isDStar && nJets==1) icountRecoPPR++;
-
- //DStar invariant mass
- invMDStar = GetInvariantMassDStar(fLorentzTrack3,fLorentzTrack4);
-
- //conversion from phi TLorentzVerctor to phi aliroot
- Double_t kconvert = 0;
- Double_t phiDStar = (fLorentzTrack3 + fLorentzTrack4).Phi();
- kconvert = phiDStar;
- if(phiDStar<0) kconvert = phiDStar + 2*(TMath::Pi());
- phiDStar = kconvert;
-
- // dphi between jet and D*
- dPhi = phiJet - phiDStar;
-
- //Just for plotting pourposal
- if(dPhi<=-1.58) dPhi = TMath::Abs(dPhi);
- if(dPhi>4.72){
- dPhi = dPhi-2*(TMath::Pi());
- }
-
- //Alternative cutting procedure
- //the cut on cosThetaStar is to reduce near collinear
- //background from jet fragmentation
-
- Bool_t tCut = EvaluateCutOnPiD0pt(vtx,aodTrack);
-
- if(ftopologicalCut && tCut) continue;
- if(ftopologicalCut && TMath::Abs(cosThetaStar)>cuts[1]) continue;
-
- if(invM >= 1.829 && invM <= 1.901){ // ~4 sigma cut on D0 mass
-
- if(isDStar && nJets==1) {
- fDStarMass->Fill(invMDStar);
- fTrueDiff->Fill(invMDStar-invM);
- }
-
- if(nJets==1) { // avoid double counting
- fDStar->Fill(invMDStar);
- fDiff->Fill(invMDStar-invM); // M(Kpipi)-M(Kpi)
- }
-
- // now the dphi signal and the fragmentation function
- if((invMDStar-invM)<=0.148 && (invMDStar-invM)>=0.142) {
-
- //fill candidates D* and soft pion reco pt
- if(nJets==1) fPtPion->Fill(aodTrack->Pt());
- if(nJets==1) fRECOPtDStar->Fill((fLorentzTrack3 + fLorentzTrack4).Pt());
- fPhi ->Fill(dPhi);
-
- if(dPhi>=-0.5 && dPhi<=0.5){ // evaluate in the near side
- Double_t dStarMom = (fLorentzTrack3 + fLorentzTrack4).P();
- Double_t zFrag = (TMath::Cos(dPhi)*dStarMom)/ejet;
- fResZ->Fill(TMath::Abs(zFrag));
- }
- }
- }
-
- // evaluate side band background
- SideBandBackground(invM, invMDStar, ejet, dPhi, nJets);
-
- invM = 0;
- invMDStar = 0;
-
- } // end PPR cuts
- } // end ITS cluster
- } // end acceptance
- } // D0 cand
- } // softpion
- } // tracks
- } // jets
-
- AliDebug(2, Form("Found %i Reco particles that are D0!!",icountReco));
-
- fCountReco+= fiDstar;
- fCountRecoAcc+= icountRecoAcc;
- fCountRecoITSClusters+= icountRecoITSClusters;
- fCountRecoPPR+= icountRecoPPR;
+ fPhi ->Fill(dPhi);
- PostData(1,fOutput);
+ Double_t jetCone = 0.4;
+ if(dPhi>=-jetCone && dPhi<=jetCone){ // evaluate in the near side inside UA1 radius
+ Double_t zFrag = (TMath::Cos(dPhi)*dStarMom)/EGjet;
+ fResZ->Fill(TMath::Abs(zFrag));
+ }
+ }
+ }
+ // evaluate side band background
+ SideBandBackground(finvM, finvMDStar, dStarMom, EGjet, dPhi);
+
+ } // D* cand
+
+AliDebug(2, Form("Found %i Reco particles that are D*!!",icountReco));
+
+PostData(1,fOutput);
}
//________________________________________ terminate ___________________________
// a query. It always runs on the client, it can be used to present
// the results graphically or save the results to file.
+ Info("Terminate"," terminate");
AliAnalysisTaskSE::Terminate();
-
- AliInfo(Form("Found %i of that MC D*->D0pi(D0->kpi) are in acceptance, in %d events",fCountDStar,fEvents));
-
- AliInfo(Form("Found %i RECO particles after cuts that are DStar, in %d events",fCountReco,fEvents));
- AliInfo(Form("Among the above, found %i reco D* that are decaying in K+pi+pi and are in the requested acceptance, in %d events",fCountRecoAcc,fEvents));
- AliInfo(Form("Among the above, found %i reco D* that are decaying in K+pi+pi and have at least %d clusters in ITS, in %d events",fCountRecoITSClusters,fMinITSClusters,fEvents));
- AliInfo(Form("Among the above, found %i reco D* that are decaying in K+pi+pi and satisfy D* cuts, in %d events",fCountRecoPPR,fEvents));
fOutput = dynamic_cast<TList*> (GetOutputData(1));
if (!fOutput) {
return;
}
- fcharmpt = dynamic_cast<TH1F*>(fOutput->FindObject("fcharmpt"));
- fdstarE = dynamic_cast<TH1F*>(fOutput->FindObject("fdstarE"));
- fdstarpt = dynamic_cast<TH1F*>(fOutput->FindObject("fdstarpt"));
fDStarMass = dynamic_cast<TH1F*>(fOutput->FindObject("fDStarMass"));
fTrueDiff = dynamic_cast<TH1F*>(fOutput->FindObject("fTrueDiff"));
fInvMass = dynamic_cast<TH1F*>(fOutput->FindObject("fInvMass"));
fDiffSideBand = dynamic_cast<TH1F*>(fOutput->FindObject("fDiffSideBand"));
ftrigger = dynamic_cast<TH1F*>(fOutput->FindObject("ftrigger"));
fRECOPtDStar = dynamic_cast<TH1F*>(fOutput->FindObject("fRECOPtDStar"));
+ fRECOPtBkg = dynamic_cast<TH1F*>(fOutput->FindObject("fRECOPtBkg"));
fEjet = dynamic_cast<TH1F*>(fOutput->FindObject("fEjet"));
fPhijet = dynamic_cast<TH1F*>(fOutput->FindObject("fPhijet"));
fEtaJet = dynamic_cast<TH1F*>(fOutput->FindObject("fEtaJet"));
fResZ = dynamic_cast<TH1F*>(fOutput->FindObject("fResZ"));
fResZBkg = dynamic_cast<TH1F*>(fOutput->FindObject("fResZBkg"));
fPhiBkg = dynamic_cast<TH1F*>(fOutput->FindObject("fPhiBkg"));
- fD0ptvsSoftPtSignal = dynamic_cast<TH2F*>(fOutput->FindObject("fD0ptvsSoftPtSignal"));
- fD0ptvsSoftPt = dynamic_cast<TH2F*>(fOutput->FindObject("fD0ptvsSoftPt"));
-
+ theMCFF = dynamic_cast<TH1F*>(fOutput->FindObject("theMCFF"));
+ fDphiD0Dstar = dynamic_cast<TH1F*>(fOutput->FindObject("fDphiD0Dstar"));
+ fPtJet = dynamic_cast<TH1F*>(fOutput->FindObject("fPtJet"));
+
}
//___________________________________________________________________________
void AliAnalysisTaskSEDStarJets::UserCreateOutputObjects() {
- // output
-
+ // output
Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
//slot #1
return;
}
-
-//_______________________________D0 invariant mass________________________________
-
-Double_t AliAnalysisTaskSEDStarJets::GetInvariantMass(TLorentzVector LorentzTrack1, TLorentzVector LorentzTrack2){
-
- return TMath::Abs((LorentzTrack1+LorentzTrack2).M()); // invariant mass of two tracks
-}
-//______________________________D* invariant mass_________________________________
-
-Double_t AliAnalysisTaskSEDStarJets::GetInvariantMassDStar(TLorentzVector LorentzTrack4, TLorentzVector LorentzTrack3){
-
- return TMath::Abs((LorentzTrack4+LorentzTrack3).M()); // invariant mass of two tracks
-}
-//_________________________________D* in MC _______________________________________
-
-Bool_t AliAnalysisTaskSEDStarJets::DstarInMC(AliAODMCParticle* const mcPart, TClonesArray* mcArray){
- // D* in MC
-
- Bool_t isOk = kFALSE;
-
- if(TMath::Abs(mcPart->GetPdgCode())!=413) return isOk;
-
- // getting the daughters
- Int_t daughter0 = mcPart->GetDaughter(0);
- Int_t daughter1 = mcPart->GetDaughter(1);
-
- AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughter0,daughter1));
- if (daughter0 == 0 || daughter1 == 0) {
- AliDebug(2, "Error! the D* MC doesn't have correct daughters!!");
- return isOk;
- }
-
- if (TMath::Abs(daughter1 - daughter0) != 1) { // should be everytime true - see PDGdatabooklet
- AliDebug(2, "The D* MC doesn't come from a 2-prong decay, skipping!!");
- return isOk;
- }
-
- AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter0));
- AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter1));
- if (!mcPartDaughter0 || !mcPartDaughter1) {
- AliWarning("At least one Daughter Particle not found in tree, skipping");
- return isOk;
- }
-
- if (!(TMath::Abs(mcPartDaughter0->GetPdgCode())==421 &&
- TMath::Abs(mcPartDaughter1->GetPdgCode())==211) &&
- !(TMath::Abs(mcPartDaughter0->GetPdgCode())==211 &&
- TMath::Abs(mcPartDaughter1->GetPdgCode())==421)) {
- AliDebug(2, "The D* MC doesn't come from a Kpi decay, skipping!!");
- return isOk;
- }
-
- Double_t vtx2daughter0[3] = {0,0,0}; // secondary vertex from daughter 0
- Double_t vtx2daughter1[3] = {0,0,0}; // secondary vertex from daughter 1
-
- // getting vertex from daughters
- mcPartDaughter0->XvYvZv(vtx2daughter0); // cm
- mcPartDaughter1->XvYvZv(vtx2daughter1); //cm
-
- // check if the secondary vertex is the same for both
- if (vtx2daughter0[0] != vtx2daughter1[0] && vtx2daughter0[1] != vtx2daughter1[1] && vtx2daughter0[2] != vtx2daughter1[2]) {
- AliError("Daughters have different secondary vertex, skipping the track");
- return isOk;
- }
-
- AliAODMCParticle* neutralDaugh = mcPartDaughter0;
-
- if (!EvaluateIfD0toKpi(neutralDaugh,mcArray)) {
- AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
- return isOk;
- }
-
- isOk = kTRUE;
-
- return isOk;
-
-}
-
-Bool_t AliAnalysisTaskSEDStarJets::EvaluateIfD0toKpi(AliAODMCParticle* neutralDaugh, TClonesArray* mcArray)const{
-
- //
- // chack wether D0 is decaing into kpi
- //
-
- Bool_t isHadronic = kFALSE;
-
- Int_t daughterD00 = neutralDaugh->GetDaughter(0);
- Int_t daughterD01 = neutralDaugh->GetDaughter(1);
-
- AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughterD00,daughterD01));
- if (daughterD00 == 0 || daughterD01 == 0) {
- AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
- return isHadronic;
- }
-
- if (TMath::Abs(daughterD01 - daughterD00) != 1) { // should be everytime true - see PDGdatabooklet
- AliDebug(2, "The D0 MC doesn't come from a 2-prong decay, skipping!!");
- return isHadronic;
- }
-
- AliAODMCParticle* mcPartDaughterD00 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughterD00));
- AliAODMCParticle* mcPartDaughterD01 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughterD01));
- if (!mcPartDaughterD00 || !mcPartDaughterD01) {
- AliWarning("At least one Daughter Particle not found in tree, skipping");
- return isHadronic;
- }
-
- if (!(TMath::Abs(mcPartDaughterD00->GetPdgCode())==321 &&
- TMath::Abs(mcPartDaughterD01->GetPdgCode())==211) &&
- !(TMath::Abs(mcPartDaughterD00->GetPdgCode())==211 &&
- TMath::Abs(mcPartDaughterD01->GetPdgCode())==321)) {
- AliDebug(2, "The D* MC doesn't come from a Kpi decay, skipping!!");
- return isHadronic;
- }
-
- isHadronic = kTRUE;
-
-
- return isHadronic;
-
-}
-
-
//___________________________________ hiostograms _______________________________________
Bool_t AliAnalysisTaskSEDStarJets::DefineHistoFroAnalysis(){
fPhi->GetYaxis()->SetTitle("Entries");
fOutput->Add(fPhi);
+ fDphiD0Dstar = new TH1F("phiD0Dstar","Delta phi between D0 and DStar ",1000,-6.5,6.5);
+ fOutput->Add(fDphiD0Dstar);
+
fPhiBkg = new TH1F("phiBkg","Delta phi between Jet axis and DStar background ",25,-1.57,4.72);
fPhiBkg->SetStats(kTRUE);
fPhiBkg->GetXaxis()->SetTitle("#Delta #phi (rad)");
fPhiBkg->GetYaxis()->SetTitle("Entries");
fOutput->Add(fPhiBkg);
- fRECOPtDStar = new TH1F("RECODStar1","RECO DStar pt distribution",600,0,15);
+ fRECOPtDStar = new TH1F("RECODStar1","RECO DStar pt distribution",30,0,30);
fRECOPtDStar->SetStats(kTRUE);
fRECOPtDStar->SetLineColor(2);
fRECOPtDStar->GetXaxis()->SetTitle("GeV/c");
fRECOPtDStar->GetYaxis()->SetTitle("Entries");
fOutput->Add(fRECOPtDStar);
+
+ fRECOPtBkg = new TH1F("RECOptBkg","RECO pt distribution side bands",30,0,30);
+ fRECOPtBkg->SetStats(kTRUE);
+ fRECOPtBkg->SetLineColor(2);
+ fRECOPtBkg->GetXaxis()->SetTitle("GeV/c");
+ fRECOPtBkg->GetYaxis()->SetTitle("Entries");
+ fOutput->Add(fRECOPtBkg);
- fPtPion = new TH1F("pionpt","Primary pions candidates pt ",500,0,5);
+ fPtPion = new TH1F("pionpt","Primary pions candidates pt ",500,0,10);
fPtPion->SetStats(kTRUE);
fPtPion->GetXaxis()->SetTitle("GeV/c");
fPtPion->GetYaxis()->SetTitle("Entries");
fOutput->Add(fPtPion);
-
- fcharmpt = new TH1F("charmpt","MC Charm pt distribution", 10000,0,5000);
- fdstarE = new TH1F("dstare", "MC D* energy distribution", 10000,0,5000);
- fdstarpt = new TH1F("dstarpt","MC D* pt distribution", 10000,0,5000);
- fOutput->Add(fcharmpt);
- fOutput->Add(fdstarE);
- fOutput->Add(fdstarpt);
-
+
// jet related fistograms
fEjet = new TH1F("ejet", "UA1 algorithm jet energy distribution",1000,0,500);
fPhijet = new TH1F("Phijet","UA1 algorithm jet #phi distribution", 200,-7,7);
- fEtaJet = new TH1F("Etajet","UA1 algorithm jet #eta distribution", 200,-7,7);
+ fEtaJet = new TH1F("Etajet","UA1 algorithm jet #eta distribution", 200,-7,7);
+ fPtJet = new TH1F("PtJet", "UA1 algorithm jet Pt distribution",1000,0,500);
fOutput->Add(fEjet);
fOutput->Add(fPhijet);
fOutput->Add(fEtaJet);
-
+ fOutput->Add(fPtJet);
+
+ theMCFF = new TH1F("FragFuncMC","Fragmentation function in MC for FC ",100,0,10);
fResZ = new TH1F("FragFunc","Fragmentation function ",50,0,1);
fResZBkg = new TH1F("FragFuncBkg","Fragmentation function background",50,0,1);
+ fOutput->Add(theMCFF);
fOutput->Add(fResZ);
fOutput->Add(fResZBkg);
- fD0ptvsSoftPt = new TH2F("D0piRec","Candidates (background + signal)",100,0,3,100,0,15);
- fD0ptvsSoftPtSignal = new TH2F("D0PiSignal","Signal",100,0,3,100,0,15);
- fOutput->Add(fD0ptvsSoftPt);
- fOutput->Add(fD0ptvsSoftPtSignal);
-
- return kTRUE;
-
-}
-
-//______________________________Phase space cut alternative to PPR _______________________________________
-
-Bool_t AliAnalysisTaskSEDStarJets::EvaluateCutOnPiD0pt(AliAODRecoDecayHF2Prong* const vtx, AliAODTrack* const aodTrack) {
-
- // The soft pion pt and DO pt are strongly correlated. It can be shown that ~ 95% of the signal in constrained
- // into a narrow band defined by 10 < D0pt/softPt < 18. This cut can be used with a relaxed CosThetaStar cut
- // to reconstruct the D*.
-
- Double_t softPt = 0;
- Double_t d0ptReco = 0;
-
- softPt = aodTrack->Pt();
- d0ptReco = vtx->Pt();
- fD0ptvsSoftPt->Fill(softPt,d0ptReco);
-
- if(softPt>0){
- Double_t ratio = d0ptReco/softPt;
- if( ratio <=10. || ratio>=18. ) return kFALSE;
- }
-
- fD0ptvsSoftPtSignal->Fill(softPt,d0ptReco);
- return kTRUE;
+ return kTRUE;
}
//______________________________ side band background for D*___________________________________
-void AliAnalysisTaskSEDStarJets::SideBandBackground(Double_t invM, Double_t invMDStar, Double_t ejet, Double_t dPhi, Int_t nJets){
+void AliAnalysisTaskSEDStarJets::SideBandBackground(Double_t invM, Double_t invMDStar, Double_t dStarMomBkg, Double_t EGjet, Double_t dPhi){
// D* side band background method. Two side bands, in M(Kpi) are taken at ~6 sigmas
// (expected detector resolution) on the left and right frm the D0 mass. Each band
// has a width of ~5 sigmas. Two band needed for opening angle considerations
- if((invM>=1.763 && invM<=1.811) || (invM>=1.919 && invM<=1.963)){
+ if((invM>=1.7 && invM<=1.8) || (invM>=1.92 && invM<=2.19)){
+ fDiffSideBand->Fill(invMDStar-invM); // M(Kpipi)-M(Kpi) side band background
+ if ((invMDStar-invM)<=0.14732 && (invMDStar-invM)>=0.14352) {
+ fPhiBkg->Fill(dPhi);
+ fRECOPtBkg->Fill(dStarMomBkg);
+ if(dPhi>=-0.4 && dPhi<=0.4){ // evaluate in the near side
+ Double_t zFragBkg = (TMath::Cos(dPhi)*dStarMomBkg)/EGjet;
+ fResZBkg->Fill(TMath::Abs(zFragBkg));
+ }
+ }
+ }
+}
+
+//_____________________________________________________________________________________________-
+double AliAnalysisTaskSEDStarJets::FillMCFF(AliAODMCParticle* mcPart, TClonesArray* mcArray, Int_t mcLabel){
+ //
+ // GS from MC
+ // UA1 jet algorithm reproduced in MC
+ //
+ Double_t zMC2 =-999;
+
+ Double_t leading =0;
+ Double_t PartE = 0;
+ Double_t PhiLeading = -999;
+ Double_t EtaLeading = -999;
+ Double_t PtLeading = -999;
+ Int_t counter =-999;
+
+ //find leading particle
+ for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) {
+ AliAODMCParticle* Part = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
+ if (!Part) {
+ AliWarning("MC Particle not found in tree, skipping");
+ continue;
+ }
+
+ // remove quarks and the leading particle (it will be counted later)
+ if(iPart == mcLabel) continue;
+ if(iPart <= 8) continue;
- if(nJets==1)fDiffSideBand->Fill(invMDStar-invM); // M(Kpipi)-M(Kpi) side band background
+ //remove resonances not directly detected in detector
+ Int_t PDGCode = Part->GetPdgCode();
+
+ // be sure the particle reach the detector
+ Double_t x = Part->Xv();
+ Double_t y = Part->Yv();
+ Double_t z = Part->Zv();
- if ((invMDStar-invM)<=0.148 && (invMDStar-invM)>=0.142) {
- fPhiBkg->Fill(dPhi);
-
- if(dPhi>=-0.5 && dPhi<=0.5){ // evaluate in the near side
-
- Double_t dStarMomBkg = (fLorentzTrack3 + fLorentzTrack4).P();
- Double_t zFragBkg = (TMath::Cos(dPhi)*dStarMomBkg)/ejet;
- fResZBkg->Fill(TMath::Abs(zFragBkg));
-
+ if(TMath::Abs(PDGCode)== 2212 && x<3 && y<3) continue;
+ if(TMath::Abs(x)>30 || TMath::Abs(y)>30 || TMath::Abs(z)>30 ) continue;
+ if(TMath::Abs(PDGCode)!=211 && TMath::Abs(PDGCode)!=321 && TMath::Abs(PDGCode)!=11 && TMath::Abs(PDGCode)!=13 && TMath::Abs(PDGCode)!=2212) continue;
+
+ Int_t daug0 = -999;
+ Double_t xd =-999;
+ Double_t yd =-999;
+ Double_t zd =-999;
+
+ daug0 = Part->GetDaughter(0);
+
+ if(daug0>=0){
+ AliAODMCParticle* tdaug = dynamic_cast<AliAODMCParticle*>(mcArray->At(daug0));
+ if(tdaug){
+ xd = tdaug->Xv();
+ yd = tdaug->Yv();
+ zd = tdaug->Zv();
}
}
+ if(TMath::Abs(xd)<3 || TMath::Abs(yd)<3) continue;
+
+ Bool_t AliceAcc = (TMath::Abs(Part->Eta()) <= 0.9);
+ if(!AliceAcc) continue;
+
+ PartE = Part->E();
+
+ if(PartE>leading){
+ leading = Part->E();
+ PhiLeading = Part->Phi();
+ EtaLeading = Part->Eta();
+ PtLeading = Part->Pt();
+ counter = iPart;
+ }
+
}
+
+ Double_t jetEnergy = 0;
+
+ //reconstruct the jet
+ for (Int_t iiPart=0; iiPart<mcArray->GetEntriesFast(); iiPart++) {
+ AliAODMCParticle* tPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iiPart));
+ if (!tPart) {
+ AliWarning("MC Particle not found in tree, skipping");
+ continue;
+ }
+ // remove quarks and the leading particle (it will be counted later)
+ if(iiPart == counter) continue; // do not count again the leading particle
+ if(iiPart == mcLabel) continue;
+ if(iiPart <= 8) continue;
+
+ //remove resonances not directly detected in detector
+ Int_t PDGCode = tPart->GetPdgCode();
+
+ // be sure the particle reach the detector
+ Double_t x = tPart->Xv();
+ Double_t y = tPart->Yv();
+ Double_t z = tPart->Zv();
+
+ if(TMath::Abs(PDGCode)== 2212 && (x<3 && y<3)) continue;
+ if(TMath::Abs(x)>30 || TMath::Abs(y)>30 || TMath::Abs(z)>30 ) continue; // has to be generated at least in the silicon tracker or beam pipe
+ if(TMath::Abs(PDGCode)!=211 && TMath::Abs(PDGCode)!=321 && TMath::Abs(PDGCode)!=11 && TMath::Abs(PDGCode)!=13 && TMath::Abs(PDGCode)!=2212) continue;
+
+
+ Int_t daug0 = -999;
+ Double_t xd =-999;
+ Double_t yd =-999;
+ Double_t zd =-999;
+
+ daug0 = tPart->GetDaughter(0);
+
+ if(daug0>=0){
+ AliAODMCParticle* tdaug = dynamic_cast<AliAODMCParticle*>(mcArray->At(daug0));
+ if(tdaug){
+ xd = tdaug->Xv();
+ yd = tdaug->Yv();
+ zd = tdaug->Zv();
+ }
+ }
+ if(TMath::Abs(xd)<3 && TMath::Abs(yd)<3) continue;
+ //remove particles not in ALICE acceptance
+ if(tPart->Pt()<0.07) continue;
+ Bool_t AliceAcc = (TMath::Abs(tPart->Eta()) <= 0.9);
+ if(!AliceAcc) continue;
+
+ Double_t EtaMCp = tPart->Eta();
+ Double_t PhiMCp = tPart->Phi();
+
+ Double_t DphiMClead = PhiLeading-PhiMCp;
+
+ if(DphiMClead<=-(TMath::Pi())/2) DphiMClead = DphiMClead+2*(TMath::Pi());
+ if(DphiMClead>(3*(TMath::Pi()))/2) DphiMClead = DphiMClead-2*(TMath::Pi());
+
+ Double_t deta = (EtaLeading-EtaMCp);
+ //cone radius
+ Double_t radius = TMath::Sqrt((DphiMClead*DphiMClead)+(deta*deta));
+
+ if(radius>0.4) continue; // in the jet cone
+ if(tPart->Charge()==0) continue; // only charged fraction
+
+ jetEnergy = jetEnergy+(tPart->E());
+ }
+
+ jetEnergy = jetEnergy + leading;
+
+ // delta phi D*, jet axis
+ Double_t dPhi = PhiLeading - (mcPart->Phi());
+
+ //plot right dphi
+ if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi());
+ if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi());
+
+ zMC2 = (TMath::Cos(dPhi)*(mcPart->P()))/jetEnergy;
+
+ return zMC2;
}
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff