PWGHF/hfe/AliHFEsmearDCA.cxx

   1 /**************************************************************************
   2 * Copyright(c) 1998-1999, 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 // Get improved dca info. by ITS upgrade implemented by the ALICE Heavy Flavour Electron Group
  17 //
  18 // Authors:
  19 //   MinJung Kweon <minjung@physi.uni-heidelberg.de>
  20 //
  21 #include <TBits.h>
  22 #include <TClass.h>
  23 #include <TH1F.h>
  24 #include <TH2F.h>
  25 #include <TList.h>
  26 #include <TString.h>
  27 #include <TMath.h>
  28 #include <TGraph.h>
  29 #include <TFile.h>
  30
  31 #include "AliAODTrack.h"
  32 #include "AliAODPid.h"
  33 #include "AliESDEvent.h"
  34 #include "AliESDVertex.h"
  35 #include "AliESDtrack.h"
  36 #include "AliLog.h"
  37 #include "AliMCParticle.h"
  38 #include "AliVEvent.h"
  39 #include "AliVTrack.h"
  40 #include "AliVParticle.h"
  41 #include "AliVertexerTracks.h"
  42 #include "AliVVertex.h"
  43 #include "AliExternalTrackParam.h"
  44 #include "AliMCEvent.h"
  45
  46 #include "AliHFEsmearDCA.h"
  47
  48 ClassImp(AliHFEsmearDCA)
  49
  50
  51 //______________________________________________________
  52 AliHFEsmearDCA::AliHFEsmearDCA(const Char_t */*name*/, const char *resfileCurURI, const char *resfileUpgURI, const Char_t */*title*/):
  53    fEvent(NULL),
  54    fMCEvent(NULL),
  55    fD0ZResPCur  (0),
  56    fD0ZResKCur  (0),
  57    fD0ZResPiCur (0),
  58    fD0ZResECur  (0),
  59    fD0RPResPCur (0),
  60    fD0RPResKCur (0),
  61    fD0RPResPiCur(0),
  62    fD0RPResECur (0),
  63    fPt1ResPCur  (0),
  64    fPt1ResKCur  (0),
  65    fPt1ResPiCur (0),
  66    fPt1ResECur  (0),
  67    fD0ZResPUpg  (0),
  68    fD0ZResKUpg  (0),
  69    fD0ZResPiUpg (0),
  70    fD0ZResEUpg  (0),
  71    fD0RPResPUpg (0),
  72    fD0RPResKUpg (0),
  73    fD0RPResPiUpg(0),
  74    fD0RPResEUpg (0),
  75    fPt1ResPUpg  (0),
  76    fPt1ResKUpg  (0),
  77    fPt1ResPiUpg (0),
  78    fPt1ResEUpg  (0)
  79 {
  80   //
  81   // Constructor to be used to create the task.
  82   // The the URIs specify the resolution files to be used.
  83   // They are expected to contain TGraphs with the resolutions
  84   // for the current and the upgraded ITS (see code for details).
  85   // One may also specify for how many tracks debug information
  86   // is written to the output.
  87   //
  88   TFile *resfileCur=TFile::Open(resfileCurURI);
  89   fD0RPResPCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResP" )));
  90   fD0RPResKCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResK" )));
  91   fD0RPResPiCur=new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResPi")));
  92   fD0RPResECur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResE")));
  93   fD0ZResPCur  =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResP"  )));
  94   fD0ZResKCur  =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResK"  )));
  95   fD0ZResPiCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResPi" )));
  96   fD0ZResECur  =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResE")));
  97   fPt1ResPCur  =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResP"  )));
  98   fPt1ResKCur  =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResK"  )));
  99   fPt1ResPiCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResPi" )));
 100   fPt1ResECur  =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResE" )));
 101   fD0RPResPCur ->SetName("D0RPResPCur" );
 102   fD0RPResKCur ->SetName("D0RPResKCur" );
 103   fD0RPResPiCur->SetName("D0RPResPiCur");
 104   fD0RPResECur ->SetName("D0RPResECur");
 105   fD0ZResPCur  ->SetName("D0ZResPCur"  );
 106   fD0ZResKCur  ->SetName("D0ZResKCur"  );
 107   fD0ZResPiCur ->SetName("D0ZResPiCur" );
 108   fD0ZResECur  ->SetName("D0ZResECur" );
 109   fPt1ResPCur  ->SetName("Pt1ResPCur"  );
 110   fPt1ResKCur  ->SetName("Pt1ResKCur"  );
 111   fPt1ResPiCur ->SetName("Pt1ResPiCur" );
 112   fPt1ResECur  ->SetName("Pt1ResECur" );
 113   delete resfileCur;
 114   TFile *resfileUpg=TFile::Open(resfileUpgURI );
 115   fD0RPResPUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResP" )));
 116   fD0RPResKUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResK" )));
 117   fD0RPResPiUpg=new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResPi")));
 118   fD0RPResEUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResE")));
 119   fD0ZResPUpg  =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResP"  )));
 120   fD0ZResKUpg  =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResK"  )));
 121   fD0ZResPiUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResPi" )));
 122   fD0ZResEUpg  =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResE" )));
 123   fPt1ResPUpg  =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResP"  )));
 124   fPt1ResKUpg  =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResK"  )));
 125   fPt1ResPiUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResPi" )));
 126   fPt1ResEUpg  =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResE" )));
 127   fD0RPResPUpg ->SetName("D0RPResPUpg" );
 128   fD0RPResKUpg ->SetName("D0RPResKUpg" );
 129   fD0RPResPiUpg->SetName("D0RPResPiUpg");
 130   fD0RPResEUpg ->SetName("D0RPResEUpg");
 131   fD0ZResPUpg  ->SetName("D0ZResPUpg"  );
 132   fD0ZResKUpg  ->SetName("D0ZResKUpg"  );
 133   fD0ZResPiUpg ->SetName("D0ZResPiUpg" );
 134   fD0ZResEUpg  ->SetName("D0ZResEUpg" );
 135   fPt1ResPUpg  ->SetName("Pt1ResPUpg"  );
 136   fPt1ResKUpg  ->SetName("Pt1ResKUpg"  );
 137   fPt1ResPiUpg ->SetName("Pt1ResPiUpg" );
 138   fPt1ResEUpg  ->SetName("Pt1ResEUpg" );
 139   delete resfileUpg;
 140
 141 }
 142
 143 //______________________________________________________
 144 AliHFEsmearDCA::AliHFEsmearDCA(const AliHFEsmearDCA &c):
 145    TObject(c),
 146    fEvent(c.fEvent),
 147    fMCEvent(c.fMCEvent),
 148    fD0ZResPCur  (c.fD0ZResPCur),
 149    fD0ZResKCur  (c.fD0ZResKCur),
 150    fD0ZResPiCur (c.fD0ZResPiCur),
 151    fD0ZResECur  (c.fD0ZResECur),
 152    fD0RPResPCur (c.fD0RPResPCur),
 153    fD0RPResKCur (c.fD0RPResKCur),
 154    fD0RPResPiCur(c.fD0RPResPiCur),
 155    fD0RPResECur (c.fD0RPResECur),
 156    fPt1ResPCur  (c.fPt1ResPCur),
 157    fPt1ResKCur  (c.fPt1ResKCur),
 158    fPt1ResPiCur (c.fPt1ResPiCur),
 159    fPt1ResECur  (c.fPt1ResECur),
 160    fD0ZResPUpg  (c.fD0ZResPUpg),
 161    fD0ZResKUpg  (c.fD0ZResKUpg),
 162    fD0ZResPiUpg (c.fD0ZResPiUpg),
 163    fD0ZResEUpg  (c.fD0ZResEUpg),
 164    fD0RPResPUpg (c.fD0RPResPUpg),
 165    fD0RPResKUpg (c.fD0RPResKUpg),
 166    fD0RPResPiUpg(c.fD0RPResPiUpg),
 167    fD0RPResEUpg (c.fD0RPResEUpg),
 168    fPt1ResPUpg  (c.fPt1ResPUpg),
 169    fPt1ResKUpg  (c.fPt1ResKUpg),
 170    fPt1ResPiUpg (c.fPt1ResPiUpg),
 171    fPt1ResEUpg  (c.fPt1ResEUpg)
 172 {
 173   //
 174   // Copy constructor
 175   // Performs a deep copy
 176   //
 177 }
 178
 179 //______________________________________________________
 180 AliHFEsmearDCA::~AliHFEsmearDCA(){
 181   //
 182   // Destructor
 183   //
 184 }
 185
 186 //______________________________________________________
 187 void AliHFEsmearDCA::SetRecEventInfo(const TObject *event){
 188   //
 189   // Set Virtual event an make a copy
 190   //
 191   if (!event) {
 192     AliError("Pointer to AliVEvent !");
 193     return;
 194   }
 195   TString className(event->ClassName());
 196   if (! (className.CompareTo("AliESDEvent")==0 || className.CompareTo("AliAODEvent")==0)) {
 197     AliError("argument must point to an AliESDEvent or AliAODEvent !");
 198     return ;
 199   }
 200   fEvent = (AliVEvent*) event;
 201
 202 }
 203
 204 //______________________________________________________
 205 void AliHFEsmearDCA::GetImproveITSImpactParameters(AliVTrack *track, Double_t &dcaxyn, Double_t &dcaxyo, Double_t &dcaxySign, Double_t &dcaxySigo, Double_t &dcazn, Double_t &dcazo, Double_t &dcazSign, Double_t &dcazSigo){
 206   //
 207   // Get HFE impact parameter (with recalculated primary vertex)
 208   //
 209   if(!fEvent){
 210     AliDebug(1, "No Input event available\n");
 211     return;
 212   }
 213   const Double_t kBeampiperadius=3.;
 214   const Double_t kBeampiperadiusnew=1.9;
 215   TString type = track->IsA()->GetName();
 216   Double_t dcan[2]={-999.,-999.};
 217   Double_t covn[3]={-999.,-999.,-999.};
 218   Double_t dcao[2]={-999.,-999.};
 219   Double_t covo[3]={-999.,-999.,-999.};
 220
 221   // Get reconstructed track parameters
 222
 223   Double_t xyz[3],pxpypz[3],cv[21];
 224   track->GetXYZ(xyz);
 225   pxpypz[0]=track->Px();
 226   pxpypz[1]=track->Py();
 227   pxpypz[2]=track->Pz();
 228   track->GetCovarianceXYZPxPyPz(cv);
 229   Short_t sign = (Short_t)track->Charge();
 230
 231   //AliExternalTrackParam et(track); //it doesn't work, i don't know why yet! so, I use the line below
 232   AliExternalTrackParam et(xyz,pxpypz,cv,sign);
 233   Double_t *param=const_cast<Double_t*>(et.GetParameter());
 234
 235   const AliMCParticle *mc=static_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(track->GetLabel())));
 236   Double_t mcx[3];
 237   Double_t mcp[3];
 238   Double_t mccv[36]={0.};
 239   Short_t  mcc;
 240   mc->XvYvZv(mcx);
 241   mc->PxPyPz(mcp);
 242   mcc=mc->Charge();
 243   AliExternalTrackParam mct(mcx,mcp,mccv,mcc);
 244   const Double_t *parammc=mct.GetParameter();
 245   const AliVVertex *tmpvtx = fMCEvent->GetPrimaryVertex();
 246
 247   AliVertex vtx(mcx,1.,1);
 248   AliVertex *tmpvtx2 = (AliVertex *)fMCEvent->GetPrimaryVertex();
 249
 250   mct.PropagateToDCA(tmpvtx2,track->GetBz(),10.); //mjtest
 251   // Correct reference points and frames according to MC
 252   // TODO: B-Field correct?
 253   // TODO: failing propagation....
 254   et.PropagateToDCA(tmpvtx2,track->GetBz(),10.);
 255   et.Rotate(mct.GetAlpha());
 256
 257
 258   // Select appropriate smearing functions
 259   Double_t ptmc=TMath::Abs(mc->Pt());
 260   Double_t sd0rpn=0.;
 261   Double_t sd0zn =0.;
 262   Double_t spt1n =0.;
 263   Double_t sd0rpo=0.;
 264   Double_t sd0zo =0.;
 265   Double_t spt1o =0.;
 266
 267   switch (mc->Particle()->GetPdgCode()) {
 268   case 2212: case -2212:
 269     sd0rpo=EvalGraph(fD0RPResPCur,ptmc);
 270     sd0zo =EvalGraph(fD0ZResPCur ,ptmc);
 271     spt1o =EvalGraph(fPt1ResPCur ,ptmc);
 272     sd0rpn=EvalGraph(fD0RPResPUpg,ptmc);
 273     sd0zn =EvalGraph(fD0ZResPUpg ,ptmc);
 274     spt1n =EvalGraph(fPt1ResPUpg ,ptmc);
 275     break;
 276   case 321: case -321:
 277     sd0rpo=EvalGraph(fD0RPResKCur,ptmc);
 278     sd0zo =EvalGraph(fD0ZResKCur ,ptmc);
 279     spt1o =EvalGraph(fPt1ResKCur ,ptmc);
 280     sd0rpn=EvalGraph(fD0RPResKUpg,ptmc);
 281     sd0zn =EvalGraph(fD0ZResKUpg ,ptmc);
 282     spt1n =EvalGraph(fPt1ResKUpg ,ptmc);
 283     break;
 284   case 211: case -211:
 285     sd0rpo=EvalGraph(fD0RPResPiCur,ptmc);
 286     sd0zo =EvalGraph(fD0ZResPiCur ,ptmc);
 287     spt1o =EvalGraph(fPt1ResPiCur ,ptmc);
 288     sd0rpn=EvalGraph(fD0RPResPiUpg,ptmc);
 289     sd0zn =EvalGraph(fD0ZResPiUpg ,ptmc);
 290     spt1n =EvalGraph(fPt1ResPiUpg ,ptmc);
 291     break;
 292   case 11: case -11:
 293     sd0rpo=EvalGraph(fD0RPResECur,ptmc);
 294     sd0zo =EvalGraph(fD0ZResECur ,ptmc);
 295     spt1o =EvalGraph(fPt1ResECur ,ptmc);
 296     sd0rpn=EvalGraph(fD0RPResEUpg,ptmc);
 297     sd0zn =EvalGraph(fD0ZResEUpg ,ptmc);
 298     spt1n =EvalGraph(fPt1ResEUpg ,ptmc);
 299     break;
 300   default:
 301     return;
 302   }
 303
 304   //mj special to check resolution smeared by 10 %
 305   sd0rpn = sd0rpo * 1.1;
 306   sd0zn = sd0zo * 1.1;
 307   spt1n = spt1o;
 308
 309   // Use the same units (i.e. cm and GeV/c)! TODO: pt!
 310   sd0rpo*=1.e-4;
 311   sd0zo *=1.e-4;
 312   sd0rpn*=1.e-4;
 313   sd0zn *=1.e-4;
 314
 315   // Apply the smearing
 316   Double_t d0zo  =param  [1];
 317   Double_t d0zmc =parammc[1];
 318   Double_t d0rpo =param  [0];
 319   Double_t d0rpmc=parammc[0];
 320   Double_t pt1o  =param  [4];
 321   Double_t pt1mc =parammc[4];
 322   Double_t dd0zo =d0zo-d0zmc;
 323   Double_t dd0zn =dd0zo *(sd0zo >0. ? (sd0zn /sd0zo ) : 1.);
 324   Double_t d0zn  =d0zmc+dd0zn;
 325   Double_t dd0rpo=d0rpo-d0rpmc;
 326   Double_t dd0rpn=dd0rpo*(sd0rpo>0. ? (sd0rpn/sd0rpo) : 1.);
 327   Double_t d0rpn =d0rpmc+dd0rpn;
 328   Double_t dpt1o =pt1o-pt1mc;
 329   Double_t dpt1n =dpt1o *(spt1o >0. ? (spt1n /spt1o ) : 1.);
 330   Double_t pt1n  =pt1mc+dpt1n;
 331
 332   param[0]=d0rpn;
 333   param[1]=d0zn ;
 334   param[4]=pt1n ;
 335
 336   //printf("d0rpo= %lf d0rpmc= %lf dd0rpo= %lf dd0rpn= %lf\n",d0rpo,d0rpmc,dd0rpo,dd0rpn);
 337   // Copy the smeared parameters to the AOD track
 338   Double_t x[3];
 339   Double_t p[3];
 340   et.GetXYZ(x);
 341   et.GetPxPyPz(p);
 342 //  printf("x[0]= %lf x[1]= %lf x[2]= %lf \n",x[0],x[1],x[2]);
 343 //  track->SetPosition(x,kFALSE);
 344 //  track->SetP(p,kTRUE);
 345
 346   AliESDEvent *esdevent = dynamic_cast<AliESDEvent *>(fEvent);
 347   if(!esdevent) return;
 348   //const AliVVertex *vtxESDSkip = esdevent->GetPrimaryVertex();
 349    // recalculate primary vertex for peri. and pp
 350    AliVertexerTracks vertexer(fEvent->GetMagneticField());
 351    vertexer.SetITSMode();
 352    vertexer.SetMinClusters(4);
 353    Int_t skipped[2];
 354    skipped[0] = track->GetID();
 355    vertexer.SetSkipTracks(1,skipped);
 356    vertexer.SetConstraintOn();
 357  //----diamond constraint-----------------------------
 358    Float_t diamondcovxy[3];
 359    esdevent->GetDiamondCovXY(diamondcovxy);
 360    Double_t pos[3]={esdevent->GetDiamondX(),esdevent->GetDiamondY(),0.};
 361    Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
 362    AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
 363    vertexer.SetVtxStart(diamond);
 364    delete diamond; diamond=NULL;
 365  //----------------------------------------------------
 366    const AliVVertex *vtxESDSkip = (const AliVVertex *) vertexer.FindPrimaryVertex(fEvent);
 367    //vtxESDSkip = vertexer.FindPrimaryVertex(fEvent);
 368
 369    // Getting the DCA
 370    // Propagation always done on a working copy to not disturb the track params of the original track
 371    AliESDtrack *esdtrack = NULL;
 372    if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
 373      // Case ESD track: take copy constructor
 374      AliESDtrack *tmptrack = dynamic_cast<AliESDtrack *>(track);
 375      if(tmptrack) esdtrack = new AliESDtrack(*tmptrack);
 376    } else {
 377     // Case AOD track: take different constructor
 378     esdtrack = new AliESDtrack(track);
 379    }
 380    if(esdtrack){
 381     et.PropagateToDCA(tmpvtx, fEvent->GetMagneticField(), kBeampiperadiusnew, dcan, covn);
 382     //et.PropagateToDCA(vtxESDSkip, fEvent->GetMagneticField(), kBeampiperadiusnew, dcan, covn);
 383     esdtrack->PropagateToDCA(vtxESDSkip, fEvent->GetMagneticField(), kBeampiperadius, dcao, covo);
 384     dcaxyn = dcan[0];
 385     dcaxyo = dcao[0];
 386     dcazn = dcan[1];
 387     dcazo = dcao[1];
 388     Double_t resfactor=1., resfactorz=1.;
 389     if(sd0rpo) resfactor = sd0rpn/sd0rpo;
 390     if(sd0zo) resfactorz = sd0zn/sd0zo;
 391     if(covo[0]) dcaxySign = dcan[0]/(TMath::Sqrt(covo[0])*resfactor);
 392     if(covo[0]) dcaxySigo = dcao[0]/TMath::Sqrt(covo[0]);
 393     if(covo[2]) dcazSign = dcan[1]/(TMath::Sqrt(covo[2])*resfactorz);
 394     if(covo[2]) dcazSigo = dcao[1]/TMath::Sqrt(covo[2]);
 395     //if(dcaxyo) printf("pt= %lf resolusion xy= %lf dcaxyo= %lf dcaxyn=  %lf ratio= %lf\n",pt1mc,resfactor,dcaxyo,dcaxyn,dcaxyn/dcaxyo);
 396
 397     if(!covo[0]){
 398      dcaxySigo = -49.;
 399      dcaxySign = -49.;
 400     }
 401     if(!covo[2]){
 402      dcazSigo = -49.;
 403      dcazSign = -49.;
 404     }
 405    }
 406    delete esdtrack;
 407    //delete vtxESDSkip;
 408
 409 }
 410
 411 //______________________________________________________
 412 Double_t AliHFEsmearDCA::EvalGraph(const TGraph *graph,Double_t x) const {
 413   //
 414   // Evaluates a TGraph without linear extrapolation. Instead the last
 415   // valid point of the graph is used when out of range.
 416   // The function assumes an ascending order of X.
 417   //
 418
 419   // TODO: find a pretty solution for this:
 420   Int_t    n   =graph->GetN();
 421   Double_t xmin=graph->GetX()[0  ];
 422   Double_t xmax=graph->GetX()[n-1];
 423   if (x<xmin) return graph->Eval(xmin);
 424   if (x>xmax) return graph->Eval(xmax);
 425   return graph->Eval(x);
 426 }