PWGCF/Correlations/DPhi/DiHadronPID/AliAODEventCutsDiHadronPID.cxx

   1 /*************************************************************************
   2 * Copyright(c) 1998-2008, 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 //  Event Cut class for the DiHadronPID analysis.
  18 // -----------------------------------------------------------------------
  19 //  Author: Misha Veldhoen (misha.veldhoen@cern.ch)
  20
  21 #include "AliAODEventCutsDiHadronPID.h"
  22
  23 #include <iostream>
  24 using namespace std;
  25
  26 #include "TList.h"
  27 #include "TMath.h"
  28 #include "TH1F.h"
  29 #include "TH2F.h"
  30 #include "TNamed.h"
  31 #include "TIterator.h"
  32 #include "AliAODHeader.h"
  33 #include "AliAODVertex.h"
  34 #include "AliCentrality.h"
  35 #include "AliAnalysisManager.h"
  36 #include "AliInputEventHandler.h"
  37
  38 ClassImp(AliAODEventCutsDiHadronPID);
  39
  40 // -----------------------------------------------------------------------
  41 AliAODEventCutsDiHadronPID::AliAODEventCutsDiHadronPID():
  42         TNamed(),
  43         fIsPbPb(kTRUE),
  44         fIsMC(kFALSE),
  45         fTrigger(AliVEvent::kMB),
  46         fMinCentrality(5.),
  47         fMaxCentrality(0.),
  48         fCentralityEstimator("V0M"),
  49         fMaxVertexZ(10.),
  50         fMinRefMult(0),
  51         fTestTrigger(kFALSE),
  52         fTestCentrality(kFALSE),
  53         fTestContributorsOrSPDVertex(kFALSE),
  54         fTestVertexZ(kFALSE),
  55         fTestMinRefMult(kFALSE),
  56         fSelectedEventQAHistos(0x0),
  57         fAllEventQAHistos(0x0),
  58         fHistTrigger(0x0),
  59         fHistRefMultiplicity(0x0),
  60         fHistCentrality(0x0),
  61         fHistCentralityQuality(0x0),
  62         fHistVertexZ(0x0),
  63         fDebug(0)
  64
  65 {
  66
  67         //
  68         // Default Constructor
  69         //
  70
  71         cout<<"AliAODEventCutsDiHadronPID Default Constructor Called."<<endl;
  72         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
  73
  74 }
  75
  76 // -----------------------------------------------------------------------
  77 AliAODEventCutsDiHadronPID::AliAODEventCutsDiHadronPID(const char* name):
  78         TNamed(name,"AOD Event Cuts"),
  79         fIsPbPb(kTRUE),
  80         fIsMC(kFALSE),
  81         fTrigger(AliVEvent::kMB),
  82         fMinCentrality(5.),
  83         fMaxCentrality(0.),
  84         fCentralityEstimator("V0M"),
  85         fMaxVertexZ(10.),
  86         fMinRefMult(0),
  87         fTestTrigger(kFALSE),
  88         fTestCentrality(kFALSE),
  89         fTestContributorsOrSPDVertex(kFALSE),
  90         fTestVertexZ(kFALSE),
  91         fTestMinRefMult(kFALSE),
  92         fSelectedEventQAHistos(0x0),
  93         fAllEventQAHistos(0x0),
  94         fHistTrigger(0x0),
  95         fHistRefMultiplicity(0x0),
  96         fHistCentrality(0x0),
  97         fHistCentralityQuality(0x0),
  98         fHistVertexZ(0x0),
  99         fDebug(0)
 100
 101 {
 102
 103         //
 104         // Named Constructor
 105         //
 106
 107         cout<<"AliAODEventCutsDiHadronPID Named Constructor Called."<<endl;
 108         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 109
 110 }
 111
 112 // -----------------------------------------------------------------------
 113 AliAODEventCutsDiHadronPID::~AliAODEventCutsDiHadronPID() {
 114
 115         //
 116         // Destructor
 117         //
 118
 119         cout<<"AliAODEventCutsDiHadronPID Destructor Called."<<endl;
 120         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 121
 122         if (fSelectedEventQAHistos) delete fSelectedEventQAHistos;
 123         fSelectedEventQAHistos = 0x0;
 124         if (fAllEventQAHistos) delete fAllEventQAHistos;
 125         fAllEventQAHistos = 0x0;
 126
 127 }
 128
 129 // -----------------------------------------------------------------------
 130 Long64_t AliAODEventCutsDiHadronPID::Merge(TCollection* list) {
 131
 132         //
 133         // Merger.
 134         //
 135
 136         cout<<"AliAODEventCutsDiHadronPID Merger Called."<<endl;
 137         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 138
 139         if (!list) return 0;
 140         if (list->IsEmpty()) return 1;
 141
 142         if (!fSelectedEventQAHistos||!fAllEventQAHistos) {
 143                 cout<<"AliAODEventCutsDiHadronPID::Merge() - Warning, current object's histograms are missing... Generating."<<endl;
 144                 CreateHistos();
 145         }
 146
 147         TIterator* iter = list->MakeIterator();
 148         TObject* obj;
 149
 150         // List of collections
 151         TList collection_fSelectedEventQAHistos;
 152         TList collection_fAllEventQAHistos;
 153
 154         Int_t count = 0;
 155
 156         while ((obj = iter->Next())) {
 157         AliAODEventCutsDiHadronPID* entry = dynamic_cast<AliAODEventCutsDiHadronPID*> (obj);
 158         if (entry == 0) continue;
 159
 160         // Check if the object to be merged really has the same name! (FIXME!)
 161
 162         // Getting the lists from obj.
 163         TList* list_fSelectedEventQAHistos = entry->GetListOfSelectedEventQAHistos();
 164         TList* list_fAllEventQAHistos = entry->GetListOfAllEventQAHistos();
 165
 166         // Adding the retrieved lists to the collection.
 167         if (list_fSelectedEventQAHistos) collection_fSelectedEventQAHistos.Add(list_fSelectedEventQAHistos);
 168         if (list_fAllEventQAHistos) collection_fAllEventQAHistos.Add(list_fAllEventQAHistos);
 169
 170         count++;
 171     }
 172
 173     // Merging. Note that we require the original list to exist.
 174     //  * Assume that if the collection happens to be empty, then nothing will happen.
 175     //  * All other variables are taken from the original object.
 176     if (fSelectedEventQAHistos) fSelectedEventQAHistos->Merge(&collection_fSelectedEventQAHistos);
 177     if (fAllEventQAHistos) fAllEventQAHistos->Merge(&collection_fAllEventQAHistos);
 178
 179     delete iter;
 180
 181         return count+1;
 182
 183 }
 184
 185 // -----------------------------------------------------------------------
 186 void AliAODEventCutsDiHadronPID::CreateHistos() {
 187
 188         cout<<"AliAODEventCutsDiHadronPID - Creating histograms"<<endl;
 189         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 190
 191         // Create list of Event related QA histograms (selected events).
 192         fSelectedEventQAHistos = new TList();
 193         fSelectedEventQAHistos->SetName("SelectedEventQAHistos");
 194         fSelectedEventQAHistos->SetOwner(kTRUE);
 195
 196         // The same, but for all events.
 197         fAllEventQAHistos = new TList();
 198         fAllEventQAHistos->SetName("AllEventQAHistos");
 199         fAllEventQAHistos->SetOwner(kTRUE);
 200
 201         // Creating arrays of pointers to the QA histos.
 202         fHistTrigger = new TH1F*[2];
 203         fHistRefMultiplicity = new TH1F*[2];
 204         fHistCentrality = new TH1F*[2];
 205         fHistCentralityQuality = new TH1F*[2];
 206         fHistVertexZ = new TH1F*[2];
 207
 208         const char* HistType[2] = {"Selected","All"};
 209
 210         for (Int_t iHistType = 0; iHistType < 2; iHistType++) {
 211
 212                 // Trigger Histogram.
 213                 fHistTrigger[iHistType] = new TH1F(Form("fHistTrigger%s",HistType[iHistType]),"Trigger;;Count",5,-0.5,4.5);
 214                 fHistTrigger[iHistType]->GetXaxis()->SetBinLabel(1,"kMB");
 215                 fHistTrigger[iHistType]->GetXaxis()->SetBinLabel(2,"kCentral");         // Trigger only defined for period LHC11h.
 216                 fHistTrigger[iHistType]->GetXaxis()->SetBinLabel(3,"kSemiCentral");     // Trigger only defined for period LHC11h.
 217                 fHistTrigger[iHistType]->GetXaxis()->SetBinLabel(4,"kINT7");
 218                 fHistTrigger[iHistType]->GetXaxis()->SetBinLabel(5,"Other");
 219                 if (iHistType == 0) fSelectedEventQAHistos->Add(fHistTrigger[iHistType]);
 220                 else fAllEventQAHistos->Add(fHistTrigger[iHistType]);
 221
 222                 // Ref Multiplicity Histogram.
 223                 if (fIsPbPb) {
 224                         fHistRefMultiplicity[iHistType] = new TH1F(Form("fHistRefMultiplicity%s",HistType[iHistType]),"Reference Multiplicity;N_{tracks};Count",100,0.,10000.);
 225                 } else {
 226                         fHistRefMultiplicity[iHistType] = new TH1F(Form("fHistRefMultiplicity%s",HistType[iHistType]),"Reference Multiplicity;N_{tracks};Count",100,0.,100.);
 227                 }
 228                 if (iHistType == 0) fSelectedEventQAHistos->Add(fHistRefMultiplicity[iHistType]);
 229                 else fAllEventQAHistos->Add(fHistRefMultiplicity[iHistType]);
 230
 231                 // Centrality Histogram.
 232                 fHistCentrality[iHistType] = new TH1F(Form("fHistCentrality%s",HistType[iHistType]),"Centrality;Centrality;Count",20,0,100);
 233                 if (iHistType == 0) fSelectedEventQAHistos->Add(fHistCentrality[iHistType]);
 234                 else fAllEventQAHistos->Add(fHistCentrality[iHistType]);
 235
 236                 // Centrality Quality.
 237                 fHistCentralityQuality[iHistType] = new TH1F(Form("fHistCentralityQuality%s",HistType[iHistType]),"Centrality Quality;Quality;Count",2,-0.5,1.5);
 238                 fHistCentralityQuality[iHistType]->GetXaxis()->SetBinLabel(1,"0");
 239                 fHistCentralityQuality[iHistType]->GetXaxis()->SetBinLabel(2,"Other");
 240                 if (iHistType == 0) fSelectedEventQAHistos->Add(fHistCentralityQuality[iHistType]);
 241                 else fAllEventQAHistos->Add(fHistCentralityQuality[iHistType]);
 242
 243                 // VertexZ Histogram.
 244                 fHistVertexZ[iHistType] = new TH1F(Form("fHistVertexZ%s",HistType[iHistType]),"VertexZ;z (cm);Count",60,-15.,15.);
 245                 if (iHistType == 0) fSelectedEventQAHistos->Add(fHistVertexZ[iHistType]);
 246                 else fAllEventQAHistos->Add(fHistVertexZ[iHistType]);
 247
 248         }
 249
 250 }
 251
 252 // -----------------------------------------------------------------------
 253 Bool_t AliAODEventCutsDiHadronPID::IsSelected(AliAODEvent* event) {
 254
 255         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 256
 257         if (!event) return kFALSE;
 258
 259         if (!fAllEventQAHistos||!fSelectedEventQAHistos) {cout<<"AliAODEventCutsDiHadronPID - Histograms were not created, you should have called CreateHistos()..."<<endl;}
 260
 261         // Input the event handler.
 262         AliInputEventHandler* InputHandler = (AliInputEventHandler*)((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
 263         if (!InputHandler) return kFALSE;
 264
 265         Bool_t select = kTRUE;
 266
 267         // Test Trigger.
 268         UInt_t trigger = InputHandler->IsEventSelected();
 269         Int_t triggerclass = -1;                // 0 = kMB, 1 = kCentral, 2 = kSemiCentral, 3 = Other.
 270     if (fTestTrigger) {
 271
 272         // Find Trigger class.
 273         if (trigger & AliVEvent::kMB) {triggerclass = 0;}
 274         else if (trigger & AliVEvent::kCentral) {triggerclass = 1;}
 275         else if (trigger & AliVEvent::kSemiCentral) {triggerclass = 2;}
 276         else if (trigger & AliVEvent::kINT7) {triggerclass = 3;}
 277         else {triggerclass = 4;}
 278
 279                 if (!(trigger & fTrigger)) {select = kFALSE;}   // Event not selected if not matched with the any of the desired triggers.
 280         }
 281
 282         AliCentrality* CurrentCentrality = 0x0;
 283         Int_t CurrentCentralityQuality = -999;
 284         Float_t percentile = -999.;
 285
 286         if (fIsPbPb) {
 287
 288                 // Get the centrality object.
 289             CurrentCentrality = event->GetCentrality();
 290             //cout<<"Centrality object: "<<CurrentCentrality<<endl;
 291             if (!CurrentCentrality) {select = kFALSE; return select;}
 292
 293             // Check the quality of the centrality estimation.
 294             // If 0 then quality is OK, c.f. TOF/PbPb276/macros/TOFmatchEff.C
 295             CurrentCentralityQuality = CurrentCentrality->GetQuality();
 296             //cout<<"Centrality: "<<CurrentCentrality->GetCentralityPercentile(fCentralityEstimator.Data())<<" Quality: "<<CurrentCentrality->GetQuality()<<endl;
 297             if (CurrentCentralityQuality) select = kFALSE;
 298
 299                 // Test Centrality.
 300             percentile = CurrentCentrality->GetCentralityPercentile(fCentralityEstimator.Data());
 301                 if (fTestCentrality) {
 302                 if ((percentile < fMaxCentrality)||(percentile > fMinCentrality)) select = kFALSE;
 303                 }
 304
 305         }
 306
 307         // Get the primary vertex.
 308         AliAODVertex* CurrentPrimaryVertex = event->GetPrimaryVertex();
 309     if (!CurrentPrimaryVertex) {select = kFALSE; return select;}
 310
 311         // Get the Vertex_z value (improved by Peter Christiansen).
 312         Double_t vtxz = -999.;
 313         if (fIsPbPb) {vtxz = CurrentPrimaryVertex->GetZ();}
 314         else if (CurrentPrimaryVertex->GetNContributors() > 0) {
 315                 vtxz = CurrentPrimaryVertex->GetZ();
 316         }
 317
 318         // Test Vertex Z.
 319         if (fTestVertexZ) {
 320         if (TMath::Abs(vtxz) > fMaxVertexZ) select = kFALSE;
 321         }
 322
 323         // Test number of contributors.
 324         if (fTestContributorsOrSPDVertex) {
 325                 Int_t nContributors = CurrentPrimaryVertex->GetNContributors();
 326                 if (nContributors < 1) {
 327                         if (CurrentPrimaryVertex->GetType() != AliAODVertex::kMainSPD ) {select = kFALSE;}
 328                 }
 329         }
 330
 331         // Get the event header.
 332         AliAODHeader* CurrentHeader = dynamic_cast<AliAODHeader*>(event->GetHeader());
 333         if(!CurrentHeader) AliFatal("Not a standard AOD");
 334
 335         // Test minimum reference multiplicity.
 336         Int_t CurrentRefMultiplicity = CurrentHeader->GetRefMultiplicity();
 337         if (fTestMinRefMult) {
 338                 if (CurrentRefMultiplicity < fMinRefMult) select = kFALSE;
 339         }
 340
 341         // Fill the histograms for selected events.
 342         if (select) {
 343                 fHistTrigger[0]->Fill(triggerclass);
 344                 fHistRefMultiplicity[0]->Fill(CurrentHeader->GetRefMultiplicity());
 345                 if (fIsPbPb) fHistCentrality[0]->Fill(percentile);
 346                 if (fIsPbPb) fHistCentralityQuality[0]->Fill(CurrentCentralityQuality);
 347                 fHistVertexZ[0]->Fill(vtxz);
 348         }
 349
 350         // Fill the histograms for all events.
 351         fHistTrigger[1]->Fill(triggerclass);
 352         fHistRefMultiplicity[1]->Fill(CurrentHeader->GetRefMultiplicity());
 353         if (fIsPbPb) fHistCentrality[1]->Fill(percentile);
 354         if (fIsPbPb) fHistCentralityQuality[1]->Fill(CurrentCentralityQuality);
 355         fHistVertexZ[1]->Fill(vtxz);
 356
 357         //cout<<"Event Selected: "<<select<<endl;
 358         return select;
 359
 360 }
 361
 362 // -----------------------------------------------------------------------
 363 void AliAODEventCutsDiHadronPID::PrintCuts() {
 364
 365         // NOT IMPLEMENTED.
 366         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 367
 368         return;
 369
 370 }