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