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