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