end-of-line normalization
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / AliAnalysisTaskSEImproveITS.cxx
index e81bdc0..5521ece 100644 (file)
-/*************************************************************************\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);
+