PWGDQ/dielectron/AliDielectronEventCuts.cxx

   1 /*************************************************************************
   2 * Copyright(c) 1998-2009, 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 //                Dielectron EventCuts                                  //
  18 //                                                                       //
  19 //                                                                       //
  20 /*
  21 Detailed description
  22
  23
  24 */
  25 //                                                                       //
  26 ///////////////////////////////////////////////////////////////////////////
  27
  28
  29 #include <AliTriggerAnalysis.h>
  30 #include <AliESDVertex.h>
  31 #include <AliAODVertex.h>
  32 #include <AliESDEvent.h>
  33 #include <AliAODEvent.h>
  34 #include <AliMultiplicity.h>
  35 #include <AliCentrality.h>
  36 #include <AliESDVZERO.h>
  37 #include <AliAODVZERO.h>
  38
  39 #include "AliDielectronVarManager.h"
  40 #include "AliDielectronEventCuts.h"
  41
  42 ClassImp(AliDielectronEventCuts)
  43
  44 const char* AliDielectronEventCuts::fgkVtxNames[AliDielectronEventCuts::kVtxTracksOrSPD+1] = {"Tracks", "SPD", "TPC", "Any", "TracksOrSPD"};
  45
  46 AliDielectronEventCuts::AliDielectronEventCuts() :
  47   AliAnalysisCuts(),
  48   fUsedVars(new TBits(AliDielectronVarManager::kNMaxValues)),
  49   fRun(),
  50   fVtxZmin(0.),
  51   fVtxZmax(0.),
  52   fRequireVtx(kFALSE),
  53   fMinVtxContributors(0),
  54   fMultITSTPC(kFALSE),
  55   fCentMin(1.),
  56   fCentMax(0.),
  57   fVtxType(kVtxTracks),
  58   fRequire13sel(kFALSE),
  59   fUtils(),
  60   fRequireV0and(0),
  61   fTriggerAnalysis(0x0),
  62   fkVertex(0x0),
  63   fkVertexAOD(0x0),
  64   fparMean(0x0),
  65   fparSigma(0x0),
  66   fcutSigma(3.),
  67   fparMinVtxContributors(0x0),
  68   fparMaxVtxContributors(0x0)
  69 {
  70   //
  71   // Default Constructor
  72   //
  73   for (Int_t icut=0; icut<5; ++icut){
  74     fCorrCutMin[icut]=0x0;
  75     fCorrCutMax[icut]=0x0;
  76   }
  77 }
  78
  79 //______________________________________________
  80 AliDielectronEventCuts::AliDielectronEventCuts(const char* name, const char* title) :
  81   AliAnalysisCuts(name, title),
  82   fUsedVars(new TBits(AliDielectronVarManager::kNMaxValues)),
  83   fRun(),
  84   fVtxZmin(0.),
  85   fVtxZmax(0.),
  86   fRequireVtx(kFALSE),
  87   fMinVtxContributors(0),
  88   fMultITSTPC(kFALSE),
  89   fCentMin(1.),
  90   fCentMax(0.),
  91   fVtxType(kVtxTracks),
  92   fRequire13sel(kFALSE),
  93   fUtils(),
  94   fRequireV0and(0),
  95   fTriggerAnalysis(0x0),
  96   fkVertex(0x0),
  97   fkVertexAOD(0x0),
  98   fparMean(0x0),
  99   fparSigma(0x0),
 100   fcutSigma(3.),
 101   fparMinVtxContributors(0x0),
 102   fparMaxVtxContributors(0x0)
 103 {
 104   //
 105   // Named Constructor
 106   //
 107   for (Int_t icut=0; icut<5; ++icut){
 108     fCorrCutMin[icut]=0x0;
 109     fCorrCutMax[icut]=0x0;
 110   }
 111 }
 112
 113 //______________________________________________
 114 AliDielectronEventCuts::~AliDielectronEventCuts()
 115 {
 116   //
 117   // Default Destructor
 118   //
 119   if (fUsedVars) delete fUsedVars;
 120   if (fTriggerAnalysis) delete fTriggerAnalysis;
 121 }
 122
 123 //______________________________________________
 124 Bool_t AliDielectronEventCuts::IsSelected(TObject* event)
 125 {
 126   //
 127   // check the cuts
 128   //
 129
 130   if(event->IsA() == AliESDEvent::Class()) return IsSelectedESD(event);
 131   else if(event->IsA() == AliAODEvent::Class()) return IsSelectedAOD(event);
 132   else return kFALSE;
 133 }
 134 //____________________________________________________________________
 135 Bool_t AliDielectronEventCuts::IsSelectedESD(TObject* event)
 136 {
 137   //
 138   // check the cuts
 139   //
 140
 141   AliESDEvent *ev=dynamic_cast<AliESDEvent*>(event);
 142   if (!ev) return kFALSE;
 143
 144   if (fCentMin<fCentMax){
 145     AliCentrality *centrality=ev->GetCentrality();
 146     Double_t centralityF=-1;
 147     if (centrality) centralityF = centrality->GetCentralityPercentile("V0M");
 148     if (centralityF<fCentMin || centralityF>=fCentMax) return kFALSE;
 149   }
 150
 151   fkVertex=0x0;
 152
 153   switch(fVtxType){
 154   case kVtxTracks:
 155   case kVtxTracksOrSPD:
 156     fkVertex=ev->GetPrimaryVertexTracks();
 157     break;
 158   case kVtxSPD:    fkVertex=ev->GetPrimaryVertexSPD(); break;
 159   case kVtxTPC:    fkVertex=ev->GetPrimaryVertexTPC(); break;
 160   case kVtxAny:    fkVertex=ev->GetPrimaryVertex(); break;
 161   }
 162
 163   if ((fRequireVtx||fVtxZmin<fVtxZmax||fMinVtxContributors>0)&&!fkVertex) return kFALSE;
 164
 165
 166   if (fMinVtxContributors>0){
 167     Int_t nCtrb = fkVertex->GetNContributors();
 168     if (nCtrb<fMinVtxContributors){
 169       if (fVtxType==kVtxTracksOrSPD){
 170         fkVertex=ev->GetPrimaryVertexSPD();
 171         nCtrb = fkVertex->GetNContributors();
 172         if (nCtrb<fMinVtxContributors) return kFALSE;
 173       } else {
 174         return kFALSE;
 175       }
 176     }
 177   }
 178
 179   if (fVtxZmin<fVtxZmax){
 180     Double_t zvtx=fkVertex->GetZv();
 181     if (zvtx<fVtxZmin||zvtx>fVtxZmax) return kFALSE;
 182   }
 183
 184   if(fRequire13sel){
 185     if(!fUtils.IsVertexSelected2013pA(ev)) return kFALSE;
 186     if(fUtils.IsFirstEventInChunk(ev)) return kFALSE;
 187   }
 188
 189   if (fRequireV0and){
 190     if (!fTriggerAnalysis) fTriggerAnalysis=new AliTriggerAnalysis;
 191     Bool_t v0AND = kFALSE;
 192     if (fRequireV0and==1){
 193       Bool_t v0A       = fTriggerAnalysis->IsOfflineTriggerFired(ev, AliTriggerAnalysis::kV0A);
 194       Bool_t v0C       = fTriggerAnalysis->IsOfflineTriggerFired(ev, AliTriggerAnalysis::kV0C);
 195       v0AND = v0A && v0C;
 196     }
 197
 198     if (fRequireV0and==2){
 199       Bool_t v0AHW     = (fTriggerAnalysis->V0Trigger(ev, AliTriggerAnalysis::kASide, kTRUE) == AliTriggerAnalysis::kV0BB);
 200       Bool_t v0CHW     = (fTriggerAnalysis->V0Trigger(ev, AliTriggerAnalysis::kCSide, kTRUE) == AliTriggerAnalysis::kV0BB);
 201       v0AND = v0AHW && v0CHW;
 202     }
 203
 204     if (!v0AND) return kFALSE;
 205   }
 206
 207   if (fMultITSTPC){
 208     const AliESDVertex *vtxESDTPC=ev->GetPrimaryVertexTPC();
 209     const AliMultiplicity *multESD = ev->GetMultiplicity();
 210     if ( vtxESDTPC && multESD && vtxESDTPC->GetNContributors() < (-10.+0.25*multESD->GetNumberOfITSClusters(0)) )
 211       return kFALSE;
 212   }
 213
 214   if(fparMean && fparSigma) {
 215     Double_t nTrks  = ev->GetNumberOfTracks();
 216     Double_t multV0 = 0.0;
 217     for(Int_t j=0; j<64; j++) multV0 += ev->GetVZEROData()->GetMultiplicity(j);
 218     Double_t mV0 = fparMean->Eval(nTrks);
 219     Double_t sV0 = fparSigma->Eval(nTrks);
 220     if(multV0 > mV0+fcutSigma*sV0 || multV0 < mV0-fcutSigma*sV0) return kFALSE;
 221   }
 222
 223   // cut on the number of vertex contributors using TPC versus global vertex
 224   if(fparMinVtxContributors && fparMaxVtxContributors) {
 225     const AliESDVertex *vtxTPC = ev->GetPrimaryVertexTPC();
 226     const AliESDVertex *vtxGbl = ev->GetPrimaryVertex();
 227     Double_t nContribTPC = (vtxTPC ? vtxTPC->GetNContributors() : 0);
 228     Double_t nContribGbl = (vtxGbl ? vtxGbl->GetNContributors() : 0);
 229     Double_t minCut = fparMinVtxContributors->Eval(nContribGbl);
 230     Double_t maxCut = fparMaxVtxContributors->Eval(nContribGbl);
 231     if(nContribTPC > maxCut || nContribTPC < minCut) return kFALSE;
 232   }
 233
 234
 235   return kTRUE;
 236 }
 237 //______________________________________________
 238 Bool_t AliDielectronEventCuts::IsSelectedAOD(TObject* event)
 239 {
 240   //
 241   // check the cuts
 242   //
 243
 244   AliAODEvent *ev=dynamic_cast<AliAODEvent*>(event);
 245   if (!ev) return kFALSE;
 246
 247   //Fill values
 248   Double_t values[AliDielectronVarManager::kNMaxValues];
 249   if(fUsedVars->CountBits()) {
 250     AliDielectronVarManager::SetFillMap(fUsedVars);
 251     AliDielectronVarManager::Fill(ev,values);
 252
 253     // correlation cuts
 254     for(Int_t i=0; i<5; i++) {
 255       if(fCorrCutMin[i]) {
 256         Double_t varx = values[fCorrCutMin[i]->GetXaxis()->GetUniqueID()];
 257         Double_t vary = values[fCorrCutMin[i]->GetYaxis()->GetUniqueID()];
 258         Double_t min  = ((TF1*)fCorrCutMin[i]->GetListOfFunctions()->At(0))->Eval(varx);
 259         //      printf("coor cut %d: varx %f -> eval %f > %f \n",i,varx,min,vary);
 260         if(vary<min) return kFALSE;
 261       }
 262       if(fCorrCutMax[i]) {
 263         Double_t varx = values[fCorrCutMax[i]->GetXaxis()->GetUniqueID()];
 264         Double_t vary = values[fCorrCutMax[i]->GetYaxis()->GetUniqueID()];
 265         Double_t max  = ((TF1*)fCorrCutMax[i]->GetListOfFunctions()->At(0))->Eval(varx);
 266         if(vary>max) return kFALSE;
 267       }
 268     }
 269   }
 270
 271   // run rejection
 272   Int_t run = ev->GetRunNumber();
 273   if(fRun.GetNrows()) {
 274     for(Int_t irun=0; irun<fRun.GetNrows(); irun++) {
 275       if(fRun(irun)==run) return kFALSE;
 276     }
 277   }
 278
 279   if (fCentMin<fCentMax){
 280     AliCentrality *centrality=ev->GetCentrality();
 281     Double_t centralityF=-1;
 282     if (centrality) centralityF = centrality->GetCentralityPercentile("V0M");
 283     if (centralityF<fCentMin || centralityF>=fCentMax) return kFALSE;
 284   }
 285
 286   fkVertexAOD=0x0;
 287
 288   switch(fVtxType){
 289   case kVtxTracks:         fkVertexAOD=0x0;                       break;
 290   case kVtxTPC:            fkVertexAOD=AliDielectronVarManager::GetVertex(ev, AliAODVertex::kMainTPC);   break;
 291   case kVtxSPD:
 292   case kVtxTracksOrSPD:    fkVertexAOD=ev->GetPrimaryVertexSPD(); break;
 293   case kVtxAny:            fkVertexAOD=ev->GetPrimaryVertex();    break;
 294   }
 295
 296   if ((fRequireVtx||fVtxZmin<fVtxZmax||fMinVtxContributors>0)&&!fkVertexAOD) return kFALSE;
 297
 298   if (fMinVtxContributors>0){
 299     Int_t nCtrb = fkVertexAOD->GetNContributors();
 300     if (nCtrb<fMinVtxContributors){
 301       // if (fVtxType==kVtxTracksOrSPD){
 302       //   fkVertexAOD=ev->GetVertex(AliAODVertex::kPrimary);
 303       //   nCtrb = fkVertexAOD->GetNContributors();
 304       //   if (nCtrb<fMinVtxContributors) return kFALSE;
 305       //      } else {
 306       return kFALSE;
 307       //}
 308     }
 309   }
 310
 311
 312   if (fVtxZmin<fVtxZmax){
 313     Double_t zvtx=fkVertexAOD->GetZ();
 314     if (zvtx<fVtxZmin||zvtx>fVtxZmax) return kFALSE;
 315   }
 316
 317   if(fRequire13sel){
 318     if(!fUtils.IsVertexSelected2013pA(ev)) return kFALSE;
 319 //     if(fUtils.IsFirstEventInChunk(ev)) return kFALSE;
 320   }
 321
 322   /*
 323   if (fRequireV0and){
 324     //    if (!fTriggerAnalysis) fTriggerAnalysis=new AliTriggerAnalysis;
 325     Bool_t v0AND = kFALSE;
 326     if (fRequireV0and==1){
 327       Bool_t v0A       = fTriggerAnalysis->IsOfflineTriggerFired(ev, AliTriggerAnalysis::kV0A);
 328       Bool_t v0A       = header->GetOfflineTrigger(); //TODO
 329       Bool_t v0C       = fTriggerAnalysis->IsOfflineTriggerFired(ev, AliTriggerAnalysis::kV0C);
 330       v0AND = v0A && v0C;
 331     }
 332
 333     if (fRequireV0and==2){
 334       Bool_t v0AHW     = (fTriggerAnalysis->V0Trigger(ev, AliTriggerAnalysis::kASide, kTRUE) == AliTriggerAnalysis::kV0BB);
 335       Bool_t v0CHW     = (fTriggerAnalysis->V0Trigger(ev, AliTriggerAnalysis::kCSide, kTRUE) == AliTriggerAnalysis::kV0BB);
 336       v0AND = v0AHW && v0CHW;
 337     }
 338
 339     if (!v0AND) return kFALSE;
 340   }
 341   */
 342   /*  if (fMultITSTPC){
 343     const AliESDVertex *vtxESDTPC=ev->GetPrimaryVertexTPC();
 344     const AliMultiplicity *multESD = ev->GetMultiplicity();
 345     if ( vtxESDTPC && multESD && vtxESDTPC->GetNContributors() < (-10.+0.25*multESD->GetNumberOfITSClusters(0)) )
 346       return kFALSE;
 347   }
 348   */
 349
 350   // correlation cut Ntrks vs. multV0
 351   if(fparMean && fparSigma) {
 352     Double_t nTrks  = ev->GetNumberOfTracks();
 353     Double_t multV0 = 0.0;
 354     for(Int_t j=0; j<64; j++) multV0 += ev->GetVZEROData()->GetMultiplicity(j);
 355     Double_t mV0 = fparMean->Eval(nTrks);
 356     Double_t sV0 = fparSigma->Eval(nTrks);
 357     if(multV0 > mV0+fcutSigma*sV0 || multV0 < mV0-fcutSigma*sV0) return kFALSE;
 358   }
 359
 360   // cut on the number of vertex contributors using TPC versus global vertex
 361   if(fparMinVtxContributors && fparMaxVtxContributors) {
 362     const AliAODVertex *vtxTPC = ev->GetVertex(AliAODVertex::kMainTPC);
 363     const AliAODVertex *vtxGbl = ev->GetPrimaryVertex();
 364     Double_t nContribTPC = (vtxTPC ? vtxTPC->GetNContributors() : 0);
 365     Double_t nContribGbl = (vtxGbl ? vtxGbl->GetNContributors() : 0);
 366     Double_t minCut = fparMinVtxContributors->Eval(nContribGbl);
 367     Double_t maxCut = fparMaxVtxContributors->Eval(nContribGbl);
 368     if(nContribTPC > maxCut || nContribTPC < minCut) return kFALSE;
 369   }
 370
 371   return kTRUE;
 372 }
 373
 374 //______________________________________________
 375 void AliDielectronEventCuts::SetMinCorrCutFunction(TF1 *fun, UInt_t varx, UInt_t vary)
 376 {
 377   //
 378   // add correlation cut using a TF1
 379   //
 380   fUsedVars->SetBitNumber(varx,kTRUE);
 381   fUsedVars->SetBitNumber(vary,kTRUE);
 382   // store variables
 383   fun->GetXaxis()->SetUniqueID(varx);
 384   fun->GetYaxis()->SetUniqueID(vary);
 385
 386   Int_t i=0;
 387   for(i=0; i<5; i++) {
 388     if(fCorrCutMin[i]) continue;
 389     else {
 390       TString key=GetName(); key+=Form("Min%d",i);
 391       // clone temporare histogram since otherwise it will not be streamed to file!
 392       fCorrCutMin[i] = (TH1*)fun->GetHistogram()->Clone(key.Data());
 393       fCorrCutMin[i]->GetListOfFunctions()->AddAt(fun,0);
 394       break;
 395     }
 396   }
 397   //printf("-----> corr cut added to %d %p \n",i,fCorrCutMin[i]);
 398   //  fCorrCutMin[i]->Print();
 399   //fCorrCutMin[i]->GetListOfFunctions()->ls();
 400 }
 401
 402 //______________________________________________
 403 void AliDielectronEventCuts::SetMaxCorrCutFunction(TF1 *fun, UInt_t varx, UInt_t vary)
 404 {
 405   //
 406   // add correlation cut using a TF1
 407   //
 408   fUsedVars->SetBitNumber(varx,kTRUE);
 409   fUsedVars->SetBitNumber(vary,kTRUE);
 410   // store variables
 411   fun->GetXaxis()->SetUniqueID(varx);
 412   fun->GetYaxis()->SetUniqueID(vary);
 413
 414   Int_t i=0;
 415   for(i=0; i<5; i++) {
 416     if(fCorrCutMax[i]) continue;
 417     else {
 418       TString key=GetName(); key+=Form("Max%d",i);
 419       // clone temporare histogram since otherwise it will not be streamed to file!
 420       fCorrCutMax[i] = (TH1*)fun->GetHistogram()->Clone(key.Data());
 421       fCorrCutMax[i]->GetListOfFunctions()->AddAt(fun,0);
 422       break;
 423     }
 424   }
 425   //printf("-----> corr cut added to %d %p \n",i,fCorrCutMax[i]);
 426   //  fCorrCutMax[i]->Print();
 427   //  fCorrCutMax[i]->GetListOfFunctions()->ls();
 428 }
 429
 430 //________________________________________________________________________
 431 void AliDielectronEventCuts::Print(const Option_t* /*option*/) const
 432 {
 433   //
 434   // Print cuts and the range
 435   //
 436   printf("cut ranges for '%s'\n",GetTitle());
 437   printf("All Cuts have to be fulfilled\n");
 438
 439   Int_t iCut=0;
 440   if(fRequireVtx) {
 441     printf("Cut %02d: vertex required \n",iCut);                                   iCut++; }
 442   printf("Cut %02d: vertex type: %s \n", iCut, fgkVtxNames[fVtxType]);             iCut++;
 443   if(fMinVtxContributors) {
 444     printf("Cut %02d: vertex contributors >= %d \n", iCut, fMinVtxContributors);   iCut++; }
 445   if(fVtxZmin<fVtxZmax) {
 446     printf("Cut %02d: %f < %s < %f\n",   iCut, fVtxZmin, "Zvtx", fVtxZmax);        iCut++; }
 447   if(fCentMin<fCentMax) {
 448     printf("Cut %02d: %f < %s < %f\n",   iCut, fCentMin, "V0centrality", fCentMax);iCut++; }
 449   if(fMultITSTPC) {
 450     printf("Cut %02d: cut on multiplcity ITS vs. TPC \n", iCut);                   iCut++; }
 451   if(fparMean&&fparSigma) {
 452     printf("Cut %02d: multplicity vs. #tracks correlation +-%.1f sigma inclusion \n", iCut, fcutSigma); iCut++; }
 453   if(fRequire13sel){
 454     printf("Cut %02d: vertex and event selection for 2013 pPb data taking required \n",iCut);   iCut++; }
 455   if(fRequireV0and) {
 456     printf("Cut %02d: require V0and type: %c \n", iCut, fRequireV0and);            iCut++; }
 457
 458 }
 459