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