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 "AliAnalysisTaskSEImproveITS.h"
  34
  35 //
  36 // Implementation of the "hybrid-approach" for ITS upgrade studies.
  37 // The tastk smears the track parameters according to estimations
  38 // from single-track upgrade studies. Afterwards it recalculates
  39 // the parameters of the reconstructed decays.
  40 //
  41 // WARNING: This will affect all tasks in a train after this one
  42 // (which is typically desired, though).
  43 //
  44
  45 AliAnalysisTaskSEImproveITS::AliAnalysisTaskSEImproveITS()
  46   :AliAnalysisTaskSE(),
  47    fD0ZResPCur  (0),
  48    fD0ZResKCur  (0),
  49    fD0ZResPiCur (0),
  50    fD0RPResPCur (0),
  51    fD0RPResKCur (0),
  52    fD0RPResPiCur(0),
  53    fPt1ResPCur  (0),
  54    fPt1ResKCur  (0),
  55    fPt1ResPiCur (0),
  56    fD0ZResPUpg  (0),
  57    fD0ZResKUpg  (0),
  58    fD0ZResPiUpg (0),
  59    fD0RPResPUpg (0),
  60    fD0RPResKUpg (0),
  61    fD0RPResPiUpg(0),
  62    fPt1ResPUpg  (0),
  63    fPt1ResKUpg  (0),
  64    fPt1ResPiUpg (0),
  65    fRunInVertexing(kFALSE),
  66    fDebugOutput (0),
  67    fDebugNtuple (0),
  68    fDebugVars   (0),
  69    fNDebug      (0)
  70 {
  71   //
  72   // Default constructor.
  73   //
  74 }
  75
  76 AliAnalysisTaskSEImproveITS::AliAnalysisTaskSEImproveITS(const char *name,
  77                            const char *resfileCurURI,
  78                            const char *resfileUpgURI,
  79                            Bool_t isRunInVertexing,
  80                            Int_t ndebug)
  81   :AliAnalysisTaskSE(name),
  82    fD0ZResPCur  (0),
  83    fD0ZResKCur  (0),
  84    fD0ZResPiCur (0),
  85    fD0RPResPCur (0),
  86    fD0RPResKCur (0),
  87    fD0RPResPiCur(0),
  88    fPt1ResPCur  (0),
  89    fPt1ResKCur  (0),
  90    fPt1ResPiCur (0),
  91    fD0ZResPUpg  (0),
  92    fD0ZResKUpg  (0),
  93    fD0ZResPiUpg (0),
  94    fD0RPResPUpg (0),
  95    fD0RPResKUpg (0),
  96    fD0RPResPiUpg(0),
  97    fPt1ResPUpg  (0),
  98    fPt1ResKUpg  (0),
  99    fPt1ResPiUpg (0),
 100    fRunInVertexing(isRunInVertexing),
 101    fDebugOutput (0),
 102    fDebugNtuple (0),
 103    fDebugVars   (0),
 104    fNDebug      (ndebug)
 105 {
 106   //
 107   // Constructor to be used to create the task.
 108   // The the URIs specify the resolution files to be used.
 109   // They are expected to contain TGraphs with the resolutions
 110   // for the current and the upgraded ITS (see code for details).
 111   // One may also specify for how many tracks debug information
 112   // is written to the output.
 113   //
 114   TFile *resfileCur=TFile::Open(resfileCurURI);
 115   fD0RPResPCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResP" )));
 116   fD0RPResKCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResK" )));
 117   fD0RPResPiCur=new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResPi")));
 118   fD0ZResPCur  =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResP"  )));
 119   fD0ZResKCur  =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResK"  )));
 120   fD0ZResPiCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResPi" )));
 121   fPt1ResPCur  =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResP"  )));
 122   fPt1ResKCur  =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResK"  )));
 123   fPt1ResPiCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResPi" )));
 124   fD0RPResPCur ->SetName("D0RPResPCur" );
 125   fD0RPResKCur ->SetName("D0RPResKCur" );
 126   fD0RPResPiCur->SetName("D0RPResPiCur");
 127   fD0ZResPCur  ->SetName("D0ZResPCur"  );
 128   fD0ZResKCur  ->SetName("D0ZResKCur"  );
 129   fD0ZResPiCur ->SetName("D0ZResPiCur" );
 130   fPt1ResPCur  ->SetName("Pt1ResPCur"  );
 131   fPt1ResKCur  ->SetName("Pt1ResKCur"  );
 132   fPt1ResPiCur ->SetName("Pt1ResPiCur" );
 133   delete resfileCur;
 134   TFile *resfileUpg=TFile::Open(resfileUpgURI );
 135   fD0RPResPUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResP" )));
 136   fD0RPResKUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResK" )));
 137   fD0RPResPiUpg=new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResPi")));
 138   fD0ZResPUpg  =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResP"  )));
 139   fD0ZResKUpg  =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResK"  )));
 140   fD0ZResPiUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResPi" )));
 141   fPt1ResPUpg  =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResP"  )));
 142   fPt1ResKUpg  =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResK"  )));
 143   fPt1ResPiUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResPi" )));
 144   fD0RPResPUpg ->SetName("D0RPResPUpg" );
 145   fD0RPResKUpg ->SetName("D0RPResKUpg" );
 146   fD0RPResPiUpg->SetName("D0RPResPiUpg");
 147   fD0ZResPUpg  ->SetName("D0ZResPUpg"  );
 148   fD0ZResKUpg  ->SetName("D0ZResKUpg"  );
 149   fD0ZResPiUpg ->SetName("D0ZResPiUpg" );
 150   fPt1ResPUpg  ->SetName("Pt1ResPUpg"  );
 151   fPt1ResKUpg  ->SetName("Pt1ResKUpg"  );
 152   fPt1ResPiUpg ->SetName("Pt1ResPiUpg" );
 153   delete resfileUpg;
 154
 155   DefineOutput(1,TNtuple::Class());
 156 }
 157
 158 AliAnalysisTaskSEImproveITS::~AliAnalysisTaskSEImproveITS() {
 159   //
 160   // Destructor.
 161   //
 162   delete fDebugOutput;
 163 }
 164
 165 void AliAnalysisTaskSEImproveITS::UserCreateOutputObjects() {
 166   //
 167   // Creation of user output objects.
 168   //
 169   fDebugOutput=new TList();
 170   fDebugOutput->SetOwner();
 171   fDebugNtuple=new TNtuple("fDebugNtuple","Smearing","pdg:ptmc:d0rpo:d0zo:pt1o:sd0rpo:sd0zo:spt1o:d0rpn:d0zn:pt1n:sd0rpn:sd0zn:spt1n:d0rpmc:d0zmc:pt1mc");
 172   fDebugVars=new Float_t[fDebugNtuple->GetNvar()];
 173
 174   fDebugOutput->Add(fDebugNtuple );
 175
 176   fDebugOutput->Add(fD0RPResPCur );
 177   fDebugOutput->Add(fD0RPResKCur );
 178   fDebugOutput->Add(fD0RPResPiCur);
 179   fDebugOutput->Add(fD0ZResPCur  );
 180   fDebugOutput->Add(fD0ZResKCur  );
 181   fDebugOutput->Add(fD0ZResPiCur );
 182   fDebugOutput->Add(fPt1ResPCur  );
 183   fDebugOutput->Add(fPt1ResKCur  );
 184   fDebugOutput->Add(fPt1ResPiCur );
 185   fDebugOutput->Add(fD0RPResPUpg );
 186   fDebugOutput->Add(fD0RPResKUpg );
 187   fDebugOutput->Add(fD0RPResPiUpg);
 188   fDebugOutput->Add(fD0ZResPUpg  );
 189   fDebugOutput->Add(fD0ZResKUpg  );
 190   fDebugOutput->Add(fD0ZResPiUpg );
 191   fDebugOutput->Add(fPt1ResPUpg  );
 192   fDebugOutput->Add(fPt1ResKUpg  );
 193   fDebugOutput->Add(fPt1ResPiUpg );
 194
 195   PostData(1,fDebugOutput);
 196 }
 197
 198 void AliAnalysisTaskSEImproveITS::UserExec(Option_t*) {
 199   //
 200   // The event loop
 201   //
 202   AliAODEvent *ev=0x0;
 203   if(!fRunInVertexing) {
 204     ev=dynamic_cast<AliAODEvent*>(InputEvent());
 205   } else {
 206     if(AODEvent() && IsStandardAOD()) ev = dynamic_cast<AliAODEvent*> (AODEvent());
 207   }
 208   if(!ev) return;
 209   Double_t bz=ev->GetMagneticField();
 210
 211
 212
 213
 214   // Smear all tracks
 215   TClonesArray *mcs=static_cast<TClonesArray*>(ev->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
 216   if (!mcs) return;
 217   for(Int_t itrack=0;itrack<ev->GetNumberOfTracks();++itrack)
 218     SmearTrack(ev->GetTrack(itrack),mcs);
 219
 220   // TODO: recalculated primary vertex
 221   AliVVertex *primaryVertex=ev->GetPrimaryVertex();
 222
 223   // Recalculate all candidates
 224   // D0->Kpi
 225   TClonesArray *array2Prong=static_cast<TClonesArray*>(ev->GetList()->FindObject("D0toKpi"));
 226   if (array2Prong) {
 227       for (Int_t icand=0;icand<array2Prong->GetEntries();++icand) {
 228       AliAODRecoDecayHF2Prong *decay=static_cast<AliAODRecoDecayHF2Prong*>(array2Prong->At(icand));
 229
 230       // recalculate vertices
 231       AliVVertex *oldSecondaryVertex=decay->GetSecondaryVtx();
 232
 233
 234       AliExternalTrackParam et1; et1.CopyFromVTrack(static_cast<AliAODTrack*>(decay->GetDaughter(0)));
 235       AliExternalTrackParam et2; et2.CopyFromVTrack(static_cast<AliAODTrack*>(decay->GetDaughter(1)));
 236
 237       TObjArray ta12;
 238
 239       ta12.Add(&et1); ta12.Add(&et2);
 240       AliESDVertex *v12 =RecalculateVertex(oldSecondaryVertex,&ta12 ,bz);
 241
 242
 243       // update secondary vertex
 244       Double_t pos[3];
 245       v12->GetXYZ(pos);
 246
 247       decay->GetSecondaryVtx()->SetPosition(pos[0],pos[1],pos[2]);
 248       decay->GetSecondaryVtx()->SetChi2perNDF(v12->GetChi2toNDF());
 249
 250       //!!!!TODO: covariance matrix
 251
 252       // update d0
 253       Double_t d0z0[2],covd0z0[3];
 254       Double_t d0[2],d0err[2];
 255       et1.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
 256       d0[0]=d0z0[0];
 257       d0err[0] = TMath::Sqrt(covd0z0[0]);
 258       et2.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
 259       d0[1]=d0z0[0];
 260       d0err[1] = TMath::Sqrt(covd0z0[0]);
 261       decay->Setd0Prongs(2,d0);
 262       decay->Setd0errProngs(2,d0err);
 263       //
 264
 265
 266       Double_t xdummy=0.,ydummy=0.;
 267       Double_t dca;
 268       dca=et1.GetDCA(&et2,bz,xdummy,ydummy);
 269       decay->SetDCA(dca);
 270
 271
 272       delete v12;
 273
 274       Double_t px[2],py[2],pz[2];
 275       for (Int_t i=0;i<2;++i) {
 276         const AliAODTrack *t=static_cast<AliAODTrack*>(decay->GetDaughter(i));
 277         px[i]=t->Px();
 278         py[i]=t->Py();
 279         pz[i]=t->Pz();
 280       }
 281       decay->SetPxPyPzProngs(2,px,py,pz);
 282     }
 283   }
 284
 285
 286
 287   // Three prong
 288   TClonesArray *array3Prong=static_cast<TClonesArray*>(ev->GetList()->FindObject("Charm3Prong"));
 289   if (array3Prong) {
 290     for (Int_t icand=0;icand<array3Prong->GetEntries();++icand) {
 291       AliAODRecoDecayHF3Prong *decay=static_cast<AliAODRecoDecayHF3Prong*>(array3Prong->At(icand));
 292
 293       // recalculate vertices
 294       AliVVertex *oldSecondaryVertex=decay->GetSecondaryVtx();
 295       AliExternalTrackParam et1; et1.CopyFromVTrack(static_cast<AliAODTrack*>(decay->GetDaughter(0)));
 296       AliExternalTrackParam et2; et2.CopyFromVTrack(static_cast<AliAODTrack*>(decay->GetDaughter(1)));
 297       AliExternalTrackParam et3; et3.CopyFromVTrack(static_cast<AliAODTrack*>(decay->GetDaughter(2)));
 298       TObjArray ta123,ta12,ta23;
 299       ta123.Add(&et1);ta123.Add(&et2);ta123.Add(&et3);
 300       ta12. Add(&et1);ta12 .Add(&et2);
 301                       ta23 .Add(&et2);ta23 .Add(&et3);
 302       AliESDVertex *v123=RecalculateVertex(oldSecondaryVertex,&ta123,bz);
 303       AliESDVertex *v12 =RecalculateVertex(oldSecondaryVertex,&ta12 ,bz);
 304       AliESDVertex *v23 =RecalculateVertex(oldSecondaryVertex,&ta23 ,bz);
 305
 306       // update secondary vertex
 307       Double_t pos[3];
 308       v123->GetXYZ(pos);
 309       decay->GetSecondaryVtx()->SetPosition(pos[0],pos[1],pos[2]);
 310       decay->GetSecondaryVtx()->SetChi2perNDF(v123->GetChi2toNDF());
 311       //TODO: covariance matrix
 312
 313       // update d0 for all progs
 314       Double_t d0z0[2],covd0z0[3];
 315       Double_t d0[3],d0err[3];
 316       et1.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
 317       d0[0]=d0z0[0];
 318       d0err[0] = TMath::Sqrt(covd0z0[0]);
 319       et2.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
 320       d0[1]=d0z0[0];
 321       d0err[1] = TMath::Sqrt(covd0z0[0]);
 322       et3.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
 323       d0[2]=d0z0[0];
 324       d0err[2] = TMath::Sqrt(covd0z0[0]);
 325       decay->Setd0Prongs   (3,d0   );
 326       decay->Setd0errProngs(3,d0err);
 327       // TODO: setter missing
 328
 329       // update dca for prong combinations
 330       Double_t xdummy=0.,ydummy=0.;
 331       Double_t dca[3];
 332       dca[0]=et1.GetDCA(&et2,bz,xdummy,ydummy);
 333       dca[1]=et3.GetDCA(&et2,bz,xdummy,ydummy);
 334       dca[2]=et1.GetDCA(&et3,bz,xdummy,ydummy);
 335       decay->SetDCAs(3,dca);
 336
 337       // update dist12 and dist23
 338       primaryVertex->GetXYZ(pos);
 339       decay->SetDist12toPrim(TMath::Sqrt((v12->GetX()-pos[0])*(v12->GetX()-pos[0])
 340                                         +(v12->GetY()-pos[1])*(v12->GetY()-pos[1])
 341                                         +(v12->GetZ()-pos[2])*(v12->GetZ()-pos[2])));
 342       decay->SetDist23toPrim(TMath::Sqrt((v23->GetX()-pos[0])*(v23->GetX()-pos[0])
 343                                         +(v23->GetY()-pos[1])*(v23->GetY()-pos[1])
 344                                         +(v23->GetZ()-pos[2])*(v23->GetZ()-pos[2])));
 345
 346       delete v123;delete v12;delete v23;
 347
 348       Double_t px[3],py[3],pz[3];
 349       for (Int_t i=0;i<3;++i) {
 350         const AliAODTrack *t=static_cast<AliAODTrack*>(decay->GetDaughter(i));
 351         px[i]=t->Px();
 352         py[i]=t->Py();
 353         pz[i]=t->Pz();
 354       }
 355       decay->SetPxPyPzProngs(3,px,py,pz);
 356     }
 357   }
 358 }
 359
 360 void AliAnalysisTaskSEImproveITS::SmearTrack(AliAODTrack *track,const TClonesArray *mcs) {
 361   // Early exit, if this track has nothing in common with the ITS
 362
 363   if (!(track->HasPointOnITSLayer(0) || track->HasPointOnITSLayer(1)))
 364     return;
 365
 366   // Get reconstructed track parameters
 367   AliExternalTrackParam et; et.CopyFromVTrack(track);
 368   Double_t *param=const_cast<Double_t*>(et.GetParameter());
 369 //TODO:  Double_t *covar=const_cast<Double_t*>(et.GetCovariance());
 370
 371   // Get MC info
 372   Int_t imc=track->GetLabel();
 373   if (imc<=0) return;
 374   const AliAODMCParticle *mc=static_cast<AliAODMCParticle*>(mcs->At(imc));
 375   Double_t mcx[3];
 376   Double_t mcp[3];
 377   Double_t mccv[36]={0.};
 378   Short_t  mcc;
 379   mc->XvYvZv(mcx);
 380   mc->PxPyPz(mcp);
 381   mcc=mc->Charge();
 382   AliExternalTrackParam mct(mcx,mcp,mccv,mcc);
 383   const Double_t *parammc=mct.GetParameter();
 384 //TODO:  const Double_t *covermc=mct.GetCovariance();
 385   AliVertex vtx(mcx,1.,1);
 386
 387   // Correct reference points and frames according to MC
 388   // TODO: B-Field correct?
 389   // TODO: failing propagation....
 390   et.PropagateToDCA(&vtx,track->GetBz(),10.);
 391   et.Rotate(mct.GetAlpha());
 392
 393   // Select appropriate smearing functions
 394   Double_t ptmc=TMath::Abs(mc->Pt());
 395   Double_t sd0rpn=0.;
 396   Double_t sd0zn =0.;
 397   Double_t spt1n =0.;
 398   Double_t sd0rpo=0.;
 399   Double_t sd0zo =0.;
 400   Double_t spt1o =0.;
 401   switch (mc->GetPdgCode()) {
 402   case 2212: case -2212:
 403     sd0rpo=EvalGraph(fD0RPResPCur,ptmc);
 404     sd0zo =EvalGraph(fD0ZResPCur ,ptmc);
 405     spt1o =EvalGraph(fPt1ResPCur ,ptmc);
 406     sd0rpn=EvalGraph(fD0RPResPUpg,ptmc);
 407     sd0zn =EvalGraph(fD0ZResPUpg ,ptmc);
 408     spt1n =EvalGraph(fPt1ResPUpg ,ptmc);
 409     break;
 410   case 321: case -321:
 411     sd0rpo=EvalGraph(fD0RPResKCur,ptmc);
 412     sd0zo =EvalGraph(fD0ZResKCur ,ptmc);
 413     spt1o =EvalGraph(fPt1ResKCur ,ptmc);
 414     sd0rpn=EvalGraph(fD0RPResKUpg,ptmc);
 415     sd0zn =EvalGraph(fD0ZResKUpg ,ptmc);
 416     spt1n =EvalGraph(fPt1ResKUpg ,ptmc);
 417     break;
 418   case 211: case -211:
 419     sd0rpo=EvalGraph(fD0RPResPiCur,ptmc);
 420     sd0zo =EvalGraph(fD0ZResPiCur ,ptmc);
 421     spt1o =EvalGraph(fPt1ResPiCur ,ptmc);
 422     sd0rpn=EvalGraph(fD0RPResPiUpg,ptmc);
 423     sd0zn =EvalGraph(fD0ZResPiUpg ,ptmc);
 424     spt1n =EvalGraph(fPt1ResPiUpg ,ptmc);
 425     break;
 426   default:
 427     return;
 428   }
 429
 430   // Use the same units (i.e. cm and GeV/c)! TODO: pt!
 431   sd0rpo*=1.e-4;
 432   sd0zo *=1.e-4;
 433   sd0rpn*=1.e-4;
 434   sd0zn *=1.e-4;
 435
 436   // Apply the smearing
 437   Double_t d0zo  =param  [1];
 438   Double_t d0zmc =parammc[1];
 439   Double_t d0rpo =param  [0];
 440   Double_t d0rpmc=parammc[0];
 441   Double_t pt1o  =param  [4];
 442   Double_t pt1mc =parammc[4];
 443   Double_t dd0zo =d0zo-d0zmc;
 444   Double_t dd0zn =dd0zo *(sd0zo >0. ? (sd0zn /sd0zo ) : 1.);
 445   Double_t d0zn  =d0zmc+dd0zn;
 446   Double_t dd0rpo=d0rpo-d0rpmc;
 447   Double_t dd0rpn=dd0rpo*(sd0rpo>0. ? (sd0rpn/sd0rpo) : 1.);
 448   Double_t d0rpn =d0rpmc+dd0rpn;
 449   Double_t dpt1o =pt1o-pt1mc;
 450   Double_t dpt1n =dpt1o *(spt1o >0. ? (spt1n /spt1o ) : 1.);
 451   Double_t pt1n  =pt1mc+dpt1n;
 452   param[0]=d0rpn;
 453   param[1]=d0zn ;
 454   param[4]=pt1n ;
 455
 456   // Copy the smeared parameters to the AOD track
 457   Double_t x[3];
 458   Double_t p[3];
 459   et.GetXYZ(x);
 460   et.GetPxPyPz(p);
 461   track->SetPosition(x,kFALSE);
 462   track->SetP(p,kTRUE);
 463
 464   // write out debug infos
 465   if (fDebugNtuple->GetEntries()<fNDebug) {
 466     Int_t idbg=0;
 467     fDebugVars[idbg++]=mc->GetPdgCode();
 468     fDebugVars[idbg++]=ptmc  ;
 469     fDebugVars[idbg++]=d0rpo ;
 470     fDebugVars[idbg++]=d0zo  ;
 471     fDebugVars[idbg++]=pt1o  ;
 472     fDebugVars[idbg++]=sd0rpo;
 473     fDebugVars[idbg++]=sd0zo ;
 474     fDebugVars[idbg++]=spt1o ;
 475     fDebugVars[idbg++]=d0rpn ;
 476     fDebugVars[idbg++]=d0zn  ;
 477     fDebugVars[idbg++]=pt1n  ;
 478     fDebugVars[idbg++]=sd0rpn;
 479     fDebugVars[idbg++]=sd0zn ;
 480     fDebugVars[idbg++]=spt1n ;
 481     fDebugVars[idbg++]=d0rpmc;
 482     fDebugVars[idbg++]=d0zmc ;
 483     fDebugVars[idbg++]=pt1mc ;
 484     fDebugNtuple->Fill(fDebugVars);
 485     PostData(1,fDebugOutput);
 486   }
 487 }
 488
 489 AliESDVertex* AliAnalysisTaskSEImproveITS::RecalculateVertex(const AliVVertex *old,TObjArray *tracks,Double_t bField) {
 490   //
 491   // Helper function to recalculate a vertex.
 492   //
 493
 494   static UShort_t ids[]={1,2,3}; //TODO: unsave...
 495   AliVertexerTracks vertexer(bField);
 496   vertexer.SetVtxStart(old->GetX(),old->GetY(),old->GetZ());
 497   AliESDVertex *vertex=vertexer.VertexForSelectedTracks(tracks,ids);
 498   return vertex;
 499 }
 500
 501 Double_t AliAnalysisTaskSEImproveITS::EvalGraph(const TGraph *graph,Double_t x) const {
 502   //
 503   // Evaluates a TGraph without linear extrapolation. Instead the last
 504   // valid point of the graph is used when out of range.
 505   // The function assumes an ascending order of X.
 506   //
 507
 508   // TODO: find a pretty solution for this:
 509   Int_t    n   =graph->GetN();
 510   Double_t xmin=graph->GetX()[0  ];
 511   Double_t xmax=graph->GetX()[n-1];
 512   if (x<xmin) return graph->Eval(xmin);
 513   if (x>xmax) return graph->Eval(xmax);
 514   return graph->Eval(x);
 515 }
 516
 517 ClassImp(AliAnalysisTaskSEImproveITS);
 518