#include "AliRDHFCutsD0toKpi.h"
#include "AliRDHFCutsDStartoKpipi.h"
#include "AliRhoParameter.h"
+#include "AliParticleContainer.h"
ClassImp(AliAnalysisTaskFlavourJetCorrelations)
fMinMass(),
fMaxMass(),
fJetArrName(0),
+fTrackArrName(0),
fCandArrName(0),
fLeadingJetOnly(kFALSE),
fJetRadius(0),
fSwitchOnOutOfConeAxis(0),
fSwitchOnSparses(1),
fNAxesBigSparse(9),
+fJetCont(0),
+fTrackCont(0),
+fClusterCont(0),
fhstat(),
fhPtJetTrks(),
fhPhiJetTrks(),
fhNJetPerEv(),
fhdeltaRJetTracks(),
fhsJet(),
+fhImpPar(),
fhNDPerEvNoJet(),
fhptDPerEvNoJet(),
fhNJetPerEvNoD(),
fhPtJetPerEvNoD(),
+fhVtxX(),
+fhVtxY(),
+fhVtxZ(),
fhsDstandalone(),
fhInvMassptD(),
fhDiffSideBand(),
fhDeltaRptDptj(),
fhDeltaRptDptjB(),
fhsDphiz()
+
{
//
// Default ctor
fMinMass(),
fMaxMass(),
fJetArrName(0),
+fTrackArrName(0),
fCandArrName(0),
fLeadingJetOnly(kFALSE),
fJetRadius(0),
fSwitchOnOutOfConeAxis(0),
fSwitchOnSparses(1),
fNAxesBigSparse(9),
+fJetCont(0),
+fTrackCont(0),
+fClusterCont(0),
fhstat(),
fhPtJetTrks(),
fhPhiJetTrks(),
fhNJetPerEv(),
fhdeltaRJetTracks(),
fhsJet(),
+fhImpPar(),
fhNDPerEvNoJet(),
fhptDPerEvNoJet(),
fhNJetPerEvNoD(),
fhPtJetPerEvNoD(),
+fhVtxX(),
+fhVtxY(),
+fhVtxZ(),
fhsDstandalone(),
fhInvMassptD(),
fhDiffSideBand(),
// output
Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
AliAnalysisTaskEmcal::UserCreateOutputObjects();
+
+ fJetCont = GetJetContainer(0);
+ if(fJetCont){
+ fTrackCont = fJetCont->GetParticleContainer();
+ fClusterCont = fJetCont->GetClusterContainer();
+ }
+
+
// define histograms
// the TList fOutput is already defined in AliAnalysisTaskEmcal::UserCreateOutputObjects()
DefineHistoForAnalysis();
AliAODEvent *aodFromExt = ext->GetAOD();
arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject(fBranchName.Data());
}
- } else {
+ } else if(aodEvent){
arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject(fBranchName.Data());
}
}
//retrieve jets
- fTrackArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("PicoTracks"));
+ //this is a duplication of fTrackCont, but is is used in the loop of line 598 and changing it needs a thorough test
+ fTrackArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTrackArrName));
//clusArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("CaloClustersCorr"));
- //jetArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetArrName));
- fJetRadius=GetJetContainer(0)->GetJetRadius();
+ //fJetArray = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetArrName));
+ //fJetContainer=GetJetContainer(0);
+ //if(!fJetContainer) {
+ // AliError("Jet Container 0 not found");
+ // return kFALSE;
+ //}
+ fJetRadius=fJetCont->GetJetRadius();
if(!fTrackArr){
AliInfo(Form("Could not find the track array\n"));
if(!iseventselected) return kFALSE;
fhstat->Fill(1);
-
+
+ const AliVVertex *vertex = aodEvent->GetPrimaryVertex();
+ Double_t vtx[3];
+ vertex->GetXYZ(vtx);
+ fhVtxX->Fill(vtx[0]);
+ fhVtxY->Fill(vtx[1]);
+ fhVtxZ->Fill(vtx[2]);
+ //Printf(">>>>>>>>.>>>Primary vertex %.3f,%.3f,%.3f",vtx[0], vtx[1],vtx[2]);
//retrieve charm candidates selected
Int_t candidates = 0;
- Int_t njets=GetJetContainer()->GetNJets();
-
+ Int_t njets=fJetCont->GetNJets();
+ //Printf("N jets in this event %d",njets);
if(!fJetOnlyMode) {
candidates = fCandidateArray->GetEntriesFast();
//If there's no background subtraction rhoval=0 and momentum is simply not took into account
AliRhoParameter *rho = 0;
Double_t rhoval = 0;
- TString sname("Rho");
- if (!sname.IsNull()) {
- rho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(sname));
+ if (!fJetCont->GetRhoName().IsNull()) {
+ rho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fJetCont->GetRhoName()));
if(rho) rhoval = rho->GetVal();
}
Double_t leadingJet =0;
Double_t pointJ[6];
- Int_t ntrarr=fTrackArr->GetEntriesFast();
+ Int_t ntrarr=fTrackCont->GetNParticles();
fhNtrArr->Fill(ntrarr);
for(Int_t i=0;i<ntrarr;i++){
- AliVTrack *jtrack=static_cast<AliVTrack*>(fTrackArr->At(i));
+ AliAODTrack *jtrack=static_cast<AliAODTrack*>(fTrackCont->GetParticle(i));
//reusing histograms
+ if(!jtrack || jtrack->IsMuonTrack()) continue; // added check on IsMuonTrack because in some cases the DCA() gives crazy YAtDCA values that issue floating point exception
fhPtJetTrks->Fill(jtrack->Pt());
fhPhiJetTrks->Fill(jtrack->Phi());
fhEtaJetTrks->Fill(jtrack->Eta());
fhEjetTrks->Fill(jtrack->E());
+ Double_t vDCAglobalxy;
+ //Double_t vDCAglobalz;
+ if(jtrack->IsGlobalConstrained()){
+ // retrieve impact parameter
+ Double_t pos[3];
+
+ jtrack->GetXYZ(pos);
+
+ vDCAglobalxy = pos[0] - vtx[0]; //should be impact parameter in transverse plane
+ //vDCAglobalz = pos[1] - vtx[1]; //should be impact parameter in z direction
+ } else {
+ vDCAglobalxy=jtrack->DCA();
+ //vDCAglobalz=jtrack->ZAtDCA();
+ }
+ fhImpPar->Fill(vDCAglobalxy,jtrack->Pt());
+ //Printf("Track position %.3f,%.3f,%.3f",pos[0],pos[1],pos[2]);
}
//option to use only the leading jet
if(fLeadingJetOnly){
for (Int_t iJetsL = 0; iJetsL<njets; iJetsL++) {
- AliEmcalJet* jetL = (AliEmcalJet*)GetJetFromArray(iJetsL);
+ AliEmcalJet* jetL = (AliEmcalJet*)fJetCont->GetJet(iJetsL);
ptjet = jetL->Pt() - jetL->Area()*rhoval; //It takes into account the background subtraction
if(ptjet>leadingJet ) leadingJet = ptjet;
fPmissing[2]=0;
//Printf("Jet N %d",iJets);
- AliEmcalJet* jet = (AliEmcalJet*)GetJetFromArray(iJets);
+ AliEmcalJet* jet = (AliEmcalJet*)fJetCont->GetJet(iJets);
//Printf("Corr task Accept Jet");
if(!AcceptJet(jet)) {
fhstat->Fill(5);
//It corrects the jet momentum due to D daughters out of the jet
RecalculateMomentum(pjet,fPmissing);
fPtJet=TMath::Sqrt(pjet[0]*pjet[0]+pjet[1]*pjet[1]); //recalculated jet pt
-
+ if((jet->Pt()-jet->Area()*rhoval)<0) fPtJet = -fPtJet;
//debugging histograms
Double_t pmissing=TMath::Sqrt(fPmissing[0]*fPmissing[0]+fPmissing[1]*fPmissing[1]+fPmissing[2]*fPmissing[2]); //recalculated jet total momentum
//It corrects the jet momentum due to D daughters out of the jet
RecalculateMomentum(pjet,fPmissing);
fPtJet=TMath::Sqrt(pjet[0]*pjet[0]+pjet[1]*pjet[1]); //recalculated jet pt
-
+ if((jet->Pt()-jet->Area()*rhoval)<0) fPtJet = -fPtJet;
SideBandBackground(sbcand,jet);
}
//It corrects the jet momentum due to D daughters out of the jet
RecalculateMomentum(pjet,fPmissing);
fPtJet=TMath::Sqrt(pjet[0]*pjet[0]+pjet[1]*pjet[1]); //recalculated jet pt
-
+ if((jet->Pt()-jet->Area()*rhoval)<0) fPtJet = -fPtJet;
MCBackground(charmbg,jet);
}
}
//It corrects the each component of the jet momentum for Z calculation
AliRhoParameter *rho = 0;
Double_t rhoval = 0;
- TString sname("Rho");
- if (!sname.IsNull()) {
- rho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(sname));
+ if (!fJetCont->GetRhoName().IsNull()) {
+ rho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fJetCont->GetRhoName()));
if(rho){
rhoval = rho->GetVal();
pj[0] = jet->Px() - jet->Area()*(rhoval*TMath::Cos(jet->AreaPhi()));
fhstat->SetNdivisions(1);
fOutput->Add(fhstat);
- const Int_t nbinsmass=200;
- const Int_t nbinsptjet=500;
+ const Int_t nbinsmass=300;
+ const Int_t nbinspttrack=500;
+ Int_t nbinsptjet=500;
+ if(!fJetCont->GetRhoName().IsNull()) nbinsptjet=400;
const Int_t nbinsptD=100;
const Int_t nbinsz=100;
const Int_t nbinsphi=200;
const Int_t nbinsSpsContrib=100;
const Int_t nbinsSpsA=100;
- const Float_t ptjetlims[2]={0.,200.};
+ const Float_t pttracklims[2]={0.,200.};
+ Float_t ptjetlims[2]={0.,200.};
+ if(!fJetCont->GetRhoName().IsNull()) ptjetlims[0]=-200.;
const Float_t ptDlims[2]={0.,50.};
const Float_t zlims[2]={0.,1.2};
const Float_t philims[2]={0.,6.3};
fhPhiJetTrks->Sumw2();
fhEtaJetTrks = new TH1F("hEtaJetTrks","Jet tracks #eta distribution; #eta", nbinseta,etalims[0],etalims[1]);
fhEtaJetTrks->Sumw2();
- fhPtJetTrks = new TH1F("hPtJetTrks", "Jet tracks Pt distribution; p_{T} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
+ fhPtJetTrks = new TH1F("hPtJetTrks", "Jet tracks Pt distribution; p_{T} (GeV/c)",nbinspttrack,pttracklims[0],pttracklims[1]);
fhPtJetTrks->Sumw2();
fhEjet = new TH1F("hEjet", "Jet energy distribution;Energy (GeV)",500,0,200);
fhNJetPerEv=new TH1F("hNJetPerEv","Number of jets used per event; number of jets/ev",10,-0.5,9.5);
fhNJetPerEv->Sumw2();
+ fhImpPar = new TH2F("hImpPar", "Impact parameter of daughter tracks; Getd0Prong();#it{p}^{daugh}_{T} (GeV/c)",200, -0.1,0.1,nbinsptD, ptDlims[0], ptDlims[1]); //same range of pt of the D, but pt daughter used
+ //NB at the moment fillet with jet track imp par!!!
+ fhVtxX = new TH1F("hVtxX", "Vertex X;vtx_x",100,-0.5,0.5);
+ fhVtxY = new TH1F("hVtxY", "Vertex Y;vtx_y",100,-0.5,0.5);
+ fhVtxZ = new TH1F("hVtxZ", "Vertex Z;vtx_z",100,-10,10);
+ /*
+ if(fUseMCInfo){
+ //understand how to do this. At the moment fhImpPar is filled with the impact parameter of jet tracks even if it is written differently in the fhImpPar definition
+ fhImpParB = new TH2F("hImpParB", "Impact parameter of daughter tracks (Background); Getd0Prong();#it{p}^{daugh}_{T} (GeV/c)",200, -0.1,0.1,nbinsptD, ptDlims[0], ptDlims[1]); //same range of pt of the D, but pt daughter used
+ fOutput->Add(fhImpParB);
+
+ }
+ */
fOutput->Add(fhEjetTrks);
fOutput->Add(fhPhiJetTrks);
fOutput->Add(fhdeltaRJetTracks);
fOutput->Add(fhNtrArr);
fOutput->Add(fhNJetPerEv);
-
+ fOutput->Add(fhImpPar);
+ fOutput->Add(fhVtxX);
+ fOutput->Add(fhVtxY);
+ fOutput->Add(fhVtxZ);
if(fJetOnlyMode && fSwitchOnSparses){
//thnsparse for jets
const Int_t nAxis=6;
}
//background (side bands for the Dstar and like sign for D0)
- fJetRadius=GetJetContainer(0)->GetJetRadius();
+ AliJetContainer *jetCont=GetJetContainer(0);
+ if(!jetCont){
+ Printf("Container 0 not found, try with name %s", fJetArrName.Data());
+ jetCont=GetJetContainer(fJetArrName);
+ if(!jetCont) Printf("NOT FOUND AGAIN");
+ return kFALSE;
+ }
+ Printf("CONTAINER NAME IS %s", jetCont->GetArrayName().Data());
+ //fJetRadius=jetCont->GetJetRadius();
fhInvMassptD = new TH2F("hInvMassptD",Form("D (Delta R < %.1f) invariant mass distribution p_{T}^{j} > threshold",fJetRadius),nbinsmass,fMinMass,fMaxMass,nbinsptD,ptDlims[0],ptDlims[1]);
fhInvMassptD->SetStats(kTRUE);
fhInvMassptD->GetXaxis()->SetTitle("mass (GeV)");
point[6]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
}
+ if(!point){
+ AliError(Form("Numer of THnSparse entries %d not valid", fNAxesBigSparse));
+ return;
+ }
-
- Printf("Candidate in FillHistogramsD0JetCorr IsA %s", (candidate->IsA())->GetName());
+ //Printf("Candidate in FillHistogramsD0JetCorr IsA %s", (candidate->IsA())->GetName());
Int_t isselected=fCuts->IsSelected(candidate,AliRDHFCuts::kAll,aodEvent);
if(isselected==1 || isselected==3) {
if(fNAxesBigSparse==6 || fNAxesBigSparse==7) point[3]=masses[1];
if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) fhsDphiz->Fill(point,1.);
}
-
+ delete[] point;
}
//_______________________________________________________________________________
point[6]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
}
+ if(!point){
+ AliError(Form("Numer of THnSparse entries %d not valid", fNAxesBigSparse));
+ return;
+ }
+
//if(fIsDInJet) hPtJetWithD->Fill(ptj,deltamass,ptD);
FillMassHistograms(deltamass, ptD);
if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) fhsDphiz->Fill(point,1.);
-
+ delete[] point;
}
//_______________________________________________________________________________
point[6]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
}
+ if(!point){
+ AliError(Form("Numer of THnSparse entries %d not valid", fNAxesBigSparse));
+ return;
+ }
if(fNAxesBigSparse==9) point[4]=pdgmass;
//if(fIsDInJet) {
// hPtJetWithD->Fill(ptjet,pdgmass,ptD); // candidates within a jet
//}
-
+ delete[] point;
}
//_______________________________________________________________________________
point[5]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
point[6]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
}
+ if(!point){
+ AliError(Form("Numer of THnSparse entries %d not valid", fNAxesBigSparse));
+ return;
+ }
if(fCandidateType==kDstartoKpipi){
AliAODRecoCascadeHF* dstarbg = (AliAODRecoCascadeHF*)candbg;
}
+ delete[] point;
}
//_______________________________________________________________________________
Double_t phi1=p1->Phi(),eta1=p1->Eta();
//It subtracts the backgroud of jets if it was asked for it.
- TString sname("Rho");
- if (!sname.IsNull()) {
+
+ if (!fJetCont->GetRhoName().IsNull()) {
AliRhoParameter *rho = 0;
- rho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(sname));
+ rho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fJetCont->GetRhoName()));
if(rho){
Double_t pj[3];
Bool_t okpj=p1->PxPyPz(pj);
Int_t nchtrk=thejet->GetNumberOfTracks();
for(Int_t itk=0;itk<nchtrk;itk++){
AliVTrack* tkjet=((AliPicoTrack*)thejet->TrackAt(itk,fTrackArr))->GetTrack();
+ if(!tkjet) continue;
Int_t idtkjet=tkjet->GetID();
if(fillH) fhIDjetTracks->Fill(idtkjet);
for(Int_t id=0;id<ndaugh;id++){
//type 0 : DeltaR < jet radius, don't check daughters (and don't correct pt jet)
//type 1 : DeltaR < jet radius and check for all daughters among jet tracks, don't correct ptjet
//type 2 (default) : DeltaR < jet radius and check for all daughters among jet tracks, if not present, correct the ptjet
+ //type 3 (under development) : DeltaR < jet radius and check for all daughters among jet tracks, if not present, correct the ptjet usign the pt-scheme
fPmissing[0]=0;
fPmissing[1]=0;
Bool_t testDeltaR=kFALSE;
if(DeltaR(thejet,charm)<fJetRadius) testDeltaR=kTRUE;
+ //for type 3 this check should be redone after the modification of the jet axis
Int_t* daughOutOfJet;
AliAODTrack** charmDaugh;
fhzDinjet->Fill(Z(charm,thejet));
}
- if(!testDaugh && testDeltaR && fTypeDInJet==2){
+
+ TVector3 thejetv; //(x=0,y=0,z=0)
+ thejetv.SetPtEtaPhi(thejet->Pt(),thejet->Eta(),thejet->Phi());
+ TVector3 newjet(thejetv);
+ TVector3 correction; //(x=0,y=0,z=0)
+ TVector3 charmcand; //(x=0,y=0,z=0)
+ charmcand.SetPtEtaPhi(charm->Pt(),charm->Eta(),charm->Phi());
+
+ if(!testDaugh && testDeltaR && fTypeDInJet>=2){
Int_t ndaugh=3;
if(fCandidateType==kD0toKpi) ndaugh=2;
if(id==2)fhControlDInJ->Fill(5);
fhDRdaughOut->Fill(DeltaR(thejet, charmDaugh[id]));
}
- fPmissing[0]+=charmDaugh[id]->Px();
- fPmissing[1]+=charmDaugh[id]->Py();
- fPmissing[2]+=charmDaugh[id]->Pz();
- }
-
+ if(fTypeDInJet==2){
+
+ newjet.SetX(newjet(0)+charmDaugh[id]->Px());
+ newjet.SetY(newjet(1)+charmDaugh[id]->Py());
+ newjet.SetZ(newjet(2)+charmDaugh[id]->Pz());
+ }
+ if(fTypeDInJet==3){
+
+ Double_t ptdaug = charmDaugh[id]->Pt();
+ Double_t ptjet = newjet.Perp();
+ Double_t ptn = ptjet+ptdaug;
+ Double_t phidaug = charmDaugh[id]->Phi();
+ Double_t phijet = newjet.Phi();
+ Double_t phin = (phijet*ptjet+phidaug*ptdaug)/(ptjet+ptdaug);
+
+ Double_t etadaug = charmDaugh[id]->Eta();
+ Double_t etajet = newjet.Eta();
+ Double_t etan = (etajet*ptjet+etadaug*ptdaug)/(ptjet+ptdaug);
+
+ newjet.SetPtEtaPhi(ptn,etan,phin);
+
+ }
+ }
}
- //now the D in within the jet
+ correction=newjet-thejetv;
+ fPmissing[0]=correction(0);
+ fPmissing[1]=correction(1);
+ fPmissing[2]=correction(2);
+
+ //now the D is within the jet
testDaugh=kTRUE;
-
+ if(fTypeDInJet==3){ //recalc R(D,jet)
+ Double_t dRnew=newjet.DeltaR(charmcand);
+ if(dRnew<fJetRadius) testDeltaR=kTRUE;
+ else testDeltaR=kFALSE;
+ }
}
-
+ delete[] daughOutOfJet;
delete[] charmDaugh;
Bool_t result=0;
case 1:
result=testDeltaR && testDaugh;
break;
- case 2:
- result=testDeltaR && testDaugh;
+ case 2: //this case defines fPmissing
+ result=testDeltaR && testDaugh;
+ break;
+ case 3: //this case defines fPmissing and recalculate R(jet,D) with the updated jet axis
+ result=testDeltaR && testDaugh;
break;
+
default:
AliInfo("Selection type not specified, use 1");
result=testDeltaR && testDaugh;
break;
}
- return result;
+ return result;
}
Bool_t AliAnalysisTaskFlavourJetCorrelations::InEMCalAcceptance(AliVParticle *vpart){
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff