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