PWG3/vertexingHF/AliRDHFCuts.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2010, 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 cuts on AOD reconstructed heavy-flavour decay
  19 //
  20 // Author: A.Dainese, andrea.dainese@pd.infn.it
  21 /////////////////////////////////////////////////////////////
  22 #include <Riostream.h>
  23
  24 #include "AliVEvent.h"
  25 #include "AliESDEvent.h"
  26 #include "AliAODEvent.h"
  27 #include "AliVVertex.h"
  28 #include "AliESDVertex.h"
  29 #include "AliAODVertex.h"
  30 #include "AliESDtrack.h"
  31 #include "AliAODTrack.h"
  32 #include "AliESDtrackCuts.h"
  33 #include "AliAODRecoDecayHF.h"
  34 #include "AliRDHFCuts.h"
  35
  36 ClassImp(AliRDHFCuts)
  37
  38 //--------------------------------------------------------------------------
  39 AliRDHFCuts::AliRDHFCuts(const Char_t* name, const Char_t* title) :
  40 AliAnalysisCuts(name,title),
  41 fMinVtxType(3),
  42 fMinVtxContr(1),
  43 fMaxVtxRedChi2(1e6),
  44 fMinSPDMultiplicity(0),
  45 fTriggerMask(0),
  46 fTrackCuts(0),
  47 fnPtBins(1),
  48 fnPtBinLimits(1),
  49 fPtBinLimits(0),
  50 fnVars(1),
  51 fVarNames(0),
  52 fnVarsForOpt(0),
  53 fVarsForOpt(0),
  54 fGlobalIndex(1),
  55 fCutsRD(0),
  56 fIsUpperCut(0),
  57 fUsePID(kFALSE),
  58 fPidHF(0),
  59 fWhyRejection(0),
  60 fRemoveDaughtersFromPrimary(kFALSE)
  61 {
  62   //
  63   // Default Constructor
  64   //
  65 }
  66 //--------------------------------------------------------------------------
  67 AliRDHFCuts::AliRDHFCuts(const AliRDHFCuts &source) :
  68   AliAnalysisCuts(source),
  69   fMinVtxType(source.fMinVtxType),
  70   fMinVtxContr(source.fMinVtxContr),
  71   fMaxVtxRedChi2(source.fMaxVtxRedChi2),
  72   fMinSPDMultiplicity(source.fMinSPDMultiplicity),
  73   fTriggerMask(source.fTriggerMask),
  74   fTrackCuts(0),
  75   fnPtBins(source.fnPtBins),
  76   fnPtBinLimits(source.fnPtBinLimits),
  77   fPtBinLimits(0),
  78   fnVars(source.fnVars),
  79   fVarNames(0),
  80   fnVarsForOpt(source.fnVarsForOpt),
  81   fVarsForOpt(0),
  82   fGlobalIndex(source.fGlobalIndex),
  83   fCutsRD(0),
  84   fIsUpperCut(0),
  85   fUsePID(source.fUsePID),
  86   fPidHF(0),
  87   fWhyRejection(source.fWhyRejection),
  88   fRemoveDaughtersFromPrimary(source.fRemoveDaughtersFromPrimary)
  89 {
  90   //
  91   // Copy constructor
  92   //
  93   cout<<"Copy constructor"<<endl;
  94   if(source.GetTrackCuts()) AddTrackCuts(source.GetTrackCuts());
  95   if(source.fPtBinLimits) SetPtBins(source.fnPtBinLimits,source.fPtBinLimits);
  96   if(source.fVarNames) SetVarNames(source.fnVars,source.fVarNames,source.fIsUpperCut);
  97   if(source.fCutsRD) SetCuts(source.fGlobalIndex,source.fCutsRD);
  98   if(source.fVarsForOpt) SetVarsForOpt(source.fnVarsForOpt,source.fVarsForOpt);
  99   if(source.fPidHF) SetPidHF(source.fPidHF);
 100   PrintAll();
 101
 102 }
 103 //--------------------------------------------------------------------------
 104 AliRDHFCuts &AliRDHFCuts::operator=(const AliRDHFCuts &source)
 105 {
 106   //
 107   // assignment operator
 108   //
 109   if(&source == this) return *this;
 110
 111   AliAnalysisCuts::operator=(source);
 112
 113   fMinVtxType=source.fMinVtxType;
 114   fMinVtxContr=source.fMinVtxContr;
 115   fMaxVtxRedChi2=source.fMaxVtxRedChi2;
 116   fMinSPDMultiplicity=source.fMinSPDMultiplicity;
 117   fTriggerMask=source.fTriggerMask;
 118   fnPtBins=source.fnPtBins;
 119   fnVars=source.fnVars;
 120   fGlobalIndex=source.fGlobalIndex;
 121   fnVarsForOpt=source.fnVarsForOpt;
 122   fUsePID=source.fUsePID;
 123   SetPidHF(source.GetPidHF());
 124   fWhyRejection=source.fWhyRejection;
 125   fRemoveDaughtersFromPrimary=source.fRemoveDaughtersFromPrimary;
 126
 127   if(source.GetTrackCuts()) AddTrackCuts(source.GetTrackCuts());
 128   if(source.fPtBinLimits) SetPtBins(source.fnPtBinLimits,source.fPtBinLimits);
 129   if(source.fVarNames) SetVarNames(source.fnVars,source.fVarNames,source.fIsUpperCut);
 130   if(source.fCutsRD) SetCuts(source.fGlobalIndex,source.fCutsRD);
 131   if(source.fVarsForOpt) SetVarsForOpt(source.fnVarsForOpt,source.fVarsForOpt);
 132   PrintAll();
 133
 134   return *this;
 135 }
 136 //--------------------------------------------------------------------------
 137 AliRDHFCuts::~AliRDHFCuts() {
 138   //
 139   // Default Destructor
 140   //
 141   if(fPtBinLimits) {delete [] fPtBinLimits; fPtBinLimits=0;}
 142   if(fVarNames) {delete [] fVarNames; fVarNames=0;}
 143   if(fVarsForOpt) {delete [] fVarsForOpt; fVarsForOpt=0;}
 144   if(fCutsRD) {
 145     delete [] fCutsRD;
 146     fCutsRD=0;
 147   }
 148   if(fIsUpperCut) {delete [] fIsUpperCut; fIsUpperCut=0;}
 149   if(fPidHF){
 150     delete fPidHF;
 151     fPidHF=0;
 152   }
 153 }
 154 //---------------------------------------------------------------------------
 155 Bool_t AliRDHFCuts::IsEventSelected(AliVEvent *event) const {
 156   //
 157   // Event selection
 158   //
 159   //if(fTriggerMask && event->GetTriggerMask()!=fTriggerMask) return kFALSE;
 160
 161   // multiplicity cuts no implemented yet
 162
 163   const AliVVertex *vertex = event->GetPrimaryVertex();
 164
 165   if(!vertex) return kFALSE;
 166
 167   TString title=vertex->GetTitle();
 168   if(title.Contains("Z") && fMinVtxType>1) return kFALSE;
 169   if(title.Contains("3D") && fMinVtxType>2) return kFALSE;
 170
 171   if(vertex->GetNContributors()<fMinVtxContr) return kFALSE;
 172
 173   return kTRUE;
 174 }
 175 //---------------------------------------------------------------------------
 176 Bool_t AliRDHFCuts::AreDaughtersSelected(AliAODRecoDecayHF *d) const {
 177   //
 178   // Daughter tracks selection
 179   //
 180   if(!fTrackCuts) return kTRUE;
 181
 182   Int_t ndaughters = d->GetNDaughters();
 183   AliAODVertex *vAOD = d->GetPrimaryVtx();
 184   Double_t pos[3],cov[6];
 185   vAOD->GetXYZ(pos);
 186   vAOD->GetCovarianceMatrix(cov);
 187   const AliESDVertex vESD(pos,cov,100.,100);
 188
 189   Bool_t retval=kTRUE;
 190
 191   for(Int_t idg=0; idg<ndaughters; idg++) {
 192     AliAODTrack *dgTrack = (AliAODTrack*)d->GetDaughter(idg);
 193     if(!dgTrack) retval = kFALSE;
 194     //printf("charge %d\n",dgTrack->Charge());
 195     if(dgTrack->Charge()==0) continue; // it's not a track, but a V0
 196
 197     if(!IsDaughterSelected(dgTrack,&vESD,fTrackCuts)) retval = kFALSE;
 198   }
 199
 200   return retval;
 201 }
 202 //---------------------------------------------------------------------------
 203 Bool_t AliRDHFCuts::IsDaughterSelected(AliAODTrack *track,const AliESDVertex *primary,AliESDtrackCuts *cuts) const {
 204   //
 205   // Convert to ESDtrack, relate to vertex and check cuts
 206   //
 207   if(!cuts) return kTRUE;
 208
 209   Bool_t retval=kTRUE;
 210
 211   // convert to ESD track here
 212   AliESDtrack esdTrack(track);
 213   // needed to calculate the impact parameters
 214   esdTrack.RelateToVertex(primary,0.,3.);
 215   if(!cuts->IsSelected(&esdTrack)) retval = kFALSE;
 216
 217   return retval;
 218 }
 219 //---------------------------------------------------------------------------
 220 void AliRDHFCuts::SetPtBins(Int_t nPtBinLimits,Float_t *ptBinLimits) {
 221   // Set the pt bins
 222
 223   if(fPtBinLimits) {
 224     delete [] fPtBinLimits;
 225     fPtBinLimits = NULL;
 226     printf("Changing the pt bins\n");
 227   }
 228
 229   if(nPtBinLimits != fnPtBins+1){
 230     cout<<"Warning: ptBinLimits dimention "<<nPtBinLimits<<" != nPtBins+1 ("<<fnPtBins+1<<")\nSetting nPtBins to "<<nPtBinLimits-1<<endl;
 231     SetNPtBins(nPtBinLimits-1);
 232   }
 233
 234   fnPtBinLimits = nPtBinLimits;
 235   SetGlobalIndex();
 236   cout<<"Changing also Global Index -> "<<fGlobalIndex<<endl;
 237   fPtBinLimits = new Float_t[fnPtBinLimits];
 238   for(Int_t ib=0; ib<nPtBinLimits; ib++) fPtBinLimits[ib]=ptBinLimits[ib];
 239
 240   return;
 241 }
 242 //---------------------------------------------------------------------------
 243 void AliRDHFCuts::SetVarNames(Int_t nVars,TString *varNames,Bool_t *isUpperCut){
 244   // Set the variable names
 245
 246   if(fVarNames) {
 247     delete [] fVarNames;
 248     fVarNames = NULL;
 249     printf("Changing the variable names\n");
 250   }
 251   if(nVars!=fnVars){
 252     printf("Wrong number of variables: it has to be %d\n",fnVars);
 253     return;
 254   }
 255   //fnVars=nVars;
 256   fVarNames = new TString[nVars];
 257   fIsUpperCut = new Bool_t[nVars];
 258   for(Int_t iv=0; iv<nVars; iv++) {
 259     fVarNames[iv] = varNames[iv];
 260     fIsUpperCut[iv] = isUpperCut[iv];
 261   }
 262
 263   return;
 264 }
 265 //---------------------------------------------------------------------------
 266 void AliRDHFCuts::SetVarsForOpt(Int_t nVars,Bool_t *forOpt) {
 267   // Set the variables to be used for cuts optimization
 268
 269   if(fVarsForOpt) {
 270     delete [] fVarsForOpt;
 271     fVarsForOpt = NULL;
 272     printf("Changing the variables for cut optimization\n");
 273   }
 274
 275   if(nVars==0){//!=fnVars) {
 276     printf("%d not accepted as number of variables: it has to be %d\n",nVars,fnVars);
 277     return;
 278   }
 279
 280   fnVarsForOpt = 0;
 281   fVarsForOpt = new Bool_t[fnVars];
 282   for(Int_t iv=0; iv<fnVars; iv++) {
 283     fVarsForOpt[iv]=forOpt[iv];
 284     if(fVarsForOpt[iv]) fnVarsForOpt++;
 285   }
 286
 287   return;
 288 }
 289 //---------------------------------------------------------------------------
 290 void AliRDHFCuts::SetCuts(Int_t nVars,Int_t nPtBins,Float_t **cutsRD) {
 291   //
 292   // store the cuts
 293   //
 294   if(nVars!=fnVars) {
 295     printf("Wrong number of variables: it has to be %d\n",fnVars);
 296     return;
 297   }
 298   if(nPtBins!=fnPtBins) {
 299     printf("Wrong number of pt bins: it has to be %d\n",fnPtBins);
 300     return;
 301   }
 302
 303   if(!fCutsRD)  fCutsRD = new Float_t[fGlobalIndex];
 304
 305
 306   for(Int_t iv=0; iv<fnVars; iv++) {
 307
 308     for(Int_t ib=0; ib<fnPtBins; ib++) {
 309
 310       //check
 311       if(GetGlobalIndex(iv,ib)>=fGlobalIndex) {
 312         cout<<"Overflow, exit..."<<endl;
 313         return;
 314       }
 315
 316       fCutsRD[GetGlobalIndex(iv,ib)] = cutsRD[iv][ib];
 317
 318     }
 319   }
 320   return;
 321 }
 322 //---------------------------------------------------------------------------
 323 void AliRDHFCuts::SetCuts(Int_t glIndex,Float_t* cutsRDGlob){
 324   //
 325   // store the cuts
 326   //
 327   if(glIndex != fGlobalIndex){
 328     cout<<"Wrong array size: it has to be "<<fGlobalIndex<<endl;
 329     return;
 330   }
 331   if(!fCutsRD)  fCutsRD = new Float_t[fGlobalIndex];
 332
 333   for(Int_t iGl=0;iGl<fGlobalIndex;iGl++){
 334     fCutsRD[iGl] = cutsRDGlob[iGl];
 335   }
 336   return;
 337 }
 338 //---------------------------------------------------------------------------
 339 void AliRDHFCuts::PrintAll() const {
 340   //
 341   // print all cuts values
 342   //
 343   if(fVarNames){
 344     cout<<"Array of variables"<<endl;
 345     for(Int_t iv=0;iv<fnVars;iv++){
 346       cout<<fVarNames[iv]<<"\t";
 347     }
 348     cout<<endl;
 349   }
 350   if(fVarsForOpt){
 351     cout<<"Array of optimization"<<endl;
 352     for(Int_t iv=0;iv<fnVars;iv++){
 353       cout<<fVarsForOpt[iv]<<"\t";
 354     }
 355     cout<<endl;
 356   }
 357   if(fIsUpperCut){
 358     cout<<"Array of upper/lower cut"<<endl;
 359    for(Int_t iv=0;iv<fnVars;iv++){
 360      cout<<fIsUpperCut[iv]<<"\t";
 361    }
 362    cout<<endl;
 363   }
 364   if(fPtBinLimits){
 365     cout<<"Array of ptbin limits"<<endl;
 366     for(Int_t ib=0;ib<fnPtBinLimits;ib++){
 367       cout<<fPtBinLimits[ib]<<"\t";
 368     }
 369     cout<<endl;
 370   }
 371   if(fCutsRD){
 372     cout<<"Matrix of cuts"<<endl;
 373    for(Int_t iv=0;iv<fnVars;iv++){
 374      for(Int_t ib=0;ib<fnPtBins;ib++){
 375        cout<<"fCutsRD["<<iv<<"]["<<ib<<"] = "<<fCutsRD[GetGlobalIndex(iv,ib)]<<"\t";
 376      }
 377      cout<<endl;
 378    }
 379    cout<<endl;
 380   }
 381   return;
 382 }
 383 //---------------------------------------------------------------------------
 384 void AliRDHFCuts::GetCuts(Float_t**& cutsRD) const{
 385   //
 386   // get the cuts
 387   //
 388
 389   //cout<<"Give back a "<<fnVars<<"x"<<fnPtBins<<" matrix."<<endl;
 390
 391
 392   Int_t iv,ib;
 393   if(!cutsRD) {
 394     //cout<<"Initialization..."<<endl;
 395     cutsRD=new Float_t*[fnVars];
 396     for(iv=0; iv<fnVars; iv++) {
 397       cutsRD[iv] = new Float_t[fnPtBins];
 398     }
 399   }
 400
 401   for(Int_t iGlobal=0; iGlobal<fGlobalIndex; iGlobal++) {
 402     GetVarPtIndex(iGlobal,iv,ib);
 403     cutsRD[iv][ib] = fCutsRD[iGlobal];
 404   }
 405
 406   return;
 407 }
 408
 409 //---------------------------------------------------------------------------
 410 Int_t AliRDHFCuts::GetGlobalIndex(Int_t iVar,Int_t iPtBin) const{
 411   //
 412   // give the global index from variable and pt bin
 413   //
 414   return iPtBin*fnVars+iVar;
 415 }
 416
 417 //---------------------------------------------------------------------------
 418 void AliRDHFCuts::GetVarPtIndex(Int_t iGlob, Int_t& iVar, Int_t& iPtBin) const {
 419   //
 420   //give the index of the variable and of the pt bin from the global index
 421   //
 422   iPtBin=(Int_t)iGlob/fnVars;
 423   iVar=iGlob%fnVars;
 424
 425   return;
 426 }
 427
 428 //---------------------------------------------------------------------------
 429 Int_t AliRDHFCuts::PtBin(Double_t pt) const {
 430   //
 431   //give the pt bin where the pt lies.
 432   //
 433   Int_t ptbin=-1;
 434   for (Int_t i=0;i<fnPtBins;i++){
 435     if(pt<fPtBinLimits[i+1]) {
 436       ptbin=i;
 437       break;
 438     }
 439   }
 440   return ptbin;
 441 }
 442 //-------------------------------------------------------------------
 443 Float_t AliRDHFCuts::GetCutValue(Int_t iVar,Int_t iPtBin) const {
 444   //
 445   // Give the value of cut set for the variable iVar and the pt bin iPtBin
 446   //
 447   if(!fCutsRD){
 448     cout<<"Cuts not iniziaisez yet"<<endl;
 449     return 0;
 450   }
 451   return fCutsRD[GetGlobalIndex(iVar,iPtBin)];
 452 }
 453