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