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 {
  92   // Constructor
  93   fOutput=new TList();
  94   fOutput->SetOwner();
  95   fOutput->SetName("fOutput");
  96
  97   fCalib=new TList();
  98   fCalib->SetOwner();
  99   fCalib->SetName("fCalib");
 100
 101   fQvecIntList=new TList();
 102   fQvecIntList->SetOwner();
 103   fQvecIntList->SetName("fQvecIntList");
 104
 105   TH1I *fHistoCuts = new TH1I("fHistoCuts", "Event Cuts", kNVtxCuts, -0.5, kNVtxCuts - 0.5);
 106   TH1F *fHistoVtxBefSel = new TH1F("fHistoVtxBefSel", "Vtx distr before event selection;z (cm)",500,-15,15);
 107   TH1F *fHistoVtxAftSel = new TH1F("fHistoVtxAftSel", "Vtx distr after event selection;z (cm)",500,-15,15);
 108   TH1F *fHistoEtaBefSel = new TH1F("fHistoEtaBefSel", "Eta distr before event selection;eta",500,-2,2);
 109   TH1F *fHistoEtaAftSel = new TH1F("fHistoEtaAftSel", "Eta distr after event selection;eta",500,-2,2);
 110   TH1F *fHistoNChAftSel = new TH1F("fHistoNChAftSel", "NCh distr after event selection;Nch",2000,-0.5,1999.5);
 111   TH2F *fHistoQVector = new TH2F("fHistoQVector", "QVector with VZERO distribution;centrality;Q vector from EP task",20,0,100,100,0,10);
 112   TH2F *fHistoEP = new TH2F("fHistoEP", "EP with VZERO distribution;centrality;Psi_{EP} from EP task",20,0,100,100,-2,2);
 113   TH2F *fPsiACor = new TH2F("fPsiACor", "EP with VZERO A distribution;centrality;Psi_{EP} VZERO-A",20,0,100,100,-2,2);
 114   TH2F *fPsiCCor = new TH2F("fPsiCCor", "EP with VZERO C distribution;centrality;Psi_{EP} VZERO-C",20,0,100,100,-2,2);
 115   TH2F *fQVecACor = new TH2F("fQVecACor", "QVec VZERO A;centrality;Qvector VZERO-A",20,0,100,100,0,10);
 116   TH2F *fQVecCCor = new TH2F("fQVecCCor", "QVec VZERO C;centrality;Qvector VZERO-C",20,0,100,100,0,10);
 117   TH2F *fV0M = new TH2F("fV0M", "V0 Multiplicity, before correction;V0 sector",64,-.5,63.5,500,0,1000);
 118   TH2F *fV0MCor = new TH2F("fV0MCor", "V0 Multiplicity, after correction;V0 sector",64,-.5,63.5,500,0,1000);
 119
 120   fSplineArrayV0A = new TObjArray();
 121   fSplineArrayV0A->SetOwner();
 122   fSplineArrayV0C = new TObjArray();
 123   fSplineArrayV0C->SetOwner();
 124
 125   fOutput->Add(fHistoCuts);
 126   fOutput->Add(fHistoVtxBefSel);
 127   fOutput->Add(fHistoVtxAftSel);
 128   fOutput->Add(fHistoEtaBefSel);
 129   fOutput->Add(fHistoEtaAftSel);
 130   fOutput->Add(fHistoNChAftSel);
 131   fOutput->Add(fHistoQVector);
 132   fOutput->Add(fHistoEP);
 133   fOutput->Add(fPsiACor);
 134   fOutput->Add(fPsiCCor);
 135   fOutput->Add(fQVecACor);
 136   fOutput->Add(fQVecCCor);
 137   fOutput->Add(fV0M);
 138   fOutput->Add(fV0MCor);
 139
 140   for (Int_t i = 0; i<10; i++){
 141     fMeanQxa2[i] = -1;
 142     fMeanQya2[i] = -1;
 143     fMeanQxc2[i] = -1;
 144     fMeanQyc2[i] = -1;
 145   }
 146 }
 147
 148 //______________________________________________________
 149 Bool_t AliSpectraAODEventCuts::IsSelected(AliAODEvent * aod,AliSpectraAODTrackCuts     *trackcuts)
 150 {
 151   // Returns true if Event Cuts are selected and applied
 152   fAOD = aod;
 153   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??
 154   // FIXME: all those references by name are slow.
 155   ((TH1I*)fOutput->FindObject("fHistoCuts"))->Fill(kProcessedEvents);
 156   Bool_t IsPhysSel = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & fSelectBit);
 157   if(!IsPhysSel)return IsPhysSel;
 158   ((TH1I*)fOutput->FindObject("fHistoCuts"))->Fill(kPhysSelEvents);
 159   //loop on tracks, before event selection, filling QA histos
 160   AliAODVertex * vertex = fAOD->GetPrimaryVertex();//FIXME vertex is recreated
 161   if(vertex)((TH1F*)fOutput->FindObject("fHistoVtxBefSel"))->Fill(vertex->GetZ());
 162   fIsSelected =kFALSE;
 163   if(CheckVtxRange() && CheckCentralityCut() && CheckMultiplicityCut()){ //selection on vertex and Centrality
 164     fIsSelected=CheckQVectorCut(); // QVector is calculated only if the centrality and vertex are correct (performance)
 165   }
 166   if(fIsSelected){
 167     ((TH1I*)fOutput->FindObject("fHistoCuts"))->Fill(kAcceptedEvents);
 168     if(vertex)((TH1F*)fOutput->FindObject("fHistoVtxAftSel"))->Fill(vertex->GetZ());
 169   }
 170   Int_t Nch=0;
 171   for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) {
 172     AliAODTrack* track = fAOD->GetTrack(iTracks);
 173     if (!fTrackCuts->IsSelected(track,kFALSE)) continue;
 174     ((TH1F*)fOutput->FindObject("fHistoEtaBefSel"))->Fill(track->Eta());
 175     if(fIsSelected){
 176       ((TH1F*)fOutput->FindObject("fHistoEtaAftSel"))->Fill(track->Eta());
 177       Nch++;
 178     }
 179   }
 180   if(fIsSelected)((TH1F*)fOutput->FindObject("fHistoNChAftSel"))->Fill(Nch);
 181   return fIsSelected;
 182 }
 183
 184 //______________________________________________________
 185 Bool_t AliSpectraAODEventCuts::CheckVtxRange()
 186 {
 187   // reject events outside of range
 188   AliAODVertex * vertex = fAOD->GetPrimaryVertex();
 189   //when moving to 2011 wìone has to add a cut using SPD vertex.
 190   //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
 191   if (!vertex)
 192     {
 193       ((TH1I*)fOutput->FindObject("fHistoCuts"))->Fill(kVtxNoEvent);
 194       return kFALSE;
 195     }
 196   if (vertex->GetZ() > fVertexCutMin && vertex->GetZ() < fVertexCutMax)
 197     {
 198       return kTRUE;
 199     }
 200   ((TH1I*)fOutput->FindObject("fHistoCuts"))->Fill(kVtxRange);
 201   return kFALSE;
 202 }
 203
 204 //______________________________________________________
 205 Bool_t AliSpectraAODEventCuts::CheckCentralityCut()
 206 {
 207   // Check centrality cut
 208   fCent=-999.;
 209   fCent=fAOD->GetCentrality()->GetCentralityPercentile(fCentralityMethod.Data());
 210   if ( (fCent <= fCentralityCutMax)  &&  (fCent >= fCentralityCutMin) )  return kTRUE;
 211   ((TH1I*)fOutput->FindObject("fHistoCuts"))->Fill(kVtxCentral);
 212   return kFALSE;
 213 }
 214
 215 //______________________________________________________
 216 Bool_t AliSpectraAODEventCuts::CheckMultiplicityCut()
 217 {
 218   // Check multiplicity cut
 219   // FIXME: why this is not tracket in the track stats histos?
 220   Int_t Ncharged=0;
 221   for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++){
 222     AliAODTrack* track = fAOD->GetTrack(iTracks);
 223     if (!fTrackCuts->IsSelected(track,kFALSE)) continue;
 224     Ncharged++;
 225   }
 226   if(Ncharged>fMultiplicityCutMin && Ncharged<fMultiplicityCutMax)return kTRUE;
 227
 228   return kFALSE;
 229 }
 230
 231 //______________________________________________________
 232 Bool_t AliSpectraAODEventCuts::CheckQVectorCut()
 233 {
 234   Double_t qxEPVZERO = -999., qyEPVZERO = -999.;
 235   Double_t qVZERO = -999.;
 236   Double_t psi=-999.;
 237
 238   if(fIsLHC10h){
 239     qVZERO=CalculateQVectorLHC10h();
 240     psi=fPsiV0A;
 241   }else{
 242     psi=fAOD->GetEventplane()->CalculateVZEROEventPlane(fAOD,10,2,qxEPVZERO,qyEPVZERO);//FIXME we can a flag for 2010 and 2011
 243     qVZERO= TMath::Sqrt(qxEPVZERO*qxEPVZERO + qyEPVZERO*qyEPVZERO);
 244   }
 245
 246   //cut on q vector
 247   if(qVZERO<fQVectorCutMin || qVZERO>fQVectorCutMax)return kFALSE;
 248   Double_t cent=fAOD->GetCentrality()->GetCentralityPercentile(fCentralityMethod.Data());
 249   ((TH2F*)fOutput->FindObject("fHistoQVector"))->Fill(cent,qVZERO);
 250   ((TH2F*)fOutput->FindObject("fHistoEP"))->Fill(cent,psi);
 251   ((TH1I*)fOutput->FindObject("fHistoCuts"))->Fill(kQVector);
 252
 253
 254   return kTRUE;
 255 }
 256
 257 //______________________________________________________
 258 Double_t AliSpectraAODEventCuts::CalculateQVectorLHC10h(){
 259
 260   Int_t run = fAOD->GetRunNumber();
 261   if(run != fRun){
 262     // Load the calibrations run dependent
 263     if(OpenInfoCalbration(run))fRun=run;
 264     else{
 265       fqV0C=-999.;
 266       fqV0A=-999.;
 267       return -999.;
 268     }
 269   }
 270
 271   //V0 info
 272   Double_t Qxa2 = 0, Qya2 = 0;
 273   Double_t Qxc2 = 0, Qyc2 = 0;
 274   Double_t sumMc = 0, sumMa = 0;
 275
 276   AliAODVZERO* aodV0 = fAOD->GetVZEROData();
 277
 278   for (Int_t iv0 = 0; iv0 < 64; iv0++) {
 279
 280     Double_t phiV0 = TMath::PiOver4()*(0.5 + iv0 % 8);
 281
 282     Float_t multv0 = aodV0->GetMultiplicity(iv0);
 283     ((TH2F*)fOutput->FindObject("fV0M"))->Fill(iv0,multv0);
 284
 285     if (iv0 < 32){
 286
 287       Double_t multCorC = -10;
 288
 289       if (iv0 < 8)
 290         multCorC = multv0*fV0Cpol1/fMultV0->GetBinContent(iv0+1);
 291       else if (iv0 >= 8 && iv0 < 16)
 292         multCorC = multv0*fV0Cpol2/fMultV0->GetBinContent(iv0+1);
 293       else if (iv0 >= 16 && iv0 < 24)
 294         multCorC = multv0*fV0Cpol3/fMultV0->GetBinContent(iv0+1);
 295       else if (iv0 >= 24 && iv0 < 32)
 296         multCorC = multv0*fV0Cpol4/fMultV0->GetBinContent(iv0+1);
 297
 298       if (multCorC < 0){
 299         cout<<"Problem with multiplicity in V0C"<<endl;
 300         fqV0C=-999.;
 301         fqV0A=-999.;
 302         return -999.;
 303       }
 304
 305       Qxc2 += TMath::Cos(2*phiV0) * multCorC;
 306       Qyc2 += TMath::Sin(2*phiV0) * multCorC;
 307
 308       sumMc += multCorC;
 309       ((TH2F*)fOutput->FindObject("fV0MCor"))->Fill(iv0,multCorC);
 310
 311     } else {
 312
 313       Double_t multCorA = -10;
 314
 315       if (iv0 >= 32 && iv0 < 40)
 316         multCorA = multv0*fV0Apol1/fMultV0->GetBinContent(iv0+1);
 317       else if (iv0 >= 40 && iv0 < 48)
 318         multCorA = multv0*fV0Apol2/fMultV0->GetBinContent(iv0+1);
 319       else if (iv0 >= 48 && iv0 < 56)
 320         multCorA = multv0*fV0Apol3/fMultV0->GetBinContent(iv0+1);
 321       else if (iv0 >= 56 && iv0 < 64)
 322         multCorA = multv0*fV0Apol4/fMultV0->GetBinContent(iv0+1);
 323
 324       if (multCorA < 0){
 325         cout<<"Problem with multiplicity in V0A"<<endl;
 326         fqV0C=-999.;
 327         fqV0A=-999.;
 328         return -999.;
 329       }
 330
 331       Qxa2 += TMath::Cos(2*phiV0) * multCorA;
 332       Qya2 += TMath::Sin(2*phiV0) * multCorA;
 333
 334       sumMa += multCorA;
 335       ((TH2F*)fOutput->FindObject("fV0MCor"))->Fill(iv0,multCorA);
 336     }
 337   }
 338
 339   Short_t centrV0  = GetCentrCode(fAOD);
 340
 341   Double_t Qxamean2 = 0.;
 342   Double_t Qyamean2 = 0.;
 343   Double_t Qxcmean2 = 0.;
 344   Double_t Qycmean2 = 0.;
 345
 346   if(centrV0!=-1){
 347     Qxamean2 = fMeanQxa2[centrV0];
 348     Qyamean2 = fMeanQya2[centrV0];
 349     Qxcmean2 = fMeanQxc2[centrV0];
 350     Qycmean2 = fMeanQyc2[centrV0];
 351   }
 352
 353   Double_t QxaCor2 = Qxa2 - Qxamean2*sumMa;
 354   Double_t QyaCor2 = Qya2 - Qyamean2*sumMa;
 355   Double_t QxcCor2 = Qxc2 - Qxcmean2*sumMc;
 356   Double_t QycCor2 = Qyc2 - Qycmean2*sumMc;
 357
 358   fPsiV0A = TMath::ATan2(QyaCor2, QxaCor2)/2.;
 359   fPsiV0C = TMath::ATan2(QycCor2, QxcCor2)/2.;
 360
 361   ((TH2F*)fOutput->FindObject("fPsiACor"))->Fill((Float_t)fAOD->GetCentrality()->GetCentralityPercentile("V0M"), fPsiV0A);
 362   ((TH2F*)fOutput->FindObject("fPsiCCor"))->Fill((Float_t)fAOD->GetCentrality()->GetCentralityPercentile("V0M"), fPsiV0C);
 363
 364   fqV0A = TMath::Sqrt((QxaCor2*QxaCor2 + QyaCor2*QyaCor2)/sumMa);
 365   fqV0C = TMath::Sqrt((QxcCor2*QxcCor2 + QycCor2*QycCor2)/sumMc);
 366   fqV0Ax = QxaCor2*TMath::Sqrt(1./sumMa);
 367   fqV0Cx = QxcCor2*TMath::Sqrt(1./sumMc);
 368   fqV0Ay = QyaCor2*TMath::Sqrt(1./sumMa);
 369   fqV0Cy = QycCor2*TMath::Sqrt(1./sumMc);
 370
 371   ((TH2F*)fOutput->FindObject("fQVecACor"))->Fill((Float_t)fAOD->GetCentrality()->GetCentralityPercentile("V0M"), fqV0A);
 372   ((TH2F*)fOutput->FindObject("fQVecCCor"))->Fill((Float_t)fAOD->GetCentrality()->GetCentralityPercentile("V0M"), fqV0C);
 373
 374   return fqV0A; //FIXME we have to combine VZERO-A and C
 375 }
 376
 377 //______________________________________________________
 378 Short_t AliSpectraAODEventCuts::GetCentrCode(AliVEvent* ev)
 379 {
 380
 381   Short_t centrCode = -1;
 382
 383   AliCentrality* centrality = 0;
 384   AliAODEvent* aod = (AliAODEvent*)ev;
 385   centrality = aod->GetHeader()->GetCentralityP();
 386
 387   Float_t centV0 = centrality->GetCentralityPercentile("V0M");
 388   Float_t centTrk = centrality->GetCentralityPercentile("TRK");
 389
 390
 391   if (TMath::Abs(centV0 - centTrk) < 5.0 && centV0 <= 80 && centV0 > 0){
 392
 393      if ((centV0 > 0) && (centV0 <= 5.0))
 394        centrCode = 0;
 395      else if ((centV0 > 5.0) && (centV0 <= 10.0))
 396        centrCode = 1;
 397      else if ((centV0 > 10.0) && (centV0 <= 20.0))
 398        centrCode = 2;
 399      else if ((centV0 > 20.0) && (centV0 <= 30.0))
 400        centrCode = 3;
 401      else if ((centV0 > 30.0) && (centV0 <= 40.0))
 402        centrCode = 4;
 403      else if ((centV0 > 40.0) && (centV0 <= 50.0))
 404        centrCode = 5;
 405      else if ((centV0 > 50.0) && (centV0 <= 60.0))
 406        centrCode = 6;
 407      else if ((centV0 > 60.0) && (centV0 <= 70.0))
 408        centrCode = 7;
 409      else if ((centV0 > 70.0) && (centV0 <= 80.0))
 410        centrCode = 8;
 411   }
 412
 413   return centrCode;
 414
 415 }
 416
 417 //______________________________________________________
 418 void AliSpectraAODEventCuts::PrintCuts()
 419 {
 420   // print info about event cuts
 421   cout << "Event Stats" << endl;
 422   cout << " > Trigger Selection: " << fSelectBit << endl;
 423   cout << " > Centrality estimator: " << fCentralityMethod << endl;
 424   cout << " > Number of accepted events: " << NumberOfEvents() << endl;
 425   cout << " > Number of processed events: " << NumberOfProcessedEvents() << endl;
 426   cout << " > Number of PhysSel events: " <<  NumberOfPhysSelEvents()  << endl;
 427   cout << " > Vertex out of range: " << ((TH1I*)fOutput->FindObject("fHistoCuts"))->GetBinContent(kVtxRange + 1) << endl;
 428   cout << " > Events cut by centrality: " << ((TH1I*)fOutput->FindObject("fHistoCuts"))->GetBinContent(kVtxCentral + 1) << endl;
 429   cout << " > Events without vertex: " << ((TH1I*)fOutput->FindObject("fHistoCuts"))->GetBinContent(kVtxNoEvent + 1) << endl;
 430   cout << " > QVector cut: " << ((TH1I*)fOutput->FindObject("fHistoCuts"))->GetBinContent(kQVector + 1) << endl;
 431   cout << " > Track type used for the QVector calculation: " << fTrackBits << endl;
 432   cout << " > QRange: [" << fQVectorCutMin <<"," <<fQVectorCutMax<<"]"<< endl;
 433   cout << " > Vertex: [" << fVertexCutMin <<"," <<fVertexCutMax<<"]"<< endl;
 434   cout << " > Multiplicity: [" << fMultiplicityCutMin <<"," <<fMultiplicityCutMax<<"]"<< endl;
 435   cout << " > Centrality: [" << fCentralityCutMin <<"," <<fCentralityCutMax<<"]"<< endl;
 436 }
 437
 438 //_____________________________________________________________________________
 439 Bool_t AliSpectraAODEventCuts::OpenInfoCalbration(Int_t run)
 440 {
 441
 442   AliOADBContainer* cont = (AliOADBContainer*) fCalib->FindObject("hMultV0BefCorr");
 443   if(!cont){
 444     printf("OADB object hMultV0BefCorr is not available in the file\n");
 445     return 0;
 446   }
 447
 448   if(!(cont->GetObject(run))){
 449     printf("OADB object hMultV0BefCorr is not available for run %i\n",run);
 450     return 0;
 451   }
 452   fMultV0 = ((TH2F*) cont->GetObject(run))->ProfileX();
 453
 454   TF1* fpolc1 = new TF1("fpolc1","pol0", 0, 7);
 455   fMultV0->Fit(fpolc1, "RN");
 456   fV0Cpol1 = fpolc1->GetParameter(0);
 457
 458   TF1* fpolc2 = new TF1("fpolc2","pol0", 8, 15);
 459   fMultV0->Fit(fpolc2, "RN");
 460   fV0Cpol2 = fpolc2->GetParameter(0);
 461
 462   TF1* fpolc3 = new TF1("fpolc3","pol0", 16, 23);
 463   fMultV0->Fit(fpolc3, "RN");
 464   fV0Cpol3 = fpolc3->GetParameter(0);
 465
 466   TF1* fpolc4 = new TF1("fpolc4","pol0", 24, 31);
 467   fMultV0->Fit(fpolc4, "RN");
 468   fV0Cpol4 = fpolc4->GetParameter(0);
 469
 470   TF1* fpola1 = new TF1("fpola1","pol0", 32, 39);
 471   fMultV0->Fit(fpola1, "RN");
 472   fV0Apol1 = fpola1->GetParameter(0);
 473
 474   TF1* fpola2 = new TF1("fpola2","pol0", 40, 47);
 475   fMultV0->Fit(fpola2, "RN");
 476   fV0Apol2 = fpola2->GetParameter(0);
 477
 478   TF1* fpola3 = new TF1("fpola3","pol0", 48, 55);
 479   fMultV0->Fit(fpola3, "RN");
 480   fV0Apol3 = fpola3->GetParameter(0);
 481
 482   TF1* fpola4 = new TF1("fpola4","pol0", 56, 63);
 483   fMultV0->Fit(fpola4, "RN");
 484   fV0Apol4 = fpola4->GetParameter(0);
 485
 486   for(Int_t i=0; i < 10; i++){
 487
 488     char nameQxa2[100];
 489     snprintf(nameQxa2,100, "hQxa2m_%i", i);
 490
 491     char nameQya2[100];
 492     snprintf(nameQya2,100, "hQya2m_%i", i);
 493
 494     char nameQxc2[100];
 495     snprintf(nameQxc2,100, "hQxc2m_%i", i);
 496
 497     char nameQyc2[100];
 498     snprintf(nameQyc2,100, "hQyc2m_%i", i);
 499
 500     AliOADBContainer* contQxa2 = (AliOADBContainer*) fCalib->FindObject(nameQxa2);
 501     if(!contQxa2){
 502       printf("OADB object %s is not available in the file\n", nameQxa2);
 503       return 0;
 504     }
 505
 506     if(!(contQxa2->GetObject(run))){
 507       printf("OADB object %s is not available for run %i\n", nameQxa2, run);
 508       return 0;
 509     }
 510
 511     fMeanQxa2[i] = ((TH1F*) contQxa2->GetObject(run))->GetMean();
 512
 513
 514     AliOADBContainer* contQya2 = (AliOADBContainer*) fCalib->FindObject(nameQya2);
 515     if(!contQya2){
 516       printf("OADB object %s is not available in the file\n", nameQya2);
 517       return 0;
 518     }
 519
 520     if(!(contQya2->GetObject(run))){
 521       printf("OADB object %s is not available for run %i\n", nameQya2, run);
 522       return 0;
 523     }
 524
 525     fMeanQya2[i] = ((TH1F*) contQya2->GetObject(run))->GetMean();
 526
 527
 528     AliOADBContainer* contQxc2 = (AliOADBContainer*) fCalib->FindObject(nameQxc2);
 529     if(!contQxc2){
 530       printf("OADB object %s is not available in the file\n", nameQxc2);
 531       return 0;
 532     }
 533
 534     if(!(contQxc2->GetObject(run))){
 535       printf("OADB object %s is not available for run %i\n", nameQxc2, run);
 536       return 0;
 537     }
 538
 539     fMeanQxc2[i] = ((TH1F*) contQxc2->GetObject(run))->GetMean();
 540
 541
 542     AliOADBContainer* contQyc2 = (AliOADBContainer*) fCalib->FindObject(nameQyc2);
 543     if(!contQyc2){
 544       printf("OADB object %s is not available in the file\n", nameQyc2);
 545       return 0;
 546     }
 547
 548     if(!(contQyc2->GetObject(run))){
 549       printf("OADB object %s is not available for run %i\n", nameQyc2, run);
 550       return 0;
 551     }
 552
 553     fMeanQyc2[i] = ((TH1F*) contQyc2->GetObject(run))->GetMean();
 554
 555   }
 556   return 1;
 557 }
 558
 559 //______________________________________________________
 560
 561
 562 Long64_t AliSpectraAODEventCuts::Merge(TCollection* list)
 563 {
 564   // Merge a list of AliSpectraAODEventCuts objects with this.
 565   // Returns the number of merged objects (including this).
 566
 567   AliInfo("Merging");
 568
 569   if (!list)
 570     return 0;
 571
 572   if (list->IsEmpty())
 573     return 1;
 574
 575   TIterator* iter = list->MakeIterator();
 576   TObject* obj;
 577
 578   // collections of all histograms
 579   TList collections;
 580
 581   Int_t count = 0;
 582
 583   while ((obj = iter->Next())) {
 584     AliSpectraAODEventCuts* entry = dynamic_cast<AliSpectraAODEventCuts*> (obj);
 585     if (entry == 0)
 586       continue;
 587
 588     TList * l = entry->GetOutputList();
 589     collections.Add(l);
 590     count++;
 591   }
 592
 593   fOutput->Merge(&collections);
 594
 595   delete iter;
 596
 597   return count+1;
 598 }
 599
 600 //______________________________________________________
 601 Double_t AliSpectraAODEventCuts::GetQvecPercentile(Int_t v0side){
 602
 603   //check Qvec and Centrality consistency
 604   if(fCent>90.) return -999.;
 605   if(v0side==0/*V0A*/){ if(fqV0A== -999.) return -999.; }
 606   if(v0side==1/*V0C*/){ if(fqV0C== -999.) return -999.; }
 607
 608   fQvecIntegral = 0x0;
 609
 610   if(v0side==0/*V0A*/){ fQvecIntegral = (TH2D*)fQvecIntList->FindObject("VZEROA"); }
 611   if(v0side==1/*V0C*/){ fQvecIntegral = (TH2D*)fQvecIntList->FindObject("VZEROC"); }
 612
 613   Double_t ic = (Int_t)fCent; //fQvecIntegral: 1% centrality bin
 614
 615   TH1D *h1D = (TH1D*)fQvecIntegral->ProjectionY("h1D",ic+1,ic+1);
 616
 617   TSpline *spline = 0x0;
 618
 619   if(v0side==0/*V0A*/){
 620     if( CheckSplineArray(fSplineArrayV0A) ) {
 621       spline = (TSpline*)fSplineArrayV0A->At(ic);
 622       //cout<<"Qvec V0A - ic: "<<ic<<" - Found TSpline..."<<endl;
 623     } else {
 624       spline = new TSpline3(h1D,"sp3");
 625       fSplineArrayV0A->AddAtAndExpand(spline,ic);
 626       //cout<<"Qvec V0A - ic: "<<ic<<" - TSpline not found. Creating..."<<endl;
 627     }
 628   }
 629   else if(v0side==1/*V0C*/){
 630     if( CheckSplineArray(fSplineArrayV0C) ) {
 631       spline = (TSpline*)fSplineArrayV0C->At(ic);
 632     } else {
 633       spline = new TSpline3(h1D,"sp3");
 634       fSplineArrayV0C->AddAtAndExpand(spline,ic);
 635     }
 636   }
 637
 638   Double_t percentile =  -999.;
 639   if(v0side==0/*V0A*/){ percentile = 100*spline->Eval(fqV0A); }
 640   if(v0side==1/*V0C*/){ percentile = 100*spline->Eval(fqV0C); }
 641
 642   if(percentile>100. || percentile<0.) return -999.;
 643
 644   return percentile;
 645 }
 646
 647 //______________________________________________________
 648 Bool_t AliSpectraAODEventCuts::CheckSplineArray(TObjArray * splarr){
 649   //check TSpline array consistency
 650
 651   Int_t n = (Int_t)fCent;
 652
 653   if(splarr->GetSize()<n) return kFALSE;
 654
 655   if(!splarr->At(n)) return kFALSE;
 656
 657   return kTRUE;
 658
 659 }