PWGHF/hfe/AliHFEextraEventCuts.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, 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 // Cut on the Event at reconstructed level: for the moment
  16 // the requirements on the number of charged tracks and on
  17 // the vertex position and resolution are implemented
  18 // The argument of IsSelected member function (passed object) is cast into
  19 // an AliESDEvent. In the future may be modified to use AliVEvent interface
  20 // and include more cut variables.
  21 // The class derives from AliCFCutBase
  22 // Author:S.Arcelli Silvia.Arcelli@cern.ch
  23 //
  24 //
  25 #include "TH1F.h"
  26 #include "TList.h"
  27 #include "TBits.h"
  28 #include "AliLog.h"
  29 #include "AliAnalysisUtils.h"
  30 #include "AliESDEvent.h"
  31 #include "AliESDVertex.h"
  32 #include "AliAODEvent.h"
  33 #include "AliAODVertex.h"
  34 #include "AliHFEextraEventCuts.h"
  35 ClassImp(AliHFEextraEventCuts)
  36 //____________________________________________________________________
  37 AliHFEextraEventCuts::AliHFEextraEventCuts() :
  38   AliCFCutBase(),
  39   fAnalysisUtils(NULL),
  40   fRequireVtxCuts(kFALSE),
  41   fVtxZMax(1.e99),
  42   fVtxZMin(-1.e99),
  43   fVtxNCtrbMin(0),
  44   fVtxMixed(0),
  45   fVtxSPD(0),
  46   fCheckCorrelationSPDVtx(0),
  47   fVtxResolution(kFALSE),
  48   fPApileupCut(kFALSE),
  49   fBitMap(0x0)
  50 {
  51   //
  52   //ctor
  53   //
  54 }
  55
  56 //____________________________________________________________________
  57 AliHFEextraEventCuts::AliHFEextraEventCuts(Char_t* name, Char_t* title) :
  58   AliCFCutBase(name,title),
  59   fAnalysisUtils(NULL),
  60   fRequireVtxCuts(kFALSE),
  61   fVtxZMax(1.e99),
  62   fVtxZMin(-1.e99),
  63   fVtxNCtrbMin(0),
  64   fVtxMixed(0),
  65   fVtxSPD(0),
  66   fCheckCorrelationSPDVtx(0),
  67   fVtxResolution(kFALSE),
  68   fPApileupCut(kFALSE),
  69   fBitMap(0x0)
  70  {
  71   //
  72   //ctor
  73   //
  74   fBitMap=new TBits(0);
  75   Initialise();
  76   fAnalysisUtils = new AliAnalysisUtils;
  77  }
  78
  79 //____________________________________________________________________
  80 AliHFEextraEventCuts::AliHFEextraEventCuts(const AliHFEextraEventCuts& c) :
  81   AliCFCutBase(c),
  82   fAnalysisUtils(NULL),
  83   fRequireVtxCuts(c.fRequireVtxCuts),
  84   fVtxZMax(c.fVtxZMax),
  85   fVtxZMin(c.fVtxZMin),
  86   fVtxNCtrbMin(c.fVtxNCtrbMin),
  87   fVtxMixed(c.fVtxMixed),
  88   fVtxSPD(c.fVtxSPD),
  89   fCheckCorrelationSPDVtx(c.fCheckCorrelationSPDVtx),
  90   fVtxResolution(c.fVtxResolution),
  91   fPApileupCut(c.fPApileupCut),
  92   fBitMap(c.fBitMap)
  93 {
  94   //
  95   //copy constructor
  96   //
  97   for (Int_t i=0; i<c.kNCuts; i++){
  98     for (Int_t j=0; j<c.kNStepQA; j++){
  99       if(c.fhQA[i][j]) fhQA[i][j] = (TH1F*)c.fhQA[i][j]->Clone();
 100     }
 101   }
 102   fAnalysisUtils = new AliAnalysisUtils;
 103 }
 104
 105 //____________________________________________________________________
 106 AliHFEextraEventCuts::~AliHFEextraEventCuts() {
 107   //
 108   //dtor
 109   //
 110
 111   for (Int_t i=0; i<kNCuts; i++){
 112     for (Int_t j=0; j<kNStepQA; j++){
 113       if(fhQA[i][j]) delete fhQA[i][j];
 114     }
 115   }
 116
 117   if(fBitMap) delete fBitMap;
 118   if(fAnalysisUtils) delete fAnalysisUtils;
 119 }
 120 //_____________________________________________________________________________
 121 void AliHFEextraEventCuts::Initialise()
 122 {
 123
 124   //
 125   //initialization
 126   //
 127
 128   //
 129   // sets pointers to histos to zero
 130   //
 131
 132   for(Int_t i=0; i<kNCuts; i++){
 133     for(Int_t j =0; j<kNStepQA; j++){
 134       fhQA[i][j]=0x0;
 135     }
 136   }
 137 }
 138
 139 //____________________________________________________________________
 140 AliHFEextraEventCuts& AliHFEextraEventCuts::operator=(const AliHFEextraEventCuts& c)
 141 {
 142   //
 143   // Assignment operator
 144   //
 145   if (this != &c) {
 146     AliCFCutBase::operator=(c) ;
 147     fRequireVtxCuts=c.fRequireVtxCuts;
 148     fVtxZMax=c.fVtxZMax;
 149     fVtxZMin=c.fVtxZMin;
 150     fVtxNCtrbMin=c.fVtxNCtrbMin;
 151     fVtxMixed=c.fVtxMixed;
 152     fVtxSPD=c.fVtxSPD;
 153     fCheckCorrelationSPDVtx=c.fCheckCorrelationSPDVtx;
 154     fVtxResolution = c.fVtxResolution;
 155     fBitMap=c.fBitMap;
 156     if(!fAnalysisUtils) fAnalysisUtils = new AliAnalysisUtils;
 157     fPApileupCut = c.fPApileupCut;
 158   }
 159
 160   for (Int_t i=0; i<c.kNCuts; i++){
 161     for (Int_t j=0; j<c.kNStepQA; j++){
 162       if(c.fhQA[i][j]) fhQA[i][j] = (TH1F*)c.fhQA[i][j]->Clone();
 163     }
 164   }
 165
 166
 167   return *this ;
 168 }
 169
 170 //____________________________________________________________________
 171 Bool_t AliHFEextraEventCuts::IsSelected(TObject* obj) {
 172   //
 173   //Check if the requested cuts are passed
 174   //
 175
 176
 177   SelectionBitMap(obj);
 178
 179   if (fIsQAOn) FillHistograms(obj,0);
 180   Bool_t isSelected = kTRUE;
 181
 182   for (UInt_t icut=0; icut<fBitMap->GetNbits();icut++)
 183         if(!fBitMap->TestBitNumber(icut)) isSelected = kFALSE;
 184
 185   if (!isSelected) return kFALSE ;
 186   if (fIsQAOn) FillHistograms(obj,1);
 187   return kTRUE;
 188
 189 }
 190 //____________________________________________________________________
 191 void AliHFEextraEventCuts::SelectionBitMap(TObject* obj) {
 192   //
 193   //cut on the number of charged tracks and on the event vertex.
 194   //so far specific to AliESDEvents
 195   //
 196
 197   //Check if the requested cuts are passed and return a bitmap
 198
 199   AliVEvent *inputEvent = dynamic_cast<AliVEvent *>(obj);
 200   if(!inputEvent){
 201     AliDebug(1, "Not a virtual event");
 202     for(Int_t j=0;j<kNCuts;j++)fBitMap->SetBitNumber(j,kFALSE);
 203     return;
 204   }
 205   const AliVVertex *vtxTracks = GetPrimaryVertexTracks(inputEvent),
 206                    *vtxSPD = GetPrimaryVertexSPD(inputEvent),
 207                    *vtxPrim(NULL);
 208
 209   //first assume the event will be accepted:
 210   for(Int_t j=0;j<kNCuts;j++)fBitMap->SetBitNumber(j,kTRUE);
 211
 212   if(fVtxMixed){
 213     // Use mixed vertex: Prefer vertex with tracks, in case not available use SPD vertex
 214     if(vtxTracks && vtxTracks->GetNContributors() > 0) vtxPrim = vtxTracks;
 215     else if(vtxSPD && vtxSPD->GetNContributors() > 0) vtxPrim = vtxSPD;
 216   } else if(fVtxSPD){
 217     if(vtxSPD && vtxSPD->GetNContributors () > 0) vtxPrim = vtxSPD;
 218   } else {
 219     if(vtxTracks && vtxTracks->GetNContributors() > 0) vtxPrim = vtxTracks;
 220   }
 221   if(!vtxPrim){
 222     // No primary vertex: Reject event
 223           for(Int_t j=1;j<kNCuts;j++)fBitMap->SetBitNumber(j,kFALSE);
 224           AliWarning("Cannot get vertex, skipping event");
 225           return;
 226   }
 227   // Check standard vertex cuts using the primary vertex
 228   if(fVtxZMin > - 1000. && fVtxZMax < 1000.){
 229     // Primary vertex z cut required
 230     // Pick up the position and uncertainties
 231     Double_t vtxPos[3];
 232     vtxPos[0] = vtxPrim->GetX();
 233     vtxPos[1] = vtxPrim->GetY();
 234     vtxPos[2] = vtxPrim->GetZ();
 235     if (vtxPos[2]>fVtxZMax || vtxPos[2]<fVtxZMin)
 236             fBitMap->SetBitNumber(kVtxPosZ,kFALSE);
 237   }
 238   if(fVtxNCtrbMin > 0){
 239     // cut required if the cut value is set to something larger than 0
 240     // same effect as setting the min. cut value to 0
 241     Int_t nCtrb = vtxPrim->GetNContributors();
 242     if (nCtrb<fVtxNCtrbMin)
 243             fBitMap->SetBitNumber(kVtxNCtrb,kFALSE);
 244   }
 245
 246   // check vertex correlation cut
 247   if(fCheckCorrelationSPDVtx){
 248     if(vtxTracks && vtxTracks->GetNContributors() && vtxSPD && vtxSPD->GetNContributors()){
 249       // Both vertices available, check correlation
 250       if(TMath::Abs(vtxTracks->GetZ() - vtxSPD->GetZ()) >= 0.5) fBitMap->SetBitNumber(kCorrelation, kFALSE);
 251     } else {
 252       // No correlation available: set the cut to false
 253       fBitMap->SetBitNumber(kCorrelation, kFALSE);
 254     }
 255   }
 256   // check cut on the SPD vertex resolution
 257   if(fVtxResolution){
 258     if(vtxSPD){
 259       TString vtxTyp = vtxSPD->GetTitle();
 260       Double_t cov[6]={0};
 261       vtxSPD->GetCovarianceMatrix(cov);
 262       Double_t zRes = TMath::Sqrt(cov[5]);
 263       if (vtxTyp.Contains("vertexer:Z") && (zRes>0.25)) fBitMap->SetBitNumber(kResolution, kFALSE);
 264     } else{
 265       fBitMap->SetBitNumber(kResolution, kFALSE);
 266     }
 267   }
 268   // check pA pileup cut
 269   if(fPApileupCut){
 270     if(fAnalysisUtils->IsPileUpEvent(inputEvent)) fBitMap->SetBitNumber(kpApileup,kFALSE);
 271   }
 272 }
 273
 274 //_____________________________________________________________________________
 275 void AliHFEextraEventCuts::FillHistograms(TObject* obj, Bool_t b)
 276 {
 277   //
 278   // fill the QA histograms
 279   //
 280
 281   if(!fIsQAOn) return;
 282   // cast TObject into VEvent
 283   AliVEvent* inputEvent = dynamic_cast<AliESDEvent *>(obj);
 284   if (!inputEvent) return;
 285
 286   // index = 0: fill histograms before cuts
 287   // index = 1: fill histograms after cuts
 288   Int_t index = -1;
 289   index = ((b) ? 1 : 0);
 290
 291   // Obtain vertices
 292   const AliVVertex *vtxTracks = GetPrimaryVertexTracks(inputEvent),
 293                    *vtxSPD = GetPrimaryVertexSPD(inputEvent),
 294                    *vtxPrim(NULL);
 295
 296   //look at vertex parameters:
 297   if(fVtxMixed){
 298     if(vtxTracks && vtxTracks->GetNContributors() > 0) vtxPrim = vtxTracks;
 299     else if(vtxSPD && vtxSPD->GetNContributors()) vtxPrim = vtxSPD;
 300   }
 301   else {
 302     vtxPrim = vtxTracks;
 303   }
 304   if(!vtxPrim)return;
 305   // vertex position and uncertainties
 306   fhQA[kVtxPosZ] [index]->Fill(vtxPrim->GetZ());
 307   fhQA[kVtxNCtrb][index]->Fill(vtxPrim->GetNContributors());
 308   // SPD Vertex Position Correlation
 309   if(vtxTracks && vtxSPD){
 310     fhQA[kCorrelation][index]->Fill(vtxTracks->GetZ()-vtxSPD->GetZ());
 311   }
 312   if(vtxSPD){
 313     Double_t cov[6]={0};
 314     vtxSPD->GetCovarianceMatrix(cov);
 315     fhQA[kResolution][index]->Fill(TMath::Sqrt(cov[5]));
 316   }
 317   double isPileup =0 ;
 318   if(fAnalysisUtils->IsPileUpEvent(inputEvent)) isPileup = 1,
 319   fhQA[kpApileup][index]->Fill(isPileup);
 320 }
 321
 322 //____________________________________________________________________
 323 void AliHFEextraEventCuts::SetHistogramBins(Int_t index, Int_t nbins, Double_t *bins)
 324 {
 325   //
 326   //setting x-axis bin limits of QA histogram fhQA[index]
 327   //
 328
 329   for(Int_t i=0;i<kNStepQA;i++){
 330     if(!fhQA[index][i]){AliWarning("non-existing histogram!");
 331     return;
 332     }
 333     fhQA[index][i]->GetXaxis()->Set(nbins,bins);
 334   }
 335 }
 336 //____________________________________________________________________
 337 void AliHFEextraEventCuts::SetHistogramBins(Int_t index, Int_t nbins, Double_t xmin, Double_t xmax)
 338 {
 339   //
 340   //setting x-axis bins and range of QA histogram fhQA[index]
 341   //
 342
 343   for(Int_t i=0;i<kNStepQA;i++){
 344     if(!fhQA[index][i]){AliWarning("non-existing histogram!");
 345     return;
 346     }
 347     fhQA[index][i]->GetXaxis()->Set(nbins,xmin,xmax);
 348   }
 349 }
 350
 351 //_____________________________________________________________________________
 352  void AliHFEextraEventCuts::DefineHistograms() {
 353   //
 354   // histograms for cut variables
 355   //
 356   Int_t color = 2;
 357
 358   if(!fIsQAOn) {
 359     AliInfo(Form("No QA histos requested, Please first set the QA flag on!"));
 360     return;
 361   }
 362
 363   // book QA histograms
 364   Char_t str[5];
 365   for (Int_t i=0; i<kNStepQA; i++) {
 366     if (i==0) snprintf(str,5," ");
 367     else snprintf(str,5,"_cut");
 368
 369     fhQA[kVtxPosZ][i]   = new  TH1F(Form("%s_Vtx_Pos_Z%s",GetName(),str),               "",200,-50.,50.);
 370     fhQA[kVtxNCtrb][i]  = new  TH1F(Form("%s_Vtx_N_Ctrb%s",GetName(),str),              "",1000,0.,1000.);
 371     fhQA[kCorrelation][i] = new TH1F(Form("%s_SPDCorrelation_%s",GetName(),str), "",200, -10., 10);
 372     fhQA[kResolution][i] = new TH1F(Form("%s_SPDResolution_%s",GetName(), str), "", 100, 0., 1.);
 373     fhQA[kpApileup][i] = new TH1F(Form("%s_pApileup_%s",GetName(),str), "", 2, 0., 2.);
 374
 375     fhQA[kVtxPosZ][i]   ->SetXTitle("Vertex Position Z (cm)");
 376     fhQA[kVtxNCtrb][i]  ->SetXTitle("Number of contributors");
 377   }
 378
 379   for(Int_t i=0; i<kNCuts; i++) fhQA[i][1]->SetLineColor(color);
 380
 381 }
 382
 383 //_____________________________________________________________________________
 384 void AliHFEextraEventCuts::AddQAHistograms(TList *qaList) {
 385   //
 386   // saves the histograms in a TList
 387   //
 388
 389   DefineHistograms();
 390
 391   for (Int_t j=0; j<kNStepQA; j++) {
 392     for(Int_t i=0; i<kNCuts; i++)
 393         qaList->Add(fhQA[i][j]);
 394   }
 395 }
 396
 397 //_____________________________________________________________________________
 398 const AliVVertex *AliHFEextraEventCuts::GetPrimaryVertexSPD(const AliVEvent * const inputEvent){
 399   //
 400   // Get SPD vertex from event
 401   //
 402   const AliVVertex *spdvtx(NULL);
 403   const AliESDEvent *esd(NULL);
 404   const AliAODEvent *aod(NULL);
 405   if((esd = dynamic_cast<const AliESDEvent *>(inputEvent))){
 406     spdvtx = esd->GetPrimaryVertexSPD();
 407   } else if((aod = dynamic_cast<const AliAODEvent *>(inputEvent))){
 408     spdvtx = aod->GetPrimaryVertexSPD();
 409   }
 410   return spdvtx;
 411 }
 412
 413 //_____________________________________________________________________________
 414 const AliVVertex *AliHFEextraEventCuts::GetPrimaryVertexTracks(const AliVEvent *const inputEvent){
 415   //
 416   // Get Primary Vertex from tracks
 417   //
 418   const AliVVertex *trkvtx(NULL);
 419   const AliESDEvent *esd(NULL);
 420   const AliAODEvent *aod(NULL);
 421   if((esd = dynamic_cast<const AliESDEvent *>(inputEvent))){
 422     trkvtx = esd->GetPrimaryVertexTracks();
 423   } else if((aod = dynamic_cast<const AliAODEvent *>(inputEvent))){
 424     const AliVVertex *vtxTmp = aod->GetPrimaryVertex();
 425     // check whether the primary vertex is the vertex from tracks
 426     TString vtxTtl = vtxTmp->GetTitle();
 427     if(vtxTtl.Contains("VertexerTracks")){
 428       trkvtx = vtxTmp;
 429     }
 430   }
 431   return trkvtx;
 432 }
 433