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