PWG3/vertexingHF/AliAODRecoDecayHF.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2006, 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 /////////////////////////////////////////////////////////////
  17 //
  18 // Base class for AOD reconstructed heavy-flavour decay
  19 //
  20 // Author: A.Dainese, andrea.dainese@lnl.infn.it
  21 /////////////////////////////////////////////////////////////
  22
  23 #include <TDatabasePDG.h>
  24 #include <TVector3.h>
  25 #include <TRandom.h>
  26 #include "AliAODRecoDecay.h"
  27 #include "AliAODRecoDecayHF.h"
  28 #include "AliAODEvent.h"
  29 #include "AliVertexerTracks.h"
  30 #include "AliExternalTrackParam.h"
  31 #include "AliKFVertex.h"
  32 #include "AliVVertex.h"
  33 #include "AliESDVertex.h"
  34
  35 ClassImp(AliAODRecoDecayHF)
  36
  37 //--------------------------------------------------------------------------
  38 AliAODRecoDecayHF::AliAODRecoDecayHF() :
  39   AliAODRecoDecay(),
  40   fOwnPrimaryVtx(0x0),
  41   fEventPrimaryVtx(),
  42   fListOfCuts(),
  43   fd0err(0x0),
  44   fProngID(0x0)
  45 {
  46   //
  47   // Default Constructor
  48   //
  49 }
  50 //--------------------------------------------------------------------------
  51 AliAODRecoDecayHF::AliAODRecoDecayHF(AliAODVertex *vtx2,Int_t nprongs,Short_t charge,
  52                                      Double_t *px,Double_t *py,Double_t *pz,
  53                                      Double_t *d0,Double_t *d0err) :
  54   AliAODRecoDecay(vtx2,nprongs,charge,px,py,pz,d0),
  55   fOwnPrimaryVtx(0x0),
  56   fEventPrimaryVtx(),
  57   fListOfCuts(),
  58   fd0err(0x0),
  59   fProngID(0x0)
  60 {
  61   //
  62   // Constructor with AliAODVertex for decay vertex
  63   //
  64   fd0err = new Double_t[GetNProngs()];
  65   for(Int_t i=0; i<GetNProngs(); i++) fd0err[i] = d0err[i];
  66 }
  67 //--------------------------------------------------------------------------
  68 AliAODRecoDecayHF::AliAODRecoDecayHF(AliAODVertex *vtx2,Int_t nprongs,Short_t charge,
  69                                      Double_t *d0,Double_t *d0err) :
  70   AliAODRecoDecay(vtx2,nprongs,charge,d0),
  71   fOwnPrimaryVtx(0x0),
  72   fEventPrimaryVtx(),
  73   fListOfCuts(),
  74   fd0err(0x0),
  75   fProngID(0x0)
  76 {
  77   //
  78   // Constructor with AliAODVertex for decay vertex and without prongs momenta
  79   //
  80   fd0err = new Double_t[GetNProngs()];
  81   for(Int_t i=0; i<GetNProngs(); i++) fd0err[i] = d0err[i];
  82 }
  83 //--------------------------------------------------------------------------
  84 AliAODRecoDecayHF::AliAODRecoDecayHF(Double_t vtx1[3],Double_t vtx2[3],
  85                                      Int_t nprongs,Short_t charge,
  86                                      Double_t *px,Double_t *py,Double_t *pz,
  87                                      Double_t *d0) :
  88   AliAODRecoDecay(0x0,nprongs,charge,px,py,pz,d0),
  89   fOwnPrimaryVtx(0x0),
  90   fEventPrimaryVtx(),
  91   fListOfCuts(),
  92   fd0err(0x0),
  93   fProngID(0x0)
  94 {
  95   //
  96   // Constructor that can used for a "MC" object
  97   //
  98
  99   fOwnPrimaryVtx = new AliAODVertex(vtx1);
 100
 101   AliAODVertex *vtx = new AliAODVertex(vtx2);
 102   SetOwnSecondaryVtx(vtx);
 103
 104 }
 105 //--------------------------------------------------------------------------
 106 AliAODRecoDecayHF::AliAODRecoDecayHF(const AliAODRecoDecayHF &source) :
 107   AliAODRecoDecay(source),
 108   fOwnPrimaryVtx(0x0),
 109   fEventPrimaryVtx(source.fEventPrimaryVtx),
 110   fListOfCuts(source.fListOfCuts),
 111   fd0err(0x0),
 112   fProngID(0x0)
 113 {
 114   //
 115   // Copy constructor
 116   //
 117   if(source.GetOwnPrimaryVtx()) fOwnPrimaryVtx = new AliAODVertex(*(source.GetOwnPrimaryVtx()));
 118
 119   if(source.GetNProngs()>0) {
 120     fd0err = new Double_t[GetNProngs()];
 121     memcpy(fd0err,source.fd0err,GetNProngs()*sizeof(Double_t));
 122     if(source.fProngID) {
 123       fProngID = new UShort_t[GetNProngs()];
 124       memcpy(fProngID,source.fProngID,GetNProngs()*sizeof(UShort_t));
 125     }
 126   }
 127 }
 128 //--------------------------------------------------------------------------
 129 AliAODRecoDecayHF &AliAODRecoDecayHF::operator=(const AliAODRecoDecayHF &source)
 130 {
 131   //
 132   // assignment operator
 133   //
 134   if(&source == this) return *this;
 135
 136   AliAODRecoDecay::operator=(source);
 137
 138   fEventPrimaryVtx = source.fEventPrimaryVtx;
 139   fListOfCuts = source.fListOfCuts;
 140
 141   if(source.GetOwnPrimaryVtx()) fOwnPrimaryVtx = new AliAODVertex(*(source.GetOwnPrimaryVtx()));
 142
 143   if(source.GetNProngs()>0) {
 144     fd0err = new Double_t[GetNProngs()];
 145     memcpy(fd0err,source.fd0err,GetNProngs()*sizeof(Double_t));
 146     if(source.fProngID) {
 147       fProngID = new UShort_t[GetNProngs()];
 148       memcpy(fProngID,source.fProngID,GetNProngs()*sizeof(UShort_t));
 149     }
 150   }
 151   return *this;
 152 }
 153 //--------------------------------------------------------------------------
 154 AliAODRecoDecayHF::~AliAODRecoDecayHF() {
 155   //
 156   // Default Destructor
 157   //
 158   if(fOwnPrimaryVtx) delete fOwnPrimaryVtx;
 159   if(fd0err) delete [] fd0err;
 160   if(fProngID) delete [] fProngID;
 161 }
 162 //---------------------------------------------------------------------------
 163 AliKFParticle *AliAODRecoDecayHF::ApplyVertexingKF(Int_t *iprongs,Int_t nprongs,Int_t *pdgs,Bool_t topoCostraint, Double_t bzkG, Double_t *mass) const {
 164   //
 165   // Applies the KF vertexer
 166   // Int_t iprongs[nprongs] = indices of the prongs to be used from the vertexer
 167   // Int_t pdgs[nprongs] = pdgs assigned to the prongs, needed to define the AliKFParticle
 168   // Bool_t topoCostraint = if kTRUE, the topological constraint is applied
 169   // Double_t bzkG = magnetic field
 170   // Double_t mass[2] = {mass, sigma} for the mass constraint (if mass[0]>0 the constraint is applied).
 171   //
 172
 173   AliKFParticle::SetField(bzkG);
 174   AliKFParticle *vertexKF=0;
 175
 176   AliKFVertex copyKF;
 177   Int_t nt=0,ntcheck=0;
 178
 179   Double_t pos[3]={0.,0.,0.};
 180   if(!fOwnPrimaryVtx) {
 181     printf("AliAODRecoDecayHF::ApplyVertexingKF(): cannot apply because primary vertex is not found\n");
 182     return vertexKF;
 183   }
 184   fOwnPrimaryVtx->GetXYZ(pos);
 185   Int_t contr=fOwnPrimaryVtx->GetNContributors();
 186   Double_t covmatrix[6]={0.,0.,0.,0.,0.,0.};
 187   fOwnPrimaryVtx->GetCovarianceMatrix(covmatrix);
 188   Double_t chi2=fOwnPrimaryVtx->GetChi2();
 189   AliESDVertex primaryVtx2(pos,covmatrix,chi2,contr,"Vertex");
 190
 191
 192   if(topoCostraint){
 193    copyKF=AliKFVertex(primaryVtx2);
 194    nt=primaryVtx2.GetNContributors();
 195    ntcheck=nt;
 196   }
 197
 198   vertexKF = new AliKFParticle();
 199   for(Int_t i= 0;i<nprongs;i++){
 200     Int_t ipr=iprongs[i];
 201     AliAODTrack *aodTrack = (AliAODTrack*)GetDaughter(ipr);
 202     if(!aodTrack) {
 203       printf("AliAODRecoDecayHF::ApplyVertexingKF(): no daughters available\n");
 204       delete vertexKF; vertexKF=NULL;
 205       return vertexKF;
 206     }
 207     AliKFParticle daughterKF(*aodTrack,pdgs[i]);
 208     vertexKF->AddDaughter(daughterKF);
 209
 210     if(topoCostraint && nt>0){
 211       //Int_t index=(Int_t)GetProngID(ipr);
 212       if(!aodTrack->GetUsedForPrimVtxFit()) continue;
 213       copyKF -= daughterKF;
 214       ntcheck--;
 215     }
 216   }
 217
 218   if(topoCostraint){
 219     if(ntcheck>0) {
 220       copyKF += (*vertexKF);
 221       vertexKF->SetProductionVertex(copyKF);
 222     }
 223  }
 224
 225   if(mass[0]>0.){
 226     vertexKF->SetMassConstraint(mass[0],mass[1]);
 227   }
 228
 229   return vertexKF;
 230 }
 231 //---------------------------------------------------------------------------
 232 AliAODVertex* AliAODRecoDecayHF::RemoveDaughtersFromPrimaryVtx(AliAODEvent *aod) {
 233   //
 234   // This method returns a primary vertex without the daughter tracks of the
 235   // candidate and it recalculates the impact parameters and errors.
 236   //
 237   // The output vertex is created with "new". The user has to
 238   // set it to the candidate with SetOwnPrimaryVtx(), unset it at the end
 239   // of processing with UnsetOwnPrimaryVtx() and delete it.
 240   // If a NULL pointer is returned, the removal failed (too few tracks left).
 241   //
 242   // For the moment, the primary vertex is recalculated from scratch without
 243   // the daughter tracks.
 244   //
 245
 246   AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
 247   if(!vtxAOD) return 0;
 248   TString title=vtxAOD->GetTitle();
 249   if(!title.Contains("VertexerTracks")) return 0;
 250
 251
 252
 253   AliVertexerTracks *vertexer = new AliVertexerTracks(aod->GetMagneticField());
 254
 255   Int_t ndg = GetNDaughters();
 256
 257   vertexer->SetITSMode();
 258   vertexer->SetMinClusters(4);
 259   vertexer->SetConstraintOff();
 260
 261   if(title.Contains("WithConstraint")) {
 262     Float_t diamondcovxy[3];
 263     aod->GetDiamondCovXY(diamondcovxy);
 264     Double_t pos[3]={aod->GetDiamondX(),aod->GetDiamondY(),0.};
 265     Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
 266     AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
 267     vertexer->SetVtxStart(diamond);
 268     delete diamond; diamond=NULL;
 269   }
 270
 271   Int_t skipped[10];
 272   Int_t nTrksToSkip=0,id;
 273   AliAODTrack *t = 0;
 274   for(Int_t i=0; i<ndg; i++) {
 275     t = (AliAODTrack*)GetDaughter(i);
 276     id = (Int_t)t->GetID();
 277     if(id<0) continue;
 278     skipped[nTrksToSkip++] = id;
 279   }
 280   vertexer->SetSkipTracks(nTrksToSkip,skipped);
 281   AliESDVertex *vtxESDNew = vertexer->FindPrimaryVertex(aod);
 282
 283   delete vertexer; vertexer=NULL;
 284
 285   if(!vtxESDNew) return 0;
 286   if(vtxESDNew->GetNContributors()<=0) {
 287     delete vtxESDNew; vtxESDNew=NULL;
 288     return 0;
 289   }
 290
 291   // convert to AliAODVertex
 292   Double_t pos[3],cov[6],chi2perNDF;
 293   vtxESDNew->GetXYZ(pos); // position
 294   vtxESDNew->GetCovMatrix(cov); //covariance matrix
 295   chi2perNDF = vtxESDNew->GetChi2toNDF();
 296   delete vtxESDNew; vtxESDNew=NULL;
 297
 298   AliAODVertex *vtxAODNew = new AliAODVertex(pos,cov,chi2perNDF);
 299
 300   // now recalculate the daughters impact parameters
 301   AliExternalTrackParam *etp = 0;
 302   Double_t dz[2],covdz[3];
 303   for(Int_t i=0; i<ndg; i++) {
 304     t = (AliAODTrack*)GetDaughter(i);
 305     etp = new AliExternalTrackParam(t);
 306     if(etp->PropagateToDCA(vtxAODNew,aod->GetMagneticField(),3.,dz,covdz)) {
 307       fd0[i]    = dz[0];
 308       fd0err[i] = TMath::Sqrt(covdz[0]);
 309     }
 310     delete etp; etp=NULL;
 311   }
 312
 313   return vtxAODNew;
 314 }
 315
 316
 317
 318
 319 //-----------------------------------------------------------------------------------
 320 void AliAODRecoDecayHF::Misalign(TString misal) {
 321   //
 322   // Method to smear the impact parameter of the duaghter tracks
 323   // and the sec. vtx position accordinlgy
 324   // Useful to study effect of misalignment.
 325   // The starting point are parameterisations of the impact parameter resolution
 326   // from MC and data
 327   // Errors on d0 and vtx are not recalculated (to be done)
 328   //
 329   if(misal=="null")return;
 330   Double_t pard0rphiMC[3]={36.7,36.,1.25};// d0(pt)=[0]+[1]/(pt^[2]); in micron, conversion to cm is done later
 331   Double_t pard0rphimisal[3]={0,0,0};
 332   Double_t pard0zMC[3]={85.,130.,0.7};// d0(pt)=[0]+[1]/(pt^[2]); in micron, conversion to cm is done later
 333   Double_t pard0zmisal[3]={0,0,0};
 334   if(misal=="data") {
 335     //use this to reproduce data d0(pt) trend for pions
 336     pard0rphimisal[0]=37.;
 337     pard0rphimisal[1]=37.5;
 338     pard0rphimisal[2]=1.25;
 339     pard0zmisal[0]=96.;
 340     pard0zmisal[1]=131.;
 341     pard0zmisal[2]=0.7;
 342   }
 343   else if(misal=="resB") {
 344     // temporary values: asymptotic value larger by a factor 1.2 w.r.t. MC
 345     pard0rphimisal[0]=44.4;
 346     pard0rphimisal[1]=37.5;
 347     pard0rphimisal[2]=1.25;
 348     pard0zmisal[0]=115.2;
 349     pard0zmisal[1]=131.;
 350     pard0zmisal[2]=0.7;
 351   }
 352   else if(misal=="resC") {
 353     // temporary values: slightly larger asymptotic value, larger values at low pt
 354     pard0rphimisal[0]=40.;
 355     pard0rphimisal[1]=40.;
 356     pard0rphimisal[2]=1.3;
 357     pard0zmisal[0]=125.;
 358     pard0zmisal[1]=131.;
 359     pard0zmisal[2]=0.85;
 360   }
 361   else printf("AliAODRecoDecayHF::Misalign():  wrong misalign type specified \n");
 362
 363
 364   AliAODVertex *evVtx=0x0,*secVtx=0x0;
 365   Double_t evVtxPos[3]={-9999.,-9999.,-9999.},secVtxPos[3]={-9999.,9999.,9999.};
 366   if(fOwnPrimaryVtx)fOwnPrimaryVtx->GetXYZ(evVtxPos);
 367   else {
 368     evVtx=(AliAODVertex*)(fEventPrimaryVtx.GetObject());
 369     evVtx->GetXYZ(evVtxPos);
 370   }
 371   secVtx=(AliAODVertex*)GetSecondaryVtx();
 372   secVtx->GetXYZ(secVtxPos);
 373
 374   TVector3 v2v1(secVtxPos[0]-evVtxPos[0],secVtxPos[1]-evVtxPos[1],0.);
 375
 376   Double_t sigmarphinull,sigmarphimisal,sigmarphiadd;
 377   Double_t sigmaznull,sigmazmisal,sigmazadd;
 378   Double_t deltad0rphi[10],deltad0z[10];
 379
 380   // loop on the two prongs
 381   for(Int_t i=0; i<fNProngs; i++) {
 382     sigmarphinull = pard0rphiMC[0]+pard0rphiMC[1]/TMath::Power(PtProng(i),pard0rphiMC[2]);
 383     sigmarphimisal = pard0rphimisal[0]+pard0rphimisal[1]/TMath::Power(PtProng(i),pard0rphimisal[2]);
 384     if(sigmarphimisal>sigmarphinull) {
 385       sigmarphiadd = TMath::Sqrt(sigmarphimisal*sigmarphimisal-
 386                                  sigmarphinull*sigmarphinull);
 387       deltad0rphi[i] = gRandom->Gaus(0.,sigmarphiadd);
 388     } else {
 389       deltad0rphi[i] = 0.;
 390     }
 391
 392     sigmaznull =  pard0zMC[0]+pard0zMC[1]/TMath::Power(PtProng(i),pard0zMC[2]);
 393     sigmazmisal = pard0zmisal[0]+pard0zmisal[1]/TMath::Power(PtProng(i),pard0zmisal[2]);
 394     if(sigmazmisal>sigmaznull) {
 395       sigmazadd = TMath::Sqrt(sigmazmisal*sigmazmisal-
 396                               sigmaznull*sigmaznull);
 397       deltad0z[i] = gRandom->Gaus(0.,sigmazadd);
 398     } else {
 399       deltad0z[i] = 0.;
 400     }
 401
 402     TVector3 pxy(fPx[i],fPy[i],0.);
 403     TVector3 pxycrossv2v1=pxy.Cross(v2v1);
 404     if( pxycrossv2v1.Z()*fd0[i] > 0 ) {
 405       secVtxPos[0]+=1.e-4*deltad0rphi[i]*(-fPy[i])/PtProng(i);// e-4: conversion to cm
 406       secVtxPos[1]+=1.e-4*deltad0rphi[i]*(+fPx[i])/PtProng(i);
 407     } else {
 408       secVtxPos[0]+=1.e-4*deltad0rphi[i]*(+fPy[i])/PtProng(i);
 409       secVtxPos[1]+=1.e-4*deltad0rphi[i]*(-fPx[i])/PtProng(i);
 410     }
 411
 412     // change d0rphi
 413     fd0[i] += 1.e-4*deltad0rphi[i]; // e-4: conversion to cm
 414     // change secondary vertex z
 415     secVtxPos[2]+=0.5e-4*deltad0z[i];
 416   }
 417   secVtx->SetX(secVtxPos[0]);
 418   secVtx->SetY(secVtxPos[1]);
 419   secVtx->SetZ(secVtxPos[2]);
 420
 421   return;
 422 }