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"
  19 #include "TFile.h"
  20 #include "TList.h"
  21 #include "TMath.h"
  22 #include "TProfile.h"
  23 #include "TVector2.h"
  24 #include "TH1D.h"
  25 #include "TH2D.h"
  26
  27 #include "AliFlowCommonConstants.h"
  28 #include "AliFlowEventSimple.h"
  29 #include "AliFlowVector.h"
  30 #include "AliFlowTrackSimple.h"
  31 #include "AliFlowCommonHist.h"
  32 #include "AliFlowCommonHistResults.h"
  33 #include "AliFlowAnalysisWithScalarProduct.h"
  34
  35 //////////////////////////////////////////////////////////////////////////////
  36 // AliFlowAnalysisWithScalarProduct:
  37 // Description:
  38 // Maker to analyze Flow with the Scalar product method.
  39 //
  40 // authors: N. van der Kolk (kolk@nikhef.nl), A. Bilandzic (anteb@nikhef.nl)
  41 //////////////////////////////////////////////////////////////////////////////
  42
  43 ClassImp(AliFlowAnalysisWithScalarProduct)
  44
  45   //-----------------------------------------------------------------------
  46   AliFlowAnalysisWithScalarProduct::AliFlowAnalysisWithScalarProduct():
  47    fEventNumber(0),
  48    fDebug(kFALSE),
  49    fApplyCorrectionForNUA(kFALSE),
  50    fRelDiffMsub(1.),
  51    fWeightsList(NULL),
  52    fUsePhiWeights(kFALSE),
  53    fPhiWeightsSub0(NULL),
  54    fPhiWeightsSub1(NULL),
  55    fHistProFlags(NULL),
  56    fHistProUQetaRP(NULL),
  57    fHistProUQetaPOI(NULL),
  58    fHistProUQetaAllEventsPOI(NULL),
  59    fHistProUQPtRP(NULL),
  60    fHistProUQPtPOI(NULL),
  61    fHistProUQPtAllEventsPOI(NULL),
  62    fHistProQNorm(NULL),
  63    fHistProQaQb(NULL),
  64    fHistProQaQbNorm(NULL),
  65    fHistProQaQbReImNorm(NULL),
  66    fHistProNonIsotropicTermsQ(NULL),
  67    fHistSumOfLinearWeights(NULL),
  68    fHistSumOfQuadraticWeights(NULL),
  69    fHistProUQQaQbPtRP(NULL),
  70    fHistProUQQaQbEtaRP(NULL),
  71    fHistProUQQaQbPtPOI(NULL),
  72    fHistProUQQaQbEtaPOI(NULL),
  73    fCommonHistsSP(NULL),
  74    fCommonHistsResSP(NULL),
  75    fCommonHistsmuQ(NULL),
  76    fHistQNorm(NULL),
  77    fHistQaQb(NULL),
  78    fHistQaQbNorm(NULL),
  79    fHistQNormvsQaQbNorm(NULL),
  80    fHistQaQbCos(NULL),
  81    fHistResolution(NULL),
  82    fHistQaNorm(NULL),
  83    fHistQaNormvsMa(NULL),
  84    fHistQbNorm(NULL),
  85    fHistQbNormvsMb(NULL),
  86    fHistMavsMb(NULL),
  87    fHistList(NULL)
  88 {
  89   // Constructor.
  90   fWeightsList = new TList();
  91   fHistList = new TList();
  92
  93   // Initialize arrays:
  94   for(Int_t i=0;i<3;i++)
  95   {
  96    fHistSumOfWeightsPtRP[i] = NULL;
  97    fHistSumOfWeightsEtaRP[i] = NULL;
  98    fHistSumOfWeightsPtPOI[i] = NULL;
  99    fHistSumOfWeightsEtaPOI[i] = NULL;
 100   }
 101   for(Int_t rp=0;rp<2;rp++)
 102   {
 103    for(Int_t pe=0;pe<2;pe++)
 104    {
 105     for(Int_t sc=0;sc<2;sc++)
 106     {
 107      fHistProNonIsotropicTermsU[rp][pe][sc] = NULL;
 108     }
 109    }
 110   }
 111 }
 112  //-----------------------------------------------------------------------
 113  AliFlowAnalysisWithScalarProduct::~AliFlowAnalysisWithScalarProduct()
 114  {
 115    //destructor
 116    delete fWeightsList;
 117    delete fHistList;
 118  }
 119
 120 //-----------------------------------------------------------------------
 121 void AliFlowAnalysisWithScalarProduct::WriteHistograms(TString* outputFileName)
 122 {
 123  //store the final results in output .root file
 124
 125   TFile *output = new TFile(outputFileName->Data(),"RECREATE");
 126   //output->WriteObject(fHistList, "cobjSP","SingleKey");
 127   fHistList->SetName("cobjSP");
 128   fHistList->SetOwner(kTRUE);
 129   fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
 130   delete output;
 131 }
 132
 133 //-----------------------------------------------------------------------
 134 void AliFlowAnalysisWithScalarProduct::WriteHistograms(TString outputFileName)
 135 {
 136  //store the final results in output .root file
 137
 138   TFile *output = new TFile(outputFileName.Data(),"RECREATE");
 139   //output->WriteObject(fHistList, "cobjSP","SingleKey");
 140   fHistList->SetName("cobjSP");
 141   fHistList->SetOwner(kTRUE);
 142   fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
 143   delete output;
 144 }
 145
 146 //-----------------------------------------------------------------------
 147 void AliFlowAnalysisWithScalarProduct::WriteHistograms(TDirectoryFile *outputFileName)
 148 {
 149  //store the final results in output .root file
 150  fHistList->SetName("cobjSP");
 151  fHistList->SetOwner(kTRUE);
 152  outputFileName->Add(fHistList);
 153  outputFileName->Write(outputFileName->GetName(), TObject::kSingleKey);
 154 }
 155
 156 //-----------------------------------------------------------------------
 157 void AliFlowAnalysisWithScalarProduct::Init() {
 158
 159   //Define all histograms
 160   cout<<"---Analysis with the Scalar Product Method--- Init"<<endl;
 161
 162   //save old value and prevent histograms from being added to directory
 163   //to avoid name clashes in case multiple analaysis objects are used
 164   //in an analysis
 165
 166   Bool_t oldHistAddStatus = TH1::AddDirectoryStatus();
 167   TH1::AddDirectory(kFALSE);
 168
 169   Int_t iNbinsPt   = AliFlowCommonConstants::GetMaster()->GetNbinsPt();
 170   Double_t dPtMin  = AliFlowCommonConstants::GetMaster()->GetPtMin();
 171   Double_t dPtMax  = AliFlowCommonConstants::GetMaster()->GetPtMax();
 172   Int_t iNbinsEta  = AliFlowCommonConstants::GetMaster()->GetNbinsEta();
 173   Double_t dEtaMin = AliFlowCommonConstants::GetMaster()->GetEtaMin();
 174   Double_t dEtaMax = AliFlowCommonConstants::GetMaster()->GetEtaMax();
 175
 176   fHistProFlags = new TProfile("FlowPro_Flags_SP","Flow_Flags_SP",1,0,1,"s");
 177   fHistProFlags->GetXaxis()->SetBinLabel(1,"fApplyCorrectionForNUA");
 178   fHistList->Add(fHistProFlags);
 179
 180   fHistProUQetaRP = new TProfile("FlowPro_UQetaRP_SP","Flow_UQetaRP_SP",iNbinsEta,dEtaMin,dEtaMax,"s");
 181   fHistProUQetaRP->SetXTitle("{eta}");
 182   fHistProUQetaRP->SetYTitle("<uQ>");
 183   fHistList->Add(fHistProUQetaRP);
 184
 185   fHistProUQetaPOI = new TProfile("FlowPro_UQetaPOI_SP","Flow_UQetaPOI_SP",iNbinsEta,dEtaMin,dEtaMax,"s");
 186   fHistProUQetaPOI->SetXTitle("{eta}");
 187   fHistProUQetaPOI->SetYTitle("<uQ>");
 188   fHistList->Add(fHistProUQetaPOI);
 189
 190   fHistProUQetaAllEventsPOI = new TProfile("FlowPro_UQetaAllEventsPOI_SP","FlowPro_UQetaAllEventsPOI_SP",iNbinsEta,dEtaMin,dEtaMax);
 191   fHistProUQetaAllEventsPOI->SetXTitle("{eta}");
 192   fHistProUQetaAllEventsPOI->SetYTitle("<uQ>");
 193   fHistList->Add(fHistProUQetaAllEventsPOI);
 194
 195   fHistProUQPtRP = new TProfile("FlowPro_UQPtRP_SP","Flow_UQPtRP_SP",iNbinsPt,dPtMin,dPtMax,"s");
 196   fHistProUQPtRP->SetXTitle("p_t (GeV)");
 197   fHistProUQPtRP->SetYTitle("<uQ>");
 198   fHistList->Add(fHistProUQPtRP);
 199
 200   fHistProUQPtPOI = new TProfile("FlowPro_UQPtPOI_SP","Flow_UQPtPOI_SP",iNbinsPt,dPtMin,dPtMax,"s");
 201   fHistProUQPtPOI->SetXTitle("p_t (GeV)");
 202   fHistProUQPtPOI->SetYTitle("<uQ>");
 203   fHistList->Add(fHistProUQPtPOI);
 204
 205   fHistProUQPtAllEventsPOI = new TProfile("FlowPro_UQPtAllEventsPOI_SP","FlowPro_UQPtAllEventsPOI_SP",iNbinsPt,dPtMin,dPtMax);
 206   fHistProUQPtAllEventsPOI->SetXTitle("p_t (GeV)");
 207   fHistProUQPtAllEventsPOI->SetYTitle("<uQ>");
 208   fHistList->Add(fHistProUQPtAllEventsPOI);
 209
 210   fHistProQNorm = new TProfile("FlowPro_QNorm_SP","FlowPro_QNorm_SP", 1, 0.5, 1.5,"s");
 211   fHistProQNorm ->SetYTitle("<|Qa+Qb|>");
 212   fHistList->Add(fHistProQNorm);
 213
 214   fHistProQaQb = new TProfile("FlowPro_QaQb_SP","FlowPro_QaQb_SP", 1, 0.5, 1.5,"s");
 215   fHistProQaQb->SetYTitle("<QaQb>");
 216   fHistList->Add(fHistProQaQb);
 217
 218   fHistProQaQbNorm = new TProfile("FlowPro_QaQbNorm_SP","FlowPro_QaQbNorm_SP", 1, 0.5, 1.5,"s");
 219   fHistProQaQbNorm->SetYTitle("<QaQb/MaMb>");
 220   fHistList->Add(fHistProQaQbNorm);
 221
 222   fHistProQaQbReImNorm = new TProfile("FlowPro_QaQbReImNorm_SP","FlowPro_QaQbReImNorm_SP", 4, 0.5, 4.5,"s");
 223   fHistProQaQbReImNorm->GetXaxis()->SetBinLabel(1,"#LT#LTsin(#phi_{a})#GT#GT");
 224   fHistProQaQbReImNorm->GetXaxis()->SetBinLabel(2,"#LT#LTcos(#phi_{a})#GT#GT");
 225   fHistProQaQbReImNorm->GetXaxis()->SetBinLabel(3,"#LT#LTsin(#phi_{b})#GT#GT");
 226   fHistProQaQbReImNorm->GetXaxis()->SetBinLabel(4,"#LT#LTcos(#phi_{b})#GT#GT");
 227   fHistList->Add(fHistProQaQbReImNorm);
 228
 229   fHistProNonIsotropicTermsQ = new TProfile("FlowPro_NonIsotropicTermsQ_SP","FlowPro_NonIsotropicTermsQ_SP", 2, 0.5, 2.5,"s");
 230   fHistProNonIsotropicTermsQ->GetXaxis()->SetBinLabel(1,"#LT#LTsin(#phi_{a+b})#GT#GT");
 231   fHistProNonIsotropicTermsQ->GetXaxis()->SetBinLabel(2,"#LT#LTcos(#phi_{a+b})#GT#GT");
 232   fHistList->Add(fHistProNonIsotropicTermsQ);
 233
 234   TString rpPoi[2] = {"RP","POI"};
 235   TString ptEta[2] = {"Pt","Eta"};
 236   TString sinCos[2] = {"sin","cos"};
 237   Int_t nBinsPtEta[2] = {iNbinsPt,iNbinsEta};
 238   Double_t minPtEta[2] = {dPtMin,dEtaMin};
 239   Double_t maxPtEta[2] = {dPtMax,dEtaMax};
 240   for(Int_t rp=0;rp<2;rp++)
 241   {
 242    for(Int_t pe=0;pe<2;pe++)
 243    {
 244     for(Int_t sc=0;sc<2;sc++)
 245     {
 246      fHistProNonIsotropicTermsU[rp][pe][sc] = new TProfile(Form("FlowPro_NonIsotropicTerms_%s_%s_%s_SP",rpPoi[rp].Data(),ptEta[pe].Data(),sinCos[sc].Data()),Form("FlowPro_NonIsotropicTerms_%s_%s_%s_SP",rpPoi[rp].Data(),ptEta[pe].Data(),sinCos[sc].Data()),nBinsPtEta[pe],minPtEta[pe],maxPtEta[pe]);
 247      fHistList->Add(fHistProNonIsotropicTermsU[rp][pe][sc]);
 248     }
 249    }
 250   }
 251
 252   fHistSumOfLinearWeights = new TH1D("Flow_SumOfLinearWeights_SP","Flow_SumOfLinearWeights_SP",1,-0.5, 0.5);
 253   fHistSumOfLinearWeights -> SetYTitle("sum (*)");
 254   fHistSumOfLinearWeights -> SetXTitle("sum (Ma*Mb)");
 255   fHistList->Add(fHistSumOfLinearWeights);
 256
 257   fHistSumOfQuadraticWeights = new TH1D("Flow_SumOfQuadraticWeights_SP","Flow_SumOfQuadraticWeights_SP",1,-0.5, 0.5);
 258   fHistSumOfQuadraticWeights -> SetYTitle("sum (*)");
 259   fHistSumOfQuadraticWeights -> SetXTitle("sum (Ma*Mb)^2");
 260   fHistList->Add(fHistSumOfQuadraticWeights);
 261
 262   fHistProUQQaQbPtRP = new TProfile("FlowPro_UQQaQbPtRP_SP","FlowPro_UQQaQbPtRP_SP",iNbinsPt,dPtMin,dPtMax);
 263   fHistProUQQaQbPtRP -> SetYTitle("<*>");
 264   fHistProUQQaQbPtRP -> SetXTitle("<Qu QaQb>");
 265   fHistList->Add(fHistProUQQaQbPtRP);
 266
 267   fHistProUQQaQbEtaRP = new TProfile("FlowPro_UQQaQbEtaRP_SP","FlowPro_UQQaQbEtaRP_SP",iNbinsEta,dEtaMin,dEtaMax);
 268   fHistProUQQaQbEtaRP -> SetYTitle("<*>");
 269   fHistProUQQaQbEtaRP -> SetXTitle("<Qu QaQb>");
 270   fHistList->Add(fHistProUQQaQbEtaRP);
 271
 272   fHistProUQQaQbPtPOI = new TProfile("FlowPro_UQQaQbPtPOI_SP","FlowPro_UQQaQbPtPOI_SP",iNbinsPt,dPtMin,dPtMax);
 273   fHistProUQQaQbPtPOI -> SetYTitle("<*>");
 274   fHistProUQQaQbPtPOI -> SetXTitle("<Qu QaQb>");
 275   fHistList->Add(fHistProUQQaQbPtPOI);
 276
 277   fHistProUQQaQbEtaPOI = new TProfile("FlowPro_UQQaQbEtaPOI_SP","FlowPro_UQQaQbEtaPOI_SP",iNbinsEta,dEtaMin,dEtaMax);
 278   fHistProUQQaQbEtaPOI -> SetYTitle("<*>");
 279   fHistProUQQaQbEtaPOI -> SetXTitle("<Qu QaQb>");
 280   fHistList->Add(fHistProUQQaQbEtaPOI);
 281
 282   TString weightFlag[3] = {"w_Qu_","w_Qu^2_","w_QuQaQb_"};
 283   for(Int_t i=0;i<3;i++)
 284   {
 285    fHistSumOfWeightsPtRP[i] = new TH1D(Form("Flow_SumOfWeights%sPtRP_SP",weightFlag[i].Data()),
 286                               Form("Flow_SumOfWeights%sPtRP_SP",weightFlag[i].Data()),iNbinsPt,dPtMin,dPtMax);
 287    fHistSumOfWeightsPtRP[i] -> SetYTitle("sum (*)");
 288    fHistSumOfWeightsPtRP[i] -> SetXTitle("p_{T}");
 289    fHistList->Add(fHistSumOfWeightsPtRP[i]);
 290
 291    fHistSumOfWeightsEtaRP[i] = new TH1D(Form("Flow_SumOfWeights%sEtaRP_SP",weightFlag[i].Data()),
 292                                Form("Flow_SumOfWeights%sEtaRP_SP",weightFlag[i].Data()),iNbinsEta,dEtaMin,dEtaMax);
 293    fHistSumOfWeightsEtaRP[i] -> SetYTitle("sum (*)");
 294    fHistSumOfWeightsEtaRP[i] -> SetXTitle("#eta");
 295    fHistList->Add(fHistSumOfWeightsEtaRP[i]);
 296
 297    fHistSumOfWeightsPtPOI[i] = new TH1D(Form("Flow_SumOfWeights%sPtPOI_SP",weightFlag[i].Data()),
 298                                Form("Flow_SumOfWeights%sPtPOI_SP",weightFlag[i].Data()),iNbinsPt,dPtMin,dPtMax);
 299    fHistSumOfWeightsPtPOI[i] -> SetYTitle("sum (*)");
 300    fHistSumOfWeightsPtPOI[i] -> SetXTitle("p_{T}");
 301    fHistList->Add(fHistSumOfWeightsPtPOI[i]);
 302
 303    fHistSumOfWeightsEtaPOI[i] = new TH1D(Form("Flow_SumOfWeights%sEtaPOI_SP",weightFlag[i].Data()),
 304                                 Form("Flow_SumOfWeights%sEtaPOI_SP",weightFlag[i].Data()),iNbinsEta,dEtaMin,dEtaMax);
 305    fHistSumOfWeightsEtaPOI[i] -> SetYTitle("sum (*)");
 306    fHistSumOfWeightsEtaPOI[i] -> SetXTitle("#eta");
 307    fHistList->Add(fHistSumOfWeightsEtaPOI[i]);
 308   }
 309
 310   fCommonHistsSP = new AliFlowCommonHist("AliFlowCommonHistSP");
 311   fHistList->Add(fCommonHistsSP);
 312   fCommonHistsResSP = new AliFlowCommonHistResults("AliFlowCommonHistResultsSP");
 313   fHistList->Add(fCommonHistsResSP);
 314   fCommonHistsmuQ = new AliFlowCommonHist("AliFlowCommonHistmuQ");
 315   fHistList->Add(fCommonHistsmuQ);
 316
 317   fHistQNorm = new TH1D("Flow_QNorm_SP","Flow_QNorm_SP",110,0.,1.1);
 318   fHistQNorm -> SetYTitle("dN/d(|(Qa+Qb)/(Ma+Mb)|)");
 319   fHistQNorm -> SetXTitle("|(Qa+Qb)/(Ma+Mb)|");
 320   fHistList->Add(fHistQNorm);
 321
 322   fHistQaQb = new TH1D("Flow_QaQb_SP","Flow_QaQb_SP",200,-100.,100.);
 323   fHistQaQb -> SetYTitle("dN/dQaQb");
 324   fHistQaQb -> SetXTitle("QaQb");
 325   fHistList->Add(fHistQaQb);
 326
 327   fHistQaQbNorm = new TH1D("Flow_QaQbNorm_SP","Flow_QaQbNorm_SP",44,-1.1,1.1);
 328   fHistQaQbNorm -> SetYTitle("dN/d(QaQb/MaMb)");
 329   fHistQaQbNorm -> SetXTitle("QaQb/MaMb");
 330   fHistList->Add(fHistQaQbNorm);
 331
 332   fHistQNormvsQaQbNorm = new TH2D("Flow_QNormvsQaQbNorm_SP","Flow_QNormvsQaQbNorm_SP",88,-1.1,1.1,22,0.,1.1);
 333   fHistQNormvsQaQbNorm -> SetYTitle("|Q/Mq|");
 334   fHistQNormvsQaQbNorm -> SetXTitle("QaQb/MaMb");
 335   fHistList->Add(fHistQNormvsQaQbNorm);
 336
 337   fHistQaQbCos = new TH1D("Flow_QaQbCos_SP","Flow_QaQbCos_SP",63,0.,TMath::Pi());
 338   fHistQaQbCos -> SetYTitle("dN/d(#phi)");
 339   fHistQaQbCos -> SetXTitle("#phi");
 340   fHistList->Add(fHistQaQbCos);
 341
 342   fHistResolution = new TH1D("Flow_resolution_SP","Flow_resolution_SP",100,-1.0,1.0);
 343   fHistResolution -> SetYTitle("dN/d(cos(2(#phi_a - #phi_b))");
 344   fHistResolution -> SetXTitle("cos(2*(#phi_a - #phi_b))");
 345   fHistList->Add(fHistResolution);
 346
 347   fHistQaNorm = new TH1D("Flow_QaNorm_SP","Flow_QaNorm_SP",22,0.,1.1);
 348   fHistQaNorm -> SetYTitle("dN/d(|Qa/Ma|)");
 349   fHistQaNorm -> SetXTitle("|Qa/Ma|");
 350   fHistList->Add(fHistQaNorm);
 351
 352   fHistQaNormvsMa = new TH2D("Flow_QaNormvsMa_SP","Flow_QaNormvsMa_SP",100,0.,100.,22,0.,1.1);
 353   fHistQaNormvsMa -> SetYTitle("|Qa/Ma|");
 354   fHistQaNormvsMa -> SetXTitle("Ma");
 355   fHistList->Add(fHistQaNormvsMa);
 356
 357   fHistQbNorm = new TH1D("Flow_QbNorm_SP","Flow_QbNorm_SP",22,0.,1.1);
 358   fHistQbNorm -> SetYTitle("dN/d(|Qb/Mb|)");
 359   fHistQbNorm -> SetXTitle("|Qb/Mb|");
 360   fHistList->Add(fHistQbNorm);
 361
 362   fHistQbNormvsMb = new TH2D("Flow_QbNormvsMb_SP","Flow_QbNormvsMb_SP",100,0.,100.,22,0.,1.1);
 363   fHistQbNormvsMb -> SetYTitle("|Qb/Mb|");
 364   fHistQbNormvsMb -> SetXTitle("|Mb|");
 365   fHistList->Add(fHistQbNormvsMb);
 366
 367   fHistMavsMb = new TH2D("Flow_MavsMb_SP","Flow_MavsMb_SP",100,0.,100.,100,0.,100.);
 368   fHistMavsMb -> SetYTitle("Ma");
 369   fHistMavsMb -> SetXTitle("Mb");
 370   fHistList->Add(fHistMavsMb);
 371
 372
 373   //weights
 374   if(fUsePhiWeights) {
 375     if(!fWeightsList) {
 376       cout<<"WARNING: fWeightsList is NULL in the Scalar Product method."<<endl;
 377       exit(0);
 378     }
 379     if(fWeightsList->FindObject("phi_weights_sub0"))  {
 380       fPhiWeightsSub0 = dynamic_cast<TH1F*>
 381         (fWeightsList->FindObject("phi_weights_sub0"));
 382       fHistList->Add(fPhiWeightsSub0);
 383     } else {
 384       cout<<"WARNING: histogram with phi weights is not accessible in Scalar Product"<<endl;
 385       exit(0);
 386     }
 387     if(fWeightsList->FindObject("phi_weights_sub1"))  {
 388       fPhiWeightsSub1 = dynamic_cast<TH1F*>
 389         (fWeightsList->FindObject("phi_weights_sub1"));
 390       fHistList->Add(fPhiWeightsSub1);
 391     } else {
 392       cout<<"WARNING: histogram with phi weights is not accessible in Scalar Product"<<endl;
 393       exit(0);
 394     }
 395
 396   } // end of if(fUsePhiWeights)
 397
 398   fEventNumber = 0;  //set number of events to zero
 399
 400   //store all boolean flags needed in Finish():
 401   this->StoreFlags();
 402
 403   TH1::AddDirectory(oldHistAddStatus);
 404 }
 405
 406 //-----------------------------------------------------------------------
 407 void AliFlowAnalysisWithScalarProduct::Make(AliFlowEventSimple* anEvent) {
 408
 409
 410   if (anEvent) {
 411
 412   //Calculate muQ (for comparing pp and PbPb)
 413   FillmuQ(anEvent);
 414
 415   //Calculate flow based on  <QaQb/MaMb> = <v^2>
 416   FillSP(anEvent);
 417
 418   }
 419 }
 420
 421 //-----------------------------------------------------------------------
 422 void AliFlowAnalysisWithScalarProduct::FillSP(AliFlowEventSimple* anEvent) {
 423
 424   //Calculate flow based on  <QaQb/MaMb> = <v^2>
 425
 426   //Fill histograms
 427   if (anEvent) {
 428
 429     //get Q vectors for the eta-subevents
 430     AliFlowVector* vQarray = new AliFlowVector[2];
 431     if (fUsePhiWeights) {
 432       anEvent->Get2Qsub(vQarray,2,fWeightsList,kTRUE);
 433     } else {
 434       anEvent->Get2Qsub(vQarray);
 435     }
 436     //subevent a
 437     AliFlowVector vQa = vQarray[0];
 438     //subevent b
 439     AliFlowVector vQb = vQarray[1];
 440
 441     //For calculating v2 only events should be taken where both subevents are not empty
 442     //check that the subevents are not empty:
 443     Double_t dMa = vQa.GetMult();
 444     Double_t dMb = vQb.GetMult();
 445     if (dMa > 0. && dMb > 0.) {
 446
 447       //request that the subevent multiplicities are not too different
 448       //fRelDiffMsub can be set from the configuration macro
 449       Double_t dRelDiff = TMath::Abs((dMa - dMb)/(dMa + dMb));
 450       if (dRelDiff < fRelDiffMsub) {
 451
 452         //fill control histograms
 453         if (fUsePhiWeights) {
 454           fCommonHistsSP->FillControlHistograms(anEvent,fWeightsList,kTRUE);
 455         } else {
 456           fCommonHistsSP->FillControlHistograms(anEvent);
 457         }
 458
 459         //fill some SP control histograms
 460         fHistProQaQb -> Fill(1.,vQa*vQb,1.); //Fill with weight 1 -> Weight with MaMb????
 461         fHistQaQbCos ->Fill(TMath::ACos((vQa/vQa.Mod())*(vQb/vQb.Mod())));  //Acos(Qa*Qb) = angle
 462         fHistResolution -> Fill(TMath::Cos( vQa.Phi()- vQb.Phi() ));  //vQa.Phi() returns 2*phi
 463         fHistQaQb -> Fill(vQa*vQb);
 464         fHistMavsMb -> Fill(dMb,dMa);
 465
 466         //get total Q vector = the sum of subevent a and subevent b
 467         AliFlowVector vQ = vQa + vQb;
 468
 469         //needed to correct for non-uniform acceptance:
 470         fHistProNonIsotropicTermsQ->Fill(1.,vQ.Y()/(dMa+dMb),dMa+dMb);
 471         fHistProNonIsotropicTermsQ->Fill(2.,vQ.X()/(dMa+dMb),dMa+dMb);
 472
 473         //weight the Q vectors for the subevents by the multiplicity
 474         //Note: Weight Q only in the particle loop when it is clear
 475         //if it should be (m-1) or M
 476         Double_t dQXa = vQa.X()/dMa;
 477         Double_t dQYa = vQa.Y()/dMa;
 478         vQa.Set(dQXa,dQYa);
 479
 480         Double_t dQXb = vQb.X()/dMb;
 481         Double_t dQYb = vQb.Y()/dMb;
 482         vQb.Set(dQXb,dQYb);
 483
 484         //scalar product of the two subevents
 485         Double_t dQaQb = (vQa*vQb);
 486         fHistProQaQbNorm -> Fill(1.,dQaQb,dMa*dMb);  //Fill (QaQb/MaMb) with weight (MaMb).
 487         //needed for the error calculation:
 488         fHistSumOfLinearWeights -> Fill(0.,dMa*dMb);
 489         fHistSumOfQuadraticWeights -> Fill(0.,pow(dMa*dMb,2.));
 490         //needed for correcting non-uniform acceptance:
 491         fHistProQaQbReImNorm->Fill(1.,dQYa,dMa); // to get <<sin(phi_a)>>
 492         fHistProQaQbReImNorm->Fill(2.,dQXa,dMa); // to get <<cos(phi_a)>>
 493         fHistProQaQbReImNorm->Fill(3.,dQYb,dMb); // to get <<sin(phi_b)>>
 494         fHistProQaQbReImNorm->Fill(4.,dQXb,dMb); // to get <<cos(phi_b)>>
 495
 496         //fill some SP control histograms
 497         fHistQaQbNorm ->Fill(vQa*vQb);
 498         fHistQaNorm ->Fill(vQa.Mod());
 499         fHistQaNormvsMa->Fill(dMa,vQa.Mod());
 500         fHistQbNorm ->Fill(vQb.Mod());
 501         fHistQbNormvsMb->Fill(dMb,vQb.Mod());
 502
 503         //loop over the tracks of the event
 504         AliFlowTrackSimple*   pTrack = NULL;
 505         Int_t iNumberOfTracks = anEvent->NumberOfTracks();
 506         Double_t dMq =  vQ.GetMult();
 507
 508         for (Int_t i=0;i<iNumberOfTracks;i++)
 509           {
 510             pTrack = anEvent->GetTrack(i) ;
 511             if (pTrack){
 512               Double_t dPhi = pTrack->Phi();
 513
 514               //calculate vU
 515               TVector2 vU;
 516               //do not need to use weight for v as the length will be made 1
 517               Double_t dUX = TMath::Cos(2*dPhi);
 518               Double_t dUY = TMath::Sin(2*dPhi);
 519               vU.Set(dUX,dUY);
 520               Double_t dModulus = vU.Mod();
 521               if (dModulus > 0.) vU.Set(dUX/dModulus,dUY/dModulus);  // make length 1
 522               else cerr<<"dModulus is zero!"<<endl;
 523
 524               //redefine the Q vector and devide by its multiplicity
 525               TVector2 vQm;
 526               Double_t dQmX = 0.;
 527               Double_t dQmY = 0.;
 528               //subtract particle from the flowvector if used to define it
 529               if (pTrack->InSubevent(0) || pTrack->InSubevent(1)) {
 530                 //set default phi weight to 1
 531                 Double_t dW = 1.;
 532                 //if phi weights are used
 533                 if(fUsePhiWeights && fPhiWeightsSub0 && fPhiWeightsSub1) {
 534                   //value of the center of the phi bin
 535                   Double_t dPhiCenter = 0.;
 536                   if (pTrack->InSubevent(0) ) {
 537                     Int_t iNBinsPhiSub0 = fPhiWeightsSub0->GetNbinsX();
 538                     Int_t phiBin = 1+(Int_t)(TMath::Floor(dPhi*iNBinsPhiSub0/TMath::TwoPi()));
 539                     dW = fPhiWeightsSub0->GetBinContent(phiBin);
 540                     dPhiCenter = fPhiWeightsSub0->GetBinCenter(phiBin);
 541                     dQmX = (vQ.X() - dW*(pTrack->Weight())* TMath::Cos(2*dPhiCenter) )/(dMq-1);
 542                     dQmY = (vQ.Y() - dW*(pTrack->Weight())* TMath::Sin(2*dPhiCenter) )/(dMq-1);
 543
 544                     vQm.Set(dQmX,dQmY);
 545                   }
 546
 547                   else if ( pTrack->InSubevent(1)) {
 548                     Int_t iNBinsPhiSub1 = fPhiWeightsSub1->GetNbinsX();
 549                     Int_t phiBin = 1+(Int_t)(TMath::Floor(dPhi*iNBinsPhiSub1/TMath::TwoPi()));
 550                     dW = fPhiWeightsSub1->GetBinContent(phiBin);
 551                     dPhiCenter = fPhiWeightsSub1->GetBinCenter(phiBin);
 552                     dQmX = (vQ.X() - dW*(pTrack->Weight())* TMath::Cos(2*dPhiCenter) )/(dMq-1);
 553                     dQmY = (vQ.Y() - dW*(pTrack->Weight())* TMath::Sin(2*dPhiCenter) )/(dMq-1);
 554
 555                     vQm.Set(dQmX,dQmY);
 556                   }
 557                   //bin = 1 + value*nbins/range
 558                   //TMath::Floor rounds to the lower integer
 559                 }
 560                 // if no phi weights are used
 561                 else {
 562                   dQmX = (vQ.X() - (pTrack->Weight())*dUX)/(dMq-1);
 563                   dQmY = (vQ.Y() - (pTrack->Weight())*dUY)/(dMq-1);
 564                   vQm.Set(dQmX,dQmY);
 565                 }
 566
 567                 //dUQ = scalar product of vU and vQm
 568                 Double_t dUQ = (vU * vQm);
 569                 Double_t dPt = pTrack->Pt();
 570                 Double_t dEta = pTrack->Eta();
 571
 572                 //fill the profile histograms
 573                 if (pTrack->InRPSelection()) {
 574                   fHistProUQetaRP -> Fill(dEta,dUQ,(dMq-1)); //Fill (Qu/(Mq-1)) with weight (Mq-1)
 575                   //needed for the error calculation:
 576                   fHistProUQQaQbEtaRP -> Fill(dEta,dUQ*dQaQb,(dMq-1)*dMa*dMb); //Fill [Qu/(Mq-1)]*[QaQb/MaMb] with weight (Mq-1)MaMb
 577                   fHistProUQPtRP -> Fill(dPt,dUQ,(dMq-1));                     //Fill (Qu/(Mq-1)) with weight (Mq-1)
 578                   fHistProUQQaQbPtRP -> Fill(dPt,dUQ*dQaQb,(dMq-1)*dMa*dMb);   //Fill [Qu/(Mq-1)]*[QaQb/MaMb] with weight (Mq-1)MaMb
 579
 580                   fHistSumOfWeightsEtaRP[0]->Fill(dEta,(dMq-1));        // sum of Mq-1
 581                   fHistSumOfWeightsEtaRP[1]->Fill(dEta,pow((dMq-1),2.));// sum of (Mq-1)^2
 582                   fHistSumOfWeightsEtaRP[2]->Fill(dEta,(dMq-1)*dMa*dMb);// sum of (Mq-1)*MaMb
 583                   fHistSumOfWeightsPtRP[0]->Fill(dPt,(dMq-1));          // sum of Mq-1
 584                   fHistSumOfWeightsPtRP[1]->Fill(dPt,pow((dMq-1),2.));  // sum of (Mq-1)^2
 585                   fHistSumOfWeightsPtRP[2]->Fill(dPt,(dMq-1)*dMa*dMb);  // sum of (Mq-1)*MaMb
 586                   //nonisotropic terms:
 587                   fHistProNonIsotropicTermsU[0][0][0]->Fill(dPt,dUY,1.);
 588                   fHistProNonIsotropicTermsU[0][0][1]->Fill(dPt,dUX,1.);
 589                   fHistProNonIsotropicTermsU[0][1][0]->Fill(dEta,dUY,1.);
 590                   fHistProNonIsotropicTermsU[0][1][1]->Fill(dEta,dUX,1.);
 591                 }
 592                 if (pTrack->InPOISelection()) {
 593                   fHistProUQetaPOI -> Fill(dEta,dUQ,(dMq-1));//Fill (Qu/(Mq-1)) with weight (Mq-1)
 594                   //needed for the error calculation:
 595                   fHistProUQQaQbEtaPOI -> Fill(dEta,dUQ*dQaQb,(dMq-1)*dMa*dMb); //Fill [Qu/(Mq-1)]*[QaQb/MaMb] with weight (Mq-1)MaMb
 596                   fHistProUQPtPOI -> Fill(dPt,dUQ,(dMq-1));                     //Fill (Qu/(Mq-1)) with weight (Mq-1)
 597                   fHistProUQQaQbPtPOI -> Fill(dPt,dUQ*dQaQb,(dMq-1)*dMa*dMb);   //Fill [Qu/(Mq-1)]*[QaQb/MaMb] with weight (Mq-1)MaMb
 598
 599                   fHistSumOfWeightsEtaPOI[0]->Fill(dEta,(dMq-1));        // sum of Mq-1
 600                   fHistSumOfWeightsEtaPOI[1]->Fill(dEta,pow((dMq-1),2.));// sum of (Mq-1)^2
 601                   fHistSumOfWeightsEtaPOI[2]->Fill(dEta,(dMq-1)*dMa*dMb);// sum of (Mq-1)*MaMb
 602                   fHistSumOfWeightsPtPOI[0]->Fill(dPt,(dMq-1));          // sum of Mq-1
 603                   fHistSumOfWeightsPtPOI[1]->Fill(dPt,pow((dMq-1),2.)); // sum of (Mq-1)^2
 604                   fHistSumOfWeightsPtPOI[2]->Fill(dPt,(dMq-1)*dMa*dMb); // sum of (Mq-1)*MaMb
 605                   //nonisotropic terms:
 606                   fHistProNonIsotropicTermsU[1][0][0]->Fill(dPt,dUY,1.);
 607                   fHistProNonIsotropicTermsU[1][0][1]->Fill(dPt,dUX,1.);
 608                   fHistProNonIsotropicTermsU[1][1][0]->Fill(dEta,dUY,1.);
 609                   fHistProNonIsotropicTermsU[1][1][1]->Fill(dEta,dUX,1.);
 610                 }
 611
 612               } else { //do not subtract the particle from the flowvector
 613                 dQmX = vQ.X()/dMq;
 614                 dQmY = vQ.Y()/dMq;
 615                 vQm.Set(dQmX,dQmY);
 616
 617                 //fill histograms with vQm
 618                 fHistProQNorm->Fill(1.,vQm.Mod(),dMq);
 619                 fHistQNorm->Fill(vQm.Mod());
 620                 fHistQNormvsQaQbNorm->Fill(vQa*vQb ,vQm.Mod());
 621
 622                 //dUQ = scalar product of vU and vQm
 623                 Double_t dUQ = (vU * vQm);
 624                 Double_t dPt = pTrack->Pt();
 625                 Double_t dEta = pTrack->Eta();
 626
 627                 //fill the profile histograms
 628                 if (pTrack->InRPSelection()) {
 629                   fHistProUQetaRP -> Fill(dEta,dUQ,dMq);                   //Fill (Qu/Mq) with weight Mq
 630                   //needed for the error calculation:
 631                   fHistProUQQaQbEtaRP -> Fill(dEta,dUQ*dQaQb,dMq*dMa*dMb); //Fill [Qu/Mq]*[QaQb/MaMb] with weight Mq*MaMb
 632                   fHistProUQPtRP -> Fill(dPt,dUQ,dMq);                     //Fill (Qu/Mq) with weight Mq
 633                   fHistProUQQaQbPtRP -> Fill(dPt,dUQ*dQaQb,dMq*dMa*dMb);   //Fill [Qu/Mq]*[QaQb/MaMb] with weight Mq*MaMb
 634
 635                   fHistSumOfWeightsEtaRP[0]->Fill(dEta,dMq);        // sum of Mq
 636                   fHistSumOfWeightsEtaRP[1]->Fill(dEta,pow(dMq,2.));// sum of Mq^2
 637                   fHistSumOfWeightsEtaRP[2]->Fill(dEta,dMq*dMa*dMb);// sum of Mq*MaMb
 638                   fHistSumOfWeightsPtRP[0]->Fill(dPt,dMq);          // sum of Mq
 639                   fHistSumOfWeightsPtRP[1]->Fill(dPt,pow(dMq,2.));  // sum of Mq^2
 640                   fHistSumOfWeightsPtRP[2]->Fill(dPt,dMq*dMa*dMb);  // sum of Mq*MaMb
 641                   //nonisotropic terms:
 642                   fHistProNonIsotropicTermsU[0][0][0]->Fill(dPt,dUY,1.);
 643                   fHistProNonIsotropicTermsU[0][0][1]->Fill(dPt,dUX,1.);
 644                   fHistProNonIsotropicTermsU[0][1][0]->Fill(dEta,dUY,1.);
 645                   fHistProNonIsotropicTermsU[0][1][1]->Fill(dEta,dUX,1.);
 646                 }
 647                 if (pTrack->InPOISelection()) {
 648                   fHistProUQetaPOI -> Fill(dEta,dUQ,dMq); //Fill (Qu/Mq) with weight Mq
 649                   //needed for the error calculation:
 650                   fHistProUQQaQbEtaPOI -> Fill(dEta,dUQ*dQaQb,dMq*dMa*dMb); //Fill [Qu/Mq]*[QaQb/MaMb] with weight Mq*MaMb
 651                   fHistProUQPtPOI -> Fill(dPt,dUQ,dMq);                     //Fill (Qu/Mq) with weight Mq
 652                   fHistProUQQaQbPtPOI -> Fill(dPt,dUQ*dQaQb,dMq*dMa*dMb);   //Fill [Qu/Mq]*[QaQb/MaMb] with weight Mq*MaMb
 653
 654                   fHistSumOfWeightsEtaPOI[0]->Fill(dEta,dMq);        // sum of Mq
 655                   fHistSumOfWeightsEtaPOI[1]->Fill(dEta,pow(dMq,2.));// sum of Mq^2
 656                   fHistSumOfWeightsEtaPOI[2]->Fill(dEta,dMq*dMa*dMb);// sum of Mq*MaMb
 657                   fHistSumOfWeightsPtPOI[0]->Fill(dPt,dMq);          // sum of Mq
 658                   fHistSumOfWeightsPtPOI[1]->Fill(dPt,pow(dMq,2.));  // sum of Mq^2
 659                   fHistSumOfWeightsPtPOI[2]->Fill(dPt,dMq*dMa*dMb);  // sum of Mq*MaMb
 660                   //nonisotropic terms:
 661                   fHistProNonIsotropicTermsU[1][0][0]->Fill(dPt,dUY,1.);
 662                   fHistProNonIsotropicTermsU[1][0][1]->Fill(dPt,dUX,1.);
 663                   fHistProNonIsotropicTermsU[1][1][0]->Fill(dEta,dUY,1.);
 664                   fHistProNonIsotropicTermsU[1][1][1]->Fill(dEta,dUX,1.);
 665                 }
 666               }//track not in subevents
 667
 668             }//track
 669
 670           }//loop over tracks
 671
 672         fEventNumber++;
 673
 674       } //difference Ma and Mb
 675
 676     }// subevents not empty
 677     delete [] vQarray;
 678
 679   } //event
 680
 681 }//end of FillSP()
 682
 683 //-----------------------------------------------------------------------
 684 void AliFlowAnalysisWithScalarProduct::FillmuQ(AliFlowEventSimple* anEvent) {
 685
 686   if (anEvent) {
 687
 688     //get Q vectors for the eta-subevents
 689     AliFlowVector* vQarray = new AliFlowVector[2];
 690     if (fUsePhiWeights) {
 691       anEvent->Get2Qsub(vQarray,2,fWeightsList,kTRUE);
 692     } else {
 693       anEvent->Get2Qsub(vQarray);
 694     }
 695     //subevent a
 696     AliFlowVector vQa = vQarray[0];
 697     //subevent b
 698     AliFlowVector vQb = vQarray[1];
 699
 700     //get total Q vector = the sum of subevent a and subevent b
 701     AliFlowVector vQ;
 702     if (vQa.GetMult() > 0 && vQb.GetMult() > 0) {
 703       vQ = vQa + vQb;
 704     } else if ( vQa.GetMult() > 0 && !(vQb.GetMult() > 0) ){
 705       vQ = vQa;
 706     } else if ( !(vQa.GetMult() > 0) && vQb.GetMult() > 0 ){
 707       vQ = vQb;
 708     }
 709
 710     //For calculating uQ for comparison all events should be taken also if one of the subevents is empty
 711     //check if the total Q vector is not empty
 712     Double_t dMq =  vQ.GetMult();
 713     if (dMq > 0.) {
 714
 715       //Fill control histograms
 716       if (fUsePhiWeights) {
 717         fCommonHistsmuQ->FillControlHistograms(anEvent,fWeightsList,kTRUE);
 718       } else {
 719         fCommonHistsmuQ->FillControlHistograms(anEvent);
 720       }
 721
 722       //loop over all POI tracks and fill uQ
 723       AliFlowTrackSimple*   pTrack = NULL;
 724       for (Int_t i=0;i<anEvent->NumberOfTracks();i++) {
 725         pTrack = anEvent->GetTrack(i) ;
 726         if (pTrack){
 727
 728           if (pTrack->InPOISelection()) {
 729
 730             Double_t dPhi = pTrack->Phi();
 731             //weights do not need to be used as the length of vU will be set to 1
 732
 733             //calculate vU
 734             TVector2 vU;
 735             Double_t dUX = TMath::Cos(2*dPhi);
 736             Double_t dUY = TMath::Sin(2*dPhi);
 737             vU.Set(dUX,dUY);
 738             Double_t dModulus = vU.Mod();
 739             // make length 1
 740             if (dModulus!=0.) vU.Set(dUX/dModulus,dUY/dModulus);
 741             else cerr<<"dModulus is zero!"<<endl;
 742
 743             //redefine the Q vector
 744             TVector2 vQm;
 745             Double_t dQmX = 0.;
 746             Double_t dQmY = 0.;
 747             //subtract particle from the flowvector if used to define it
 748             if (pTrack->InSubevent(0) || pTrack->InSubevent(1)) {
 749               //the number of tracks contributing to vQ must be more than 1
 750               if (dMq > 1) {
 751                 //set default phi weight to 1
 752                 Double_t dW = 1.;
 753                 //if phi weights are used
 754                 if(fUsePhiWeights && fPhiWeightsSub0 && fPhiWeightsSub1) {
 755                   //value of the center of the phi bin
 756                   Double_t dPhiCenter = 0.;
 757                   if (pTrack->InSubevent(0) ) {
 758                     Int_t iNBinsPhiSub0 = fPhiWeightsSub0->GetNbinsX();
 759                     Int_t phiBin = 1+(Int_t)(TMath::Floor(dPhi*iNBinsPhiSub0/TMath::TwoPi()));
 760                     dW = fPhiWeightsSub0->GetBinContent(phiBin);
 761                     dPhiCenter = fPhiWeightsSub0->GetBinCenter(phiBin);
 762                     dQmX = (vQ.X() - dW*(pTrack->Weight())* TMath::Cos(2*dPhiCenter) );
 763                     dQmY = (vQ.Y() - dW*(pTrack->Weight())* TMath::Sin(2*dPhiCenter) );
 764
 765                     vQm.Set(dQmX,dQmY);
 766                   }
 767
 768                   else if ( pTrack->InSubevent(1)) {
 769                     Int_t iNBinsPhiSub1 = fPhiWeightsSub1->GetNbinsX();
 770                     Int_t phiBin = 1+(Int_t)(TMath::Floor(dPhi*iNBinsPhiSub1/TMath::TwoPi()));
 771                     dW = fPhiWeightsSub1->GetBinContent(phiBin);
 772                     dPhiCenter = fPhiWeightsSub1->GetBinCenter(phiBin);
 773                     dQmX = (vQ.X() - dW*(pTrack->Weight())* TMath::Cos(2*dPhiCenter) );
 774                     dQmY = (vQ.Y() - dW*(pTrack->Weight())* TMath::Sin(2*dPhiCenter) );
 775
 776                     vQm.Set(dQmX,dQmY);
 777                   }
 778                   //bin = 1 + value*nbins/range
 779                   //TMath::Floor rounds to the lower integer
 780                 }
 781                 // if no phi weights are used
 782                 else {
 783                   dQmX = (vQ.X() - (pTrack->Weight())*dUX);
 784                   dQmY = (vQ.Y() - (pTrack->Weight())*dUY);
 785                   vQm.Set(dQmX,dQmY);
 786                 }
 787
 788                 //dUQ = scalar product of vU and vQm
 789                 Double_t dUQ = (vU * vQm);
 790                 Double_t dPt = pTrack->Pt();
 791                 Double_t dEta = pTrack->Eta();
 792                 //fill the profile histograms
 793                 fHistProUQetaAllEventsPOI -> Fill(dEta,dUQ);   //Fill (Qu)
 794                 fHistProUQPtAllEventsPOI -> Fill(dPt,dUQ);     //Fill (Qu)
 795
 796               } //dMq > 1
 797             }
 798             else { //do not subtract the particle from the flowvector
 799
 800               dQmX = vQ.X();
 801               dQmY = vQ.Y();
 802               vQm.Set(dQmX,dQmY);
 803
 804               //dUQ = scalar product of vU and vQm
 805               Double_t dUQ = (vU * vQm);
 806               Double_t dPt = pTrack->Pt();
 807               Double_t dEta = pTrack->Eta();
 808               //fill the profile histograms
 809               fHistProUQetaAllEventsPOI -> Fill(dEta,dUQ);   //Fill (Qu)
 810               fHistProUQPtAllEventsPOI -> Fill(dPt,dUQ);     //Fill (Qu)
 811
 812             }
 813
 814           } //in POI selection
 815         } //track valid
 816       } //end of loop over tracks
 817     } //Q vector is not empty
 818
 819   } //anEvent valid
 820
 821 } //end of FillmuQ
 822
 823 //--------------------------------------------------------------------
 824 void AliFlowAnalysisWithScalarProduct::GetOutputHistograms(TList *outputListHistos){
 825
 826   //get pointers to all output histograms (called before Finish())
 827
 828   if (outputListHistos) {
 829   //Get the common histograms from the output list
 830     AliFlowCommonHist *pCommonHistSP = dynamic_cast<AliFlowCommonHist*>
 831       (outputListHistos->FindObject("AliFlowCommonHistSP"));
 832     AliFlowCommonHistResults *pCommonHistResultsSP = dynamic_cast<AliFlowCommonHistResults*>
 833       (outputListHistos->FindObject("AliFlowCommonHistResultsSP"));
 834     AliFlowCommonHist *pCommonHistmuQ = dynamic_cast<AliFlowCommonHist*>
 835       (outputListHistos->FindObject("AliFlowCommonHistmuQ"));
 836
 837     TProfile* pHistProQNorm    = dynamic_cast<TProfile*>(outputListHistos->FindObject("FlowPro_QNorm_SP"));
 838     TProfile* pHistProQaQb     = dynamic_cast<TProfile*>(outputListHistos->FindObject("FlowPro_QaQb_SP"));
 839     TProfile* pHistProQaQbNorm = dynamic_cast<TProfile*>(outputListHistos->FindObject("FlowPro_QaQbNorm_SP"));
 840     TProfile* pHistProQaQbReImNorm = dynamic_cast<TProfile*>(outputListHistos->FindObject("FlowPro_QaQbReImNorm_SP"));
 841     TProfile* pHistProNonIsotropicTermsQ = dynamic_cast<TProfile*>(outputListHistos->FindObject("FlowPro_NonIsotropicTermsQ_SP"));
 842     TH1D*     pHistSumOfLinearWeights    = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_SumOfLinearWeights_SP"));
 843     TH1D*     pHistSumOfQuadraticWeights = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_SumOfQuadraticWeights_SP"));
 844
 845     TProfile* pHistProFlags    = dynamic_cast<TProfile*>(outputListHistos->FindObject("FlowPro_Flags_SP"));
 846     TProfile* pHistProUQetaRP  = dynamic_cast<TProfile*>(outputListHistos->FindObject("FlowPro_UQetaRP_SP"));
 847     TProfile* pHistProUQetaPOI = dynamic_cast<TProfile*>(outputListHistos->FindObject("FlowPro_UQetaPOI_SP"));
 848     TProfile* pHistProUQPtRP   = dynamic_cast<TProfile*>(outputListHistos->FindObject("FlowPro_UQPtRP_SP"));
 849     TProfile* pHistProUQPtPOI  = dynamic_cast<TProfile*>(outputListHistos->FindObject("FlowPro_UQPtPOI_SP"));
 850     TProfile* pHistProUQQaQbPtRP    = dynamic_cast<TProfile*>(outputListHistos->FindObject("FlowPro_UQQaQbPtRP_SP"));
 851     TProfile* pHistProUQQaQbEtaRP   = dynamic_cast<TProfile*>(outputListHistos->FindObject("FlowPro_UQQaQbEtaRP_SP"));
 852     TProfile* pHistProUQQaQbPtPOI   = dynamic_cast<TProfile*>(outputListHistos->FindObject("FlowPro_UQQaQbPtPOI_SP"));
 853     TProfile* pHistProUQQaQbEtaPOI  = dynamic_cast<TProfile*>(outputListHistos->FindObject("FlowPro_UQQaQbEtaPOI_SP"));
 854     TString weightFlag[3] = {"w_Qu_","w_Qu^2_","w_QuQaQb_"};
 855
 856
 857     TH1D* pHistSumOfWeightsPtRP[3] = {NULL};
 858     TH1D* pHistSumOfWeightsEtaRP[3] = {NULL};
 859     TH1D* pHistSumOfWeightsPtPOI[3] = {NULL};
 860     TH1D* pHistSumOfWeightsEtaPOI[3] = {NULL};
 861
 862     for(Int_t i=0;i<3;i++) {
 863       pHistSumOfWeightsPtRP[i]   = dynamic_cast<TH1D*>(outputListHistos->FindObject(Form("Flow_SumOfWeights%sPtRP_SP",weightFlag[i].Data())));
 864       pHistSumOfWeightsEtaRP[i]  = dynamic_cast<TH1D*>(outputListHistos->FindObject(Form("Flow_SumOfWeights%sEtaRP_SP",weightFlag[i].Data())));
 865       pHistSumOfWeightsPtPOI[i]  = dynamic_cast<TH1D*>(outputListHistos->FindObject(Form("Flow_SumOfWeights%sPtPOI_SP",weightFlag[i].Data())));
 866       pHistSumOfWeightsEtaPOI[i] = dynamic_cast<TH1D*>(outputListHistos->FindObject(Form("Flow_SumOfWeights%sEtaPOI_SP",weightFlag[i].Data())));
 867     }
 868
 869     TString rpPoi[2] = {"RP","POI"};
 870     TString ptEta[2] = {"Pt","Eta"};
 871     TString sinCos[2] = {"sin","cos"};
 872     TProfile *pHistProNonIsotropicTermsU[2][2][2] = {{{NULL}}};
 873     for(Int_t rp=0;rp<2;rp++) {
 874       for(Int_t pe=0;pe<2;pe++) {
 875         for(Int_t sc=0;sc<2;sc++) {
 876           pHistProNonIsotropicTermsU[rp][pe][sc] = dynamic_cast<TProfile*>(outputListHistos->FindObject(Form("FlowPro_NonIsotropicTerms_%s_%s_%s_SP",rpPoi[rp].Data(),ptEta[pe].Data(),sinCos[sc].Data())));
 877         }
 878       }
 879     }
 880
 881     TH1D*     pHistQNorm    = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_QNorm_SP"));
 882     TH1D*     pHistQaQb     = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_QaQb_SP"));
 883     TH1D*     pHistQaQbNorm = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_QaQbNorm_SP"));
 884     TH2D*     pHistQNormvsQaQbNorm = dynamic_cast<TH2D*>(outputListHistos->FindObject("Flow_QNormvsQaQbNorm_SP"));
 885     TH1D*     pHistQaQbCos  = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_QaQbCos_SP"));
 886     TH1D*     pHistResolution = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_resolution_SP"));
 887     TH1D*     pHistQaNorm   = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_QaNorm_SP"));
 888     TH2D*     pHistQaNormvsMa   = dynamic_cast<TH2D*>(outputListHistos->FindObject("Flow_QaNormvsMa_SP"));
 889     TH1D*     pHistQbNorm   = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_QbNorm_SP"));
 890     TH2D*     pHistQbNormvsMb   = dynamic_cast<TH2D*>(outputListHistos->FindObject("Flow_QbNormvsMb_SP"));
 891     TH2D*     pHistMavsMb = dynamic_cast<TH2D*>(outputListHistos->FindObject("Flow_MavsMb_SP"));
 892
 893     //pass the pointers to the task
 894     if (pCommonHistSP &&
 895         pCommonHistResultsSP &&
 896         pCommonHistmuQ &&
 897         pHistProQNorm &&
 898         pHistProQaQb &&
 899         pHistProQaQbNorm &&
 900         pHistProQaQbReImNorm &&
 901         pHistProNonIsotropicTermsQ &&
 902         pHistSumOfLinearWeights &&
 903         pHistSumOfQuadraticWeights &&
 904         pHistProFlags &&
 905         pHistProUQetaRP &&
 906         pHistProUQetaPOI &&
 907         pHistProUQPtRP &&
 908         pHistProUQPtPOI &&
 909         pHistProUQQaQbPtRP &&
 910         pHistProUQQaQbEtaRP &&
 911         pHistProUQQaQbPtPOI &&
 912         pHistProUQQaQbEtaPOI &&
 913         pHistSumOfWeightsPtRP[0] && pHistSumOfWeightsPtRP[1] && pHistSumOfWeightsPtRP[2] &&
 914         pHistSumOfWeightsEtaRP[0] && pHistSumOfWeightsEtaRP[1] && pHistSumOfWeightsEtaRP[2] &&
 915         pHistSumOfWeightsPtPOI[0] && pHistSumOfWeightsPtPOI[1] && pHistSumOfWeightsPtPOI[2] &&
 916         pHistSumOfWeightsEtaPOI[0] && pHistSumOfWeightsEtaPOI[1] && pHistSumOfWeightsEtaPOI[2] &&
 917         pHistProNonIsotropicTermsU[0][0][0] && pHistProNonIsotropicTermsU[1][0][0] && pHistProNonIsotropicTermsU[0][1][0] && pHistProNonIsotropicTermsU[0][0][1] &&
 918         pHistProNonIsotropicTermsU[1][1][0] && pHistProNonIsotropicTermsU[1][0][1] && pHistProNonIsotropicTermsU[0][1][1] && pHistProNonIsotropicTermsU[1][1][1] &&
 919         pHistQNorm &&
 920         pHistQaQb &&
 921         pHistQaQbNorm &&
 922         pHistQNormvsQaQbNorm &&
 923         pHistQaQbCos &&
 924         pHistResolution &&
 925         pHistQaNorm &&
 926         pHistQaNormvsMa &&
 927         pHistQbNorm &&
 928         pHistQbNormvsMb &&
 929         pHistMavsMb
 930         ) {
 931
 932       this -> SetCommonHistsSP(pCommonHistSP);
 933       this -> SetCommonHistsResSP(pCommonHistResultsSP);
 934       this -> SetCommonHistsmuQ(pCommonHistmuQ);
 935       this -> SetHistProQNorm(pHistProQNorm);
 936       this -> SetHistProQaQb(pHistProQaQb);
 937       this -> SetHistProQaQbNorm(pHistProQaQbNorm);
 938       this -> SetHistProQaQbReImNorm(pHistProQaQbReImNorm);
 939       this -> SetHistProNonIsotropicTermsQ(pHistProNonIsotropicTermsQ);
 940       this -> SetHistSumOfLinearWeights(pHistSumOfLinearWeights);
 941       this -> SetHistSumOfQuadraticWeights(pHistSumOfQuadraticWeights);
 942       this -> SetHistProFlags(pHistProFlags);
 943       this -> SetHistProUQetaRP(pHistProUQetaRP);
 944       this -> SetHistProUQetaPOI(pHistProUQetaPOI);
 945       this -> SetHistProUQPtRP(pHistProUQPtRP);
 946       this -> SetHistProUQPtPOI(pHistProUQPtPOI);
 947       this -> SetHistProUQQaQbPtRP(pHistProUQQaQbPtRP);
 948       this -> SetHistProUQQaQbEtaRP(pHistProUQQaQbEtaRP);
 949       this -> SetHistProUQQaQbPtPOI(pHistProUQQaQbPtPOI);
 950       this -> SetHistProUQQaQbEtaPOI(pHistProUQQaQbEtaPOI);
 951       for(Int_t i=0;i<3;i++) {
 952         if(pHistSumOfWeightsPtRP[i]) this -> SetHistSumOfWeightsPtRP(pHistSumOfWeightsPtRP[i],i);
 953         if(pHistSumOfWeightsEtaRP[i]) this -> SetHistSumOfWeightsEtaRP(pHistSumOfWeightsEtaRP[i],i);
 954         if(pHistSumOfWeightsPtPOI[i]) this -> SetHistSumOfWeightsPtPOI(pHistSumOfWeightsPtPOI[i],i);
 955         if(pHistSumOfWeightsEtaPOI[i]) this -> SetHistSumOfWeightsEtaPOI(pHistSumOfWeightsEtaPOI[i],i);
 956       }
 957       for(Int_t rp=0;rp<2;rp++)  {
 958         for(Int_t pe=0;pe<2;pe++) {
 959           for(Int_t sc=0;sc<2;sc++) {
 960             if(pHistProNonIsotropicTermsU[rp][pe][sc]) {
 961               this->SetHistProNonIsotropicTermsU(pHistProNonIsotropicTermsU[rp][pe][sc],rp,pe,sc);
 962             }
 963           }
 964         }
 965       }
 966       this -> SetHistQNorm(pHistQNorm);
 967       this -> SetHistQaQb(pHistQaQb);
 968       this -> SetHistQaQbNorm(pHistQaQbNorm);
 969       this -> SetHistQNormvsQaQbNorm(pHistQNormvsQaQbNorm);
 970       this -> SetHistQaQbCos(pHistQaQbCos);
 971       this -> SetHistResolution(pHistResolution);
 972       this -> SetHistQaNorm(pHistQaNorm);
 973       this -> SetHistQaNormvsMa(pHistQaNormvsMa);
 974       this -> SetHistQbNorm(pHistQbNorm);
 975       this -> SetHistQbNormvsMb(pHistQbNormvsMb);
 976       this -> SetHistMavsMb(pHistMavsMb);
 977
 978     } else {
 979       cout<<"WARNING: Histograms needed to run Finish() in SP are not accessable!"<<endl; }
 980
 981   } // end of if(outputListHistos)
 982 }
 983
 984 //--------------------------------------------------------------------
 985 void AliFlowAnalysisWithScalarProduct::Finish() {
 986
 987   //calculate flow and fill the AliFlowCommonHistResults
 988   if (fDebug) cout<<"AliFlowAnalysisWithScalarProduct::Finish()"<<endl;
 989
 990   //access all boolean flags needed in Finish():
 991   this->AccessFlags();
 992
 993   cout<<"*************************************"<<endl;
 994   cout<<"*************************************"<<endl;
 995   cout<<"      Integrated flow from           "<<endl;
 996   cout<<"         Scalar product              "<<endl;
 997   cout<<endl;
 998
 999   Int_t iNbinsPt  = AliFlowCommonConstants::GetMaster()->GetNbinsPt();
1000   Int_t iNbinsEta = AliFlowCommonConstants::GetMaster()->GetNbinsEta();
1001
1002   //Calculate the event plane resolution
1003   //----------------------------------
1004   Double_t dCos2phi = fHistResolution->GetMean();
1005   if (dCos2phi > 0.0){
1006     Double_t dResolution = TMath::Sqrt(2*dCos2phi);
1007     cout<<"An estimate of the event plane resolution is: "<<dResolution<<endl;
1008     cout<<endl;
1009   }
1010
1011   //Calculate reference flow (noname)
1012   //----------------------------------
1013   //weighted average over (QaQb/MaMb) with weight (MaMb)
1014   Double_t dQaQb  = fHistProQaQbNorm->GetBinContent(1);
1015   Double_t dSpreadQaQb = fHistProQaQbNorm->GetBinError(1);
1016   Double_t dEntriesQaQb = fHistProQaQbNorm->GetEntries();
1017
1018   //non-isotropic terms:
1019   Double_t dImQa = fHistProQaQbReImNorm->GetBinContent(1);
1020   Double_t dReQa = fHistProQaQbReImNorm->GetBinContent(2);
1021   Double_t dImQb = fHistProQaQbReImNorm->GetBinContent(3);
1022   Double_t dReQb = fHistProQaQbReImNorm->GetBinContent(4);
1023
1024   if(fApplyCorrectionForNUA)
1025   {
1026    dQaQb = dQaQb - dImQa*dImQb - dReQa*dReQb;
1027   }
1028
1029   if (dEntriesQaQb > 0.) {
1030     cout<<"QaQb = "<<dQaQb<<" +- "<<(dSpreadQaQb/TMath::Sqrt(dEntriesQaQb))<<endl;
1031     cout<<endl;
1032   }
1033
1034   Double_t dV = -999.;
1035   if(dQaQb>=0.)
1036   {
1037    dV = TMath::Sqrt(dQaQb);
1038   }
1039   //statistical error of dQaQb:
1040   //  statistical error = term1 * spread * term2:
1041   //  term1 = sqrt{sum_{i=1}^{N} w^2}/(sum_{i=1}^{N} w)
1042   //  term2 = 1/sqrt(1-term1^2)
1043   Double_t dSumOfLinearWeights = fHistSumOfLinearWeights->GetBinContent(1);
1044   Double_t dSumOfQuadraticWeights = fHistSumOfQuadraticWeights->GetBinContent(1);
1045   Double_t dTerm1 = 0.;
1046   Double_t dTerm2 = 0.;
1047   if(dSumOfLinearWeights) {
1048     dTerm1 = pow(dSumOfQuadraticWeights,0.5)/dSumOfLinearWeights;
1049   }
1050   if(1.-pow(dTerm1,2.) > 0.) {
1051     dTerm2 = 1./pow(1-pow(dTerm1,2.),0.5);
1052   }
1053   Double_t dStatErrorQaQb = dTerm1 * dSpreadQaQb * dTerm2;
1054   //calculate the statistical error
1055   Double_t dVerr = 0.;
1056   if(dQaQb > 0.) {
1057     dVerr = (1./(2.*pow(dQaQb,0.5)))*dStatErrorQaQb;
1058   }
1059   fCommonHistsResSP->FillIntegratedFlow(dV,dVerr);
1060   cout<<"v2(subevents) = "<<dV<<" +- "<<dVerr<<endl;
1061
1062   //Calculate differential flow and integrated flow (RP, POI)
1063   //---------------------------------------------------------
1064   //v as a function of eta for RP selection
1065   for(Int_t b=1;b<iNbinsEta+1;b++) {
1066     Double_t duQpro = fHistProUQetaRP->GetBinContent(b);
1067     if(fApplyCorrectionForNUA) {
1068       duQpro = duQpro
1069         - fHistProNonIsotropicTermsU[0][1][1]->GetBinContent(b)*fHistProNonIsotropicTermsQ->GetBinContent(2)
1070         - fHistProNonIsotropicTermsU[0][1][0]->GetBinContent(b)*fHistProNonIsotropicTermsQ->GetBinContent(1);
1071     }
1072     Double_t dv2pro = -999.;
1073     if (dV!=0.) { dv2pro = duQpro/dV; }
1074     //calculate the statistical error
1075     Double_t dv2ProErr = CalculateStatisticalError(b, dStatErrorQaQb, fHistProUQetaRP, fHistProUQQaQbEtaRP, fHistSumOfWeightsEtaRP);
1076     //fill TH1D
1077     fCommonHistsResSP->FillDifferentialFlowEtaRP(b, dv2pro, dv2ProErr);
1078   } //loop over bins b
1079
1080
1081   //v as a function of eta for POI selection
1082   for(Int_t b=1;b<iNbinsEta+1;b++) {
1083     Double_t duQpro = fHistProUQetaPOI->GetBinContent(b);
1084     if(fApplyCorrectionForNUA)  {
1085       duQpro = duQpro
1086         - fHistProNonIsotropicTermsU[1][1][1]->GetBinContent(b)*fHistProNonIsotropicTermsQ->GetBinContent(2)
1087         - fHistProNonIsotropicTermsU[1][1][0]->GetBinContent(b)*fHistProNonIsotropicTermsQ->GetBinContent(1);
1088     }
1089     Double_t dv2pro = -999.;
1090     if (dV!=0.) { dv2pro = duQpro/dV; }
1091     //calculate the statistical error
1092     Double_t dv2ProErr = CalculateStatisticalError(b, dStatErrorQaQb, fHistProUQetaPOI, fHistProUQQaQbEtaPOI, fHistSumOfWeightsEtaPOI);
1093
1094     //fill TH1D
1095     fCommonHistsResSP->FillDifferentialFlowEtaPOI(b, dv2pro, dv2ProErr);
1096   } //loop over bins b
1097
1098
1099   //v as a function of Pt for RP selection
1100   TH1F* fHistPtRP = fCommonHistsSP->GetHistPtRP(); //for calculating integrated flow
1101   Double_t dVRP = 0.;
1102   Double_t dSumRP = 0.;
1103   Double_t dErrVRP =0.;
1104
1105   for(Int_t b=1;b<iNbinsPt+1;b++) {
1106     Double_t duQpro = fHistProUQPtRP->GetBinContent(b);
1107     if(fApplyCorrectionForNUA) {
1108       duQpro = duQpro
1109         - fHistProNonIsotropicTermsU[0][0][1]->GetBinContent(b)*fHistProNonIsotropicTermsQ->GetBinContent(2)
1110         - fHistProNonIsotropicTermsU[0][0][0]->GetBinContent(b)*fHistProNonIsotropicTermsQ->GetBinContent(1);
1111     }
1112     Double_t dv2pro = -999.;
1113     if (dV!=0.) { dv2pro = duQpro/dV; }
1114     //calculate the statistical error
1115     Double_t dv2ProErr = CalculateStatisticalError(b, dStatErrorQaQb, fHistProUQPtRP, fHistProUQQaQbPtRP, fHistSumOfWeightsPtRP);
1116
1117     //fill TH1D
1118     fCommonHistsResSP->FillDifferentialFlowPtRP(b, dv2pro, dv2ProErr);
1119
1120     //calculate integrated flow for RP selection
1121     if (fHistPtRP){
1122       Double_t dYieldPt = fHistPtRP->GetBinContent(b);
1123       dVRP += dv2pro*dYieldPt;
1124       dSumRP +=dYieldPt;
1125       dErrVRP += dYieldPt*dYieldPt*dv2ProErr*dv2ProErr;
1126     } else { cout<<"fHistPtRP is NULL"<<endl; }
1127   } //loop over bins b
1128
1129   if (dSumRP != 0.) {
1130     dVRP /= dSumRP; //the pt distribution should be normalised
1131     dErrVRP /= (dSumRP*dSumRP);
1132     dErrVRP = TMath::Sqrt(dErrVRP);
1133   }
1134   fCommonHistsResSP->FillIntegratedFlowRP(dVRP,dErrVRP);
1135   cout<<"v2(RP) = "<<dVRP<<" +- "<<dErrVRP<<endl;
1136
1137
1138   //v as a function of Pt for POI selection
1139   TH1F* fHistPtPOI = fCommonHistsSP->GetHistPtPOI(); //for calculating integrated flow
1140   Double_t dVPOI = 0.;
1141   Double_t dSumPOI = 0.;
1142   Double_t dErrVPOI =0.;
1143
1144   for(Int_t b=1;b<iNbinsPt+1;b++) {
1145     Double_t duQpro = fHistProUQPtPOI->GetBinContent(b);
1146     if(fApplyCorrectionForNUA)  {
1147      duQpro = duQpro
1148        - fHistProNonIsotropicTermsU[1][0][1]->GetBinContent(b)*fHistProNonIsotropicTermsQ->GetBinContent(2)
1149        - fHistProNonIsotropicTermsU[1][0][0]->GetBinContent(b)*fHistProNonIsotropicTermsQ->GetBinContent(1);
1150     }
1151     Double_t dv2pro = -999.;
1152     if (dV!=0.) { dv2pro = duQpro/dV; }
1153     //calculate the statistical error
1154     Double_t dv2ProErr = CalculateStatisticalError(b, dStatErrorQaQb, fHistProUQPtPOI, fHistProUQQaQbPtPOI, fHistSumOfWeightsPtPOI);
1155
1156     //fill TH1D
1157     fCommonHistsResSP->FillDifferentialFlowPtPOI(b, dv2pro, dv2ProErr);
1158
1159     //calculate integrated flow for POI selection
1160     if (fHistPtPOI){
1161       Double_t dYieldPt = fHistPtPOI->GetBinContent(b);
1162       dVPOI += dv2pro*dYieldPt;
1163       dSumPOI +=dYieldPt;
1164       dErrVPOI += dYieldPt*dYieldPt*dv2ProErr*dv2ProErr;
1165     } else { cout<<"fHistPtPOI is NULL"<<endl; }
1166   } //loop over bins b
1167
1168   if (dSumPOI > 0.) {
1169     dVPOI /= dSumPOI; //the pt distribution should be normalised
1170     dErrVPOI /= (dSumPOI*dSumPOI);
1171     dErrVPOI = TMath::Sqrt(dErrVPOI);
1172   }
1173   fCommonHistsResSP->FillIntegratedFlowPOI(dVPOI,dErrVPOI);
1174   cout<<"v2(POI) = "<<dVPOI<<" +- "<<dErrVPOI<<endl;
1175
1176   cout<<endl;
1177   cout<<"*************************************"<<endl;
1178   cout<<"*************************************"<<endl;
1179
1180   //cout<<".....finished"<<endl;
1181 }
1182
1183
1184 //--------------------------------------------------------------------
1185 Double_t AliFlowAnalysisWithScalarProduct::CalculateStatisticalError(Int_t b, Double_t aStatErrorQaQb, TProfile* pHistProUQ, TProfile* pHistProUQQaQb, TH1D** pHistSumOfWeights) {
1186   //calculate the statistical error for differential flow for bin b
1187   Double_t duQproSpread = pHistProUQ->GetBinError(b);
1188   Double_t sumOfMq = pHistSumOfWeights[0]->GetBinContent(b);
1189   Double_t sumOfMqSquared = pHistSumOfWeights[1]->GetBinContent(b);
1190   Double_t dQaQb = fHistProQaQbNorm->GetBinContent(1);
1191   Double_t dTerm1 = 0.;
1192   Double_t dTerm2 = 0.;
1193   if(sumOfMq) {
1194     dTerm1 = (pow(sumOfMqSquared,0.5)/sumOfMq);
1195   }
1196   if(1.-pow(dTerm1,2.)>0.) {
1197     dTerm2 = 1./pow(1.-pow(dTerm1,2.),0.5);
1198   }
1199   Double_t duQproErr = dTerm1*duQproSpread*dTerm2;
1200   // covariances:
1201   Double_t dTerm1Cov = pHistSumOfWeights[2]->GetBinContent(b);
1202   Double_t dTerm2Cov = fHistSumOfLinearWeights->GetBinContent(1);
1203   Double_t dTerm3Cov = sumOfMq;
1204   Double_t dWeightedCovariance = 0.;
1205   if(dTerm2Cov*dTerm3Cov>0.) {
1206     Double_t dDenominator = 1.-dTerm1Cov/(dTerm2Cov*dTerm3Cov);
1207     Double_t dPrefactor = dTerm1Cov/(dTerm2Cov*dTerm3Cov);
1208     if(dDenominator!=0) {
1209       Double_t dCovariance = (pHistProUQQaQb->GetBinContent(b)-dQaQb*pHistProUQ->GetBinContent(b))/dDenominator;
1210       dWeightedCovariance = dCovariance*dPrefactor;
1211     }
1212   }
1213
1214   Double_t dv2ProErr = 0.; // final statitical error:
1215   if(dQaQb>0.) {
1216     Double_t dv2ProErrorSquared = (1./4.)*pow(dQaQb,-3.)*
1217       (pow(pHistProUQ->GetBinContent(b),2.)*pow(aStatErrorQaQb,2.)
1218        + 4.*pow(dQaQb,2.)*pow(duQproErr,2.)
1219        - 4.*dQaQb*pHistProUQ->GetBinContent(b)*dWeightedCovariance);
1220     if(dv2ProErrorSquared>0.) dv2ProErr = pow(dv2ProErrorSquared,0.5);
1221   }
1222
1223   return dv2ProErr;
1224 }
1225
1226
1227 //--------------------------------------------------------------------
1228
1229 void AliFlowAnalysisWithScalarProduct::StoreFlags()
1230 {
1231  // Store all boolean flags needed in Finish() in profile fHistProFlags.
1232
1233  // Apply correction for non-uniform acceptance or not:
1234  fHistProFlags->Fill(0.5,fApplyCorrectionForNUA);
1235
1236 }
1237
1238 //--------------------------------------------------------------------
1239
1240 void AliFlowAnalysisWithScalarProduct::AccessFlags()
1241 {
1242  // Access all boolean flags needed in Finish() from profile fHistProFlags.
1243
1244  // Apply correction for non-uniform acceptance or not:
1245  fApplyCorrectionForNUA = (Bool_t) fHistProFlags->GetBinContent(1);
1246
1247 }
1248
1249 //--------------------------------------------------------------------