PWGLF/SPECTRA/PiKaPr/TestAOD/AliSpectraAODEventCuts.cxx

   1
   2 /**************************************************************************
   3  * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
   4  *                                                                        *
   5  * Author: The ALICE Off-line Project.                                    *
   6  * Contributors are mentioned in the code where appropriate.              *
   7  *                                                                        *
   8  * Permission to use, copy, modify and distribute this software and its   *
   9  * documentation strictly for non-commercial purposes is hereby granted   *
  10  * without fee, provided that the above copyright notice appears in all   *
  11  * copies and that both the copyright notice and this permission notice   *
  12  * appear in the supporting documentation. The authors make no claims     *
  13  * about the suitability of this software for any purpose. It is          *
  14  * provided "as is" without express or implied warranty.                  *
  15  **************************************************************************/
  16
  17 //-----------------------------------------------------------------
  18 //         AliSpectraAODEventCuts class
  19 //-----------------------------------------------------------------
  20
  21 #include "TChain.h"
  22 #include "TTree.h"
  23 #include "TLegend.h"
  24 #include "TH1F.h"
  25 #include "TH2F.h"
  26 #include "TCanvas.h"
  27 #include "AliAnalysisTask.h"
  28 #include "AliAnalysisManager.h"
  29 #include "AliAODTrack.h"
  30 #include "AliAODMCParticle.h"
  31 #include "AliAODEvent.h"
  32 #include "AliAODInputHandler.h"
  33 #include "AliAnalysisTaskESDfilter.h"
  34 #include "AliAnalysisDataContainer.h"
  35 #include "AliSpectraAODEventCuts.h"
  36 #include "AliSpectraAODTrackCuts.h"
  37 #include "AliSpectraAODHistoManager.h"
  38 #include <iostream>
  39
  40 using namespace std;
  41
  42 ClassImp(AliSpectraAODEventCuts)
  43
  44 AliSpectraAODEventCuts::AliSpectraAODEventCuts(const char *name) : TNamed(name, "AOD Event Cuts"), fAOD(0),fIsMC(0), fIsSelected(0), fCentralityCutMin(0), fCentralityCutMax(0), fQVectorPosCutMin(0), fQVectorPosCutMax(0), fQVectorNegCutMin(0), fQVectorNegCutMax(0), fHistoCuts(0),fHistoVtxBefSel(0),fHistoVtxAftSel(0),fHistoEtaBefSel(0),fHistoEtaAftSel(0),fHistoNChAftSel(0),fHistoQVectorPos(0),fHistoQVectorNeg(0)
  45 {
  46   // Constructor
  47   fHistoCuts = new TH1I("fEventCuts", "Event Cuts", kNVtxCuts, -0.5, kNVtxCuts - 0.5);
  48   fHistoVtxBefSel = new TH1F("fHistoVtxBefSel", "Vtx distr before event selection",500,-15,15);
  49   fHistoVtxAftSel = new TH1F("fHistoVtxAftSel", "Vtx distr after event selection",500,-15,15);
  50   fHistoEtaBefSel = new TH1F("fHistoEtaBefSel", "Eta distr before event selection",500,-2,2);
  51   fHistoEtaAftSel = new TH1F("fHistoEtaAftSel", "Eta distr after event selection",500,-2,2);
  52   fHistoNChAftSel = new TH1F("fHistoNChAftSel", "NCh distr after event selection",2000,-0.5,1999.5);
  53   fHistoQVectorPos = new TH1F("fHistoQVectorPos", "QVectorPos distribution",100,0,10);
  54   fHistoQVectorNeg = new TH1F("fHistoQVectorNeg", "QVectorNeg distribution",100,0,10);
  55   fCentralityCutMin = 0.0;      // default value of centrality cut minimum, 0 ~ no cut
  56   fCentralityCutMax = 10000.0;  // default value of centrality cut maximum,  ~ no cut
  57   fQVectorPosCutMin=0.0;
  58   fQVectorPosCutMax=10000.0;
  59   fQVectorNegCutMin=0.0;
  60   fQVectorNegCutMax=10000.0;
  61 }
  62
  63 //______________________________________________________
  64 Bool_t AliSpectraAODEventCuts::IsSelected(AliAODEvent * aod,AliSpectraAODTrackCuts     *trackcuts)
  65 {
  66   // Returns true if Event Cuts are selected and applied
  67   fAOD = aod;
  68   fTrackCuts = trackcuts;
  69   fHistoCuts->Fill(kProcessedEvents);
  70   Bool_t IsPhysSel = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);//FIXME we can add the trigger mask here
  71   if(!IsPhysSel)return IsPhysSel;
  72   fHistoCuts->Fill(kPhysSelEvents);
  73   //loop on tracks, before event selection, filling QA histos
  74   AliAODVertex * vertex = fAOD->GetPrimaryVertex();//FIXME vertex is recreated
  75   if(vertex)fHistoVtxBefSel->Fill(vertex->GetZ());
  76   fIsSelected =kFALSE;
  77   if(CheckVtxRange() && CheckCentralityCut()){ //selection on vertex and Centrality
  78     fIsSelected=CheckQVectorCut(); // QVector is calculated only if the centrality and vertex are correct (performance)
  79   }
  80   if(fIsSelected){
  81     fHistoCuts->Fill(kAcceptedEvents);
  82     if(vertex)fHistoVtxAftSel->Fill(vertex->GetZ());
  83   }
  84   Int_t Nch=0;
  85   for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) {
  86     AliAODTrack* track = fAOD->GetTrack(iTracks);
  87     if (!fTrackCuts->IsSelected(track,kFALSE)) continue;
  88     fHistoEtaBefSel->Fill(track->Eta());
  89     if(fIsSelected){
  90       fHistoEtaAftSel->Fill(track->Eta());
  91       Nch++;
  92     }
  93   }
  94   if(fIsSelected)fHistoNChAftSel->Fill(Nch);
  95   return fIsSelected;
  96 }
  97
  98 //______________________________________________________
  99 Bool_t AliSpectraAODEventCuts::CheckVtxRange()
 100 {
 101   // reject events outside of range
 102   AliAODVertex * vertex = fAOD->GetPrimaryVertex();
 103   if (!vertex)
 104     {
 105       fHistoCuts->Fill(kVtxNoEvent);
 106       return kFALSE;
 107     }
 108   if (TMath::Abs(vertex->GetZ()) < 10)
 109     {
 110       return kTRUE;
 111     }
 112   fHistoCuts->Fill(kVtxRange);
 113   return kFALSE;
 114 }
 115
 116 //______________________________________________________
 117 Bool_t AliSpectraAODEventCuts::CheckCentralityCut()
 118 {
 119   // Check centrality cut
 120   Double_t cent=0;
 121   if(!fUseCentPatchAOD049)cent=fAOD->GetCentrality()->GetCentralityPercentile("V0M");
 122   else cent=ApplyCentralityPatchAOD049();
 123
 124   if ( (cent <= fCentralityCutMax)  &&  (cent >= fCentralityCutMin) )  return kTRUE;
 125   fHistoCuts->Fill(kVtxCentral);
 126   return kFALSE;
 127 }
 128
 129 //______________________________________________________
 130 Bool_t AliSpectraAODEventCuts::CheckQVectorCut()
 131 {
 132   // Check qvector cut
 133   /// FIXME: Q vector
 134   // //Selection on QVector, before ANY other selection on the event
 135   // //Spectra MUST be normalized wrt events AFTER the selection on Qvector
 136   Double_t Qx2EtaPos = 0, Qy2EtaPos = 0;
 137   Double_t Qx2EtaNeg = 0, Qy2EtaNeg = 0;
 138   Int_t multPos = 0;
 139   Int_t multNeg = 0;
 140   for(Int_t iT = 0; iT < fAOD->GetNumberOfTracks(); iT++) {
 141     AliAODTrack* aodTrack = fAOD->GetTrack(iT);
 142     if (!aodTrack->TestFilterBit(1)) continue; //FilterBit 1 is the standard TPC track
 143     if (aodTrack->Eta() >= 0){
 144       multPos++;
 145       Qx2EtaPos += TMath::Cos(2*aodTrack->Phi());
 146       Qy2EtaPos += TMath::Sin(2*aodTrack->Phi());
 147     } else {
 148       multNeg++;
 149       Qx2EtaNeg += TMath::Cos(2*aodTrack->Phi());
 150       Qy2EtaNeg += TMath::Sin(2*aodTrack->Phi());
 151     }
 152   }
 153   Double_t qPos=-999;
 154   if(multPos!=0)qPos= TMath::Sqrt((Qx2EtaPos*Qx2EtaPos + Qy2EtaPos*Qy2EtaPos)/multPos);
 155   Double_t qNeg=-999;
 156   if(multNeg!=0)qNeg= TMath::Sqrt((Qx2EtaNeg*Qx2EtaNeg + Qy2EtaNeg*Qy2EtaNeg)/multNeg);
 157
 158   //Printf("Cuts on Q vector: %f,%f,%f,%f",fQVectorPosCutMin,fQVectorPosCutMax,fQVectorNegCutMin,fQVectorNegCutMax);
 159   //Printf("Q vectors: %f,%f",qPos,qNeg);
 160   if(qPos<fQVectorPosCutMin || qPos>fQVectorPosCutMax || qNeg<fQVectorNegCutMin || qNeg>fQVectorNegCutMax)return kFALSE;
 161
 162
 163   fHistoCuts->Fill(kQVector);
 164   fHistoQVectorPos->Fill(qPos);
 165   fHistoQVectorNeg->Fill(qNeg);
 166   return kTRUE;
 167
 168 }
 169
 170 //______________________________________________________
 171 void AliSpectraAODEventCuts::PrintCuts()
 172 {
 173   // print info about event cuts
 174   cout << "Event Stats" << endl;
 175   cout << " > Number of accepted events: " << fHistoCuts->GetBinContent(kAcceptedEvents + 1) << endl;
 176   cout << " > Number of processed events: " << fHistoCuts->GetBinContent(kProcessedEvents + 1) << endl;
 177   cout << " > Number of PhysSel events: " << fHistoCuts->GetBinContent(kPhysSelEvents + 1) << endl;
 178   cout << " > Vertex out of range: " << fHistoCuts->GetBinContent(kVtxRange + 1) << endl;
 179   cout << " > Events cut by centrality: " << fHistoCuts->GetBinContent(kVtxCentral + 1) << endl;
 180   cout << " > Events without vertex: " << fHistoCuts->GetBinContent(kVtxNoEvent + 1) << endl;
 181   cout << " > QVector cut: " << fHistoCuts->GetBinContent(kQVector + 1) << endl;
 182 }
 183 //______________________________________________________
 184
 185 Double_t AliSpectraAODEventCuts::ApplyCentralityPatchAOD049()
 186 {
 187    //
 188    //Apply centrality patch for AOD049 outliers
 189    //
 190   // if (fCentralityType!="V0M") {
 191   //   AliWarning("Requested patch forAOD049 for wrong value (not centrality from V0).");
 192   //   return -999.0;
 193   // }
 194   AliCentrality *centrality = fAOD->GetCentrality();
 195    if (!centrality) {
 196      AliWarning("Cannot get centrality from AOD event.");
 197      return -999.0;
 198    }
 199
 200    Float_t cent = (Float_t)(centrality->GetCentralityPercentile("V0M"));
 201    /*
 202      Bool_t isSelRun = kFALSE;
 203      Int_t selRun[5] = {138364, 138826, 138828, 138836, 138871};
 204      if(cent<0){
 205      Int_t quality = centrality->GetQuality();
 206      if(quality<=1){
 207      cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
 208      } else {
 209      Int_t runnum=aodEvent->GetRunNumber();
 210      for(Int_t ir=0;ir<5;ir++){
 211      if(runnum==selRun[ir]){
 212      isSelRun=kTRUE;
 213      break;
 214      }
 215      }
 216      if((quality==8||quality==9)&&isSelRun) cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
 217      }
 218      }
 219    */
 220    if(cent>=0.0) {
 221      Float_t v0 = 0.0;
 222      AliAODEvent *aodEvent = (AliAODEvent *)fAOD;
 223      AliAODVZERO *aodV0 = (AliAODVZERO *) aodEvent->GetVZEROData();
 224      v0+=aodV0->GetMTotV0A();
 225      v0+=aodV0->GetMTotV0C();
 226      if ( (cent==0) && (v0<19500) ) {
 227        AliDebug(3, Form("Filtering issue in centrality -> cent = %5.2f",cent));
 228        return -999.0;
 229      }
 230      Float_t tkl = (Float_t)(aodEvent->GetTracklets()->GetNumberOfTracklets());
 231      Float_t val = 1.30552 +  0.147931 * v0;
 232
 233      Float_t tklSigma[101] = {176.644, 156.401, 153.789, 153.015, 142.476, 137.951, 136.127, 129.852, 127.436, 124.86,
 234                               120.788, 115.611, 113.172, 110.496, 109.127, 104.421, 102.479, 99.9766, 97.5152, 94.0654,
 235                               92.4602, 89.3364, 87.1342, 83.3497, 82.6216, 81.1084, 78.0793, 76.1234, 72.9434, 72.1334,
 236                               68.0056, 68.2755, 66.0376, 62.9666, 62.4274, 59.65, 58.3776, 56.6361, 54.5184, 53.4224,
 237                               51.932, 50.8922, 48.2848, 47.912, 46.5717, 43.4114, 43.2083, 41.3065, 40.1863, 38.5255,
 238                               37.2851, 37.5396, 34.4949, 33.8366, 31.8043, 31.7412, 30.8392, 30.0274, 28.8793, 27.6398,
 239                               26.6488, 25.0183, 25.1489, 24.4185, 22.9107, 21.2002, 21.6977, 20.1242, 20.4963, 19.0235,
 240                               19.298, 17.4103, 16.868, 15.2939, 15.2939, 16.0295, 14.186, 14.186, 15.2173, 12.9504, 12.9504,
 241                               12.9504, 15.264, 12.3674, 12.3674, 12.3674, 12.3674, 12.3674, 18.3811, 13.7544, 13.7544,
 242                               13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544
 243      };
 244
 245      if ( TMath::Abs(tkl-val) > 6.*tklSigma[(Int_t)cent] )  {
 246        AliDebug(3, Form("Outlier event in centrality -> cent = %5.2f",cent));
 247        return -999.0;
 248      }
 249    } else {
 250      //force it to be -999. whatever the negative value was
 251      cent = -999.;
 252    }
 253    return cent;
 254 }
 255
 256 //______________________________________________________
 257
 258
 259 Long64_t AliSpectraAODEventCuts::Merge(TCollection* list)
 260 {
 261   // Merge a list of AliSpectraAODEventCuts objects with this.
 262   // Returns the number of merged objects (including this).
 263
 264   //  AliInfo("Merging");
 265
 266   if (!list)
 267     return 0;
 268
 269   if (list->IsEmpty())
 270     return 1;
 271
 272   TIterator* iter = list->MakeIterator();
 273   TObject* obj;
 274
 275   // collections of all histograms
 276   TList collections;//FIXME we should use only 1 collection
 277   TList collections_histoVtxBefSel;
 278   TList collections_histoVtxAftSel;
 279   TList collections_histoEtaBefSel;
 280   TList collections_histoEtaAftSel;
 281   TList collections_histoNChAftSel;
 282   TList collections_histoQVectorPos;
 283   TList collections_histoQVectorNeg;
 284
 285   Int_t count = 0;
 286
 287   while ((obj = iter->Next())) {
 288     AliSpectraAODEventCuts* entry = dynamic_cast<AliSpectraAODEventCuts*> (obj);
 289     if (entry == 0)
 290       continue;
 291
 292     TH1I * histo = entry->GetHistoCuts();
 293     collections.Add(histo);
 294     TH1F * histo_histoVtxBefSel = entry->GetHistoVtxBefSel();
 295     collections_histoVtxBefSel.Add(histo_histoVtxBefSel);
 296     TH1F * histo_histoVtxAftSel = entry->GetHistoVtxAftSel();
 297     collections_histoVtxAftSel.Add(histo_histoVtxAftSel);
 298     TH1F * histo_histoEtaBefSel = entry->GetHistoEtaBefSel();
 299     collections_histoEtaBefSel.Add(histo_histoEtaBefSel);
 300     TH1F * histo_histoEtaAftSel = entry->GetHistoEtaAftSel();
 301     collections_histoEtaAftSel.Add(histo_histoEtaAftSel);
 302     TH1F * histo_histoNChAftSel = entry->GetHistoNChAftSel();
 303     collections_histoNChAftSel.Add(histo_histoNChAftSel);
 304     TH1F * histo_histoQVectorPos = entry->GetHistoQVectorPos();
 305     collections_histoQVectorPos.Add(histo_histoQVectorPos);
 306     TH1F * histo_histoQVectorNeg = entry->GetHistoQVectorNeg();
 307     collections_histoQVectorNeg.Add(histo_histoQVectorNeg);
 308     count++;
 309   }
 310
 311   fHistoCuts->Merge(&collections);
 312   fHistoVtxBefSel->Merge(&collections_histoVtxBefSel);
 313   fHistoVtxAftSel->Merge(&collections_histoVtxAftSel);
 314   fHistoEtaBefSel->Merge(&collections_histoEtaBefSel);
 315   fHistoEtaAftSel->Merge(&collections_histoEtaAftSel);
 316   fHistoNChAftSel->Merge(&collections_histoNChAftSel);
 317   fHistoQVectorPos->Merge(&collections_histoQVectorPos);
 318   fHistoQVectorNeg->Merge(&collections_histoQVectorNeg);
 319
 320   delete iter;
 321
 322   return count+1;
 323 }
 324