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), fHistoCuts(0),fHistoVtxBefSel(0),fHistoVtxAftSel(0),fHistoEtaBefSel(0),fHistoEtaAftSel(0),fHistoNChAftSel(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",3000,-0.5,2999.5);
  53   fCentralityCutMin = 0.0;      // default value of centrality cut minimum, 0 ~ no cut
  54   fCentralityCutMax = 10000.0;  // default value of centrality cut maximum,  ~ no cut
  55
  56 }
  57
  58 //______________________________________________________
  59 Bool_t AliSpectraAODEventCuts::IsSelected(AliAODEvent * aod,AliSpectraAODTrackCuts     *trackcuts)
  60 {
  61   // Returns true if Event Cuts are selected and applied
  62   fAOD = aod;
  63   fTrackCuts = trackcuts;
  64   fHistoCuts->Fill(kProcessedEvents);
  65   //loop on tracks, before event selection, filling QA histos
  66   AliAODVertex * vertex = fAOD->GetPrimaryVertex();//FIXME vertex is recreated
  67   if(vertex)fHistoVtxBefSel->Fill(vertex->GetZ());
  68   fIsSelected = (CheckVtxRange() && CheckCentralityCut());
  69   if(fIsSelected){
  70     fHistoCuts->Fill(kAcceptedEvents);
  71     if(vertex)fHistoVtxAftSel->Fill(vertex->GetZ());
  72   }
  73   Int_t Nch=0;
  74   for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) {
  75     AliAODTrack* track = fAOD->GetTrack(iTracks);
  76     if (!fTrackCuts->IsSelected(track)) continue;
  77     fHistoEtaBefSel->Fill(track->Eta());
  78     if(fIsSelected){
  79       fHistoEtaAftSel->Fill(track->Eta());
  80       Nch++;
  81     }
  82   }
  83   if(fIsSelected)fHistoNChAftSel->Fill(Nch);
  84   return fIsSelected;
  85 }
  86
  87 //______________________________________________________
  88 Bool_t AliSpectraAODEventCuts::CheckVtxRange()
  89 {
  90   // reject events outside of range
  91   AliAODVertex * vertex = fAOD->GetPrimaryVertex();
  92   if (!vertex)
  93     {
  94       fHistoCuts->Fill(kVtxNoEvent);
  95       return kFALSE;
  96     }
  97   if (TMath::Abs(vertex->GetZ()) < 10)
  98     {
  99       return kTRUE;
 100     }
 101   fHistoCuts->Fill(kVtxRange);
 102   return kFALSE;
 103 }
 104
 105 //______________________________________________________
 106 Bool_t AliSpectraAODEventCuts::CheckCentralityCut()
 107 {
 108   // Check centrality cut
 109   Double_t cent=0;
 110   if(fIsMC)cent=fAOD->GetCentrality()->GetCentralityPercentile("V0M");
 111   else cent=ApplyCentralityPatchAOD049();
 112
 113   if ( (cent <= fCentralityCutMax)  &&  (cent >= fCentralityCutMin) )  return kTRUE;
 114   fHistoCuts->Fill(kVtxCentral);
 115   return kFALSE;
 116 }
 117
 118 //______________________________________________________
 119 void AliSpectraAODEventCuts::PrintCuts()
 120 {
 121   // print info about event cuts
 122   cout << "Event Stats" << endl;
 123   cout << " > Number of accepted events: " << fHistoCuts->GetBinContent(kAcceptedEvents + 1) << endl;
 124   cout << " > Number of processed events: " << fHistoCuts->GetBinContent(kProcessedEvents + 1) << endl;
 125   cout << " > Vertex out of range: " << fHistoCuts->GetBinContent(kVtxRange + 1) << endl;
 126   cout << " > Events cut by centrality: " << fHistoCuts->GetBinContent(kVtxCentral + 1) << endl;
 127   cout << " > Events without vertex: " << fHistoCuts->GetBinContent(kVtxNoEvent + 1) << endl;
 128 }
 129 //______________________________________________________
 130
 131 Double_t AliSpectraAODEventCuts::ApplyCentralityPatchAOD049()
 132 {
 133    //
 134    //Apply centrality patch for AOD049 outliers
 135    //
 136   // if (fCentralityType!="V0M") {
 137   //   AliWarning("Requested patch forAOD049 for wrong value (not centrality from V0).");
 138   //   return -999.0;
 139   // }
 140   AliCentrality *centrality = fAOD->GetCentrality();
 141    if (!centrality) {
 142      AliWarning("Cannot get centrality from AOD event.");
 143      return -999.0;
 144    }
 145
 146    Float_t cent = (Float_t)(centrality->GetCentralityPercentile("V0M"));
 147    /*
 148      Bool_t isSelRun = kFALSE;
 149      Int_t selRun[5] = {138364, 138826, 138828, 138836, 138871};
 150      if(cent<0){
 151      Int_t quality = centrality->GetQuality();
 152      if(quality<=1){
 153      cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
 154      } else {
 155      Int_t runnum=aodEvent->GetRunNumber();
 156      for(Int_t ir=0;ir<5;ir++){
 157      if(runnum==selRun[ir]){
 158      isSelRun=kTRUE;
 159      break;
 160      }
 161      }
 162      if((quality==8||quality==9)&&isSelRun) cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
 163      }
 164      }
 165    */
 166    if(cent>=0.0) {
 167      Float_t v0 = 0.0;
 168      AliAODEvent *aodEvent = (AliAODEvent *)fAOD;
 169      AliAODVZERO *aodV0 = (AliAODVZERO *) aodEvent->GetVZEROData();
 170      v0+=aodV0->GetMTotV0A();
 171      v0+=aodV0->GetMTotV0C();
 172      if ( (cent==0) && (v0<19500) ) {
 173        AliDebug(3, Form("Filtering issue in centrality -> cent = %5.2f",cent));
 174        return -999.0;
 175      }
 176      Float_t tkl = (Float_t)(aodEvent->GetTracklets()->GetNumberOfTracklets());
 177      Float_t val = 1.30552 +  0.147931 * v0;
 178
 179      Float_t tklSigma[101] = {176.644, 156.401, 153.789, 153.015, 142.476, 137.951, 136.127, 129.852, 127.436, 124.86,
 180                               120.788, 115.611, 113.172, 110.496, 109.127, 104.421, 102.479, 99.9766, 97.5152, 94.0654,
 181                               92.4602, 89.3364, 87.1342, 83.3497, 82.6216, 81.1084, 78.0793, 76.1234, 72.9434, 72.1334,
 182                               68.0056, 68.2755, 66.0376, 62.9666, 62.4274, 59.65, 58.3776, 56.6361, 54.5184, 53.4224,
 183                               51.932, 50.8922, 48.2848, 47.912, 46.5717, 43.4114, 43.2083, 41.3065, 40.1863, 38.5255,
 184                               37.2851, 37.5396, 34.4949, 33.8366, 31.8043, 31.7412, 30.8392, 30.0274, 28.8793, 27.6398,
 185                               26.6488, 25.0183, 25.1489, 24.4185, 22.9107, 21.2002, 21.6977, 20.1242, 20.4963, 19.0235,
 186                               19.298, 17.4103, 16.868, 15.2939, 15.2939, 16.0295, 14.186, 14.186, 15.2173, 12.9504, 12.9504,
 187                               12.9504, 15.264, 12.3674, 12.3674, 12.3674, 12.3674, 12.3674, 18.3811, 13.7544, 13.7544,
 188                               13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544
 189      };
 190
 191      if ( TMath::Abs(tkl-val) > 6.*tklSigma[(Int_t)cent] )  {
 192        AliDebug(3, Form("Outlier event in centrality -> cent = %5.2f",cent));
 193        return -999.0;
 194      }
 195    } else {
 196      //force it to be -999. whatever the negative value was
 197      cent = -999.;
 198    }
 199    return cent;
 200 }
 201
 202 //______________________________________________________
 203
 204
 205 Long64_t AliSpectraAODEventCuts::Merge(TCollection* list)
 206 {
 207   // Merge a list of AliSpectraAODEventCuts objects with this.
 208   // Returns the number of merged objects (including this).
 209
 210   //  AliInfo("Merging");
 211
 212   if (!list)
 213     return 0;
 214
 215   if (list->IsEmpty())
 216     return 1;
 217
 218   TIterator* iter = list->MakeIterator();
 219   TObject* obj;
 220
 221   // collections of all histograms
 222   TList collections;//FIXME we should use only 1 collection
 223   TList collections_histoVtxBefSel;
 224   TList collections_histoVtxAftSel;
 225   TList collections_histoEtaBefSel;
 226   TList collections_histoEtaAftSel;
 227   TList collections_histoNChAftSel;
 228
 229   Int_t count = 0;
 230
 231   while ((obj = iter->Next())) {
 232     AliSpectraAODEventCuts* entry = dynamic_cast<AliSpectraAODEventCuts*> (obj);
 233     if (entry == 0)
 234       continue;
 235
 236     TH1I * histo = entry->GetHistoCuts();
 237     collections.Add(histo);
 238     TH1F * histo_histoVtxBefSel = entry->GetHistoVtxBefSel();
 239     collections_histoVtxBefSel.Add(histo_histoVtxBefSel);
 240     TH1F * histo_histoVtxAftSel = entry->GetHistoVtxAftSel();
 241     collections_histoVtxAftSel.Add(histo_histoVtxAftSel);
 242     TH1F * histo_histoEtaBefSel = entry->GetHistoEtaBefSel();
 243     collections_histoEtaBefSel.Add(histo_histoEtaBefSel);
 244     TH1F * histo_histoEtaAftSel = entry->GetHistoEtaAftSel();
 245     collections_histoEtaAftSel.Add(histo_histoEtaAftSel);
 246     TH1F * histo_histoNChAftSel = entry->GetHistoNChAftSel();
 247     collections_histoNChAftSel.Add(histo_histoNChAftSel);
 248     count++;
 249   }
 250
 251   fHistoCuts->Merge(&collections);
 252   fHistoVtxBefSel->Merge(&collections_histoVtxBefSel);
 253   fHistoVtxAftSel->Merge(&collections_histoVtxAftSel);
 254   fHistoEtaBefSel->Merge(&collections_histoEtaBefSel);
 255   fHistoEtaAftSel->Merge(&collections_histoEtaAftSel);
 256   fHistoNChAftSel->Merge(&collections_histoNChAftSel);
 257
 258   delete iter;
 259
 260   return count+1;
 261 }
 262
 263 /// FIXME: Q vector
 264 // //Selection on QVector, before ANY other selection on the event
 265 // //Spectra MUST be normalized wrt events AFTER the selection on Qvector
 266 // // Can we include this in fEventCuts
 267 // Double_t Qx2EtaPos = 0, Qy2EtaPos = 0;
 268 // Double_t Qx2EtaNeg = 0, Qy2EtaNeg = 0;
 269 // Int_t multPos = 0;
 270 // Int_t multNeg = 0;
 271 // for(Int_t iT = 0; iT < fAOD->GetNumberOfTracks(); iT++) {
 272 //   AliAODTrack* aodTrack = fAOD->GetTrack(iT);
 273 //   if (!fTrackCuts->IsSelected(aodTrack)) continue;
 274 //   if (aodTrack->Eta() >= 0){
 275 //     multPos++;
 276 //     Qx2EtaPos += TMath::Cos(2*aodTrack->Phi());
 277 //     Qy2EtaPos += TMath::Sin(2*aodTrack->Phi());
 278 //   } else {
 279 //     multNeg++;
 280 //     Qx2EtaNeg += TMath::Cos(2*aodTrack->Phi());
 281 //     Qy2EtaNeg += TMath::Sin(2*aodTrack->Phi());
 282 //   }
 283 // }
 284 // Double_t qPos=-999;
 285 // if(multPos!=0)qPos= TMath::Sqrt((Qx2EtaPos*Qx2EtaPos + Qy2EtaPos*Qy2EtaPos)/multPos);
 286 // Double_t qNeg=-999;
 287 // if(multNeg!=0)qNeg= TMath::Sqrt((Qx2EtaNeg*Qx2EtaNeg + Qy2EtaNeg*Qy2EtaNeg)/multNeg);
 288
 289 // if((qPos>fTrackCuts->GetQvecMin() && qPos<fTrackCuts->GetQvecMax()) || (qNeg>fTrackCuts->GetQvecMin() && qNeg<fTrackCuts->GetQvecMax())){
 290
 291 //fill q distributions vs centrality, after all event selection
 292 // fHistMan->GetqVecHistogram(kHistqVecPos)->Fill(qPos,fAOD->GetCentrality()->GetCentralityPercentile("V0M"));  // qVector distribution
 293 // fHistMan->GetqVecHistogram(kHistqVecNeg)->Fill(qNeg,fAOD->GetCentrality()->GetCentralityPercentile("V0M"));  // qVector distribution