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