PWG2/FLOW/AliFlowCommon/AliFlowAnalysisWithScalarProduct.cxx

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