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