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 "TH1I.h"
  26 #include "TH2F.h"
  27 #include "TCanvas.h"
  28 #include "TProfile.h"
  29 #include "TSpline.h"
  30 #include "AliAnalysisTask.h"
  31 #include "AliAnalysisManager.h"
  32 #include "AliAODTrack.h"
  33 #include "AliAODMCParticle.h"
  34 #include "AliAODEvent.h"
  35 #include "AliAODInputHandler.h"
  36 #include "AliAnalysisTaskESDfilter.h"
  37 #include "AliAnalysisDataContainer.h"
  38 #include "AliESDVZERO.h"
  39 #include "AliAODVZERO.h"
  40 #include "AliSpectraAODEventCuts.h"
  41 #include "AliSpectraAODTrackCuts.h"
  42 #include <iostream>
  43
  44 using namespace std;
  45
  46 ClassImp(AliSpectraAODEventCuts)
  47
  48 AliSpectraAODEventCuts::AliSpectraAODEventCuts(const char *name) :
  49   TNamed(name, "AOD Event Cuts"),
  50   fAOD(0),
  51   fSelectBit(AliVEvent::kMB),
  52   fCentralityMethod("V0M"),
  53   fTrackBits(1),
  54   fIsMC(0),
  55   fIsLHC10h(1),
  56   fTrackCuts(0),
  57   fIsSelected(0),
  58   fCentralityCutMin(0.),
  59   fCentralityCutMax(999),
  60   fQVectorCutMin(-999.),
  61   fQVectorCutMax(999.),
  62   fVertexCutMin(-10.),
  63   fVertexCutMax(10.),
  64   fMultiplicityCutMin(-999.),
  65   fMultiplicityCutMax(99999.),
  66   fqV0C(-999.),
  67   fqV0A(-999.),
  68   fqV0Cx(-999.),
  69   fqV0Ax(-999.),
  70   fqV0Cy(-999.),
  71   fqV0Ay(-999.),
  72   fPsiV0C(-999.),
  73   fPsiV0A(-999.),
  74   fCent(-999.),
  75   fOutput(0),
  76   fCalib(0),
  77   fRun(-1),
  78   fMultV0(0),
  79   fV0Cpol1(-1),
  80   fV0Cpol2(-1),
  81   fV0Cpol3(-1),
  82   fV0Cpol4(-1),
  83   fV0Apol1(-1),
  84   fV0Apol2(-1),
  85   fV0Apol3(-1),
  86   fV0Apol4(-1),
  87   fQvecIntList(0),
  88   fQvecIntegral(0),
  89   fSplineArrayV0A(0),
  90   fSplineArrayV0C(0),
  91   fQgenIntegral(0),
  92   fSplineArrayV0Agen(0),
  93   fSplineArrayV0Cgen(0),
  94   fQvecMC(0),
  95   fNch(0),
  96   fQvecCalibType(0),
  97   fV0Aeff(0)
  98 {
  99   // Constructor
 100   fOutput=new TList();
 101   fOutput->SetOwner();
 102   fOutput->SetName("fOutput");
 103
 104   fCalib=new TList();
 105   fCalib->SetOwner();
 106   fCalib->SetName("fCalib");
 107
 108   fQvecIntList=new TList();
 109   fQvecIntList->SetOwner();
 110   fQvecIntList->SetName("fQvecIntList");
 111
 112   TH1I *fHistoCuts = new TH1I("fHistoCuts", "Event Cuts", kNVtxCuts, -0.5, kNVtxCuts - 0.5);
 113   TH1F *fHistoVtxBefSel = new TH1F("fHistoVtxBefSel", "Vtx distr before event selection;z (cm)",500,-15,15);
 114   TH1F *fHistoVtxAftSel = new TH1F("fHistoVtxAftSel", "Vtx distr after event selection;z (cm)",500,-15,15);
 115   TH1F *fHistoEtaBefSel = new TH1F("fHistoEtaBefSel", "Eta distr before event selection;eta",500,-2,2);
 116   TH1F *fHistoEtaAftSel = new TH1F("fHistoEtaAftSel", "Eta distr after event selection;eta",500,-2,2);
 117   TH1F *fHistoNChAftSel = new TH1F("fHistoNChAftSel", "NCh distr after event selection;Nch",2000,-0.5,1999.5);
 118   TH2F *fHistoQVector = new TH2F("fHistoQVector", "QVector with VZERO distribution;centrality;Q vector from EP task",20,0,100,100,0,10);
 119   TH2F *fHistoEP = new TH2F("fHistoEP", "EP with VZERO distribution;centrality;Psi_{EP} from EP task",20,0,100,100,-2,2);
 120   TH2F *fPsiACor = new TH2F("fPsiACor", "EP with VZERO A distribution;centrality;Psi_{EP} VZERO-A",20,0,100,100,-2,2);
 121   TH2F *fPsiCCor = new TH2F("fPsiCCor", "EP with VZERO C distribution;centrality;Psi_{EP} VZERO-C",20,0,100,100,-2,2);
 122   TH2F *fQVecACor = new TH2F("fQVecACor", "QVec VZERO A;centrality;Qvector VZERO-A",20,0,100,100,0,10);
 123   TH2F *fQVecCCor = new TH2F("fQVecCCor", "QVec VZERO C;centrality;Qvector VZERO-C",20,0,100,100,0,10);
 124   TH2F *fV0M = new TH2F("fV0M", "V0 Multiplicity, before correction;V0 sector",64,-.5,63.5,500,0,1000);
 125   TH2F *fV0MCor = new TH2F("fV0MCor", "V0 Multiplicity, after correction;V0 sector",64,-.5,63.5,500,0,1000);
 126   TH2F *fV0Mmc = new TH2F("fV0Mmc", "V0 Multiplicity, before correction;V0 sector",64,-.5,63.5,500,0,1000);
 127
 128   fSplineArrayV0A = new TObjArray();
 129   fSplineArrayV0A->SetOwner();
 130   fSplineArrayV0C = new TObjArray();
 131   fSplineArrayV0C->SetOwner();
 132   fSplineArrayV0Agen = new TObjArray();
 133   fSplineArrayV0Agen->SetOwner();
 134   fSplineArrayV0Cgen = new TObjArray();
 135   fSplineArrayV0Cgen->SetOwner();
 136
 137   fOutput->Add(fHistoCuts);
 138   fOutput->Add(fHistoVtxBefSel);
 139   fOutput->Add(fHistoVtxAftSel);
 140   fOutput->Add(fHistoEtaBefSel);
 141   fOutput->Add(fHistoEtaAftSel);
 142   fOutput->Add(fHistoNChAftSel);
 143   fOutput->Add(fHistoQVector);
 144   fOutput->Add(fHistoEP);
 145   fOutput->Add(fPsiACor);
 146   fOutput->Add(fPsiCCor);
 147   fOutput->Add(fQVecACor);
 148   fOutput->Add(fQVecCCor);
 149   fOutput->Add(fV0M);
 150   fOutput->Add(fV0MCor);
 151   fOutput->Add(fV0Mmc);
 152
 153   for (Int_t i = 0; i<10; i++){
 154     fMeanQxa2[i] = -1;
 155     fMeanQya2[i] = -1;
 156     fMeanQxc2[i] = -1;
 157     fMeanQyc2[i] = -1;
 158   }
 159 }
 160
 161 //______________________________________________________
 162 Bool_t AliSpectraAODEventCuts::IsSelected(AliAODEvent * aod,AliSpectraAODTrackCuts     *trackcuts)
 163 {
 164   // Returns true if Event Cuts are selected and applied
 165   fAOD = aod;
 166   fTrackCuts = trackcuts; // FIXME: if track cuts is 0, do not set (use the pre-set member). Do we need to pass this here at all??
 167   // FIXME: all those references by name are slow.
 168   ((TH1I*)fOutput->FindObject("fHistoCuts"))->Fill(kProcessedEvents);
 169   Bool_t IsPhysSel = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & fSelectBit);
 170   if(!IsPhysSel)return IsPhysSel;
 171   ((TH1I*)fOutput->FindObject("fHistoCuts"))->Fill(kPhysSelEvents);
 172   //loop on tracks, before event selection, filling QA histos
 173   AliAODVertex * vertex = fAOD->GetPrimaryVertex();//FIXME vertex is recreated
 174   if(vertex)((TH1F*)fOutput->FindObject("fHistoVtxBefSel"))->Fill(vertex->GetZ());
 175   fIsSelected =kFALSE;
 176   if(CheckVtxRange() && CheckCentralityCut() && CheckMultiplicityCut()){ //selection on vertex and Centrality
 177     fIsSelected=CheckQVectorCut(); // QVector is calculated only if the centrality and vertex are correct (performance)
 178   }
 179   if(fIsSelected){
 180     ((TH1I*)fOutput->FindObject("fHistoCuts"))->Fill(kAcceptedEvents);
 181     if(vertex)((TH1F*)fOutput->FindObject("fHistoVtxAftSel"))->Fill(vertex->GetZ());
 182   }
 183   Int_t Nch=0;
 184   for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) {
 185     AliAODTrack* track = fAOD->GetTrack(iTracks);
 186     if (!fTrackCuts->IsSelected(track,kFALSE)) continue;
 187     ((TH1F*)fOutput->FindObject("fHistoEtaBefSel"))->Fill(track->Eta());
 188     if(fIsSelected){
 189       ((TH1F*)fOutput->FindObject("fHistoEtaAftSel"))->Fill(track->Eta());
 190       Nch++;
 191     }
 192   }
 193   fNch = Nch;
 194   if(fIsSelected)((TH1F*)fOutput->FindObject("fHistoNChAftSel"))->Fill(Nch);
 195   return fIsSelected;
 196 }
 197
 198 //______________________________________________________
 199 Bool_t AliSpectraAODEventCuts::CheckVtxRange()
 200 {
 201   // reject events outside of range
 202   AliAODVertex * vertex = fAOD->GetPrimaryVertex();
 203   //when moving to 2011 wìone has to add a cut using SPD vertex.
 204   //The point is that for events with |z|>20 the vertexer tracks is not working (only 2011!). One has to put a safety cut using SPD vertex large e.g. 15cm
 205   if (!vertex)
 206     {
 207       ((TH1I*)fOutput->FindObject("fHistoCuts"))->Fill(kVtxNoEvent);
 208       return kFALSE;
 209     }
 210   if (vertex->GetZ() > fVertexCutMin && vertex->GetZ() < fVertexCutMax)
 211     {
 212       return kTRUE;
 213     }
 214   ((TH1I*)fOutput->FindObject("fHistoCuts"))->Fill(kVtxRange);
 215   return kFALSE;
 216 }
 217
 218 //______________________________________________________
 219 Bool_t AliSpectraAODEventCuts::CheckCentralityCut()
 220 {
 221   // Check centrality cut
 222   fCent=-999.;
 223   fCent=fAOD->GetCentrality()->GetCentralityPercentile(fCentralityMethod.Data());
 224   if ( (fCent <= fCentralityCutMax)  &&  (fCent >= fCentralityCutMin) )  return kTRUE;
 225   ((TH1I*)fOutput->FindObject("fHistoCuts"))->Fill(kVtxCentral);
 226   return kFALSE;
 227 }
 228
 229 //______________________________________________________
 230 Bool_t AliSpectraAODEventCuts::CheckMultiplicityCut()
 231 {
 232   // Check multiplicity cut
 233   // FIXME: why this is not tracket in the track stats histos?
 234   Int_t Ncharged=0;
 235   for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++){
 236     AliAODTrack* track = fAOD->GetTrack(iTracks);
 237     if (!fTrackCuts->IsSelected(track,kFALSE)) continue;
 238     Ncharged++;
 239   }
 240   if(Ncharged>fMultiplicityCutMin && Ncharged<fMultiplicityCutMax)return kTRUE;
 241
 242   return kFALSE;
 243 }
 244
 245 //______________________________________________________
 246 Bool_t AliSpectraAODEventCuts::CheckQVectorCut()
 247 {
 248   Double_t qVZERO = -999.;
 249   Double_t psi=-999.;
 250
 251   if(fIsLHC10h){
 252     qVZERO=CalculateQVectorLHC10h();
 253     psi=fPsiV0A;
 254   }else{
 255     qVZERO=CalculateQVector();
 256     psi=fPsiV0A;
 257   }
 258
 259   //cut on q vector
 260   if(qVZERO<fQVectorCutMin || qVZERO>fQVectorCutMax)return kFALSE;
 261   Double_t cent=fAOD->GetCentrality()->GetCentralityPercentile(fCentralityMethod.Data());
 262   ((TH2F*)fOutput->FindObject("fHistoQVector"))->Fill(cent,qVZERO);
 263   ((TH2F*)fOutput->FindObject("fHistoEP"))->Fill(cent,psi);
 264   ((TH1I*)fOutput->FindObject("fHistoCuts"))->Fill(kQVector);
 265
 266
 267   return kTRUE;
 268 }
 269
 270 //______________________________________________________
 271 Double_t AliSpectraAODEventCuts::CalculateQVectorLHC10h(){
 272
 273   Int_t run = fAOD->GetRunNumber();
 274   if(run != fRun){
 275     // Load the calibrations run dependent
 276     if(OpenInfoCalbration(run))fRun=run;
 277     else{
 278       fqV0C=-999.;
 279       fqV0A=-999.;
 280       return -999.;
 281     }
 282   }
 283
 284   //V0 info
 285   Double_t Qxa2 = 0, Qya2 = 0;
 286   Double_t Qxc2 = 0, Qyc2 = 0;
 287   Double_t sumMc = 0, sumMa = 0;
 288
 289   AliAODVZERO* aodV0 = fAOD->GetVZEROData();
 290
 291   for (Int_t iv0 = 0; iv0 < 64; iv0++) {
 292
 293     Double_t phiV0 = TMath::PiOver4()*(0.5 + iv0 % 8);
 294
 295     Float_t multv0 = aodV0->GetMultiplicity(iv0);
 296     ((TH2F*)fOutput->FindObject("fV0M"))->Fill(iv0,multv0);
 297
 298     if (iv0 < 32){
 299
 300       Double_t multCorC = -10;
 301
 302       if (iv0 < 8)
 303         multCorC = multv0*fV0Cpol1/fMultV0->GetBinContent(iv0+1);
 304       else if (iv0 >= 8 && iv0 < 16)
 305         multCorC = multv0*fV0Cpol2/fMultV0->GetBinContent(iv0+1);
 306       else if (iv0 >= 16 && iv0 < 24)
 307         multCorC = multv0*fV0Cpol3/fMultV0->GetBinContent(iv0+1);
 308       else if (iv0 >= 24 && iv0 < 32)
 309         multCorC = multv0*fV0Cpol4/fMultV0->GetBinContent(iv0+1);
 310
 311       if (multCorC < 0){
 312         cout<<"Problem with multiplicity in V0C"<<endl;
 313         fqV0C=-999.;
 314         fqV0A=-999.;
 315         return -999.;
 316       }
 317
 318       Qxc2 += TMath::Cos(2*phiV0) * multCorC;
 319       Qyc2 += TMath::Sin(2*phiV0) * multCorC;
 320
 321       sumMc += multCorC;
 322       ((TH2F*)fOutput->FindObject("fV0MCor"))->Fill(iv0,multCorC);
 323
 324     } else {
 325
 326       Double_t multCorA = -10;
 327
 328       if (iv0 >= 32 && iv0 < 40)
 329         multCorA = multv0*fV0Apol1/fMultV0->GetBinContent(iv0+1);
 330       else if (iv0 >= 40 && iv0 < 48)
 331         multCorA = multv0*fV0Apol2/fMultV0->GetBinContent(iv0+1);
 332       else if (iv0 >= 48 && iv0 < 56)
 333         multCorA = multv0*fV0Apol3/fMultV0->GetBinContent(iv0+1);
 334       else if (iv0 >= 56 && iv0 < 64)
 335         multCorA = multv0*fV0Apol4/fMultV0->GetBinContent(iv0+1);
 336
 337       if (multCorA < 0){
 338         cout<<"Problem with multiplicity in V0A"<<endl;
 339         fqV0C=-999.;
 340         fqV0A=-999.;
 341         return -999.;
 342       }
 343
 344       Qxa2 += TMath::Cos(2*phiV0) * multCorA;
 345       Qya2 += TMath::Sin(2*phiV0) * multCorA;
 346
 347       sumMa += multCorA;
 348       ((TH2F*)fOutput->FindObject("fV0MCor"))->Fill(iv0,multCorA);
 349     }
 350   }
 351
 352   Short_t centrV0  = GetCentrCode(fAOD);
 353
 354   Double_t Qxamean2 = 0.;
 355   Double_t Qyamean2 = 0.;
 356   Double_t Qxcmean2 = 0.;
 357   Double_t Qycmean2 = 0.;
 358
 359   if(centrV0!=-1){
 360     Qxamean2 = fMeanQxa2[centrV0];
 361     Qyamean2 = fMeanQya2[centrV0];
 362     Qxcmean2 = fMeanQxc2[centrV0];
 363     Qycmean2 = fMeanQyc2[centrV0];
 364   }
 365
 366   Double_t QxaCor2 = Qxa2 - Qxamean2*sumMa;
 367   Double_t QyaCor2 = Qya2 - Qyamean2*sumMa;
 368   Double_t QxcCor2 = Qxc2 - Qxcmean2*sumMc;
 369   Double_t QycCor2 = Qyc2 - Qycmean2*sumMc;
 370
 371   fPsiV0A = TMath::ATan2(QyaCor2, QxaCor2)/2.;
 372   fPsiV0C = TMath::ATan2(QycCor2, QxcCor2)/2.;
 373
 374   ((TH2F*)fOutput->FindObject("fPsiACor"))->Fill((Float_t)fAOD->GetCentrality()->GetCentralityPercentile("V0M"), fPsiV0A);
 375   ((TH2F*)fOutput->FindObject("fPsiCCor"))->Fill((Float_t)fAOD->GetCentrality()->GetCentralityPercentile("V0M"), fPsiV0C);
 376
 377   fqV0A = TMath::Sqrt((QxaCor2*QxaCor2 + QyaCor2*QyaCor2)/sumMa);
 378   fqV0C = TMath::Sqrt((QxcCor2*QxcCor2 + QycCor2*QycCor2)/sumMc);
 379   fqV0Ax = QxaCor2*TMath::Sqrt(1./sumMa);
 380   fqV0Cx = QxcCor2*TMath::Sqrt(1./sumMc);
 381   fqV0Ay = QyaCor2*TMath::Sqrt(1./sumMa);
 382   fqV0Cy = QycCor2*TMath::Sqrt(1./sumMc);
 383
 384   ((TH2F*)fOutput->FindObject("fQVecACor"))->Fill((Float_t)fAOD->GetCentrality()->GetCentralityPercentile("V0M"), fqV0A);
 385   ((TH2F*)fOutput->FindObject("fQVecCCor"))->Fill((Float_t)fAOD->GetCentrality()->GetCentralityPercentile("V0M"), fqV0C);
 386
 387   return fqV0A; //FIXME we have to combine VZERO-A and C
 388 }
 389
 390 //______________________________________________________
 391 Double_t AliSpectraAODEventCuts::CalculateQVector(){
 392
 393   //V0 info
 394   Double_t Qxa2 = 0, Qya2 = 0;
 395   Double_t Qxc2 = 0, Qyc2 = 0;
 396
 397   AliAODVZERO* aodV0 = fAOD->GetVZEROData();
 398
 399   for (Int_t iv0 = 0; iv0 < 64; iv0++) {
 400
 401     Float_t multv0 = aodV0->GetMultiplicity(iv0);
 402
 403     ((TH2F*)fOutput->FindObject("fV0M"))->Fill(iv0,multv0);
 404
 405   }
 406
 407   fPsiV0A = fAOD->GetEventplane()->CalculateVZEROEventPlane(fAOD,8,2,Qxa2,Qya2); // V0A
 408   fPsiV0C = fAOD->GetEventplane()->CalculateVZEROEventPlane(fAOD,9,2,Qxc2,Qyc2); // V0C
 409
 410   ((TH2F*)fOutput->FindObject("fPsiACor"))->Fill((Float_t)fAOD->GetCentrality()->GetCentralityPercentile("V0M"), fPsiV0A);
 411   ((TH2F*)fOutput->FindObject("fPsiCCor"))->Fill((Float_t)fAOD->GetCentrality()->GetCentralityPercentile("V0M"), fPsiV0C);
 412
 413   fqV0A = TMath::Sqrt((Qxa2*Qxa2 + Qya2*Qya2));
 414   fqV0C = TMath::Sqrt((Qxc2*Qxc2 + Qyc2*Qyc2));
 415   fqV0Ax = Qxa2;
 416   fqV0Cx = Qxc2;
 417   fqV0Ay = Qya2;
 418   fqV0Cy = Qyc2;
 419
 420   ((TH2F*)fOutput->FindObject("fQVecACor"))->Fill((Float_t)fAOD->GetCentrality()->GetCentralityPercentile("V0M"), fqV0A);
 421   ((TH2F*)fOutput->FindObject("fQVecCCor"))->Fill((Float_t)fAOD->GetCentrality()->GetCentralityPercentile("V0M"), fqV0C);
 422
 423   return fqV0A; //FIXME we have to combine VZERO-A and C
 424
 425 }
 426
 427 //______________________________________________________
 428 Short_t AliSpectraAODEventCuts::GetCentrCode(AliVEvent* ev)
 429 {
 430
 431   Short_t centrCode = -1;
 432
 433   AliCentrality* centrality = 0;
 434   AliAODEvent* aod = (AliAODEvent*)ev;
 435   centrality = aod->GetHeader()->GetCentralityP();
 436
 437   Float_t centV0 = centrality->GetCentralityPercentile("V0M");
 438   Float_t centTrk = centrality->GetCentralityPercentile("TRK");
 439
 440
 441   if (TMath::Abs(centV0 - centTrk) < 5.0 && centV0 <= 80 && centV0 > 0){
 442
 443      if ((centV0 > 0) && (centV0 <= 5.0))
 444        centrCode = 0;
 445      else if ((centV0 > 5.0) && (centV0 <= 10.0))
 446        centrCode = 1;
 447      else if ((centV0 > 10.0) && (centV0 <= 20.0))
 448        centrCode = 2;
 449      else if ((centV0 > 20.0) && (centV0 <= 30.0))
 450        centrCode = 3;
 451      else if ((centV0 > 30.0) && (centV0 <= 40.0))
 452        centrCode = 4;
 453      else if ((centV0 > 40.0) && (centV0 <= 50.0))
 454        centrCode = 5;
 455      else if ((centV0 > 50.0) && (centV0 <= 60.0))
 456        centrCode = 6;
 457      else if ((centV0 > 60.0) && (centV0 <= 70.0))
 458        centrCode = 7;
 459      else if ((centV0 > 70.0) && (centV0 <= 80.0))
 460        centrCode = 8;
 461   }
 462
 463   return centrCode;
 464
 465 }
 466
 467 //______________________________________________________
 468 void AliSpectraAODEventCuts::PrintCuts()
 469 {
 470   // print info about event cuts
 471   cout << "Event Stats" << endl;
 472   cout << " > Trigger Selection: " << fSelectBit << endl;
 473   cout << " > Centrality estimator: " << fCentralityMethod << endl;
 474   cout << " > Number of accepted events: " << NumberOfEvents() << endl;
 475   cout << " > Number of processed events: " << NumberOfProcessedEvents() << endl;
 476   cout << " > Number of PhysSel events: " <<  NumberOfPhysSelEvents()  << endl;
 477   cout << " > Vertex out of range: " << ((TH1I*)fOutput->FindObject("fHistoCuts"))->GetBinContent(kVtxRange + 1) << endl;
 478   cout << " > Events cut by centrality: " << ((TH1I*)fOutput->FindObject("fHistoCuts"))->GetBinContent(kVtxCentral + 1) << endl;
 479   cout << " > Events without vertex: " << ((TH1I*)fOutput->FindObject("fHistoCuts"))->GetBinContent(kVtxNoEvent + 1) << endl;
 480   cout << " > QVector cut: " << ((TH1I*)fOutput->FindObject("fHistoCuts"))->GetBinContent(kQVector + 1) << endl;
 481   cout << " > Track type used for the QVector calculation: " << fTrackBits << endl;
 482   cout << " > QRange: [" << fQVectorCutMin <<"," <<fQVectorCutMax<<"]"<< endl;
 483   cout << " > Vertex: [" << fVertexCutMin <<"," <<fVertexCutMax<<"]"<< endl;
 484   cout << " > Multiplicity: [" << fMultiplicityCutMin <<"," <<fMultiplicityCutMax<<"]"<< endl;
 485   cout << " > Centrality: [" << fCentralityCutMin <<"," <<fCentralityCutMax<<"]"<< endl;
 486 }
 487
 488 //_____________________________________________________________________________
 489 Bool_t AliSpectraAODEventCuts::OpenInfoCalbration(Int_t run)
 490 {
 491
 492   AliOADBContainer* cont = (AliOADBContainer*) fCalib->FindObject("hMultV0BefCorr");
 493   if(!cont){
 494     printf("OADB object hMultV0BefCorr is not available in the file\n");
 495     return 0;
 496   }
 497
 498   if(!(cont->GetObject(run))){
 499     printf("OADB object hMultV0BefCorr is not available for run %i\n",run);
 500     return 0;
 501   }
 502   fMultV0 = ((TH2F*) cont->GetObject(run))->ProfileX();
 503
 504   TF1* fpolc1 = new TF1("fpolc1","pol0", 0, 7);
 505   fMultV0->Fit(fpolc1, "RN");
 506   fV0Cpol1 = fpolc1->GetParameter(0);
 507
 508   TF1* fpolc2 = new TF1("fpolc2","pol0", 8, 15);
 509   fMultV0->Fit(fpolc2, "RN");
 510   fV0Cpol2 = fpolc2->GetParameter(0);
 511
 512   TF1* fpolc3 = new TF1("fpolc3","pol0", 16, 23);
 513   fMultV0->Fit(fpolc3, "RN");
 514   fV0Cpol3 = fpolc3->GetParameter(0);
 515
 516   TF1* fpolc4 = new TF1("fpolc4","pol0", 24, 31);
 517   fMultV0->Fit(fpolc4, "RN");
 518   fV0Cpol4 = fpolc4->GetParameter(0);
 519
 520   TF1* fpola1 = new TF1("fpola1","pol0", 32, 39);
 521   fMultV0->Fit(fpola1, "RN");
 522   fV0Apol1 = fpola1->GetParameter(0);
 523
 524   TF1* fpola2 = new TF1("fpola2","pol0", 40, 47);
 525   fMultV0->Fit(fpola2, "RN");
 526   fV0Apol2 = fpola2->GetParameter(0);
 527
 528   TF1* fpola3 = new TF1("fpola3","pol0", 48, 55);
 529   fMultV0->Fit(fpola3, "RN");
 530   fV0Apol3 = fpola3->GetParameter(0);
 531
 532   TF1* fpola4 = new TF1("fpola4","pol0", 56, 63);
 533   fMultV0->Fit(fpola4, "RN");
 534   fV0Apol4 = fpola4->GetParameter(0);
 535
 536   for(Int_t i=0; i < 10; i++){
 537
 538     char nameQxa2[100];
 539     snprintf(nameQxa2,100, "hQxa2m_%i", i);
 540
 541     char nameQya2[100];
 542     snprintf(nameQya2,100, "hQya2m_%i", i);
 543
 544     char nameQxc2[100];
 545     snprintf(nameQxc2,100, "hQxc2m_%i", i);
 546
 547     char nameQyc2[100];
 548     snprintf(nameQyc2,100, "hQyc2m_%i", i);
 549
 550     AliOADBContainer* contQxa2 = (AliOADBContainer*) fCalib->FindObject(nameQxa2);
 551     if(!contQxa2){
 552       printf("OADB object %s is not available in the file\n", nameQxa2);
 553       return 0;
 554     }
 555
 556     if(!(contQxa2->GetObject(run))){
 557       printf("OADB object %s is not available for run %i\n", nameQxa2, run);
 558       return 0;
 559     }
 560
 561     fMeanQxa2[i] = ((TH1F*) contQxa2->GetObject(run))->GetMean();
 562
 563
 564     AliOADBContainer* contQya2 = (AliOADBContainer*) fCalib->FindObject(nameQya2);
 565     if(!contQya2){
 566       printf("OADB object %s is not available in the file\n", nameQya2);
 567       return 0;
 568     }
 569
 570     if(!(contQya2->GetObject(run))){
 571       printf("OADB object %s is not available for run %i\n", nameQya2, run);
 572       return 0;
 573     }
 574
 575     fMeanQya2[i] = ((TH1F*) contQya2->GetObject(run))->GetMean();
 576
 577
 578     AliOADBContainer* contQxc2 = (AliOADBContainer*) fCalib->FindObject(nameQxc2);
 579     if(!contQxc2){
 580       printf("OADB object %s is not available in the file\n", nameQxc2);
 581       return 0;
 582     }
 583
 584     if(!(contQxc2->GetObject(run))){
 585       printf("OADB object %s is not available for run %i\n", nameQxc2, run);
 586       return 0;
 587     }
 588
 589     fMeanQxc2[i] = ((TH1F*) contQxc2->GetObject(run))->GetMean();
 590
 591
 592     AliOADBContainer* contQyc2 = (AliOADBContainer*) fCalib->FindObject(nameQyc2);
 593     if(!contQyc2){
 594       printf("OADB object %s is not available in the file\n", nameQyc2);
 595       return 0;
 596     }
 597
 598     if(!(contQyc2->GetObject(run))){
 599       printf("OADB object %s is not available for run %i\n", nameQyc2, run);
 600       return 0;
 601     }
 602
 603     fMeanQyc2[i] = ((TH1F*) contQyc2->GetObject(run))->GetMean();
 604
 605   }
 606   return 1;
 607 }
 608
 609 //______________________________________________________
 610
 611
 612 Long64_t AliSpectraAODEventCuts::Merge(TCollection* list)
 613 {
 614   // Merge a list of AliSpectraAODEventCuts objects with this.
 615   // Returns the number of merged objects (including this).
 616
 617   AliInfo("Merging");
 618
 619   if (!list)
 620     return 0;
 621
 622   if (list->IsEmpty())
 623     return 1;
 624
 625   TIterator* iter = list->MakeIterator();
 626   TObject* obj;
 627
 628   // collections of all histograms
 629   TList collections;
 630
 631   Int_t count = 0;
 632
 633   while ((obj = iter->Next())) {
 634     AliSpectraAODEventCuts* entry = dynamic_cast<AliSpectraAODEventCuts*> (obj);
 635     if (entry == 0)
 636       continue;
 637
 638     TList * l = entry->GetOutputList();
 639     collections.Add(l);
 640     count++;
 641   }
 642
 643   fOutput->Merge(&collections);
 644
 645   delete iter;
 646
 647   return count+1;
 648 }
 649
 650 //______________________________________________________
 651 Double_t AliSpectraAODEventCuts::GetQvecPercentile(Int_t v0side){
 652
 653   //check Qvec and Centrality consistency
 654   if(fCent>90.) return -999.;
 655   if(v0side==0/*V0A*/){ if(fqV0A== -999.) return -999.; }
 656   if(v0side==1/*V0C*/){ if(fqV0C== -999.) return -999.; }
 657
 658   fQvecIntegral = 0x0;
 659
 660   if(v0side==0/*V0A*/){ fQvecIntegral = (TH2D*)fQvecIntList->FindObject("VZEROA"); }
 661   if(v0side==1/*V0C*/){ fQvecIntegral = (TH2D*)fQvecIntList->FindObject("VZEROC"); }
 662
 663
 664   Int_t ic = -999;
 665
 666   if(fQvecCalibType==1){
 667     if(fNch<0.) return -999.;
 668     ic = GetNchBin(fQvecIntegral);
 669   } else ic = (Int_t)fCent; //fQvecIntegral: 1% centrality bin
 670
 671   TH1D *h1D = (TH1D*)fQvecIntegral->ProjectionY("h1D",ic+1,ic+1);
 672
 673   TSpline *spline = 0x0;
 674
 675   if(v0side==0/*V0A*/){
 676     if( CheckSplineArray(fSplineArrayV0A, ic) ) {
 677       spline = (TSpline*)fSplineArrayV0A->At(ic);
 678       //cout<<"Qvec V0A - ic: "<<ic<<" - Found TSpline..."<<endl;
 679     } else {
 680       spline = new TSpline3(h1D,"sp3");
 681       fSplineArrayV0A->AddAtAndExpand(spline,ic);
 682       //cout<<"Qvec V0A - ic: "<<ic<<" - TSpline not found. Creating..."<<endl;
 683     }
 684   }
 685   else if(v0side==1/*V0C*/){
 686     if( CheckSplineArray(fSplineArrayV0C, ic) ) {
 687       spline = (TSpline*)fSplineArrayV0C->At(ic);
 688     } else {
 689       spline = new TSpline3(h1D,"sp3");
 690       fSplineArrayV0C->AddAtAndExpand(spline,ic);
 691     }
 692   }
 693
 694   Double_t percentile =  -999.;
 695   if(v0side==0/*V0A*/){ percentile = 100*spline->Eval(fqV0A); }
 696   if(v0side==1/*V0C*/){ percentile = 100*spline->Eval(fqV0C); }
 697
 698   if(percentile>100. || percentile<0.) return -999.;
 699
 700   return percentile;
 701 }
 702
 703 //______________________________________________________
 704 Double_t AliSpectraAODEventCuts::CalculateQVectorMC(Int_t v0side, Int_t type=1){
 705
 706   if(!fIsMC) return -999.;
 707
 708   // V0A efficiecy
 709   if(type==1){
 710     fV0Aeff = 0x0;
 711     fV0Aeff = (TH1F*)fQvecIntList->FindObject("vzeroa_efficiency_prim_plus_sec_over_gen");
 712     if(!fV0Aeff) return -999.;
 713   }
 714
 715   // get MC array
 716   TClonesArray *arrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 717
 718   if (!arrayMC) AliFatal("Error: MC particles branch not found!\n");
 719
 720   Double_t Qx2mc = 0., Qy2mc = 0.;
 721   Double_t mult2mc = 0;
 722
 723   Int_t nMC = arrayMC->GetEntries();
 724
 725   if(type==0){ // type==0: q-vec from tracks in vzero acceptance
 726
 727     // loop on generated
 728     for (Int_t iMC = 0; iMC < nMC; iMC++){
 729       AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iMC);
 730       if(!partMC->Charge()) continue;//Skip neutrals
 731
 732       // check vzero side
 733       if( CheckVZEROacceptance(partMC->Eta()) != v0side ) continue;
 734
 735       // Calculate Qvec components
 736       Qx2mc += TMath::Cos(2.*partMC->Phi());
 737       Qy2mc += TMath::Sin(2.*partMC->Phi());
 738       mult2mc++;
 739
 740     }// end loop on generated.
 741
 742   }//end if on type==0
 743
 744   else if(type==1){ // type==1 (default): q-vec from vzero
 745
 746     // only used in qgen_vzero
 747     Double_t multv0mc[64];
 748     for(Int_t iCh=0;iCh<64;iCh++) multv0mc[iCh]=0; // initialize multv0mc to zero
 749
 750     // loop on generated
 751     for (Int_t iMC = 0; iMC < nMC; iMC++){
 752
 753       AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iMC);
 754       if(!partMC->Charge()) continue;//Skip neutrals
 755
 756       // check vzero side
 757       if( CheckVZEROacceptance(partMC->Eta()) != v0side ) continue;
 758
 759       //get v0 channel of gen track
 760       Int_t iv0 = CheckVZEROchannel(v0side, partMC->Eta(), partMC->Phi());
 761
 762       //use the efficiecy as a weigth
 763       multv0mc[iv0] += fV0Aeff->GetFunction("f")->Eval(partMC->Pt());
 764
 765       //calculate multiplicity for each vzero channell
 766       //multv0mc[iv0]++;
 767
 768     }// end loop on generated.
 769
 770     Int_t firstCh=-1,lastCh=-1;
 771     if (v0side == 0) {
 772       firstCh = 32;
 773       lastCh = 64;
 774     }
 775     else if (v0side == 1) {
 776       firstCh = 0;
 777       lastCh = 32;
 778     }
 779     for(Int_t iCh = firstCh; iCh < lastCh; ++iCh) {
 780       Double_t phi = TMath::Pi()/8. + TMath::Pi()/4.*(iCh%8);
 781       Double_t mult = multv0mc[iCh];
 782       Qx2mc += mult*TMath::Cos(2*phi);
 783       Qy2mc += mult*TMath::Sin(2*phi);
 784       mult2mc += mult;
 785
 786       ((TH2F*)fOutput->FindObject("fV0Mmc"))->Fill(iCh,multv0mc[iCh]);
 787     }
 788
 789   }//end if on type==1
 790
 791   // return q vector
 792   fQvecMC = TMath::Sqrt((Qx2mc*Qx2mc + Qy2mc*Qy2mc)/mult2mc);
 793
 794   return fQvecMC;
 795
 796 }
 797
 798 //______________________________________________________
 799 Int_t AliSpectraAODEventCuts::CheckVZEROchannel(Int_t vzeroside, Double_t eta, Double_t phi){
 800
 801   //VZEROA eta acceptance
 802   Int_t ring = -1.;
 803
 804   if ( vzeroside == 0){
 805     if (eta > 4.5 && eta <= 5.1 ) ring = 0;
 806     if (eta > 3.9 && eta <= 4.5 ) ring = 1;
 807     if (eta > 3.4 && eta <= 3.9 ) ring = 2;
 808     if (eta > 2.8 && eta <= 3.4 ) ring = 3;
 809   } else if (vzeroside == 1){
 810     if (eta > -3.7 && eta <= -3.2 ) ring = 0;
 811     if (eta > -3.2 && eta <= -2.7 ) ring = 1;
 812     if (eta > -2.7 && eta <= -2.2 ) ring = 2;
 813     if (eta > -2.2 && eta <= -1.7 ) ring = 3;
 814   }
 815
 816
 817   //VZEROAC phi acceptance
 818   Int_t phisector= -1;
 819
 820
 821   if ( phi > 0. && phi <= TMath::Pi()/4. ) phisector = 0;
 822
 823   else if (phi > TMath::Pi()/4. && phi <= TMath::Pi()/2.) phisector = 1;
 824
 825   else if (phi > TMath::Pi()/2 && phi <= (3./4.)*TMath::Pi() ) phisector =2;
 826
 827   else if (phi > (3./4.)*TMath::Pi() && phi <= TMath::Pi() ) phisector = 3;
 828
 829   else if (phi > TMath::Pi() && phi <= (5./4.)*TMath::Pi() ) phisector = 4;
 830
 831   else if (phi > (5./4.)*TMath::Pi() && phi <= (3./2.)*TMath::Pi() ) phisector = 5;
 832
 833   else if (phi > (3./2.)*TMath::Pi() && phi <= (7./4.)*TMath::Pi() ) phisector = 6;
 834
 835   else if (phi > (7./4.)*TMath::Pi() && phi <= TMath::TwoPi() ) phisector = 7;
 836
 837   if (vzeroside ==0) return Int_t(32+(phisector+(ring*8.))); // iv0 >= 32 -> V0A
 838   if (vzeroside ==1) return Int_t(phisector+(ring*8.));  // iv0 < 32 -> V0C
 839
 840   else return -999.;
 841 }
 842
 843 //______________________________________________________
 844 Int_t AliSpectraAODEventCuts::CheckVZEROacceptance(Double_t eta){
 845
 846   // eval VZERO side - FIXME Add TPC!
 847
 848   if(eta > 2.8  && eta < 5.1) return 0; //VZEROA
 849
 850   else if (eta > -3.7  && eta < -1.7) return 1; //VZEROC
 851
 852   return -999.;
 853
 854 }
 855
 856 //______________________________________________________
 857 Double_t AliSpectraAODEventCuts::GetQvecPercentileMC(Int_t v0side, Int_t type=1){
 858
 859   //check Qvec and Centrality consistency
 860   if(fCent>90.) return -999.;
 861
 862   Double_t qvec = CalculateQVectorMC(v0side, type);
 863   if(qvec==-999.) return -999.;
 864
 865   fQgenIntegral = 0x0;
 866
 867   if(type==0){
 868     if(v0side==0/*V0A*/){ fQgenIntegral = (TH2D*)fQvecIntList->FindObject("VZEROAgen"); }
 869     if(v0side==1/*V0C*/){ fQgenIntegral = (TH2D*)fQvecIntList->FindObject("VZEROCgen"); }
 870   } else if (type==1){
 871     if(v0side==0/*V0A*/){ fQgenIntegral = (TH2D*)fQvecIntList->FindObject("VZEROAgen_vzero"); }
 872     if(v0side==1/*V0C*/){ fQgenIntegral = (TH2D*)fQvecIntList->FindObject("VZEROCgen_vzero"); }
 873   }
 874
 875
 876   //FIXME you need a check on the TH2D, add AliFatal or a return.
 877
 878   Int_t ic = -999;
 879
 880   if(fQvecCalibType==1){
 881     if(fNch<0.) return -999.;
 882     ic = GetNchBin(fQgenIntegral);
 883   } else ic = (Int_t)fCent; //fQvecIntegral: 1% centrality bin
 884
 885   TH1D *h1D = (TH1D*)fQgenIntegral->ProjectionY("h1Dgen",ic+1,ic+1);
 886
 887   TSpline *spline = 0x0;
 888
 889   if(v0side==0/*V0A*/){
 890     if( CheckSplineArray(fSplineArrayV0Agen, ic) ) {
 891       spline = (TSpline*)fSplineArrayV0Agen->At(ic);
 892       //cout<<"Qvec V0A - ic: "<<ic<<" - Found TSpline..."<<endl;
 893     } else {
 894       spline = new TSpline3(h1D,"sp3");
 895       fSplineArrayV0Agen->AddAtAndExpand(spline,ic);
 896     }
 897   }
 898   else if(v0side==1/*V0C*/){
 899     if( CheckSplineArray(fSplineArrayV0Cgen, ic) ) {
 900       spline = (TSpline*)fSplineArrayV0Cgen->At(ic);
 901     } else {
 902       spline = new TSpline3(h1D,"sp3");
 903       fSplineArrayV0Cgen->AddAtAndExpand(spline,ic);
 904     }
 905   }
 906
 907   Double_t percentile = 100*spline->Eval(qvec);
 908
 909   if(percentile>100. || percentile<0.) return -999.;
 910
 911   return percentile;
 912
 913 }
 914
 915 //______________________________________________________
 916 Bool_t AliSpectraAODEventCuts::CheckSplineArray(TObjArray * splarr, Int_t n){
 917   //check TSpline array consistency
 918
 919 //   Int_t n = (Int_t)fCent; //FIXME should be ok for icentr and nch
 920
 921   if(splarr->GetSize()<n) return kFALSE;
 922
 923   if(!splarr->At(n)) return kFALSE;
 924
 925   return kTRUE;
 926
 927 }
 928
 929 //______________________________________________________
 930 Int_t AliSpectraAODEventCuts::GetNchBin(TH2D * h){
 931
 932   Double_t xmax = h->GetXaxis()->GetXmax();
 933
 934   if(fNch>xmax) return (Int_t)h->GetNbinsX();
 935
 936   return (fNch*h->GetNbinsX())/h->GetXaxis()->GetXmax();
 937
 938 }
 939