PWG2/FLOW/AliFlowCommon/AliFlowAnalysisWithScalarProduct.cxx

   1 /*************************************************************************
   2 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15
  16 #define AliFlowAnalysisWithScalarProduct_cxx
  17
  18 #include "Riostream.h"  //needed as include
  19 #include "TFile.h"      //needed as include
  20 #include "TList.h"
  21 #include "TMath.h"
  22 #include "TProfile.h"
  23 #include "TVector2.h"
  24
  25 class TH1F;
  26
  27 #include "AliFlowCommonConstants.h"    //needed as include
  28 #include "AliFlowEventSimple.h"
  29 #include "AliFlowTrackSimple.h"
  30 #include "AliFlowCommonHist.h"
  31 #include "AliFlowCommonHistResults.h"
  32 #include "AliFlowAnalysisWithScalarProduct.h"
  33
  34 class AliFlowVector;
  35
  36 // AliFlowAnalysisWithScalarProduct:
  37 // Description:
  38 // Maker to analyze Flow with the Scalar product method.
  39 //
  40 // author: N. van der Kolk (kolk@nikhef.nl)
  41
  42 ClassImp(AliFlowAnalysisWithScalarProduct)
  43
  44   //-----------------------------------------------------------------------
  45
  46  AliFlowAnalysisWithScalarProduct::AliFlowAnalysisWithScalarProduct():
  47    fEventNumber(0),
  48    fDebug(kFALSE),
  49    fWeightsList(NULL),
  50    fUsePhiWeights(kFALSE),
  51    fPhiWeights(NULL),
  52    fHistList(NULL),
  53    fHistProUQetaRP(NULL),
  54    fHistProUQetaPOI(NULL),
  55    fHistProUQPtRP(NULL),
  56    fHistProUQPtPOI(NULL),
  57    fHistProQaQb(NULL),
  58    fHistProM(NULL),
  59    fHistM(NULL),
  60    fCommonHists(NULL),
  61    fCommonHistsRes(NULL)
  62 {
  63   // Constructor.
  64   fWeightsList = new TList();
  65   fHistList = new TList();
  66 }
  67  //-----------------------------------------------------------------------
  68
  69
  70  AliFlowAnalysisWithScalarProduct::~AliFlowAnalysisWithScalarProduct()
  71  {
  72    //destructor
  73    delete fWeightsList;
  74    delete fHistList;
  75  }
  76
  77
  78 //-----------------------------------------------------------------------
  79
  80 void AliFlowAnalysisWithScalarProduct::WriteHistograms(TString* outputFileName)
  81 {
  82  //store the final results in output .root file
  83
  84   TFile *output = new TFile(outputFileName->Data(),"RECREATE");
  85   //output->WriteObject(fHistList, "cobjSP","SingleKey");
  86   fHistList->SetName("cobjSP");
  87   fHistList->SetOwner(kTRUE);
  88   fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
  89   delete output;
  90 }
  91
  92 //-----------------------------------------------------------------------
  93
  94 void AliFlowAnalysisWithScalarProduct::WriteHistograms(TString outputFileName)
  95 {
  96  //store the final results in output .root file
  97
  98   TFile *output = new TFile(outputFileName.Data(),"RECREATE");
  99   //output->WriteObject(fHistList, "cobjSP","SingleKey");
 100   fHistList->SetName("cobjSP");
 101   fHistList->SetOwner(kTRUE);
 102   fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
 103   delete output;
 104 }
 105
 106 //-----------------------------------------------------------------------
 107
 108 void AliFlowAnalysisWithScalarProduct::WriteHistograms(TDirectoryFile *outputFileName)
 109 {
 110  //store the final results in output .root file
 111  fHistList->SetName("cobjSP");
 112  fHistList->SetOwner(kTRUE);
 113  outputFileName->Add(fHistList);
 114  outputFileName->Write(outputFileName->GetName(), TObject::kSingleKey);
 115 }
 116
 117 //-----------------------------------------------------------------------
 118 void AliFlowAnalysisWithScalarProduct::Init() {
 119
 120   //Define all histograms
 121   cout<<"---Analysis with the Scalar Product Method--- Init"<<endl;
 122
 123   Int_t iNbinsPt   = AliFlowCommonConstants::GetMaster()->GetNbinsPt();
 124   Double_t dPtMin  = AliFlowCommonConstants::GetMaster()->GetPtMin();
 125   Double_t dPtMax  = AliFlowCommonConstants::GetMaster()->GetPtMax();
 126   Int_t iNbinsEta  = AliFlowCommonConstants::GetMaster()->GetNbinsEta();
 127   Double_t dEtaMin = AliFlowCommonConstants::GetMaster()->GetEtaMin();
 128   Double_t dEtaMax = AliFlowCommonConstants::GetMaster()->GetEtaMax();
 129
 130   fHistProUQetaRP = new TProfile("Flow_UQetaRP_SP","Flow_UQetaRP_SP",iNbinsEta,dEtaMin,dEtaMax);
 131   fHistProUQetaRP->SetXTitle("{eta}");
 132   fHistProUQetaRP->SetYTitle("<uQ>");
 133   fHistList->Add(fHistProUQetaRP);
 134
 135   fHistProUQetaPOI = new TProfile("Flow_UQetaPOI_SP","Flow_UQetaPOI_SP",iNbinsEta,dEtaMin,dEtaMax);
 136   fHistProUQetaPOI->SetXTitle("{eta}");
 137   fHistProUQetaPOI->SetYTitle("<uQ>");
 138   fHistList->Add(fHistProUQetaPOI);
 139
 140   fHistProUQPtRP = new TProfile("Flow_UQPtRP_SP","Flow_UQPtRP_SP",iNbinsPt,dPtMin,dPtMax);
 141   fHistProUQPtRP->SetXTitle("p_t (GeV)");
 142   fHistProUQPtRP->SetYTitle("<uQ>");
 143   fHistList->Add(fHistProUQPtRP);
 144
 145   fHistProUQPtPOI = new TProfile("Flow_UQPtPOI_SP","Flow_UQPtPOI_SP",iNbinsPt,dPtMin,dPtMax);
 146   fHistProUQPtPOI->SetXTitle("p_t (GeV)");
 147   fHistProUQPtPOI->SetYTitle("<uQ>");
 148   fHistList->Add(fHistProUQPtPOI);
 149
 150   fHistProQaQb = new TProfile("Flow_QaQb_SP","Flow_QaQb_SP", 1, -0.5, 0.5);
 151   fHistProQaQb->SetYTitle("<QaQb>");
 152   fHistList->Add(fHistProQaQb);
 153
 154   fHistProM = new TProfile("Flow_M_SP","Flow_M_SP",2,0.5, 2.5);
 155   fHistProM -> SetYTitle("<*>");
 156   fHistProM -> SetXTitle("<M-1>, <Ma*Mb>");
 157   fHistList->Add(fHistProM);
 158
 159   fHistM = new TH1D("Flow_Msum_SP","Flow_Msum_SP",2,0.5, 2.5);
 160   fHistM -> SetYTitle("sum (*)");
 161   fHistM -> SetXTitle("sum (M-1), sum (Ma*Mb)");
 162   fHistList->Add(fHistM);
 163
 164   fCommonHists = new AliFlowCommonHist("AliFlowCommonHistSP");
 165   fHistList->Add(fCommonHists);
 166   fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsSP");
 167   fHistList->Add(fCommonHistsRes);
 168
 169   //weights
 170   if(fUsePhiWeights) {
 171     if(!fWeightsList) {
 172       cout<<"WARNING: fWeightsList is NULL in the Scalar Product method."<<endl;
 173       exit(0);
 174     }
 175     if(fWeightsList->FindObject("phi_weights"))  {
 176       fPhiWeights = dynamic_cast<TH1F*>(fWeightsList->FindObject("phi_weights"));
 177       fHistList->Add(fPhiWeights);
 178     } else {
 179       cout<<"WARNING: histogram with phi weights is not accessible in Scalar Product"<<endl;
 180       exit(0);
 181     }
 182   } // end of if(fUsePhiWeights)
 183
 184   fEventNumber = 0;  //set number of events to zero
 185 }
 186
 187 //-----------------------------------------------------------------------
 188
 189 void AliFlowAnalysisWithScalarProduct::Make(AliFlowEventSimple* anEvent) {
 190
 191   //Fill histogram
 192   if (anEvent) {
 193
 194     //fill control histograms
 195     fCommonHists->FillControlHistograms(anEvent);
 196
 197     //get Q vectors for the eta-subevents
 198     AliFlowVector* vQarray = new AliFlowVector[2];
 199     if (fUsePhiWeights) {
 200       anEvent->Get2Qsub(vQarray,2,fWeightsList,kTRUE);
 201     } else {
 202       anEvent->Get2Qsub(vQarray);
 203     }
 204     AliFlowVector vQa = vQarray[0];
 205     AliFlowVector vQb = vQarray[1];
 206     //get total Q vector
 207     AliFlowVector vQ = vQa + vQb;
 208
 209     //fill the multiplicity histograms for the prefactor
 210     Double_t dMmin1 = vQ.GetMult()-1; //M-1
 211     Double_t dMaMb = vQa.GetMult()*vQb.GetMult();     //Ma*Mb
 212     fHistM -> Fill(1,dMmin1);//sum of (M-1)
 213     fHistM -> Fill(2,vQa.GetMult()*vQb.GetMult()); //sum of (Ma*Mb)
 214     fHistProM -> Fill(1,dMmin1);      //<M-1>
 215     fHistProM -> Fill(2,vQa.GetMult()*vQb.GetMult()); //<Ma*Mb>
 216     //scalar product of the two subevents
 217     Double_t dQaQb = (vQa*vQb)*dMaMb;
 218     fHistProQaQb -> Fill(0.,dQaQb);
 219
 220     //loop over the tracks of the event
 221     AliFlowTrackSimple*   pTrack = NULL;
 222     Int_t iNumberOfTracks = anEvent->NumberOfTracks();
 223     for (Int_t i=0;i<iNumberOfTracks;i++)
 224       {
 225         pTrack = anEvent->GetTrack(i) ;
 226         if (pTrack){
 227           Double_t dPhi = pTrack->Phi();
 228           //set default weight to 1
 229           Double_t dW = 1.;
 230           //phi weight of pTrack
 231           if(fUsePhiWeights && fPhiWeights) {
 232             Int_t iNBinsPhi = fPhiWeights->GetNbinsX();
 233             dW = fPhiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*iNBinsPhi/TMath::TwoPi())));  //bin = 1 + value*nbins/range
 234           }                                                                                           //TMath::Floor rounds to the lower integer
 235           //calculate vU
 236           TVector2 vU;
 237           Double_t dUX = TMath::Cos(2*dPhi);
 238           Double_t dUY = TMath::Sin(2*dPhi);
 239           vU.Set(dUX,dUY);
 240           Double_t dModulus = vU.Mod();
 241           if (dModulus!=0.) vU.Set(dUX/dModulus,dUY/dModulus);  // make length 1
 242           else cerr<<"dModulus is zero!"<<endl;
 243
 244           TVector2 vQm = vQ;
 245           //subtract particle from the flowvector if used to define it
 246           if (pTrack->InRPSelection()) {
 247             if (pTrack->InSubevent(0) || pTrack->InSubevent(1)) {
 248               Double_t dQmX = vQm.X() - dW*dUX;
 249               Double_t dQmY = vQm.Y() - dW*dUY;
 250               vQm.Set(dQmX,dQmY);
 251             }
 252           }
 253
 254           //dUQ = scalar product of vU and vQm
 255           Double_t dUQ = (vU * vQm)*dMmin1;
 256           Double_t dPt = pTrack->Pt();
 257           Double_t dEta = pTrack->Eta();
 258           //fill the profile histograms
 259           if (pTrack->InRPSelection()) {
 260             fHistProUQetaRP -> Fill(dEta,dUQ);
 261             fHistProUQPtRP -> Fill(dPt,dUQ);
 262           }
 263           if (pTrack->InPOISelection()) {
 264             fHistProUQetaPOI -> Fill(dEta,dUQ);
 265             fHistProUQPtPOI -> Fill(dPt,dUQ);
 266           }
 267         }//track selected
 268       }//loop over tracks
 269
 270     fEventNumber++;
 271     //    cout<<"@@@@@ "<<fEventNumber<<" events processed"<<endl;
 272     delete [] vQarray;
 273   }
 274 }
 275
 276   //--------------------------------------------------------------------
 277 void AliFlowAnalysisWithScalarProduct::GetOutputHistograms(TList *outputListHistos){
 278
 279   //get pointers to all output histograms (called before Finish())
 280   if (outputListHistos) {
 281   //Get the common histograms from the output list
 282     AliFlowCommonHist *pCommonHist = dynamic_cast<AliFlowCommonHist*>
 283       (outputListHistos->FindObject("AliFlowCommonHistSP"));
 284     AliFlowCommonHistResults *pCommonHistResults = dynamic_cast<AliFlowCommonHistResults*>
 285       (outputListHistos->FindObject("AliFlowCommonHistResultsSP"));
 286     TProfile* pHistProQaQb     = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_QaQb_SP"));
 287     TProfile* pHistProM        = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_M_SP"));
 288     TProfile* pHistProUQetaRP  = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQetaRP_SP"));
 289     TProfile* pHistProUQetaPOI = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQetaPOI_SP"));
 290     TProfile* pHistProUQPtRP   = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQPtRP_SP"));
 291     TProfile* pHistProUQPtPOI  = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQPtPOI_SP"));
 292     if (pCommonHist && pCommonHistResults && pHistProQaQb && pHistProM &&
 293         pHistProUQetaRP && pHistProUQetaPOI && pHistProUQPtRP && pHistProUQPtPOI) {
 294       this -> SetCommonHists(pCommonHist);
 295       this -> SetCommonHistsRes(pCommonHistResults);
 296       this -> SetHistProQaQb(pHistProQaQb);
 297       this -> SetHistProM(pHistProM);
 298       this -> SetHistProUQetaRP(pHistProUQetaRP);
 299       this -> SetHistProUQetaPOI(pHistProUQetaPOI);
 300       this -> SetHistProUQPtRP(pHistProUQPtRP);
 301       this -> SetHistProUQPtPOI(pHistProUQPtPOI);
 302       }
 303    }
 304 }
 305
 306 //--------------------------------------------------------------------
 307 void AliFlowAnalysisWithScalarProduct::Finish() {
 308
 309   //calculate flow and fill the AliFlowCommonHistResults
 310   if (fDebug) cout<<"AliFlowAnalysisWithScalarProduct::Finish()"<<endl;
 311
 312   cout<<"*************************************"<<endl;
 313   cout<<"*************************************"<<endl;
 314   cout<<"      Integrated flow from           "<<endl;
 315   cout<<"         Scalar product              "<<endl;
 316   cout<<endl;
 317
 318   Int_t iNbinsPt  = AliFlowCommonConstants::GetMaster()->GetNbinsPt();
 319   Int_t iNbinsEta = AliFlowCommonConstants::GetMaster()->GetNbinsEta();
 320
 321   Double_t dMmin1    = fHistProM->GetBinContent(1);  //average over M-1
 322   Double_t dMmin1Err = fHistProM->GetBinError(1);    //error on average over M-1
 323   Double_t dMaMb     = fHistProM->GetBinContent(2);  //average over Ma*Mb
 324   Double_t dMaMbErr  = fHistProM->GetBinError(2);    //error on average over Ma*Mb
 325
 326   //Double_t dSumMmin1 = fHistM->GetBinContent(1);  //sum over (M-1)
 327   //Double_t dSumMaMb  = fHistM->GetBinContent(2);  //sum over (Ma*Mb)
 328
 329   Double_t dMcorrection = 0.;     //correction factor for Ma != Mb
 330   Double_t dMcorrectionErr = 0.;
 331   Double_t dMcorrectionErrRel = 0.;
 332   Double_t dMcorrectionErrRel2 = 0.;
 333
 334   if (dMaMb != 0. && dMmin1 != 0.) {
 335     dMcorrection    = dMmin1/(TMath::Sqrt(dMaMb));
 336     dMcorrectionErr = dMcorrection*(dMmin1Err/dMmin1 + dMaMbErr/(2*dMaMb));
 337     dMcorrectionErrRel = dMcorrectionErr/dMcorrection;
 338     dMcorrectionErrRel2 = dMcorrectionErrRel*dMcorrectionErrRel;
 339   }
 340
 341   Double_t dQaQbAv  = TMath::Abs(fHistProQaQb->GetBinContent(1)); //average over events //TEST TAKE ABS
 342   dQaQbAv /= dMaMb;  //normalize for weighted average
 343   Double_t dQaQbErr = fHistProQaQb->GetBinError(1);
 344   dQaQbErr /= dMaMb;  //normalize for weighted average
 345   Double_t dQaQbErrRel = 0.;
 346   if (dQaQbAv != 0.) {
 347     dQaQbErrRel = dQaQbErr/dQaQbAv; }
 348   Double_t dQaQbErrRel2 = dQaQbErrRel*dQaQbErrRel;
 349
 350   if (dQaQbAv <= 0.){
 351     //set v to -0
 352     fCommonHistsRes->FillIntegratedFlowRP(-0.,0.);
 353     fCommonHistsRes->FillIntegratedFlow(-0.,0.);
 354     cout<<"dV(RP) = -0. +- 0."<<endl;
 355     fCommonHistsRes->FillIntegratedFlowPOI(-0.,0.);
 356     cout<<"dV(POI) = -0. +- 0."<<endl;
 357   } else {
 358   Double_t dQaQbSqrt = TMath::Sqrt(dQaQbAv);  //DOES NOT WORK IF dQaQbAv IS NEGATIVE
 359     if (dMaMb>0.) { dQaQbSqrt *= dMcorrection; }
 360     else { dQaQbSqrt = 0.; }
 361     Double_t dQaQbSqrtErrRel2 = dMcorrectionErrRel2 + (1/4)*dQaQbErrRel2;
 362
 363     //v as a function of eta for RP selection
 364     for(Int_t b=0;b<iNbinsEta;b++) {
 365       Double_t duQpro = fHistProUQetaRP->GetBinContent(b);
 366       duQpro /= dMmin1;  //normalize for weighted average
 367       Double_t duQerr = fHistProUQetaRP->GetBinError(b); //copy error for now
 368       duQerr /= dMmin1;  //normalize for weighted average
 369       Double_t duQerrRel = 0.;
 370       if (duQpro != 0.) {duQerrRel = duQerr/duQpro;}
 371       Double_t duQerrRel2 = duQerrRel*duQerrRel;
 372
 373       Double_t dv2pro     = 0.;
 374       if (dQaQbSqrt!=0.) { dv2pro = duQpro/dQaQbSqrt; }
 375       Double_t dv2errRel2 = duQerrRel2 + dQaQbSqrtErrRel2;
 376       Double_t dv2errRel  = 0.;
 377       if (dv2errRel2>0.) { dv2errRel  = TMath::Sqrt(dv2errRel2); }
 378       Double_t dv2err     = dv2pro*dv2errRel;
 379       //fill TH1D
 380       fCommonHistsRes->FillDifferentialFlowEtaRP(b, dv2pro, dv2err);
 381     } //loop over bins b
 382
 383     //v as a function of eta for POI selection
 384     for(Int_t b=0;b<iNbinsEta;b++) {
 385       Double_t duQpro = fHistProUQetaPOI->GetBinContent(b);
 386       duQpro /= dMmin1;  //normalize for weighted average
 387       Double_t duQerr = fHistProUQetaPOI->GetBinError(b); //copy error for now
 388       duQerr /= dMmin1;  //normalize for weighted average
 389       Double_t duQerrRel = 0.;
 390       if (duQpro != 0.) {duQerrRel = duQerr/duQpro;}
 391       Double_t duQerrRel2 = duQerrRel*duQerrRel;
 392
 393       Double_t dv2pro     = 0.;
 394       if (dQaQbSqrt!=0.) { dv2pro = duQpro/dQaQbSqrt; }
 395       Double_t dv2errRel2 = duQerrRel2 + dQaQbSqrtErrRel2;
 396       Double_t dv2errRel  = 0.;
 397       if (dv2errRel2>0.) { dv2errRel  = TMath::Sqrt(dv2errRel2); }
 398       Double_t dv2err     = dv2pro*dv2errRel;
 399       //fill TH1D
 400       fCommonHistsRes->FillDifferentialFlowEtaPOI(b, dv2pro, dv2err);
 401     } //loop over bins b
 402
 403     //v as a function of Pt for RP selection
 404     TH1F* fHistPtRP = fCommonHists->GetHistPtRP(); //for calculating integrated flow
 405     Double_t dVRP = 0.;
 406     Double_t dSum = 0.;
 407     Double_t dErrV =0.;
 408
 409     for(Int_t b=0;b<iNbinsPt;b++) {
 410       Double_t duQpro = fHistProUQPtRP->GetBinContent(b);
 411       duQpro /= dMmin1;  //normalize for weighted average
 412       Double_t duQerr = fHistProUQPtRP->GetBinError(b); //copy error for now
 413       duQerr /= dMmin1;  //normalize for weighted average
 414       Double_t duQerrRel = 0.;
 415       if (duQpro != 0.) {duQerrRel = duQerr/duQpro;}
 416       Double_t duQerrRel2 = duQerrRel*duQerrRel;
 417
 418       Double_t dv2pro     = 0.;
 419       if (dQaQbSqrt!=0.) { dv2pro = duQpro/dQaQbSqrt; }
 420       Double_t dv2errRel2 = duQerrRel2 + dQaQbSqrtErrRel2;
 421       Double_t dv2errRel  = 0.;
 422       if (dv2errRel2>0.) { dv2errRel  = TMath::Sqrt(dv2errRel2); }
 423       Double_t dv2err     = dv2pro*dv2errRel;
 424       //fill TH1D
 425       fCommonHistsRes->FillDifferentialFlowPtRP(b, dv2pro, dv2err);
 426       //calculate integrated flow for RP selection
 427       if (fHistPtRP){
 428         Double_t dYieldPt = fHistPtRP->GetBinContent(b);
 429         dVRP += dv2pro*dYieldPt;
 430         dSum +=dYieldPt;
 431         dErrV += dYieldPt*dYieldPt*dv2err*dv2err;
 432       } else { cout<<"fHistPtRP is NULL"<<endl; }
 433     } //loop over bins b
 434
 435     if (dSum != 0.) {
 436       dVRP /= dSum; //the pt distribution should be normalised
 437       dErrV /= (dSum*dSum);
 438       dErrV = TMath::Sqrt(dErrV);
 439     }
 440     fCommonHistsRes->FillIntegratedFlowRP(dVRP,dErrV);
 441     fCommonHistsRes->FillIntegratedFlow(dVRP,dErrV);
 442
 443     cout<<"dV(RP) = "<<dVRP<<" +- "<<dErrV<<endl;
 444
 445     //v as a function of Pt for POI selection
 446     TH1F* fHistPtPOI = fCommonHists->GetHistPtPOI(); //for calculating integrated flow
 447     Double_t dVPOI = 0.;
 448     dSum = 0.;
 449     dErrV =0.;
 450
 451     for(Int_t b=0;b<iNbinsPt;b++) {
 452       Double_t duQpro = fHistProUQPtPOI->GetBinContent(b);
 453       duQpro /= dMmin1;  //normalize for weighted average
 454       Double_t duQerr = fHistProUQPtPOI->GetBinError(b); //copy error for now
 455       duQerr /= dMmin1;  //normalize for weighted average
 456       Double_t duQerrRel = 0.;
 457       if (duQpro != 0.) {duQerrRel = duQerr/duQpro;}
 458       Double_t duQerrRel2 = duQerrRel*duQerrRel;
 459
 460       Double_t dv2pro     = 0.;
 461       if (dQaQbSqrt!=0.) { dv2pro = duQpro/dQaQbSqrt; }
 462       Double_t dv2errRel2 = duQerrRel2 + dQaQbSqrtErrRel2;
 463       Double_t dv2errRel  = 0.;
 464       if (dv2errRel2>0.) { dv2errRel  = TMath::Sqrt(dv2errRel2); }
 465       Double_t dv2err     = dv2pro*dv2errRel;
 466       //fill TH1D
 467       fCommonHistsRes->FillDifferentialFlowPtPOI(b, dv2pro, dv2err);
 468
 469       //calculate integrated flow for POI selection
 470       if (fHistPtPOI){
 471         Double_t dYieldPt = fHistPtPOI->GetBinContent(b);
 472         dVPOI += dv2pro*dYieldPt;
 473         dSum +=dYieldPt;
 474         dErrV += dYieldPt*dYieldPt*dv2err*dv2err;
 475       } else { cout<<"fHistPtPOI is NULL"<<endl; }
 476     } //loop over bins b
 477
 478     if (dSum != 0.) {
 479       dVPOI /= dSum; //the pt distribution should be normalised
 480       dErrV /= (dSum*dSum);
 481       dErrV = TMath::Sqrt(dErrV);
 482     }
 483     fCommonHistsRes->FillIntegratedFlowPOI(dVPOI,dErrV);
 484
 485     cout<<"dV(POI) = "<<dVPOI<<" +- "<<dErrV<<endl;
 486     }
 487   cout<<endl;
 488   cout<<"*************************************"<<endl;
 489   cout<<"*************************************"<<endl;
 490
 491   //cout<<".....finished"<<endl;
 492  }
 493
 494