#include "AliExternalTrackParam.h"
#include "AliAODRecoDecayHF2Prong.h"
#include "AliAODRecoDecayHF3Prong.h"
+#include "AliAODRecoCascadeHF.h"
+#include "AliNeutralTrackParam.h"
#include "AliAnalysisTaskSEImproveITS.h"
//
fPt1ResPUpg (0),
fPt1ResKUpg (0),
fPt1ResPiUpg (0),
+ fD0ZResPCurSA (0),
+ fD0ZResKCurSA (0),
+ fD0ZResPiCurSA (0),
+ fD0RPResPCurSA (0),
+ fD0RPResKCurSA (0),
+ fD0RPResPiCurSA(0),
+ fPt1ResPCurSA (0),
+ fPt1ResKCurSA (0),
+ fPt1ResPiCurSA (0),
+ fD0ZResPUpgSA (0),
+ fD0ZResKUpgSA (0),
+ fD0ZResPiUpgSA (0),
+ fD0RPResPUpgSA (0),
+ fD0RPResKUpgSA (0),
+ fD0RPResPiUpgSA(0),
+ fPt1ResPUpgSA (0),
+ fPt1ResKUpgSA (0),
+ fPt1ResPiUpgSA (0),
fRunInVertexing(kFALSE),
+ fImproveTracks(kTRUE),
fDebugOutput (0),
fDebugNtuple (0),
fDebugVars (0),
fPt1ResPUpg (0),
fPt1ResKUpg (0),
fPt1ResPiUpg (0),
+ fD0ZResPCurSA (0),
+ fD0ZResKCurSA (0),
+ fD0ZResPiCurSA (0),
+ fD0RPResPCurSA (0),
+ fD0RPResKCurSA (0),
+ fD0RPResPiCurSA(0),
+ fPt1ResPCurSA (0),
+ fPt1ResKCurSA (0),
+ fPt1ResPiCurSA (0),
+ fD0ZResPUpgSA (0),
+ fD0ZResKUpgSA (0),
+ fD0ZResPiUpgSA (0),
+ fD0RPResPUpgSA (0),
+ fD0RPResKUpgSA (0),
+ fD0RPResPiUpgSA(0),
+ fPt1ResPUpgSA (0),
+ fPt1ResKUpgSA (0),
+ fPt1ResPiUpgSA (0),
fRunInVertexing(isRunInVertexing),
+ fImproveTracks(kTRUE),
fDebugOutput (0),
fDebugNtuple (0),
fDebugVars (0),
fPt1ResPCur ->SetName("Pt1ResPCur" );
fPt1ResKCur ->SetName("Pt1ResKCur" );
fPt1ResPiCur ->SetName("Pt1ResPiCur" );
+ fD0RPResPCurSA =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResPSA" )));
+ fD0RPResKCurSA =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResKSA" )));
+ fD0RPResPiCurSA=new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResPiSA")));
+ fD0ZResPCurSA =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResPSA" )));
+ fD0ZResKCurSA =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResKSA" )));
+ fD0ZResPiCurSA =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResPiSA" )));
+ fPt1ResPCurSA =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResPSA" )));
+ fPt1ResKCurSA =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResKSA" )));
+ fPt1ResPiCurSA =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResPiSA" )));
+ fD0RPResPCurSA ->SetName("D0RPResPCurSA" );
+ fD0RPResKCurSA ->SetName("D0RPResKCurSA" );
+ fD0RPResPiCurSA->SetName("D0RPResPiCurSA");
+ fD0ZResPCurSA ->SetName("D0ZResPCurSA" );
+ fD0ZResKCurSA ->SetName("D0ZResKCurSA" );
+ fD0ZResPiCurSA ->SetName("D0ZResPiCurSA" );
+ fPt1ResPCurSA ->SetName("Pt1ResPCurSA" );
+ fPt1ResKCurSA ->SetName("Pt1ResKCurSA" );
+ fPt1ResPiCurSA ->SetName("Pt1ResPiCurSA" );
delete resfileCur;
TFile *resfileUpg=TFile::Open(resfileUpgURI );
fD0RPResPUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResP" )));
fPt1ResPUpg ->SetName("Pt1ResPUpg" );
fPt1ResKUpg ->SetName("Pt1ResKUpg" );
fPt1ResPiUpg ->SetName("Pt1ResPiUpg" );
+ fD0RPResPUpgSA =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResPSA" )));
+ fD0RPResKUpgSA =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResKSA" )));
+ fD0RPResPiUpgSA=new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResPiSA")));
+ fD0ZResPUpgSA =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResPSA" )));
+ fD0ZResKUpgSA =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResKSA" )));
+ fD0ZResPiUpgSA =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResPiSA" )));
+ fPt1ResPUpgSA =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResPSA" )));
+ fPt1ResKUpgSA =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResKSA" )));
+ fPt1ResPiUpgSA =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResPiSA" )));
+ fD0RPResPUpgSA ->SetName("D0RPResPUpgSA" );
+ fD0RPResKUpgSA ->SetName("D0RPResKUpgSA" );
+ fD0RPResPiUpgSA->SetName("D0RPResPiUpgSA");
+ fD0ZResPUpgSA ->SetName("D0ZResPUpgSA" );
+ fD0ZResKUpgSA ->SetName("D0ZResKUpgSA" );
+ fD0ZResPiUpgSA ->SetName("D0ZResPiUpgSA" );
+ fPt1ResPUpgSA ->SetName("Pt1ResPUpgSA" );
+ fPt1ResKUpgSA ->SetName("Pt1ResKUpgSA" );
+ fPt1ResPiUpgSA ->SetName("Pt1ResPiUpgSA" );
delete resfileUpg;
DefineOutput(1,TNtuple::Class());
// Smear all tracks
TClonesArray *mcs=static_cast<TClonesArray*>(ev->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
if (!mcs) return;
- for(Int_t itrack=0;itrack<ev->GetNumberOfTracks();++itrack)
- SmearTrack(ev->GetTrack(itrack),mcs);
+ if (fImproveTracks) {
+ for(Int_t itrack=0;itrack<ev->GetNumberOfTracks();++itrack) {
+ SmearTrack(ev->GetTrack(itrack),mcs);
+ }
+ }
// TODO: recalculated primary vertex
AliVVertex *primaryVertex=ev->GetPrimaryVertex();
}
+ // Dstar->Kpipi
+ TClonesArray *arrayCascade=static_cast<TClonesArray*>(ev->GetList()->FindObject("Dstar"));
+
+ if (arrayCascade) {
+ for (Int_t icand=0;icand<arrayCascade->GetEntries();++icand) {
+ AliAODRecoCascadeHF *decayDstar=static_cast<AliAODRecoCascadeHF*>(arrayCascade->At(icand));
+ //Get D0 from D*
+ AliAODRecoDecayHF2Prong* decay=(AliAODRecoDecayHF2Prong*)decayDstar->Get2Prong();
+
+ // recalculate vertices
+ //AliVVertex *oldSecondaryVertex=decay->GetSecondaryVtx();
+
+ //soft pion
+ AliExternalTrackParam et3; et3.CopyFromVTrack(static_cast<AliAODTrack*>(decayDstar->GetBachelor()));
+
+ //track D0
+ AliNeutralTrackParam *trackD0 = new AliNeutralTrackParam(decay);
+
+ //!!!!TODO: covariance matrix
+
+ // update d0
+ Double_t d0z0[2],covd0z0[3];
+ Double_t d01[2],d01err[2];
+
+ //the D*
+ et3.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
+ d01[0]=d0z0[0];
+ d01err[0] = TMath::Sqrt(covd0z0[0]);
+ trackD0->PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
+ d01[1]=d0z0[0];
+ d01err[1] = TMath::Sqrt(covd0z0[0]);
+ decayDstar->Setd0Prongs(2,d01);
+ decayDstar->Setd0errProngs(2,d01err);
+
+ // delete v12;
+ delete trackD0; trackD0=NULL;
+
+ // a run for D*
+ Double_t px1[2],py1[2],pz1[2];
+ for (Int_t i=0;i<2;++i) {
+ const AliAODTrack *t1=static_cast<AliAODTrack*>(decayDstar->GetDaughter(i));
+ px1[i]=t1->Px();
+ py1[i]=t1->Py();
+ pz1[i]=t1->Pz();
+ }
+ decayDstar->SetPxPyPzProngs(2,px1,py1,pz1);
+
+ }
+ }
+
// Three prong
TClonesArray *array3Prong=static_cast<TClonesArray*>(ev->GetList()->FindObject("Charm3Prong"));
if (!(track->HasPointOnITSLayer(0) || track->HasPointOnITSLayer(1)))
return;
+ // Check if the track was already "improved" (this is done with a trick using layer 7 (ie the 8th))
+ if (TESTBIT(track->GetITSClusterMap(),7)) return;
+ //
+
+
// Get reconstructed track parameters
AliExternalTrackParam et; et.CopyFromVTrack(track);
Double_t *param=const_cast<Double_t*>(et.GetParameter());
Double_t spt1o =0.;
switch (mc->GetPdgCode()) {
case 2212: case -2212:
- sd0rpo=EvalGraph(fD0RPResPCur,ptmc);
- sd0zo =EvalGraph(fD0ZResPCur ,ptmc);
- spt1o =EvalGraph(fPt1ResPCur ,ptmc);
- sd0rpn=EvalGraph(fD0RPResPUpg,ptmc);
- sd0zn =EvalGraph(fD0ZResPUpg ,ptmc);
- spt1n =EvalGraph(fPt1ResPUpg ,ptmc);
+ sd0rpo=EvalGraph(ptmc,fD0RPResPCur,fD0RPResPCurSA);
+ sd0zo =EvalGraph(ptmc,fD0ZResPCur,fD0ZResPCurSA);
+ spt1o =EvalGraph(ptmc,fPt1ResPCur,fPt1ResPCurSA);
+ sd0rpn=EvalGraph(ptmc,fD0RPResPUpg,fD0RPResPUpgSA);
+ sd0zn =EvalGraph(ptmc,fD0ZResPUpg,fD0ZResPUpgSA);
+ spt1n =EvalGraph(ptmc,fPt1ResPUpg,fPt1ResPUpgSA);
break;
case 321: case -321:
- sd0rpo=EvalGraph(fD0RPResKCur,ptmc);
- sd0zo =EvalGraph(fD0ZResKCur ,ptmc);
- spt1o =EvalGraph(fPt1ResKCur ,ptmc);
- sd0rpn=EvalGraph(fD0RPResKUpg,ptmc);
- sd0zn =EvalGraph(fD0ZResKUpg ,ptmc);
- spt1n =EvalGraph(fPt1ResKUpg ,ptmc);
+ sd0rpo=EvalGraph(ptmc,fD0RPResKCur,fD0RPResKCurSA);
+ sd0zo =EvalGraph(ptmc,fD0ZResKCur,fD0ZResKCurSA);
+ spt1o =EvalGraph(ptmc,fPt1ResKCur,fPt1ResKCurSA);
+ sd0rpn=EvalGraph(ptmc,fD0RPResKUpg,fD0RPResKUpgSA);
+ sd0zn =EvalGraph(ptmc,fD0ZResKUpg,fD0ZResKUpgSA);
+ spt1n =EvalGraph(ptmc,fPt1ResKUpg,fPt1ResKUpgSA);
break;
case 211: case -211:
- sd0rpo=EvalGraph(fD0RPResPiCur,ptmc);
- sd0zo =EvalGraph(fD0ZResPiCur ,ptmc);
- spt1o =EvalGraph(fPt1ResPiCur ,ptmc);
- sd0rpn=EvalGraph(fD0RPResPiUpg,ptmc);
- sd0zn =EvalGraph(fD0ZResPiUpg ,ptmc);
- spt1n =EvalGraph(fPt1ResPiUpg ,ptmc);
+ sd0rpo=EvalGraph(ptmc,fD0RPResPiCur,fD0RPResPiCurSA);
+ sd0zo =EvalGraph(ptmc,fD0ZResPiCur,fD0ZResPiCurSA);
+ spt1o =EvalGraph(ptmc,fPt1ResPiCur,fPt1ResPiCurSA);
+ sd0rpn=EvalGraph(ptmc,fD0RPResPiUpg,fD0RPResPiUpgSA);
+ sd0zn =EvalGraph(ptmc,fD0ZResPiUpg,fD0ZResPiUpgSA);
+ spt1n =EvalGraph(ptmc,fPt1ResPiUpg,fPt1ResPiUpgSA);
break;
default:
return;
track->SetPosition(x,kFALSE);
track->SetP(p,kTRUE);
+
+ // Mark the track as "improved" with a trick (this is done with a trick using layer 7 (ie the 8th))
+ UChar_t itsClusterMap = track->GetITSClusterMap();
+ SETBIT(itsClusterMap,7);
+ track->SetITSClusterMap(itsClusterMap);
+ //
+
// write out debug infos
if (fDebugNtuple->GetEntries()<fNDebug) {
Int_t idbg=0;
return vertex;
}
-Double_t AliAnalysisTaskSEImproveITS::EvalGraph(const TGraph *graph,Double_t x) const {
+Double_t AliAnalysisTaskSEImproveITS::EvalGraph(Double_t x,const TGraph *graph,const TGraph *graphSA) const {
//
// Evaluates a TGraph without linear extrapolation. Instead the last
// valid point of the graph is used when out of range.
Int_t n =graph->GetN();
Double_t xmin=graph->GetX()[0 ];
Double_t xmax=graph->GetX()[n-1];
- if (x<xmin) return graph->Eval(xmin);
+ if (x<xmin) {
+ if(!graphSA) return graph->Eval(xmin);
+ Double_t xminSA=graphSA->GetX()[0];
+ if(x<xminSA) return graphSA->Eval(xminSA);
+ return graphSA->Eval(x);
+ }
if (x>xmax) return graph->Eval(xmax);
return graph->Eval(x);
}