-/*************************************************************************\r
-* Copyright(c) 1998-2011, ALICE Experiment at CERN, All rights reserved. *\r
-* *\r
-* Author: The ALICE Off-line Project. *\r
-* Contributors are mentioned in the code where appropriate. *\r
-* *\r
-* Permission to use, copy, modify and distribute this software and its *\r
-* documentation strictly for non-commercial purposes is hereby granted *\r
-* without fee, provided that the above copyright notice appears in all *\r
-* copies and that both the copyright notice and this permission notice *\r
-* appear in the supporting documentation. The authors make no claims *\r
-* about the suitability of this software for any purpose. It is *\r
-* provided "as is" without express or implied warranty. * \r
-**************************************************************************/\r
-\r
-#include <TObjArray.h>\r
-#include <TClonesArray.h>\r
-#include <TGraph.h>\r
-#include <TFile.h>\r
-#include <TList.h>\r
-#include <TNtuple.h>\r
-\r
-#include "AliVertex.h"\r
-#include "AliVVertex.h"\r
-#include "AliESDVertex.h"\r
-#include "AliVertexerTracks.h"\r
-#include "AliAODEvent.h"\r
-#include "AliAODTrack.h"\r
-#include "AliAODMCParticle.h"\r
-#include "AliExternalTrackParam.h"\r
-#include "AliAODRecoDecayHF2Prong.h"\r
-#include "AliAODRecoDecayHF3Prong.h"\r
-#include "AliAODRecoCascadeHF.h"\r
-#include "AliNeutralTrackParam.h"\r
-#include "AliAnalysisTaskSEImproveITS.h"\r
-\r
-//\r
-// Implementation of the "hybrid-approach" for ITS upgrade studies.\r
-// The tastk smears the track parameters according to estimations\r
-// from single-track upgrade studies. Afterwards it recalculates\r
-// the parameters of the reconstructed decays.\r
-//\r
-// WARNING: This will affect all tasks in a train after this one\r
-// (which is typically desired, though).\r
-//\r
-\r
-AliAnalysisTaskSEImproveITS::AliAnalysisTaskSEImproveITS() \r
- :AliAnalysisTaskSE(),\r
- fD0ZResPCur (0),\r
- fD0ZResKCur (0),\r
- fD0ZResPiCur (0),\r
- fD0RPResPCur (0),\r
- fD0RPResKCur (0),\r
- fD0RPResPiCur(0),\r
- fPt1ResPCur (0),\r
- fPt1ResKCur (0),\r
- fPt1ResPiCur (0),\r
- fD0ZResPUpg (0),\r
- fD0ZResKUpg (0),\r
- fD0ZResPiUpg (0),\r
- fD0RPResPUpg (0),\r
- fD0RPResKUpg (0),\r
- fD0RPResPiUpg(0),\r
- fPt1ResPUpg (0),\r
- fPt1ResKUpg (0),\r
- fPt1ResPiUpg (0),\r
- fD0ZResPCurSA (0),\r
- fD0ZResKCurSA (0),\r
- fD0ZResPiCurSA (0),\r
- fD0RPResPCurSA (0),\r
- fD0RPResKCurSA (0),\r
- fD0RPResPiCurSA(0),\r
- fPt1ResPCurSA (0),\r
- fPt1ResKCurSA (0),\r
- fPt1ResPiCurSA (0),\r
- fD0ZResPUpgSA (0),\r
- fD0ZResKUpgSA (0),\r
- fD0ZResPiUpgSA (0),\r
- fD0RPResPUpgSA (0),\r
- fD0RPResKUpgSA (0),\r
- fD0RPResPiUpgSA(0),\r
- fPt1ResPUpgSA (0),\r
- fPt1ResKUpgSA (0),\r
- fPt1ResPiUpgSA (0),\r
- fRunInVertexing(kFALSE),\r
- fImproveTracks(kTRUE),\r
- fDebugOutput (0),\r
- fDebugNtuple (0),\r
- fDebugVars (0), \r
- fNDebug (0)\r
-{\r
- //\r
- // Default constructor.\r
- //\r
-}\r
-\r
-AliAnalysisTaskSEImproveITS::AliAnalysisTaskSEImproveITS(const char *name,\r
- const char *resfileCurURI,\r
- const char *resfileUpgURI,\r
- Bool_t isRunInVertexing,\r
- Int_t ndebug)\r
- :AliAnalysisTaskSE(name),\r
- fD0ZResPCur (0),\r
- fD0ZResKCur (0),\r
- fD0ZResPiCur (0),\r
- fD0RPResPCur (0),\r
- fD0RPResKCur (0),\r
- fD0RPResPiCur(0),\r
- fPt1ResPCur (0),\r
- fPt1ResKCur (0),\r
- fPt1ResPiCur (0),\r
- fD0ZResPUpg (0),\r
- fD0ZResKUpg (0),\r
- fD0ZResPiUpg (0),\r
- fD0RPResPUpg (0),\r
- fD0RPResKUpg (0),\r
- fD0RPResPiUpg(0),\r
- fPt1ResPUpg (0),\r
- fPt1ResKUpg (0),\r
- fPt1ResPiUpg (0),\r
- fD0ZResPCurSA (0),\r
- fD0ZResKCurSA (0),\r
- fD0ZResPiCurSA (0),\r
- fD0RPResPCurSA (0),\r
- fD0RPResKCurSA (0),\r
- fD0RPResPiCurSA(0),\r
- fPt1ResPCurSA (0),\r
- fPt1ResKCurSA (0),\r
- fPt1ResPiCurSA (0),\r
- fD0ZResPUpgSA (0),\r
- fD0ZResKUpgSA (0),\r
- fD0ZResPiUpgSA (0),\r
- fD0RPResPUpgSA (0),\r
- fD0RPResKUpgSA (0),\r
- fD0RPResPiUpgSA(0),\r
- fPt1ResPUpgSA (0),\r
- fPt1ResKUpgSA (0),\r
- fPt1ResPiUpgSA (0),\r
- fRunInVertexing(isRunInVertexing),\r
- fImproveTracks(kTRUE),\r
- fDebugOutput (0),\r
- fDebugNtuple (0),\r
- fDebugVars (0),\r
- fNDebug (ndebug)\r
-{\r
- //\r
- // Constructor to be used to create the task.\r
- // The the URIs specify the resolution files to be used. \r
- // They are expected to contain TGraphs with the resolutions\r
- // for the current and the upgraded ITS (see code for details).\r
- // One may also specify for how many tracks debug information\r
- // is written to the output.\r
- //\r
- TFile *resfileCur=TFile::Open(resfileCurURI);\r
- fD0RPResPCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResP" )));\r
- fD0RPResKCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResK" )));\r
- fD0RPResPiCur=new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResPi")));\r
- fD0ZResPCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResP" )));\r
- fD0ZResKCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResK" )));\r
- fD0ZResPiCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResPi" )));\r
- fPt1ResPCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResP" )));\r
- fPt1ResKCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResK" )));\r
- fPt1ResPiCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResPi" )));\r
- fD0RPResPCur ->SetName("D0RPResPCur" );\r
- fD0RPResKCur ->SetName("D0RPResKCur" );\r
- fD0RPResPiCur->SetName("D0RPResPiCur");\r
- fD0ZResPCur ->SetName("D0ZResPCur" ); \r
- fD0ZResKCur ->SetName("D0ZResKCur" );\r
- fD0ZResPiCur ->SetName("D0ZResPiCur" );\r
- fPt1ResPCur ->SetName("Pt1ResPCur" );\r
- fPt1ResKCur ->SetName("Pt1ResKCur" );\r
- fPt1ResPiCur ->SetName("Pt1ResPiCur" );\r
- fD0RPResPCurSA =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResPSA" )));\r
- fD0RPResKCurSA =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResKSA" )));\r
- fD0RPResPiCurSA=new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResPiSA")));\r
- fD0ZResPCurSA =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResPSA" )));\r
- fD0ZResKCurSA =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResKSA" )));\r
- fD0ZResPiCurSA =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResPiSA" )));\r
- fPt1ResPCurSA =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResPSA" )));\r
- fPt1ResKCurSA =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResKSA" )));\r
- fPt1ResPiCurSA =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResPiSA" )));\r
- fD0RPResPCurSA ->SetName("D0RPResPCurSA" );\r
- fD0RPResKCurSA ->SetName("D0RPResKCurSA" );\r
- fD0RPResPiCurSA->SetName("D0RPResPiCurSA");\r
- fD0ZResPCurSA ->SetName("D0ZResPCurSA" ); \r
- fD0ZResKCurSA ->SetName("D0ZResKCurSA" );\r
- fD0ZResPiCurSA ->SetName("D0ZResPiCurSA" );\r
- fPt1ResPCurSA ->SetName("Pt1ResPCurSA" );\r
- fPt1ResKCurSA ->SetName("Pt1ResKCurSA" );\r
- fPt1ResPiCurSA ->SetName("Pt1ResPiCurSA" );\r
- delete resfileCur;\r
- TFile *resfileUpg=TFile::Open(resfileUpgURI );\r
- fD0RPResPUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResP" )));\r
- fD0RPResKUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResK" )));\r
- fD0RPResPiUpg=new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResPi")));\r
- fD0ZResPUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResP" )));\r
- fD0ZResKUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResK" )));\r
- fD0ZResPiUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResPi" )));\r
- fPt1ResPUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResP" )));\r
- fPt1ResKUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResK" )));\r
- fPt1ResPiUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResPi" )));\r
- fD0RPResPUpg ->SetName("D0RPResPUpg" );\r
- fD0RPResKUpg ->SetName("D0RPResKUpg" );\r
- fD0RPResPiUpg->SetName("D0RPResPiUpg");\r
- fD0ZResPUpg ->SetName("D0ZResPUpg" );\r
- fD0ZResKUpg ->SetName("D0ZResKUpg" );\r
- fD0ZResPiUpg ->SetName("D0ZResPiUpg" );\r
- fPt1ResPUpg ->SetName("Pt1ResPUpg" );\r
- fPt1ResKUpg ->SetName("Pt1ResKUpg" );\r
- fPt1ResPiUpg ->SetName("Pt1ResPiUpg" );\r
- fD0RPResPUpgSA =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResPSA" )));\r
- fD0RPResKUpgSA =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResKSA" )));\r
- fD0RPResPiUpgSA=new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResPiSA")));\r
- fD0ZResPUpgSA =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResPSA" )));\r
- fD0ZResKUpgSA =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResKSA" )));\r
- fD0ZResPiUpgSA =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResPiSA" )));\r
- fPt1ResPUpgSA =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResPSA" )));\r
- fPt1ResKUpgSA =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResKSA" )));\r
- fPt1ResPiUpgSA =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResPiSA" )));\r
- fD0RPResPUpgSA ->SetName("D0RPResPUpgSA" );\r
- fD0RPResKUpgSA ->SetName("D0RPResKUpgSA" );\r
- fD0RPResPiUpgSA->SetName("D0RPResPiUpgSA");\r
- fD0ZResPUpgSA ->SetName("D0ZResPUpgSA" );\r
- fD0ZResKUpgSA ->SetName("D0ZResKUpgSA" );\r
- fD0ZResPiUpgSA ->SetName("D0ZResPiUpgSA" );\r
- fPt1ResPUpgSA ->SetName("Pt1ResPUpgSA" );\r
- fPt1ResKUpgSA ->SetName("Pt1ResKUpgSA" );\r
- fPt1ResPiUpgSA ->SetName("Pt1ResPiUpgSA" );\r
- delete resfileUpg;\r
-\r
- DefineOutput(1,TNtuple::Class());\r
-}\r
-\r
-AliAnalysisTaskSEImproveITS::~AliAnalysisTaskSEImproveITS() {\r
- //\r
- // Destructor.\r
- //\r
- delete fDebugOutput;\r
-}\r
-\r
-void AliAnalysisTaskSEImproveITS::UserCreateOutputObjects() {\r
- //\r
- // Creation of user output objects.\r
- //\r
- fDebugOutput=new TList();\r
- fDebugOutput->SetOwner();\r
- fDebugNtuple=new TNtuple("fDebugNtuple","Smearing","pdg:ptmc:d0rpo:d0zo:pt1o:sd0rpo:sd0zo:spt1o:d0rpn:d0zn:pt1n:sd0rpn:sd0zn:spt1n:d0rpmc:d0zmc:pt1mc");\r
- fDebugVars=new Float_t[fDebugNtuple->GetNvar()];\r
- \r
- fDebugOutput->Add(fDebugNtuple );\r
-\r
- fDebugOutput->Add(fD0RPResPCur );\r
- fDebugOutput->Add(fD0RPResKCur );\r
- fDebugOutput->Add(fD0RPResPiCur);\r
- fDebugOutput->Add(fD0ZResPCur ); \r
- fDebugOutput->Add(fD0ZResKCur );\r
- fDebugOutput->Add(fD0ZResPiCur );\r
- fDebugOutput->Add(fPt1ResPCur );\r
- fDebugOutput->Add(fPt1ResKCur );\r
- fDebugOutput->Add(fPt1ResPiCur );\r
- fDebugOutput->Add(fD0RPResPUpg );\r
- fDebugOutput->Add(fD0RPResKUpg );\r
- fDebugOutput->Add(fD0RPResPiUpg);\r
- fDebugOutput->Add(fD0ZResPUpg );\r
- fDebugOutput->Add(fD0ZResKUpg );\r
- fDebugOutput->Add(fD0ZResPiUpg );\r
- fDebugOutput->Add(fPt1ResPUpg );\r
- fDebugOutput->Add(fPt1ResKUpg );\r
- fDebugOutput->Add(fPt1ResPiUpg );\r
-\r
- PostData(1,fDebugOutput);\r
-}\r
-\r
-void AliAnalysisTaskSEImproveITS::UserExec(Option_t*) {\r
- //\r
- // The event loop\r
- //\r
- AliAODEvent *ev=0x0;\r
- if(!fRunInVertexing) {\r
- ev=dynamic_cast<AliAODEvent*>(InputEvent());\r
- } else {\r
- if(AODEvent() && IsStandardAOD()) ev = dynamic_cast<AliAODEvent*> (AODEvent());\r
- } \r
- if(!ev) return;\r
- Double_t bz=ev->GetMagneticField();\r
-\r
-\r
-\r
-\r
- // Smear all tracks\r
- TClonesArray *mcs=static_cast<TClonesArray*>(ev->GetList()->FindObject(AliAODMCParticle::StdBranchName()));\r
- if (!mcs) return;\r
- if (fImproveTracks) {\r
- for(Int_t itrack=0;itrack<ev->GetNumberOfTracks();++itrack) {\r
- SmearTrack(ev->GetTrack(itrack),mcs);\r
- }\r
- }\r
-\r
- // TODO: recalculated primary vertex\r
- AliVVertex *primaryVertex=ev->GetPrimaryVertex();\r
-\r
- // Recalculate all candidates\r
- // D0->Kpi\r
- TClonesArray *array2Prong=static_cast<TClonesArray*>(ev->GetList()->FindObject("D0toKpi"));\r
- if (array2Prong) {\r
- for (Int_t icand=0;icand<array2Prong->GetEntries();++icand) {\r
- AliAODRecoDecayHF2Prong *decay=static_cast<AliAODRecoDecayHF2Prong*>(array2Prong->At(icand));\r
-\r
- // recalculate vertices\r
- AliVVertex *oldSecondaryVertex=decay->GetSecondaryVtx();\r
-\r
-\r
- AliExternalTrackParam et1; et1.CopyFromVTrack(static_cast<AliAODTrack*>(decay->GetDaughter(0)));\r
- AliExternalTrackParam et2; et2.CopyFromVTrack(static_cast<AliAODTrack*>(decay->GetDaughter(1)));\r
-\r
- TObjArray ta12;\r
- \r
- ta12.Add(&et1); ta12.Add(&et2); \r
- AliESDVertex *v12 =RecalculateVertex(oldSecondaryVertex,&ta12 ,bz);\r
- \r
-\r
- // update secondary vertex\r
- Double_t pos[3];\r
- v12->GetXYZ(pos);\r
- \r
- decay->GetSecondaryVtx()->SetPosition(pos[0],pos[1],pos[2]);\r
- decay->GetSecondaryVtx()->SetChi2perNDF(v12->GetChi2toNDF()); \r
- \r
- //!!!!TODO: covariance matrix\r
-\r
- // update d0 \r
- Double_t d0z0[2],covd0z0[3];\r
- Double_t d0[2],d0err[2];\r
- et1.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);\r
- d0[0]=d0z0[0];\r
- d0err[0] = TMath::Sqrt(covd0z0[0]);\r
- et2.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);\r
- d0[1]=d0z0[0];\r
- d0err[1] = TMath::Sqrt(covd0z0[0]); \r
- decay->Setd0Prongs(2,d0);\r
- decay->Setd0errProngs(2,d0err);\r
- // \r
-\r
-\r
- Double_t xdummy=0.,ydummy=0.;\r
- Double_t dca;\r
- dca=et1.GetDCA(&et2,bz,xdummy,ydummy);\r
- decay->SetDCA(dca);\r
-\r
- \r
- delete v12;\r
-\r
- Double_t px[2],py[2],pz[2];\r
- for (Int_t i=0;i<2;++i) {\r
- const AliAODTrack *t=static_cast<AliAODTrack*>(decay->GetDaughter(i));\r
- px[i]=t->Px();\r
- py[i]=t->Py();\r
- pz[i]=t->Pz();\r
- }\r
- decay->SetPxPyPzProngs(2,px,py,pz);\r
- }\r
- }\r
-\r
-\r
- // Dstar->Kpipi\r
- TClonesArray *arrayCascade=static_cast<TClonesArray*>(ev->GetList()->FindObject("Dstar"));\r
- \r
- if (arrayCascade) {\r
- for (Int_t icand=0;icand<arrayCascade->GetEntries();++icand) {\r
- AliAODRecoCascadeHF *decayDstar=static_cast<AliAODRecoCascadeHF*>(arrayCascade->At(icand));\r
- //Get D0 from D*\r
- AliAODRecoDecayHF2Prong* decay=(AliAODRecoDecayHF2Prong*)decayDstar->Get2Prong();\r
- \r
- // recalculate vertices\r
- //AliVVertex *oldSecondaryVertex=decay->GetSecondaryVtx();\r
-\r
- //soft pion\r
- AliExternalTrackParam et3; et3.CopyFromVTrack(static_cast<AliAODTrack*>(decayDstar->GetBachelor()));\r
- \r
- //track D0\r
- AliNeutralTrackParam *trackD0 = new AliNeutralTrackParam(decay);\r
-\r
- //!!!!TODO: covariance matrix\r
- \r
- // update d0 \r
- Double_t d0z0[2],covd0z0[3];\r
- Double_t d01[2],d01err[2];\r
-\r
- //the D*\r
- et3.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);\r
- d01[0]=d0z0[0];\r
- d01err[0] = TMath::Sqrt(covd0z0[0]); \r
- trackD0->PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);\r
- d01[1]=d0z0[0];\r
- d01err[1] = TMath::Sqrt(covd0z0[0]); \r
- decayDstar->Setd0Prongs(2,d01);\r
- decayDstar->Setd0errProngs(2,d01err);\r
- \r
- // delete v12;\r
- delete trackD0; trackD0=NULL;\r
-\r
- // a run for D*\r
- Double_t px1[2],py1[2],pz1[2];\r
- for (Int_t i=0;i<2;++i) {\r
- const AliAODTrack *t1=static_cast<AliAODTrack*>(decayDstar->GetDaughter(i));\r
- px1[i]=t1->Px();\r
- py1[i]=t1->Py();\r
- pz1[i]=t1->Pz();\r
- }\r
- decayDstar->SetPxPyPzProngs(2,px1,py1,pz1);\r
- \r
- }\r
- }\r
-\r
-\r
- // Three prong\r
- TClonesArray *array3Prong=static_cast<TClonesArray*>(ev->GetList()->FindObject("Charm3Prong"));\r
- if (array3Prong) {\r
- for (Int_t icand=0;icand<array3Prong->GetEntries();++icand) {\r
- AliAODRecoDecayHF3Prong *decay=static_cast<AliAODRecoDecayHF3Prong*>(array3Prong->At(icand));\r
-\r
- // recalculate vertices\r
- AliVVertex *oldSecondaryVertex=decay->GetSecondaryVtx();\r
- AliExternalTrackParam et1; et1.CopyFromVTrack(static_cast<AliAODTrack*>(decay->GetDaughter(0)));\r
- AliExternalTrackParam et2; et2.CopyFromVTrack(static_cast<AliAODTrack*>(decay->GetDaughter(1)));\r
- AliExternalTrackParam et3; et3.CopyFromVTrack(static_cast<AliAODTrack*>(decay->GetDaughter(2)));\r
- TObjArray ta123,ta12,ta23;\r
- ta123.Add(&et1);ta123.Add(&et2);ta123.Add(&et3);\r
- ta12. Add(&et1);ta12 .Add(&et2);\r
- ta23 .Add(&et2);ta23 .Add(&et3);\r
- AliESDVertex *v123=RecalculateVertex(oldSecondaryVertex,&ta123,bz);\r
- AliESDVertex *v12 =RecalculateVertex(oldSecondaryVertex,&ta12 ,bz);\r
- AliESDVertex *v23 =RecalculateVertex(oldSecondaryVertex,&ta23 ,bz);\r
-\r
- // update secondary vertex\r
- Double_t pos[3];\r
- v123->GetXYZ(pos);\r
- decay->GetSecondaryVtx()->SetPosition(pos[0],pos[1],pos[2]);\r
- decay->GetSecondaryVtx()->SetChi2perNDF(v123->GetChi2toNDF()); \r
- //TODO: covariance matrix\r
-\r
- // update d0 for all progs\r
- Double_t d0z0[2],covd0z0[3];\r
- Double_t d0[3],d0err[3];\r
- et1.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);\r
- d0[0]=d0z0[0];\r
- d0err[0] = TMath::Sqrt(covd0z0[0]);\r
- et2.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);\r
- d0[1]=d0z0[0];\r
- d0err[1] = TMath::Sqrt(covd0z0[0]);\r
- et3.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);\r
- d0[2]=d0z0[0];\r
- d0err[2] = TMath::Sqrt(covd0z0[0]);\r
- decay->Setd0Prongs (3,d0 );\r
- decay->Setd0errProngs(3,d0err);\r
- // TODO: setter missing\r
-\r
- // update dca for prong combinations\r
- Double_t xdummy=0.,ydummy=0.;\r
- Double_t dca[3];\r
- dca[0]=et1.GetDCA(&et2,bz,xdummy,ydummy);\r
- dca[1]=et3.GetDCA(&et2,bz,xdummy,ydummy);\r
- dca[2]=et1.GetDCA(&et3,bz,xdummy,ydummy);\r
- decay->SetDCAs(3,dca);\r
-\r
- // update dist12 and dist23\r
- primaryVertex->GetXYZ(pos);\r
- decay->SetDist12toPrim(TMath::Sqrt((v12->GetX()-pos[0])*(v12->GetX()-pos[0])\r
- +(v12->GetY()-pos[1])*(v12->GetY()-pos[1])\r
- +(v12->GetZ()-pos[2])*(v12->GetZ()-pos[2])));\r
- decay->SetDist23toPrim(TMath::Sqrt((v23->GetX()-pos[0])*(v23->GetX()-pos[0])\r
- +(v23->GetY()-pos[1])*(v23->GetY()-pos[1])\r
- +(v23->GetZ()-pos[2])*(v23->GetZ()-pos[2])));\r
- \r
- delete v123;delete v12;delete v23;\r
-\r
- Double_t px[3],py[3],pz[3];\r
- for (Int_t i=0;i<3;++i) {\r
- const AliAODTrack *t=static_cast<AliAODTrack*>(decay->GetDaughter(i));\r
- px[i]=t->Px();\r
- py[i]=t->Py();\r
- pz[i]=t->Pz();\r
- }\r
- decay->SetPxPyPzProngs(3,px,py,pz);\r
- }\r
- }\r
-}\r
-\r
-void AliAnalysisTaskSEImproveITS::SmearTrack(AliAODTrack *track,const TClonesArray *mcs) {\r
- // Early exit, if this track has nothing in common with the ITS\r
-\r
- if (!(track->HasPointOnITSLayer(0) || track->HasPointOnITSLayer(1)))\r
- return;\r
-\r
- // Check if the track was already "improved" (this is done with a trick using layer 7 (ie the 8th))\r
- if (TESTBIT(track->GetITSClusterMap(),7)) return;\r
- //\r
-\r
-\r
- // Get reconstructed track parameters\r
- AliExternalTrackParam et; et.CopyFromVTrack(track);\r
- Double_t *param=const_cast<Double_t*>(et.GetParameter());\r
-//TODO: Double_t *covar=const_cast<Double_t*>(et.GetCovariance());\r
-\r
- // Get MC info\r
- Int_t imc=track->GetLabel();\r
- if (imc<=0) return;\r
- const AliAODMCParticle *mc=static_cast<AliAODMCParticle*>(mcs->At(imc));\r
- Double_t mcx[3];\r
- Double_t mcp[3];\r
- Double_t mccv[36]={0.};\r
- Short_t mcc;\r
- mc->XvYvZv(mcx);\r
- mc->PxPyPz(mcp);\r
- mcc=mc->Charge();\r
- AliExternalTrackParam mct(mcx,mcp,mccv,mcc);\r
- const Double_t *parammc=mct.GetParameter();\r
-//TODO: const Double_t *covermc=mct.GetCovariance();\r
- AliVertex vtx(mcx,1.,1);\r
-\r
- // Correct reference points and frames according to MC\r
- // TODO: B-Field correct?\r
- // TODO: failing propagation....\r
- et.PropagateToDCA(&vtx,track->GetBz(),10.);\r
- et.Rotate(mct.GetAlpha());\r
-\r
- // Select appropriate smearing functions\r
- Double_t ptmc=TMath::Abs(mc->Pt());\r
- Double_t sd0rpn=0.;\r
- Double_t sd0zn =0.;\r
- Double_t spt1n =0.;\r
- Double_t sd0rpo=0.;\r
- Double_t sd0zo =0.;\r
- Double_t spt1o =0.;\r
- switch (mc->GetPdgCode()) {\r
- case 2212: case -2212:\r
- sd0rpo=EvalGraph(ptmc,fD0RPResPCur,fD0RPResPCurSA);\r
- sd0zo =EvalGraph(ptmc,fD0ZResPCur,fD0ZResPCurSA);\r
- spt1o =EvalGraph(ptmc,fPt1ResPCur,fPt1ResPCurSA);\r
- sd0rpn=EvalGraph(ptmc,fD0RPResPUpg,fD0RPResPUpgSA);\r
- sd0zn =EvalGraph(ptmc,fD0ZResPUpg,fD0ZResPUpgSA);\r
- spt1n =EvalGraph(ptmc,fPt1ResPUpg,fPt1ResPUpgSA);\r
- break;\r
- case 321: case -321:\r
- sd0rpo=EvalGraph(ptmc,fD0RPResKCur,fD0RPResKCurSA);\r
- sd0zo =EvalGraph(ptmc,fD0ZResKCur,fD0ZResKCurSA);\r
- spt1o =EvalGraph(ptmc,fPt1ResKCur,fPt1ResKCurSA);\r
- sd0rpn=EvalGraph(ptmc,fD0RPResKUpg,fD0RPResKUpgSA);\r
- sd0zn =EvalGraph(ptmc,fD0ZResKUpg,fD0ZResKUpgSA);\r
- spt1n =EvalGraph(ptmc,fPt1ResKUpg,fPt1ResKUpgSA);\r
- break;\r
- case 211: case -211:\r
- sd0rpo=EvalGraph(ptmc,fD0RPResPiCur,fD0RPResPiCurSA);\r
- sd0zo =EvalGraph(ptmc,fD0ZResPiCur,fD0ZResPiCurSA);\r
- spt1o =EvalGraph(ptmc,fPt1ResPiCur,fPt1ResPiCurSA);\r
- sd0rpn=EvalGraph(ptmc,fD0RPResPiUpg,fD0RPResPiUpgSA);\r
- sd0zn =EvalGraph(ptmc,fD0ZResPiUpg,fD0ZResPiUpgSA);\r
- spt1n =EvalGraph(ptmc,fPt1ResPiUpg,fPt1ResPiUpgSA);\r
- break;\r
- default:\r
- return;\r
- }\r
-\r
- // Use the same units (i.e. cm and GeV/c)! TODO: pt!\r
- sd0rpo*=1.e-4;\r
- sd0zo *=1.e-4;\r
- sd0rpn*=1.e-4;\r
- sd0zn *=1.e-4;\r
-\r
- // Apply the smearing\r
- Double_t d0zo =param [1];\r
- Double_t d0zmc =parammc[1];\r
- Double_t d0rpo =param [0];\r
- Double_t d0rpmc=parammc[0];\r
- Double_t pt1o =param [4];\r
- Double_t pt1mc =parammc[4];\r
- Double_t dd0zo =d0zo-d0zmc;\r
- Double_t dd0zn =dd0zo *(sd0zo >0. ? (sd0zn /sd0zo ) : 1.);\r
- Double_t d0zn =d0zmc+dd0zn;\r
- Double_t dd0rpo=d0rpo-d0rpmc;\r
- Double_t dd0rpn=dd0rpo*(sd0rpo>0. ? (sd0rpn/sd0rpo) : 1.);\r
- Double_t d0rpn =d0rpmc+dd0rpn;\r
- Double_t dpt1o =pt1o-pt1mc;\r
- Double_t dpt1n =dpt1o *(spt1o >0. ? (spt1n /spt1o ) : 1.);\r
- Double_t pt1n =pt1mc+dpt1n;\r
- param[0]=d0rpn;\r
- param[1]=d0zn ;\r
- param[4]=pt1n ;\r
-\r
- // Copy the smeared parameters to the AOD track\r
- Double_t x[3];\r
- Double_t p[3];\r
- et.GetXYZ(x);\r
- et.GetPxPyPz(p);\r
- track->SetPosition(x,kFALSE);\r
- track->SetP(p,kTRUE);\r
-\r
-\r
- // Mark the track as "improved" with a trick (this is done with a trick using layer 7 (ie the 8th))\r
- UChar_t itsClusterMap = track->GetITSClusterMap();\r
- SETBIT(itsClusterMap,7);\r
- track->SetITSClusterMap(itsClusterMap);\r
- //\r
-\r
- // write out debug infos\r
- if (fDebugNtuple->GetEntries()<fNDebug) {\r
- Int_t idbg=0;\r
- fDebugVars[idbg++]=mc->GetPdgCode();\r
- fDebugVars[idbg++]=ptmc ;\r
- fDebugVars[idbg++]=d0rpo ;\r
- fDebugVars[idbg++]=d0zo ;\r
- fDebugVars[idbg++]=pt1o ;\r
- fDebugVars[idbg++]=sd0rpo;\r
- fDebugVars[idbg++]=sd0zo ;\r
- fDebugVars[idbg++]=spt1o ;\r
- fDebugVars[idbg++]=d0rpn ;\r
- fDebugVars[idbg++]=d0zn ;\r
- fDebugVars[idbg++]=pt1n ;\r
- fDebugVars[idbg++]=sd0rpn;\r
- fDebugVars[idbg++]=sd0zn ;\r
- fDebugVars[idbg++]=spt1n ;\r
- fDebugVars[idbg++]=d0rpmc;\r
- fDebugVars[idbg++]=d0zmc ;\r
- fDebugVars[idbg++]=pt1mc ;\r
- fDebugNtuple->Fill(fDebugVars);\r
- PostData(1,fDebugOutput);\r
- }\r
-}\r
-\r
-AliESDVertex* AliAnalysisTaskSEImproveITS::RecalculateVertex(const AliVVertex *old,TObjArray *tracks,Double_t bField) {\r
- //\r
- // Helper function to recalculate a vertex.\r
- //\r
-\r
- static UShort_t ids[]={1,2,3}; //TODO: unsave...\r
- AliVertexerTracks vertexer(bField);\r
- vertexer.SetVtxStart(old->GetX(),old->GetY(),old->GetZ());\r
- AliESDVertex *vertex=vertexer.VertexForSelectedTracks(tracks,ids);\r
- return vertex;\r
-}\r
-\r
-Double_t AliAnalysisTaskSEImproveITS::EvalGraph(Double_t x,const TGraph *graph,const TGraph *graphSA) const {\r
- //\r
- // Evaluates a TGraph without linear extrapolation. Instead the last\r
- // valid point of the graph is used when out of range.\r
- // The function assumes an ascending order of X.\r
- //\r
-\r
- // TODO: find a pretty solution for this:\r
- Int_t n =graph->GetN();\r
- Double_t xmin=graph->GetX()[0 ];\r
- Double_t xmax=graph->GetX()[n-1];\r
- if (x<xmin) {\r
- if(!graphSA) return graph->Eval(xmin);\r
- Double_t xminSA=graphSA->GetX()[0];\r
- if(x<xminSA) return graphSA->Eval(xminSA);\r
- return graphSA->Eval(x);\r
- }\r
- if (x>xmax) return graph->Eval(xmax);\r
- return graph->Eval(x);\r
-}\r
-\r
-ClassImp(AliAnalysisTaskSEImproveITS);\r
-\r
+/*************************************************************************
+* Copyright(c) 1998-2011, 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. *
+**************************************************************************/
+
+#include <TObjArray.h>
+#include <TClonesArray.h>
+#include <TGraph.h>
+#include <TFile.h>
+#include <TList.h>
+#include <TNtuple.h>
+
+#include "AliVertex.h"
+#include "AliVVertex.h"
+#include "AliESDVertex.h"
+#include "AliVertexerTracks.h"
+#include "AliAODEvent.h"
+#include "AliAODTrack.h"
+#include "AliAODMCParticle.h"
+#include "AliExternalTrackParam.h"
+#include "AliAODRecoDecayHF2Prong.h"
+#include "AliAODRecoDecayHF3Prong.h"
+#include "AliAODRecoCascadeHF.h"
+#include "AliNeutralTrackParam.h"
+#include "AliAnalysisTaskSEImproveITS.h"
+
+//
+// Implementation of the "hybrid-approach" for ITS upgrade studies.
+// The tastk smears the track parameters according to estimations
+// from single-track upgrade studies. Afterwards it recalculates
+// the parameters of the reconstructed decays.
+//
+// WARNING: This will affect all tasks in a train after this one
+// (which is typically desired, though).
+//
+
+AliAnalysisTaskSEImproveITS::AliAnalysisTaskSEImproveITS()
+ :AliAnalysisTaskSE(),
+ fD0ZResPCur (0),
+ fD0ZResKCur (0),
+ fD0ZResPiCur (0),
+ fD0RPResPCur (0),
+ fD0RPResKCur (0),
+ fD0RPResPiCur(0),
+ fPt1ResPCur (0),
+ fPt1ResKCur (0),
+ fPt1ResPiCur (0),
+ fD0ZResPUpg (0),
+ fD0ZResKUpg (0),
+ fD0ZResPiUpg (0),
+ fD0RPResPUpg (0),
+ fD0RPResKUpg (0),
+ fD0RPResPiUpg(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(kFALSE),
+ fImproveTracks(kTRUE),
+ fDebugOutput (0),
+ fDebugNtuple (0),
+ fDebugVars (0),
+ fNDebug (0)
+{
+ //
+ // Default constructor.
+ //
+}
+
+AliAnalysisTaskSEImproveITS::AliAnalysisTaskSEImproveITS(const char *name,
+ const char *resfileCurURI,
+ const char *resfileUpgURI,
+ Bool_t isRunInVertexing,
+ Int_t ndebug)
+ :AliAnalysisTaskSE(name),
+ fD0ZResPCur (0),
+ fD0ZResKCur (0),
+ fD0ZResPiCur (0),
+ fD0RPResPCur (0),
+ fD0RPResKCur (0),
+ fD0RPResPiCur(0),
+ fPt1ResPCur (0),
+ fPt1ResKCur (0),
+ fPt1ResPiCur (0),
+ fD0ZResPUpg (0),
+ fD0ZResKUpg (0),
+ fD0ZResPiUpg (0),
+ fD0RPResPUpg (0),
+ fD0RPResKUpg (0),
+ fD0RPResPiUpg(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),
+ fNDebug (ndebug)
+{
+ //
+ // Constructor to be used to create the task.
+ // The the URIs specify the resolution files to be used.
+ // They are expected to contain TGraphs with the resolutions
+ // for the current and the upgraded ITS (see code for details).
+ // One may also specify for how many tracks debug information
+ // is written to the output.
+ //
+ TFile *resfileCur=TFile::Open(resfileCurURI);
+ fD0RPResPCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResP" )));
+ fD0RPResKCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResK" )));
+ fD0RPResPiCur=new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResPi")));
+ fD0ZResPCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResP" )));
+ fD0ZResKCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResK" )));
+ fD0ZResPiCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResPi" )));
+ fPt1ResPCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResP" )));
+ fPt1ResKCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResK" )));
+ fPt1ResPiCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResPi" )));
+ fD0RPResPCur ->SetName("D0RPResPCur" );
+ fD0RPResKCur ->SetName("D0RPResKCur" );
+ fD0RPResPiCur->SetName("D0RPResPiCur");
+ fD0ZResPCur ->SetName("D0ZResPCur" );
+ fD0ZResKCur ->SetName("D0ZResKCur" );
+ fD0ZResPiCur ->SetName("D0ZResPiCur" );
+ 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" )));
+ fD0RPResKUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResK" )));
+ fD0RPResPiUpg=new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResPi")));
+ fD0ZResPUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResP" )));
+ fD0ZResKUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResK" )));
+ fD0ZResPiUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResPi" )));
+ fPt1ResPUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResP" )));
+ fPt1ResKUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResK" )));
+ fPt1ResPiUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResPi" )));
+ fD0RPResPUpg ->SetName("D0RPResPUpg" );
+ fD0RPResKUpg ->SetName("D0RPResKUpg" );
+ fD0RPResPiUpg->SetName("D0RPResPiUpg");
+ fD0ZResPUpg ->SetName("D0ZResPUpg" );
+ fD0ZResKUpg ->SetName("D0ZResKUpg" );
+ fD0ZResPiUpg ->SetName("D0ZResPiUpg" );
+ 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());
+}
+
+AliAnalysisTaskSEImproveITS::~AliAnalysisTaskSEImproveITS() {
+ //
+ // Destructor.
+ //
+ delete fDebugOutput;
+}
+
+void AliAnalysisTaskSEImproveITS::UserCreateOutputObjects() {
+ //
+ // Creation of user output objects.
+ //
+ fDebugOutput=new TList();
+ fDebugOutput->SetOwner();
+ fDebugNtuple=new TNtuple("fDebugNtuple","Smearing","pdg:ptmc:d0rpo:d0zo:pt1o:sd0rpo:sd0zo:spt1o:d0rpn:d0zn:pt1n:sd0rpn:sd0zn:spt1n:d0rpmc:d0zmc:pt1mc");
+ fDebugVars=new Float_t[fDebugNtuple->GetNvar()];
+
+ fDebugOutput->Add(fDebugNtuple );
+
+ fDebugOutput->Add(fD0RPResPCur );
+ fDebugOutput->Add(fD0RPResKCur );
+ fDebugOutput->Add(fD0RPResPiCur);
+ fDebugOutput->Add(fD0ZResPCur );
+ fDebugOutput->Add(fD0ZResKCur );
+ fDebugOutput->Add(fD0ZResPiCur );
+ fDebugOutput->Add(fPt1ResPCur );
+ fDebugOutput->Add(fPt1ResKCur );
+ fDebugOutput->Add(fPt1ResPiCur );
+ fDebugOutput->Add(fD0RPResPUpg );
+ fDebugOutput->Add(fD0RPResKUpg );
+ fDebugOutput->Add(fD0RPResPiUpg);
+ fDebugOutput->Add(fD0ZResPUpg );
+ fDebugOutput->Add(fD0ZResKUpg );
+ fDebugOutput->Add(fD0ZResPiUpg );
+ fDebugOutput->Add(fPt1ResPUpg );
+ fDebugOutput->Add(fPt1ResKUpg );
+ fDebugOutput->Add(fPt1ResPiUpg );
+
+ PostData(1,fDebugOutput);
+}
+
+void AliAnalysisTaskSEImproveITS::UserExec(Option_t*) {
+ //
+ // The event loop
+ //
+ AliAODEvent *ev=0x0;
+ if(!fRunInVertexing) {
+ ev=dynamic_cast<AliAODEvent*>(InputEvent());
+ } else {
+ if(AODEvent() && IsStandardAOD()) ev = dynamic_cast<AliAODEvent*> (AODEvent());
+ }
+ if(!ev) return;
+ Double_t bz=ev->GetMagneticField();
+
+
+
+
+ // Smear all tracks
+ TClonesArray *mcs=static_cast<TClonesArray*>(ev->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
+ if (!mcs) return;
+ if (fImproveTracks) {
+ for(Int_t itrack=0;itrack<ev->GetNumberOfTracks();++itrack) {
+ SmearTrack(ev->GetTrack(itrack),mcs);
+ }
+ }
+
+ // TODO: recalculated primary vertex
+ AliVVertex *primaryVertex=ev->GetPrimaryVertex();
+
+ // Recalculate all candidates
+ // D0->Kpi
+ TClonesArray *array2Prong=static_cast<TClonesArray*>(ev->GetList()->FindObject("D0toKpi"));
+ if (array2Prong) {
+ for (Int_t icand=0;icand<array2Prong->GetEntries();++icand) {
+ AliAODRecoDecayHF2Prong *decay=static_cast<AliAODRecoDecayHF2Prong*>(array2Prong->At(icand));
+
+ // recalculate vertices
+ AliVVertex *oldSecondaryVertex=decay->GetSecondaryVtx();
+
+
+ AliExternalTrackParam et1; et1.CopyFromVTrack(static_cast<AliAODTrack*>(decay->GetDaughter(0)));
+ AliExternalTrackParam et2; et2.CopyFromVTrack(static_cast<AliAODTrack*>(decay->GetDaughter(1)));
+
+ TObjArray ta12;
+
+ ta12.Add(&et1); ta12.Add(&et2);
+ AliESDVertex *v12 =RecalculateVertex(oldSecondaryVertex,&ta12 ,bz);
+
+
+ // update secondary vertex
+ Double_t pos[3];
+ v12->GetXYZ(pos);
+
+ decay->GetSecondaryVtx()->SetPosition(pos[0],pos[1],pos[2]);
+ decay->GetSecondaryVtx()->SetChi2perNDF(v12->GetChi2toNDF());
+
+ //!!!!TODO: covariance matrix
+
+ // update d0
+ Double_t d0z0[2],covd0z0[3];
+ Double_t d0[2],d0err[2];
+ et1.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
+ d0[0]=d0z0[0];
+ d0err[0] = TMath::Sqrt(covd0z0[0]);
+ et2.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
+ d0[1]=d0z0[0];
+ d0err[1] = TMath::Sqrt(covd0z0[0]);
+ decay->Setd0Prongs(2,d0);
+ decay->Setd0errProngs(2,d0err);
+ //
+
+
+ Double_t xdummy=0.,ydummy=0.;
+ Double_t dca;
+ dca=et1.GetDCA(&et2,bz,xdummy,ydummy);
+ decay->SetDCA(dca);
+
+
+ delete v12;
+
+ Double_t px[2],py[2],pz[2];
+ for (Int_t i=0;i<2;++i) {
+ const AliAODTrack *t=static_cast<AliAODTrack*>(decay->GetDaughter(i));
+ px[i]=t->Px();
+ py[i]=t->Py();
+ pz[i]=t->Pz();
+ }
+ decay->SetPxPyPzProngs(2,px,py,pz);
+ }
+ }
+
+
+ // 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 (array3Prong) {
+ for (Int_t icand=0;icand<array3Prong->GetEntries();++icand) {
+ AliAODRecoDecayHF3Prong *decay=static_cast<AliAODRecoDecayHF3Prong*>(array3Prong->At(icand));
+
+ // recalculate vertices
+ AliVVertex *oldSecondaryVertex=decay->GetSecondaryVtx();
+ AliExternalTrackParam et1; et1.CopyFromVTrack(static_cast<AliAODTrack*>(decay->GetDaughter(0)));
+ AliExternalTrackParam et2; et2.CopyFromVTrack(static_cast<AliAODTrack*>(decay->GetDaughter(1)));
+ AliExternalTrackParam et3; et3.CopyFromVTrack(static_cast<AliAODTrack*>(decay->GetDaughter(2)));
+ TObjArray ta123,ta12,ta23;
+ ta123.Add(&et1);ta123.Add(&et2);ta123.Add(&et3);
+ ta12. Add(&et1);ta12 .Add(&et2);
+ ta23 .Add(&et2);ta23 .Add(&et3);
+ AliESDVertex *v123=RecalculateVertex(oldSecondaryVertex,&ta123,bz);
+ AliESDVertex *v12 =RecalculateVertex(oldSecondaryVertex,&ta12 ,bz);
+ AliESDVertex *v23 =RecalculateVertex(oldSecondaryVertex,&ta23 ,bz);
+
+ // update secondary vertex
+ Double_t pos[3];
+ v123->GetXYZ(pos);
+ decay->GetSecondaryVtx()->SetPosition(pos[0],pos[1],pos[2]);
+ decay->GetSecondaryVtx()->SetChi2perNDF(v123->GetChi2toNDF());
+ //TODO: covariance matrix
+
+ // update d0 for all progs
+ Double_t d0z0[2],covd0z0[3];
+ Double_t d0[3],d0err[3];
+ et1.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
+ d0[0]=d0z0[0];
+ d0err[0] = TMath::Sqrt(covd0z0[0]);
+ et2.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
+ d0[1]=d0z0[0];
+ d0err[1] = TMath::Sqrt(covd0z0[0]);
+ et3.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
+ d0[2]=d0z0[0];
+ d0err[2] = TMath::Sqrt(covd0z0[0]);
+ decay->Setd0Prongs (3,d0 );
+ decay->Setd0errProngs(3,d0err);
+ // TODO: setter missing
+
+ // update dca for prong combinations
+ Double_t xdummy=0.,ydummy=0.;
+ Double_t dca[3];
+ dca[0]=et1.GetDCA(&et2,bz,xdummy,ydummy);
+ dca[1]=et3.GetDCA(&et2,bz,xdummy,ydummy);
+ dca[2]=et1.GetDCA(&et3,bz,xdummy,ydummy);
+ decay->SetDCAs(3,dca);
+
+ // update dist12 and dist23
+ primaryVertex->GetXYZ(pos);
+ decay->SetDist12toPrim(TMath::Sqrt((v12->GetX()-pos[0])*(v12->GetX()-pos[0])
+ +(v12->GetY()-pos[1])*(v12->GetY()-pos[1])
+ +(v12->GetZ()-pos[2])*(v12->GetZ()-pos[2])));
+ decay->SetDist23toPrim(TMath::Sqrt((v23->GetX()-pos[0])*(v23->GetX()-pos[0])
+ +(v23->GetY()-pos[1])*(v23->GetY()-pos[1])
+ +(v23->GetZ()-pos[2])*(v23->GetZ()-pos[2])));
+
+ delete v123;delete v12;delete v23;
+
+ Double_t px[3],py[3],pz[3];
+ for (Int_t i=0;i<3;++i) {
+ const AliAODTrack *t=static_cast<AliAODTrack*>(decay->GetDaughter(i));
+ px[i]=t->Px();
+ py[i]=t->Py();
+ pz[i]=t->Pz();
+ }
+ decay->SetPxPyPzProngs(3,px,py,pz);
+ }
+ }
+}
+
+void AliAnalysisTaskSEImproveITS::SmearTrack(AliAODTrack *track,const TClonesArray *mcs) {
+ // Early exit, if this track has nothing in common with the ITS
+
+ 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());
+//TODO: Double_t *covar=const_cast<Double_t*>(et.GetCovariance());
+
+ // Get MC info
+ Int_t imc=track->GetLabel();
+ if (imc<=0) return;
+ const AliAODMCParticle *mc=static_cast<AliAODMCParticle*>(mcs->At(imc));
+ Double_t mcx[3];
+ Double_t mcp[3];
+ Double_t mccv[36]={0.};
+ Short_t mcc;
+ mc->XvYvZv(mcx);
+ mc->PxPyPz(mcp);
+ mcc=mc->Charge();
+ AliExternalTrackParam mct(mcx,mcp,mccv,mcc);
+ const Double_t *parammc=mct.GetParameter();
+//TODO: const Double_t *covermc=mct.GetCovariance();
+ AliVertex vtx(mcx,1.,1);
+
+ // Correct reference points and frames according to MC
+ // TODO: B-Field correct?
+ // TODO: failing propagation....
+ et.PropagateToDCA(&vtx,track->GetBz(),10.);
+ et.Rotate(mct.GetAlpha());
+
+ // Select appropriate smearing functions
+ Double_t ptmc=TMath::Abs(mc->Pt());
+ Double_t sd0rpn=0.;
+ Double_t sd0zn =0.;
+ Double_t spt1n =0.;
+ Double_t sd0rpo=0.;
+ Double_t sd0zo =0.;
+ Double_t spt1o =0.;
+ switch (mc->GetPdgCode()) {
+ case 2212: case -2212:
+ 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(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(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;
+ }
+
+ // Use the same units (i.e. cm and GeV/c)! TODO: pt!
+ sd0rpo*=1.e-4;
+ sd0zo *=1.e-4;
+ sd0rpn*=1.e-4;
+ sd0zn *=1.e-4;
+
+ // Apply the smearing
+ Double_t d0zo =param [1];
+ Double_t d0zmc =parammc[1];
+ Double_t d0rpo =param [0];
+ Double_t d0rpmc=parammc[0];
+ Double_t pt1o =param [4];
+ Double_t pt1mc =parammc[4];
+ Double_t dd0zo =d0zo-d0zmc;
+ Double_t dd0zn =dd0zo *(sd0zo >0. ? (sd0zn /sd0zo ) : 1.);
+ Double_t d0zn =d0zmc+dd0zn;
+ Double_t dd0rpo=d0rpo-d0rpmc;
+ Double_t dd0rpn=dd0rpo*(sd0rpo>0. ? (sd0rpn/sd0rpo) : 1.);
+ Double_t d0rpn =d0rpmc+dd0rpn;
+ Double_t dpt1o =pt1o-pt1mc;
+ Double_t dpt1n =dpt1o *(spt1o >0. ? (spt1n /spt1o ) : 1.);
+ Double_t pt1n =pt1mc+dpt1n;
+ param[0]=d0rpn;
+ param[1]=d0zn ;
+ param[4]=pt1n ;
+
+ // Copy the smeared parameters to the AOD track
+ Double_t x[3];
+ Double_t p[3];
+ et.GetXYZ(x);
+ et.GetPxPyPz(p);
+ 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;
+ fDebugVars[idbg++]=mc->GetPdgCode();
+ fDebugVars[idbg++]=ptmc ;
+ fDebugVars[idbg++]=d0rpo ;
+ fDebugVars[idbg++]=d0zo ;
+ fDebugVars[idbg++]=pt1o ;
+ fDebugVars[idbg++]=sd0rpo;
+ fDebugVars[idbg++]=sd0zo ;
+ fDebugVars[idbg++]=spt1o ;
+ fDebugVars[idbg++]=d0rpn ;
+ fDebugVars[idbg++]=d0zn ;
+ fDebugVars[idbg++]=pt1n ;
+ fDebugVars[idbg++]=sd0rpn;
+ fDebugVars[idbg++]=sd0zn ;
+ fDebugVars[idbg++]=spt1n ;
+ fDebugVars[idbg++]=d0rpmc;
+ fDebugVars[idbg++]=d0zmc ;
+ fDebugVars[idbg++]=pt1mc ;
+ fDebugNtuple->Fill(fDebugVars);
+ PostData(1,fDebugOutput);
+ }
+}
+
+AliESDVertex* AliAnalysisTaskSEImproveITS::RecalculateVertex(const AliVVertex *old,TObjArray *tracks,Double_t bField) {
+ //
+ // Helper function to recalculate a vertex.
+ //
+
+ static UShort_t ids[]={1,2,3}; //TODO: unsave...
+ AliVertexerTracks vertexer(bField);
+ vertexer.SetVtxStart(old->GetX(),old->GetY(),old->GetZ());
+ AliESDVertex *vertex=vertexer.VertexForSelectedTracks(tracks,ids);
+ return vertex;
+}
+
+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.
+ // The function assumes an ascending order of X.
+ //
+
+ // TODO: find a pretty solution for this:
+ Int_t n =graph->GetN();
+ Double_t xmin=graph->GetX()[0 ];
+ Double_t xmax=graph->GetX()[n-1];
+ 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);
+}
+
+ClassImp(AliAnalysisTaskSEImproveITS);
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff