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