PWGHF/vertexingHF/AliAnalysisTaskSEImproveITS.cxx

   1 /*************************************************************************
   2 * Copyright(c) 1998-2011, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15
  16 #include <TObjArray.h>
  17 #include <TClonesArray.h>
  18 #include <TGraph.h>
  19 #include <TFile.h>
  20 #include <TList.h>
  21 #include <TNtuple.h>
  22
  23 #include "AliVertex.h"
  24 #include "AliVVertex.h"
  25 #include "AliESDVertex.h"
  26 #include "AliVertexerTracks.h"
  27 #include "AliAODEvent.h"
  28 #include "AliAODTrack.h"
  29 #include "AliAODMCParticle.h"
  30 #include "AliExternalTrackParam.h"
  31 #include "AliAODRecoDecayHF2Prong.h"
  32 #include "AliAODRecoDecayHF3Prong.h"
  33 #include "AliAODRecoCascadeHF.h"
  34 #include "AliNeutralTrackParam.h"
  35 #include "AliAnalysisTaskSEImproveITS.h"
  36
  37 //
  38 // Implementation of the "hybrid-approach" for ITS upgrade studies.
  39 // The tastk smears the track parameters according to estimations
  40 // from single-track upgrade studies. Afterwards it recalculates
  41 // the parameters of the reconstructed decays.
  42 //
  43 // WARNING: This will affect all tasks in a train after this one
  44 // (which is typically desired, though).
  45 //
  46
  47 AliAnalysisTaskSEImproveITS::AliAnalysisTaskSEImproveITS()
  48   :AliAnalysisTaskSE(),
  49    fD0ZResPCur  (0),
  50    fD0ZResKCur  (0),
  51    fD0ZResPiCur (0),
  52    fD0RPResPCur (0),
  53    fD0RPResKCur (0),
  54    fD0RPResPiCur(0),
  55    fPt1ResPCur  (0),
  56    fPt1ResKCur  (0),
  57    fPt1ResPiCur (0),
  58    fD0ZResPUpg  (0),
  59    fD0ZResKUpg  (0),
  60    fD0ZResPiUpg (0),
  61    fD0RPResPUpg (0),
  62    fD0RPResKUpg (0),
  63    fD0RPResPiUpg(0),
  64    fPt1ResPUpg  (0),
  65    fPt1ResKUpg  (0),
  66    fPt1ResPiUpg (0),
  67    fD0ZResPCurSA  (0),
  68    fD0ZResKCurSA  (0),
  69    fD0ZResPiCurSA (0),
  70    fD0RPResPCurSA (0),
  71    fD0RPResKCurSA (0),
  72    fD0RPResPiCurSA(0),
  73    fPt1ResPCurSA  (0),
  74    fPt1ResKCurSA  (0),
  75    fPt1ResPiCurSA (0),
  76    fD0ZResPUpgSA  (0),
  77    fD0ZResKUpgSA  (0),
  78    fD0ZResPiUpgSA (0),
  79    fD0RPResPUpgSA (0),
  80    fD0RPResKUpgSA (0),
  81    fD0RPResPiUpgSA(0),
  82    fPt1ResPUpgSA  (0),
  83    fPt1ResKUpgSA  (0),
  84    fPt1ResPiUpgSA (0),
  85    fRunInVertexing(kFALSE),
  86    fImproveTracks(kTRUE),
  87    fDebugOutput (0),
  88    fDebugNtuple (0),
  89    fDebugVars   (0),
  90    fNDebug      (0)
  91 {
  92   //
  93   // Default constructor.
  94   //
  95 }
  96
  97 AliAnalysisTaskSEImproveITS::AliAnalysisTaskSEImproveITS(const char *name,
  98                            const char *resfileCurURI,
  99                            const char *resfileUpgURI,
 100                            Bool_t isRunInVertexing,
 101                            Int_t ndebug)
 102   :AliAnalysisTaskSE(name),
 103    fD0ZResPCur  (0),
 104    fD0ZResKCur  (0),
 105    fD0ZResPiCur (0),
 106    fD0RPResPCur (0),
 107    fD0RPResKCur (0),
 108    fD0RPResPiCur(0),
 109    fPt1ResPCur  (0),
 110    fPt1ResKCur  (0),
 111    fPt1ResPiCur (0),
 112    fD0ZResPUpg  (0),
 113    fD0ZResKUpg  (0),
 114    fD0ZResPiUpg (0),
 115    fD0RPResPUpg (0),
 116    fD0RPResKUpg (0),
 117    fD0RPResPiUpg(0),
 118    fPt1ResPUpg  (0),
 119    fPt1ResKUpg  (0),
 120    fPt1ResPiUpg (0),
 121    fD0ZResPCurSA  (0),
 122    fD0ZResKCurSA  (0),
 123    fD0ZResPiCurSA (0),
 124    fD0RPResPCurSA (0),
 125    fD0RPResKCurSA (0),
 126    fD0RPResPiCurSA(0),
 127    fPt1ResPCurSA  (0),
 128    fPt1ResKCurSA  (0),
 129    fPt1ResPiCurSA (0),
 130    fD0ZResPUpgSA  (0),
 131    fD0ZResKUpgSA  (0),
 132    fD0ZResPiUpgSA (0),
 133    fD0RPResPUpgSA (0),
 134    fD0RPResKUpgSA (0),
 135    fD0RPResPiUpgSA(0),
 136    fPt1ResPUpgSA  (0),
 137    fPt1ResKUpgSA  (0),
 138    fPt1ResPiUpgSA (0),
 139    fRunInVertexing(isRunInVertexing),
 140    fImproveTracks(kTRUE),
 141    fDebugOutput (0),
 142    fDebugNtuple (0),
 143    fDebugVars   (0),
 144    fNDebug      (ndebug)
 145 {
 146   //
 147   // Constructor to be used to create the task.
 148   // The the URIs specify the resolution files to be used.
 149   // They are expected to contain TGraphs with the resolutions
 150   // for the current and the upgraded ITS (see code for details).
 151   // One may also specify for how many tracks debug information
 152   // is written to the output.
 153   //
 154   TFile *resfileCur=TFile::Open(resfileCurURI);
 155   fD0RPResPCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResP" )));
 156   fD0RPResKCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResK" )));
 157   fD0RPResPiCur=new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResPi")));
 158   fD0ZResPCur  =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResP"  )));
 159   fD0ZResKCur  =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResK"  )));
 160   fD0ZResPiCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResPi" )));
 161   fPt1ResPCur  =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResP"  )));
 162   fPt1ResKCur  =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResK"  )));
 163   fPt1ResPiCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResPi" )));
 164   fD0RPResPCur ->SetName("D0RPResPCur" );
 165   fD0RPResKCur ->SetName("D0RPResKCur" );
 166   fD0RPResPiCur->SetName("D0RPResPiCur");
 167   fD0ZResPCur  ->SetName("D0ZResPCur"  );
 168   fD0ZResKCur  ->SetName("D0ZResKCur"  );
 169   fD0ZResPiCur ->SetName("D0ZResPiCur" );
 170   fPt1ResPCur  ->SetName("Pt1ResPCur"  );
 171   fPt1ResKCur  ->SetName("Pt1ResKCur"  );
 172   fPt1ResPiCur ->SetName("Pt1ResPiCur" );
 173   fD0RPResPCurSA =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResPSA" )));
 174   fD0RPResKCurSA =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResKSA" )));
 175   fD0RPResPiCurSA=new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResPiSA")));
 176   fD0ZResPCurSA  =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResPSA"  )));
 177   fD0ZResKCurSA  =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResKSA"  )));
 178   fD0ZResPiCurSA =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResPiSA" )));
 179   fPt1ResPCurSA  =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResPSA"  )));
 180   fPt1ResKCurSA  =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResKSA"  )));
 181   fPt1ResPiCurSA =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResPiSA" )));
 182   fD0RPResPCurSA ->SetName("D0RPResPCurSA" );
 183   fD0RPResKCurSA ->SetName("D0RPResKCurSA" );
 184   fD0RPResPiCurSA->SetName("D0RPResPiCurSA");
 185   fD0ZResPCurSA  ->SetName("D0ZResPCurSA"  );
 186   fD0ZResKCurSA  ->SetName("D0ZResKCurSA"  );
 187   fD0ZResPiCurSA ->SetName("D0ZResPiCurSA" );
 188   fPt1ResPCurSA  ->SetName("Pt1ResPCurSA"  );
 189   fPt1ResKCurSA  ->SetName("Pt1ResKCurSA"  );
 190   fPt1ResPiCurSA ->SetName("Pt1ResPiCurSA" );
 191   delete resfileCur;
 192   TFile *resfileUpg=TFile::Open(resfileUpgURI );
 193   fD0RPResPUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResP" )));
 194   fD0RPResKUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResK" )));
 195   fD0RPResPiUpg=new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResPi")));
 196   fD0ZResPUpg  =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResP"  )));
 197   fD0ZResKUpg  =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResK"  )));
 198   fD0ZResPiUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResPi" )));
 199   fPt1ResPUpg  =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResP"  )));
 200   fPt1ResKUpg  =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResK"  )));
 201   fPt1ResPiUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResPi" )));
 202   fD0RPResPUpg ->SetName("D0RPResPUpg" );
 203   fD0RPResKUpg ->SetName("D0RPResKUpg" );
 204   fD0RPResPiUpg->SetName("D0RPResPiUpg");
 205   fD0ZResPUpg  ->SetName("D0ZResPUpg"  );
 206   fD0ZResKUpg  ->SetName("D0ZResKUpg"  );
 207   fD0ZResPiUpg ->SetName("D0ZResPiUpg" );
 208   fPt1ResPUpg  ->SetName("Pt1ResPUpg"  );
 209   fPt1ResKUpg  ->SetName("Pt1ResKUpg"  );
 210   fPt1ResPiUpg ->SetName("Pt1ResPiUpg" );
 211   fD0RPResPUpgSA =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResPSA" )));
 212   fD0RPResKUpgSA =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResKSA" )));
 213   fD0RPResPiUpgSA=new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResPiSA")));
 214   fD0ZResPUpgSA  =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResPSA"  )));
 215   fD0ZResKUpgSA  =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResKSA"  )));
 216   fD0ZResPiUpgSA =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResPiSA" )));
 217   fPt1ResPUpgSA  =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResPSA"  )));
 218   fPt1ResKUpgSA  =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResKSA"  )));
 219   fPt1ResPiUpgSA =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResPiSA" )));
 220   fD0RPResPUpgSA ->SetName("D0RPResPUpgSA" );
 221   fD0RPResKUpgSA ->SetName("D0RPResKUpgSA" );
 222   fD0RPResPiUpgSA->SetName("D0RPResPiUpgSA");
 223   fD0ZResPUpgSA  ->SetName("D0ZResPUpgSA"  );
 224   fD0ZResKUpgSA  ->SetName("D0ZResKUpgSA"  );
 225   fD0ZResPiUpgSA ->SetName("D0ZResPiUpgSA" );
 226   fPt1ResPUpgSA  ->SetName("Pt1ResPUpgSA"  );
 227   fPt1ResKUpgSA  ->SetName("Pt1ResKUpgSA"  );
 228   fPt1ResPiUpgSA ->SetName("Pt1ResPiUpgSA" );
 229   delete resfileUpg;
 230
 231   DefineOutput(1,TNtuple::Class());
 232 }
 233
 234 AliAnalysisTaskSEImproveITS::~AliAnalysisTaskSEImproveITS() {
 235   //
 236   // Destructor.
 237   //
 238   delete fDebugOutput;
 239 }
 240
 241 void AliAnalysisTaskSEImproveITS::UserCreateOutputObjects() {
 242   //
 243   // Creation of user output objects.
 244   //
 245   fDebugOutput=new TList();
 246   fDebugOutput->SetOwner();
 247   fDebugNtuple=new TNtuple("fDebugNtuple","Smearing","pdg:ptmc:d0rpo:d0zo:pt1o:sd0rpo:sd0zo:spt1o:d0rpn:d0zn:pt1n:sd0rpn:sd0zn:spt1n:d0rpmc:d0zmc:pt1mc");
 248   fDebugVars=new Float_t[fDebugNtuple->GetNvar()];
 249
 250   fDebugOutput->Add(fDebugNtuple );
 251
 252   fDebugOutput->Add(fD0RPResPCur );
 253   fDebugOutput->Add(fD0RPResKCur );
 254   fDebugOutput->Add(fD0RPResPiCur);
 255   fDebugOutput->Add(fD0ZResPCur  );
 256   fDebugOutput->Add(fD0ZResKCur  );
 257   fDebugOutput->Add(fD0ZResPiCur );
 258   fDebugOutput->Add(fPt1ResPCur  );
 259   fDebugOutput->Add(fPt1ResKCur  );
 260   fDebugOutput->Add(fPt1ResPiCur );
 261   fDebugOutput->Add(fD0RPResPUpg );
 262   fDebugOutput->Add(fD0RPResKUpg );
 263   fDebugOutput->Add(fD0RPResPiUpg);
 264   fDebugOutput->Add(fD0ZResPUpg  );
 265   fDebugOutput->Add(fD0ZResKUpg  );
 266   fDebugOutput->Add(fD0ZResPiUpg );
 267   fDebugOutput->Add(fPt1ResPUpg  );
 268   fDebugOutput->Add(fPt1ResKUpg  );
 269   fDebugOutput->Add(fPt1ResPiUpg );
 270
 271   PostData(1,fDebugOutput);
 272 }
 273
 274 void AliAnalysisTaskSEImproveITS::UserExec(Option_t*) {
 275   //
 276   // The event loop
 277   //
 278   AliAODEvent *ev=0x0;
 279   if(!fRunInVertexing) {
 280     ev=dynamic_cast<AliAODEvent*>(InputEvent());
 281   } else {
 282     if(AODEvent() && IsStandardAOD()) ev = dynamic_cast<AliAODEvent*> (AODEvent());
 283   }
 284   if(!ev) return;
 285   Double_t bz=ev->GetMagneticField();
 286
 287
 288
 289
 290   // Smear all tracks
 291   TClonesArray *mcs=static_cast<TClonesArray*>(ev->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
 292   if (!mcs) return;
 293   if (fImproveTracks) {
 294     for(Int_t itrack=0;itrack<ev->GetNumberOfTracks();++itrack) {
 295       SmearTrack(ev->GetTrack(itrack),mcs);
 296     }
 297   }
 298
 299   // TODO: recalculated primary vertex
 300   AliVVertex *primaryVertex=ev->GetPrimaryVertex();
 301
 302   // Recalculate all candidates
 303   // D0->Kpi
 304   TClonesArray *array2Prong=static_cast<TClonesArray*>(ev->GetList()->FindObject("D0toKpi"));
 305   if (array2Prong) {
 306       for (Int_t icand=0;icand<array2Prong->GetEntries();++icand) {
 307       AliAODRecoDecayHF2Prong *decay=static_cast<AliAODRecoDecayHF2Prong*>(array2Prong->At(icand));
 308
 309       // recalculate vertices
 310       AliVVertex *oldSecondaryVertex=decay->GetSecondaryVtx();
 311
 312
 313       AliExternalTrackParam et1; et1.CopyFromVTrack(static_cast<AliAODTrack*>(decay->GetDaughter(0)));
 314       AliExternalTrackParam et2; et2.CopyFromVTrack(static_cast<AliAODTrack*>(decay->GetDaughter(1)));
 315
 316       TObjArray ta12;
 317
 318       ta12.Add(&et1); ta12.Add(&et2);
 319       AliESDVertex *v12 =RecalculateVertex(oldSecondaryVertex,&ta12 ,bz);
 320
 321
 322       // update secondary vertex
 323       Double_t pos[3];
 324       v12->GetXYZ(pos);
 325
 326       decay->GetSecondaryVtx()->SetPosition(pos[0],pos[1],pos[2]);
 327       decay->GetSecondaryVtx()->SetChi2perNDF(v12->GetChi2toNDF());
 328
 329       //!!!!TODO: covariance matrix
 330
 331       // update d0
 332       Double_t d0z0[2],covd0z0[3];
 333       Double_t d0[2],d0err[2];
 334       et1.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
 335       d0[0]=d0z0[0];
 336       d0err[0] = TMath::Sqrt(covd0z0[0]);
 337       et2.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
 338       d0[1]=d0z0[0];
 339       d0err[1] = TMath::Sqrt(covd0z0[0]);
 340       decay->Setd0Prongs(2,d0);
 341       decay->Setd0errProngs(2,d0err);
 342       //
 343
 344
 345       Double_t xdummy=0.,ydummy=0.;
 346       Double_t dca;
 347       dca=et1.GetDCA(&et2,bz,xdummy,ydummy);
 348       decay->SetDCA(dca);
 349
 350
 351       delete v12;
 352
 353       Double_t px[2],py[2],pz[2];
 354       for (Int_t i=0;i<2;++i) {
 355         const AliAODTrack *t=static_cast<AliAODTrack*>(decay->GetDaughter(i));
 356         px[i]=t->Px();
 357         py[i]=t->Py();
 358         pz[i]=t->Pz();
 359       }
 360       decay->SetPxPyPzProngs(2,px,py,pz);
 361     }
 362   }
 363
 364
 365   // Dstar->Kpipi
 366   TClonesArray *arrayCascade=static_cast<TClonesArray*>(ev->GetList()->FindObject("Dstar"));
 367
 368   if (arrayCascade) {
 369     for (Int_t icand=0;icand<arrayCascade->GetEntries();++icand) {
 370       AliAODRecoCascadeHF *decayDstar=static_cast<AliAODRecoCascadeHF*>(arrayCascade->At(icand));
 371       //Get D0 from D*
 372       AliAODRecoDecayHF2Prong* decay=(AliAODRecoDecayHF2Prong*)decayDstar->Get2Prong();
 373
 374       // recalculate vertices
 375       //AliVVertex *oldSecondaryVertex=decay->GetSecondaryVtx();
 376
 377       //soft pion
 378       AliExternalTrackParam et3; et3.CopyFromVTrack(static_cast<AliAODTrack*>(decayDstar->GetBachelor()));
 379
 380       //track D0
 381       AliNeutralTrackParam *trackD0 = new AliNeutralTrackParam(decay);
 382
 383       //!!!!TODO: covariance matrix
 384
 385       // update d0
 386       Double_t d0z0[2],covd0z0[3];
 387       Double_t d01[2],d01err[2];
 388
 389       //the D*
 390       et3.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
 391       d01[0]=d0z0[0];
 392       d01err[0] = TMath::Sqrt(covd0z0[0]);
 393       trackD0->PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
 394       d01[1]=d0z0[0];
 395       d01err[1] = TMath::Sqrt(covd0z0[0]);
 396       decayDstar->Setd0Prongs(2,d01);
 397       decayDstar->Setd0errProngs(2,d01err);
 398
 399       // delete v12;
 400       delete trackD0; trackD0=NULL;
 401
 402        // a run for D*
 403       Double_t px1[2],py1[2],pz1[2];
 404       for (Int_t i=0;i<2;++i) {
 405         const AliAODTrack *t1=static_cast<AliAODTrack*>(decayDstar->GetDaughter(i));
 406         px1[i]=t1->Px();
 407         py1[i]=t1->Py();
 408         pz1[i]=t1->Pz();
 409       }
 410       decayDstar->SetPxPyPzProngs(2,px1,py1,pz1);
 411
 412      }
 413   }
 414
 415
 416   // Three prong
 417   TClonesArray *array3Prong=static_cast<TClonesArray*>(ev->GetList()->FindObject("Charm3Prong"));
 418   if (array3Prong) {
 419     for (Int_t icand=0;icand<array3Prong->GetEntries();++icand) {
 420       AliAODRecoDecayHF3Prong *decay=static_cast<AliAODRecoDecayHF3Prong*>(array3Prong->At(icand));
 421
 422       // recalculate vertices
 423       AliVVertex *oldSecondaryVertex=decay->GetSecondaryVtx();
 424       AliExternalTrackParam et1; et1.CopyFromVTrack(static_cast<AliAODTrack*>(decay->GetDaughter(0)));
 425       AliExternalTrackParam et2; et2.CopyFromVTrack(static_cast<AliAODTrack*>(decay->GetDaughter(1)));
 426       AliExternalTrackParam et3; et3.CopyFromVTrack(static_cast<AliAODTrack*>(decay->GetDaughter(2)));
 427       TObjArray ta123,ta12,ta23;
 428       ta123.Add(&et1);ta123.Add(&et2);ta123.Add(&et3);
 429       ta12. Add(&et1);ta12 .Add(&et2);
 430                       ta23 .Add(&et2);ta23 .Add(&et3);
 431       AliESDVertex *v123=RecalculateVertex(oldSecondaryVertex,&ta123,bz);
 432       AliESDVertex *v12 =RecalculateVertex(oldSecondaryVertex,&ta12 ,bz);
 433       AliESDVertex *v23 =RecalculateVertex(oldSecondaryVertex,&ta23 ,bz);
 434
 435       // update secondary vertex
 436       Double_t pos[3];
 437       v123->GetXYZ(pos);
 438       decay->GetSecondaryVtx()->SetPosition(pos[0],pos[1],pos[2]);
 439       decay->GetSecondaryVtx()->SetChi2perNDF(v123->GetChi2toNDF());
 440       //TODO: covariance matrix
 441
 442       // update d0 for all progs
 443       Double_t d0z0[2],covd0z0[3];
 444       Double_t d0[3],d0err[3];
 445       et1.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
 446       d0[0]=d0z0[0];
 447       d0err[0] = TMath::Sqrt(covd0z0[0]);
 448       et2.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
 449       d0[1]=d0z0[0];
 450       d0err[1] = TMath::Sqrt(covd0z0[0]);
 451       et3.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
 452       d0[2]=d0z0[0];
 453       d0err[2] = TMath::Sqrt(covd0z0[0]);
 454       decay->Setd0Prongs   (3,d0   );
 455       decay->Setd0errProngs(3,d0err);
 456       // TODO: setter missing
 457
 458       // update dca for prong combinations
 459       Double_t xdummy=0.,ydummy=0.;
 460       Double_t dca[3];
 461       dca[0]=et1.GetDCA(&et2,bz,xdummy,ydummy);
 462       dca[1]=et3.GetDCA(&et2,bz,xdummy,ydummy);
 463       dca[2]=et1.GetDCA(&et3,bz,xdummy,ydummy);
 464       decay->SetDCAs(3,dca);
 465
 466       // update dist12 and dist23
 467       primaryVertex->GetXYZ(pos);
 468       decay->SetDist12toPrim(TMath::Sqrt((v12->GetX()-pos[0])*(v12->GetX()-pos[0])
 469                                         +(v12->GetY()-pos[1])*(v12->GetY()-pos[1])
 470                                         +(v12->GetZ()-pos[2])*(v12->GetZ()-pos[2])));
 471       decay->SetDist23toPrim(TMath::Sqrt((v23->GetX()-pos[0])*(v23->GetX()-pos[0])
 472                                         +(v23->GetY()-pos[1])*(v23->GetY()-pos[1])
 473                                         +(v23->GetZ()-pos[2])*(v23->GetZ()-pos[2])));
 474
 475       delete v123;delete v12;delete v23;
 476
 477       Double_t px[3],py[3],pz[3];
 478       for (Int_t i=0;i<3;++i) {
 479         const AliAODTrack *t=static_cast<AliAODTrack*>(decay->GetDaughter(i));
 480         px[i]=t->Px();
 481         py[i]=t->Py();
 482         pz[i]=t->Pz();
 483       }
 484       decay->SetPxPyPzProngs(3,px,py,pz);
 485     }
 486   }
 487 }
 488
 489 void AliAnalysisTaskSEImproveITS::SmearTrack(AliAODTrack *track,const TClonesArray *mcs) {
 490   // Early exit, if this track has nothing in common with the ITS
 491
 492   if (!(track->HasPointOnITSLayer(0) || track->HasPointOnITSLayer(1)))
 493     return;
 494
 495   // Check if the track was already "improved" (this is done with a trick using layer 7 (ie the 8th))
 496   if (TESTBIT(track->GetITSClusterMap(),7)) return;
 497   //
 498
 499
 500   // Get reconstructed track parameters
 501   AliExternalTrackParam et; et.CopyFromVTrack(track);
 502   Double_t *param=const_cast<Double_t*>(et.GetParameter());
 503 //TODO:  Double_t *covar=const_cast<Double_t*>(et.GetCovariance());
 504
 505   // Get MC info
 506   Int_t imc=track->GetLabel();
 507   if (imc<=0) return;
 508   const AliAODMCParticle *mc=static_cast<AliAODMCParticle*>(mcs->At(imc));
 509   Double_t mcx[3];
 510   Double_t mcp[3];
 511   Double_t mccv[36]={0.};
 512   Short_t  mcc;
 513   mc->XvYvZv(mcx);
 514   mc->PxPyPz(mcp);
 515   mcc=mc->Charge();
 516   AliExternalTrackParam mct(mcx,mcp,mccv,mcc);
 517   const Double_t *parammc=mct.GetParameter();
 518 //TODO:  const Double_t *covermc=mct.GetCovariance();
 519   AliVertex vtx(mcx,1.,1);
 520
 521   // Correct reference points and frames according to MC
 522   // TODO: B-Field correct?
 523   // TODO: failing propagation....
 524   et.PropagateToDCA(&vtx,track->GetBz(),10.);
 525   et.Rotate(mct.GetAlpha());
 526
 527   // Select appropriate smearing functions
 528   Double_t ptmc=TMath::Abs(mc->Pt());
 529   Double_t sd0rpn=0.;
 530   Double_t sd0zn =0.;
 531   Double_t spt1n =0.;
 532   Double_t sd0rpo=0.;
 533   Double_t sd0zo =0.;
 534   Double_t spt1o =0.;
 535   switch (mc->GetPdgCode()) {
 536   case 2212: case -2212:
 537     sd0rpo=EvalGraph(ptmc,fD0RPResPCur,fD0RPResPCurSA);
 538     sd0zo =EvalGraph(ptmc,fD0ZResPCur,fD0ZResPCurSA);
 539     spt1o =EvalGraph(ptmc,fPt1ResPCur,fPt1ResPCurSA);
 540     sd0rpn=EvalGraph(ptmc,fD0RPResPUpg,fD0RPResPUpgSA);
 541     sd0zn =EvalGraph(ptmc,fD0ZResPUpg,fD0ZResPUpgSA);
 542     spt1n =EvalGraph(ptmc,fPt1ResPUpg,fPt1ResPUpgSA);
 543     break;
 544   case 321: case -321:
 545     sd0rpo=EvalGraph(ptmc,fD0RPResKCur,fD0RPResKCurSA);
 546     sd0zo =EvalGraph(ptmc,fD0ZResKCur,fD0ZResKCurSA);
 547     spt1o =EvalGraph(ptmc,fPt1ResKCur,fPt1ResKCurSA);
 548     sd0rpn=EvalGraph(ptmc,fD0RPResKUpg,fD0RPResKUpgSA);
 549     sd0zn =EvalGraph(ptmc,fD0ZResKUpg,fD0ZResKUpgSA);
 550     spt1n =EvalGraph(ptmc,fPt1ResKUpg,fPt1ResKUpgSA);
 551     break;
 552   case 211: case -211:
 553     sd0rpo=EvalGraph(ptmc,fD0RPResPiCur,fD0RPResPiCurSA);
 554     sd0zo =EvalGraph(ptmc,fD0ZResPiCur,fD0ZResPiCurSA);
 555     spt1o =EvalGraph(ptmc,fPt1ResPiCur,fPt1ResPiCurSA);
 556     sd0rpn=EvalGraph(ptmc,fD0RPResPiUpg,fD0RPResPiUpgSA);
 557     sd0zn =EvalGraph(ptmc,fD0ZResPiUpg,fD0ZResPiUpgSA);
 558     spt1n =EvalGraph(ptmc,fPt1ResPiUpg,fPt1ResPiUpgSA);
 559     break;
 560   default:
 561     return;
 562   }
 563
 564   // Use the same units (i.e. cm and GeV/c)! TODO: pt!
 565   sd0rpo*=1.e-4;
 566   sd0zo *=1.e-4;
 567   sd0rpn*=1.e-4;
 568   sd0zn *=1.e-4;
 569
 570   // Apply the smearing
 571   Double_t d0zo  =param  [1];
 572   Double_t d0zmc =parammc[1];
 573   Double_t d0rpo =param  [0];
 574   Double_t d0rpmc=parammc[0];
 575   Double_t pt1o  =param  [4];
 576   Double_t pt1mc =parammc[4];
 577   Double_t dd0zo =d0zo-d0zmc;
 578   Double_t dd0zn =dd0zo *(sd0zo >0. ? (sd0zn /sd0zo ) : 1.);
 579   Double_t d0zn  =d0zmc+dd0zn;
 580   Double_t dd0rpo=d0rpo-d0rpmc;
 581   Double_t dd0rpn=dd0rpo*(sd0rpo>0. ? (sd0rpn/sd0rpo) : 1.);
 582   Double_t d0rpn =d0rpmc+dd0rpn;
 583   Double_t dpt1o =pt1o-pt1mc;
 584   Double_t dpt1n =dpt1o *(spt1o >0. ? (spt1n /spt1o ) : 1.);
 585   Double_t pt1n  =pt1mc+dpt1n;
 586   param[0]=d0rpn;
 587   param[1]=d0zn ;
 588   param[4]=pt1n ;
 589
 590   // Copy the smeared parameters to the AOD track
 591   Double_t x[3];
 592   Double_t p[3];
 593   et.GetXYZ(x);
 594   et.GetPxPyPz(p);
 595   track->SetPosition(x,kFALSE);
 596   track->SetP(p,kTRUE);
 597
 598
 599   // Mark the track as "improved" with a trick (this is done with a trick using layer 7 (ie the 8th))
 600   UChar_t itsClusterMap = track->GetITSClusterMap();
 601   SETBIT(itsClusterMap,7);
 602   track->SetITSClusterMap(itsClusterMap);
 603   //
 604
 605   // write out debug infos
 606   if (fDebugNtuple->GetEntries()<fNDebug) {
 607     Int_t idbg=0;
 608     fDebugVars[idbg++]=mc->GetPdgCode();
 609     fDebugVars[idbg++]=ptmc  ;
 610     fDebugVars[idbg++]=d0rpo ;
 611     fDebugVars[idbg++]=d0zo  ;
 612     fDebugVars[idbg++]=pt1o  ;
 613     fDebugVars[idbg++]=sd0rpo;
 614     fDebugVars[idbg++]=sd0zo ;
 615     fDebugVars[idbg++]=spt1o ;
 616     fDebugVars[idbg++]=d0rpn ;
 617     fDebugVars[idbg++]=d0zn  ;
 618     fDebugVars[idbg++]=pt1n  ;
 619     fDebugVars[idbg++]=sd0rpn;
 620     fDebugVars[idbg++]=sd0zn ;
 621     fDebugVars[idbg++]=spt1n ;
 622     fDebugVars[idbg++]=d0rpmc;
 623     fDebugVars[idbg++]=d0zmc ;
 624     fDebugVars[idbg++]=pt1mc ;
 625     fDebugNtuple->Fill(fDebugVars);
 626     PostData(1,fDebugOutput);
 627   }
 628 }
 629
 630 AliESDVertex* AliAnalysisTaskSEImproveITS::RecalculateVertex(const AliVVertex *old,TObjArray *tracks,Double_t bField) {
 631   //
 632   // Helper function to recalculate a vertex.
 633   //
 634
 635   static UShort_t ids[]={1,2,3}; //TODO: unsave...
 636   AliVertexerTracks vertexer(bField);
 637   vertexer.SetVtxStart(old->GetX(),old->GetY(),old->GetZ());
 638   AliESDVertex *vertex=vertexer.VertexForSelectedTracks(tracks,ids);
 639   return vertex;
 640 }
 641
 642 Double_t AliAnalysisTaskSEImproveITS::EvalGraph(Double_t x,const TGraph *graph,const TGraph *graphSA) const {
 643   //
 644   // Evaluates a TGraph without linear extrapolation. Instead the last
 645   // valid point of the graph is used when out of range.
 646   // The function assumes an ascending order of X.
 647   //
 648
 649   // TODO: find a pretty solution for this:
 650   Int_t    n   =graph->GetN();
 651   Double_t xmin=graph->GetX()[0  ];
 652   Double_t xmax=graph->GetX()[n-1];
 653   if (x<xmin) {
 654     if(!graphSA) return graph->Eval(xmin);
 655     Double_t xminSA=graphSA->GetX()[0];
 656     if(x<xminSA) return graphSA->Eval(xminSA);
 657     return graphSA->Eval(x);
 658   }
 659   if (x>xmax) return graph->Eval(xmax);
 660   return graph->Eval(x);
 661 }
 662
 663 ClassImp(AliAnalysisTaskSEImproveITS);
 664