PWG2/FLOW/AliFlowCommon/AliFlowAnalysisWithScalarProduct.cxx

   1 /*************************************************************************
   2 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15
  16 #define AliFlowAnalysisWithScalarProduct_cxx
  17
  18 #include "Riostream.h"  //needed as include
  19 #include "TFile.h"      //needed as include
  20 #include "TList.h"
  21 #include "TMath.h"
  22 #include "TProfile.h"
  23 #include "TVector2.h"
  24
  25 class TH1F;
  26
  27 #include "AliFlowCommonConstants.h"    //needed as include
  28 #include "AliFlowEventSimple.h"
  29 #include "AliFlowTrackSimple.h"
  30 #include "AliFlowCommonHist.h"
  31 #include "AliFlowCommonHistResults.h"
  32 #include "AliFlowAnalysisWithScalarProduct.h"
  33
  34 class AliFlowVector;
  35
  36 //////////////////////////////////////////////////////////////////////////////
  37 // AliFlowAnalysisWithScalarProduct:
  38 // Description:
  39 // Maker to analyze Flow with the Scalar product method.
  40 //
  41 // authors: N. van der Kolk (kolk@nikhef.nl), A. Bilandzic (anteb@nikhef.nl)
  42 //////////////////////////////////////////////////////////////////////////////
  43
  44 ClassImp(AliFlowAnalysisWithScalarProduct)
  45
  46   //-----------------------------------------------------------------------
  47   AliFlowAnalysisWithScalarProduct::AliFlowAnalysisWithScalarProduct():
  48    fEventNumber(0),
  49    fDebug(kFALSE),
  50    fWeightsList(NULL),
  51    fUsePhiWeights(kFALSE),
  52    fPhiWeights(NULL),
  53    fHistList(NULL),
  54    fHistProUQetaRP(NULL),
  55    fHistProUQetaPOI(NULL),
  56    fHistProUQPtRP(NULL),
  57    fHistProUQPtPOI(NULL),
  58    fHistProQaQb(NULL),
  59    fHistSumOfLinearWeights(NULL),
  60    fHistSumOfQuadraticWeights(NULL),
  61    fHistProUQQaQbPtRP(NULL),
  62    fHistProUQQaQbEtaRP(NULL),
  63    fHistProUQQaQbPtPOI(NULL),
  64    fHistProUQQaQbEtaPOI(NULL),
  65    fCommonHists(NULL),
  66    fCommonHistsRes(NULL)
  67 {
  68   // Constructor.
  69   fWeightsList = new TList();
  70   fHistList = new TList();
  71
  72   // Initialize arrays:
  73   for(Int_t i=0;i<3;i++)
  74   {
  75    fHistSumOfWeightsPtRP[i] = NULL;
  76    fHistSumOfWeightsEtaRP[i] = NULL;
  77    fHistSumOfWeightsPtPOI[i] = NULL;
  78    fHistSumOfWeightsEtaPOI[i] = NULL;
  79   }
  80 }
  81  //-----------------------------------------------------------------------
  82  AliFlowAnalysisWithScalarProduct::~AliFlowAnalysisWithScalarProduct()
  83  {
  84    //destructor
  85    delete fWeightsList;
  86    delete fHistList;
  87  }
  88
  89 //-----------------------------------------------------------------------
  90 void AliFlowAnalysisWithScalarProduct::WriteHistograms(TString* outputFileName)
  91 {
  92  //store the final results in output .root file
  93
  94   TFile *output = new TFile(outputFileName->Data(),"RECREATE");
  95   //output->WriteObject(fHistList, "cobjSP","SingleKey");
  96   fHistList->SetName("cobjSP");
  97   fHistList->SetOwner(kTRUE);
  98   fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
  99   delete output;
 100 }
 101
 102 //-----------------------------------------------------------------------
 103 void AliFlowAnalysisWithScalarProduct::WriteHistograms(TString outputFileName)
 104 {
 105  //store the final results in output .root file
 106
 107   TFile *output = new TFile(outputFileName.Data(),"RECREATE");
 108   //output->WriteObject(fHistList, "cobjSP","SingleKey");
 109   fHistList->SetName("cobjSP");
 110   fHistList->SetOwner(kTRUE);
 111   fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
 112   delete output;
 113 }
 114
 115 //-----------------------------------------------------------------------
 116 void AliFlowAnalysisWithScalarProduct::WriteHistograms(TDirectoryFile *outputFileName)
 117 {
 118  //store the final results in output .root file
 119  fHistList->SetName("cobjSP");
 120  fHistList->SetOwner(kTRUE);
 121  outputFileName->Add(fHistList);
 122  outputFileName->Write(outputFileName->GetName(), TObject::kSingleKey);
 123 }
 124
 125 //-----------------------------------------------------------------------
 126 void AliFlowAnalysisWithScalarProduct::Init() {
 127
 128   //Define all histograms
 129   cout<<"---Analysis with the Scalar Product Method--- Init"<<endl;
 130
 131   Int_t iNbinsPt   = AliFlowCommonConstants::GetMaster()->GetNbinsPt();
 132   Double_t dPtMin  = AliFlowCommonConstants::GetMaster()->GetPtMin();
 133   Double_t dPtMax  = AliFlowCommonConstants::GetMaster()->GetPtMax();
 134   Int_t iNbinsEta  = AliFlowCommonConstants::GetMaster()->GetNbinsEta();
 135   Double_t dEtaMin = AliFlowCommonConstants::GetMaster()->GetEtaMin();
 136   Double_t dEtaMax = AliFlowCommonConstants::GetMaster()->GetEtaMax();
 137
 138   fHistProUQetaRP = new TProfile("Flow_UQetaRP_SP","Flow_UQetaRP_SP",iNbinsEta,dEtaMin,dEtaMax,"s");
 139   fHistProUQetaRP->SetXTitle("{eta}");
 140   fHistProUQetaRP->SetYTitle("<uQ>");
 141   fHistList->Add(fHistProUQetaRP);
 142
 143   fHistProUQetaPOI = new TProfile("Flow_UQetaPOI_SP","Flow_UQetaPOI_SP",iNbinsEta,dEtaMin,dEtaMax,"s");
 144   fHistProUQetaPOI->SetXTitle("{eta}");
 145   fHistProUQetaPOI->SetYTitle("<uQ>");
 146   fHistList->Add(fHistProUQetaPOI);
 147
 148   fHistProUQPtRP = new TProfile("Flow_UQPtRP_SP","Flow_UQPtRP_SP",iNbinsPt,dPtMin,dPtMax,"s");
 149   fHistProUQPtRP->SetXTitle("p_t (GeV)");
 150   fHistProUQPtRP->SetYTitle("<uQ>");
 151   fHistList->Add(fHistProUQPtRP);
 152
 153   fHistProUQPtPOI = new TProfile("Flow_UQPtPOI_SP","Flow_UQPtPOI_SP",iNbinsPt,dPtMin,dPtMax,"s");
 154   fHistProUQPtPOI->SetXTitle("p_t (GeV)");
 155   fHistProUQPtPOI->SetYTitle("<uQ>");
 156   fHistList->Add(fHistProUQPtPOI);
 157
 158   fHistProQaQb = new TProfile("Flow_QaQb_SP","Flow_QaQb_SP", 1, 0.5, 1.5,"s");
 159   fHistProQaQb->SetYTitle("<QaQb>");
 160   fHistList->Add(fHistProQaQb);
 161
 162   fHistSumOfLinearWeights = new TH1D("Flow_SumOfLinearWeights_SP","Flow_SumOfLinearWeights_SP",1,-0.5, 0.5);
 163   fHistSumOfLinearWeights -> SetYTitle("sum (*)");
 164   fHistSumOfLinearWeights -> SetXTitle("sum (Ma*Mb)");
 165   fHistList->Add(fHistSumOfLinearWeights);
 166
 167   fHistSumOfQuadraticWeights = new TH1D("Flow_SumOfQuadraticWeights_SP","Flow_SumOfQuadraticWeights_SP",1,-0.5, 0.5);
 168   fHistSumOfQuadraticWeights -> SetYTitle("sum (*)");
 169   fHistSumOfQuadraticWeights -> SetXTitle("sum (Ma*Mb)^2");
 170   fHistList->Add(fHistSumOfQuadraticWeights);
 171
 172   fHistProUQQaQbPtRP = new TProfile("Flow_UQQaQbPtRP_SP","Flow_UQQaQbPtRP_SP",iNbinsPt,dPtMin,dPtMax);
 173   fHistProUQQaQbPtRP -> SetYTitle("<*>");
 174   fHistProUQQaQbPtRP -> SetXTitle("<Qu QaQb>");
 175   fHistList->Add(fHistProUQQaQbPtRP);
 176
 177   fHistProUQQaQbEtaRP = new TProfile("Flow_UQQaQbEtaRP_SP","Flow_UQQaQbEtaRP_SP",iNbinsEta,dEtaMin,dEtaMax);
 178   fHistProUQQaQbEtaRP -> SetYTitle("<*>");
 179   fHistProUQQaQbEtaRP -> SetXTitle("<Qu QaQb>");
 180   fHistList->Add(fHistProUQQaQbEtaRP);
 181
 182   fHistProUQQaQbPtPOI = new TProfile("Flow_UQQaQbPtPOI_SP","Flow_UQQaQbPtPOI_SP",iNbinsPt,dPtMin,dPtMax);
 183   fHistProUQQaQbPtPOI -> SetYTitle("<*>");
 184   fHistProUQQaQbPtPOI -> SetXTitle("<Qu QaQb>");
 185   fHistList->Add(fHistProUQQaQbPtPOI);
 186
 187   fHistProUQQaQbEtaPOI = new TProfile("Flow_UQQaQbEtaPOI_SP","Flow_UQQaQbEtaPOI_SP",iNbinsEta,dEtaMin,dEtaMax);
 188   fHistProUQQaQbEtaPOI -> SetYTitle("<*>");
 189   fHistProUQQaQbEtaPOI -> SetXTitle("<Qu QaQb>");
 190   fHistList->Add(fHistProUQQaQbEtaPOI);
 191
 192   TString weightFlag[3] = {"w_Qu_","w_Qu^2_","w_QuQaQb_"};
 193   for(Int_t i=0;i<3;i++)
 194   {
 195    fHistSumOfWeightsPtRP[i] = new TH1D(Form("Flow_SumOfWeights%sPtRP_SP",weightFlag[i].Data()),
 196                               Form("Flow_SumOfWeights%sPtRP_SP",weightFlag[i].Data()),iNbinsPt,dPtMin,dPtMax);
 197    fHistSumOfWeightsPtRP[i] -> SetYTitle("sum (*)");
 198    fHistSumOfWeightsPtRP[i] -> SetXTitle("p_{T}");
 199    fHistList->Add(fHistSumOfWeightsPtRP[i]);
 200
 201    fHistSumOfWeightsEtaRP[i] = new TH1D(Form("Flow_SumOfWeights%sEtaRP_SP",weightFlag[i].Data()),
 202                                Form("Flow_SumOfWeights%sEtaRP_SP",weightFlag[i].Data()),iNbinsEta,dEtaMin,dEtaMax);
 203    fHistSumOfWeightsEtaRP[i] -> SetYTitle("sum (*)");
 204    fHistSumOfWeightsEtaRP[i] -> SetXTitle("#eta");
 205    fHistList->Add(fHistSumOfWeightsEtaRP[i]);
 206
 207    fHistSumOfWeightsPtPOI[i] = new TH1D(Form("Flow_SumOfWeights%sPtPOI_SP",weightFlag[i].Data()),
 208                                Form("Flow_SumOfWeights%sPtPOI_SP",weightFlag[i].Data()),iNbinsPt,dPtMin,dPtMax);
 209    fHistSumOfWeightsPtPOI[i] -> SetYTitle("sum (*)");
 210    fHistSumOfWeightsPtPOI[i] -> SetXTitle("p_{T}");
 211    fHistList->Add(fHistSumOfWeightsPtPOI[i]);
 212
 213    fHistSumOfWeightsEtaPOI[i] = new TH1D(Form("Flow_SumOfWeights%sEtaPOI_SP",weightFlag[i].Data()),
 214                                 Form("Flow_SumOfWeights%sEtaPOI_SP",weightFlag[i].Data()),iNbinsEta,dEtaMin,dEtaMax);
 215    fHistSumOfWeightsEtaPOI[i] -> SetYTitle("sum (*)");
 216    fHistSumOfWeightsEtaPOI[i] -> SetXTitle("#eta");
 217    fHistList->Add(fHistSumOfWeightsEtaPOI[i]);
 218   }
 219
 220   fCommonHists = new AliFlowCommonHist("AliFlowCommonHistSP");
 221   fHistList->Add(fCommonHists);
 222   fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsSP");
 223   fHistList->Add(fCommonHistsRes);
 224
 225   //weights
 226   if(fUsePhiWeights) {
 227     if(!fWeightsList) {
 228       cout<<"WARNING: fWeightsList is NULL in the Scalar Product method."<<endl;
 229       exit(0);
 230     }
 231     if(fWeightsList->FindObject("phi_weights"))  {
 232       fPhiWeights = dynamic_cast<TH1F*>(fWeightsList->FindObject("phi_weights"));
 233       fHistList->Add(fPhiWeights);
 234     } else {
 235       cout<<"WARNING: histogram with phi weights is not accessible in Scalar Product"<<endl;
 236       exit(0);
 237     }
 238   } // end of if(fUsePhiWeights)
 239
 240   fEventNumber = 0;  //set number of events to zero
 241 }
 242
 243 //-----------------------------------------------------------------------
 244 void AliFlowAnalysisWithScalarProduct::Make(AliFlowEventSimple* anEvent) {
 245
 246   //Calculate flow based on  <QaQb/MaMb> = <v^2>\r
 247
 248   //Fill histograms
 249   if (anEvent) {
 250
 251     //get Q vectors for the eta-subevents
 252     AliFlowVector* vQarray = new AliFlowVector[2];
 253     if (fUsePhiWeights) {
 254       anEvent->Get2Qsub(vQarray,2,fWeightsList,kTRUE);
 255     } else {
 256       anEvent->Get2Qsub(vQarray);
 257     }
 258     //subevent a
 259     AliFlowVector vQa = vQarray[0];
 260     //subevent b
 261     AliFlowVector vQb = vQarray[1];
 262
 263     //check that the subevents are not empty:
 264     Double_t dMa = vQa.GetMult();
 265     Double_t dMb = vQb.GetMult();
 266     if (dMa != 0 && dMb != 0) {
 267
 268       //fill control histograms
 269       fCommonHists->FillControlHistograms(anEvent);
 270
 271       //get total Q vector = the sum of subevent a and subevent b
 272       AliFlowVector vQ = vQa + vQb;
 273       //weight the Q vectors for the subevents by the multiplicity
 274       //Note: Weight Q only in the particle loop when it is clear if it should be (m-1) or M
 275       Double_t dQXa = vQa.X()/dMa;
 276       Double_t dQYa = vQa.Y()/dMa;
 277       vQa.Set(dQXa,dQYa);
 278
 279       Double_t dQXb = vQb.X()/dMb;
 280       Double_t dQYb = vQb.Y()/dMb;
 281       vQb.Set(dQXb,dQYb);
 282
 283       //scalar product of the two subevents
 284       Double_t dQaQb = (vQa*vQb);
 285       fHistProQaQb -> Fill(1.,dQaQb,dMa*dMb);  //Fill (QaQb/MaMb) with weight (MaMb).
 286       //needed for the error calculation:
 287       fHistSumOfLinearWeights -> Fill(0.,dMa*dMb);
 288       fHistSumOfQuadraticWeights -> Fill(0.,pow(dMa*dMb,2.));
 289
 290       //loop over the tracks of the event
 291       AliFlowTrackSimple*   pTrack = NULL;
 292       Int_t iNumberOfTracks = anEvent->NumberOfTracks();
 293       Double_t dMq =  vQ.GetMult();
 294
 295       for (Int_t i=0;i<iNumberOfTracks;i++)
 296         {
 297           pTrack = anEvent->GetTrack(i) ;
 298           if (pTrack){
 299             Double_t dPhi = pTrack->Phi();
 300             //set default phi weight to 1
 301             Double_t dW = 1.;
 302             //phi weight of pTrack
 303             if(fUsePhiWeights && fPhiWeights) {
 304               Int_t iNBinsPhi = fPhiWeights->GetNbinsX();
 305               dW = fPhiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*iNBinsPhi/TMath::TwoPi())));
 306               //bin = 1 + value*nbins/range
 307               //TMath::Floor rounds to the lower integer
 308             }
 309
 310             //calculate vU
 311             TVector2 vU;
 312             Double_t dUX = TMath::Cos(2*dPhi);
 313             Double_t dUY = TMath::Sin(2*dPhi);
 314             vU.Set(dUX,dUY);
 315             Double_t dModulus = vU.Mod();
 316             if (dModulus!=0.) vU.Set(dUX/dModulus,dUY/dModulus);  // make length 1
 317             else cerr<<"dModulus is zero!"<<endl;
 318
 319             //redefine the Q vector and devide by its multiplicity
 320             TVector2 vQm;
 321             Double_t dQmX = 0.;
 322             Double_t dQmY = 0.;
 323             //subtract particle from the flowvector if used to define it
 324             if (pTrack->InSubevent(0) || pTrack->InSubevent(1)) {
 325               dQmX = (vQ.X() - dW*dUX)/(dMq-1);
 326               dQmY = (vQ.Y() - dW*dUY)/(dMq-1);
 327
 328               vQm.Set(dQmX,dQmY);
 329
 330               //dUQ = scalar product of vU and vQm
 331               Double_t dUQ = (vU * vQm);
 332               Double_t dPt = pTrack->Pt();
 333               Double_t dEta = pTrack->Eta();
 334               //fill the profile histograms
 335               if (pTrack->InRPSelection()) {
 336                 fHistProUQetaRP -> Fill(dEta,dUQ,(dMq-1)); //Fill (Qu/(Mq-1)) with weight (Mq-1)
 337                 //needed for the error calculation:
 338                 fHistProUQQaQbEtaRP -> Fill(dEta,dUQ*dQaQb,(dMq-1)*dMa*dMb); //Fill [Qu/(Mq-1)]*[QaQb/MaMb] with weight (Mq-1)MaMb
 339                 fHistProUQPtRP -> Fill(dPt,dUQ,(dMq-1));                     //Fill (Qu/(Mq-1)) with weight (Mq-1)
 340                 fHistProUQQaQbPtRP -> Fill(dPt,dUQ*dQaQb,(dMq-1)*dMa*dMb);   //Fill [Qu/(Mq-1)]*[QaQb/MaMb] with weight (Mq-1)MaMb
 341
 342                 fHistSumOfWeightsEtaRP[0]->Fill(dEta,(dMq-1));        // sum of Mq-1
 343                 fHistSumOfWeightsEtaRP[1]->Fill(dEta,pow((dMq-1),2.));// sum of (Mq-1)^2
 344                 fHistSumOfWeightsEtaRP[2]->Fill(dEta,(dMq-1)*dMa*dMb);// sum of (Mq-1)*MaMb
 345                 fHistSumOfWeightsPtRP[0]->Fill(dPt,(dMq-1));          // sum of Mq-1
 346                 fHistSumOfWeightsPtRP[1]->Fill(dPt,pow((dMq-1),2.));  // sum of (Mq-1)^2
 347                 fHistSumOfWeightsPtRP[2]->Fill(dPt,(dMq-1)*dMa*dMb);  // sum of (Mq-1)*MaMb
 348               }
 349               if (pTrack->InPOISelection()) {
 350                 fHistProUQetaPOI -> Fill(dEta,dUQ,(dMq-1));//Fill (Qu/(Mq-1)) with weight (Mq-1)
 351                 //needed for the error calculation:
 352                 fHistProUQQaQbEtaPOI -> Fill(dEta,dUQ*dQaQb,(dMq-1)*dMa*dMb); //Fill [Qu/(Mq-1)]*[QaQb/MaMb] with weight (Mq-1)MaMb
 353                 fHistProUQPtPOI -> Fill(dPt,dUQ,(dMq-1));                     //Fill (Qu/(Mq-1)) with weight (Mq-1)
 354                 fHistProUQQaQbPtPOI -> Fill(dPt,dUQ*dQaQb,(dMq-1)*dMa*dMb);   //Fill [Qu/(Mq-1)]*[QaQb/MaMb] with weight (Mq-1)MaMb
 355
 356                 fHistSumOfWeightsEtaPOI[0]->Fill(dEta,(dMq-1));        // sum of Mq-1
 357                 fHistSumOfWeightsEtaPOI[1]->Fill(dEta,pow((dMq-1),2.));// sum of (Mq-1)^2
 358                 fHistSumOfWeightsEtaPOI[2]->Fill(dEta,(dMq-1)*dMa*dMb);// sum of (Mq-1)*MaMb
 359                 fHistSumOfWeightsPtPOI[0]->Fill(dPt,(dMq-1));          // sum of Mq-1
 360                 fHistSumOfWeightsPtPOI[1]->Fill(dPt,pow((dMq-1),2.)); // sum of (Mq-1)^2
 361                 fHistSumOfWeightsPtPOI[2]->Fill(dPt,(dMq-1)*dMa*dMb); // sum of (Mq-1)*MaMb
 362               }
 363
 364             } else { //do not subtract the particle from the flowvector
 365               dQmX = vQ.X()/dMq;
 366               dQmY = vQ.Y()/dMq;
 367
 368               vQm.Set(dQmX,dQmY);
 369
 370               //dUQ = scalar product of vU and vQm
 371               Double_t dUQ = (vU * vQm);
 372               Double_t dPt = pTrack->Pt();
 373               Double_t dEta = pTrack->Eta();
 374               //fill the profile histograms
 375               if (pTrack->InRPSelection()) {
 376                 fHistProUQetaRP -> Fill(dEta,dUQ,dMq);                   //Fill (Qu/Mq) with weight Mq
 377                 //needed for the error calculation:
 378                 fHistProUQQaQbEtaRP -> Fill(dEta,dUQ*dQaQb,dMq*dMa*dMb); //Fill [Qu/Mq]*[QaQb/MaMb] with weight Mq*MaMb
 379                 fHistProUQPtRP -> Fill(dPt,dUQ,dMq);                     //Fill (Qu/Mq) with weight Mq
 380                 fHistProUQQaQbPtRP -> Fill(dPt,dUQ*dQaQb,dMq*dMa*dMb);   //Fill [Qu/Mq]*[QaQb/MaMb] with weight Mq*MaMb
 381
 382                 fHistSumOfWeightsEtaRP[0]->Fill(dEta,dMq);        // sum of Mq
 383                 fHistSumOfWeightsEtaRP[1]->Fill(dEta,pow(dMq,2.));// sum of Mq^2
 384                 fHistSumOfWeightsEtaRP[2]->Fill(dEta,dMq*dMa*dMb);// sum of Mq*MaMb
 385                 fHistSumOfWeightsPtRP[0]->Fill(dPt,dMq);          // sum of Mq
 386                 fHistSumOfWeightsPtRP[1]->Fill(dPt,pow(dMq,2.));  // sum of Mq^2
 387                 fHistSumOfWeightsPtRP[2]->Fill(dPt,dMq*dMa*dMb);  // sum of Mq*MaMb
 388               }
 389               if (pTrack->InPOISelection()) {
 390                 fHistProUQetaPOI -> Fill(dEta,dUQ,dMq); //Fill (Qu/Mq) with weight Mq
 391                 //needed for the error calculation:
 392                 fHistProUQQaQbEtaPOI -> Fill(dEta,dUQ*dQaQb,dMq*dMa*dMb); //Fill [Qu/Mq]*[QaQb/MaMb] with weight Mq*MaMb
 393                 fHistProUQPtPOI -> Fill(dPt,dUQ,dMq);                     //Fill (Qu/Mq) with weight Mq
 394                 fHistProUQQaQbPtPOI -> Fill(dPt,dUQ*dQaQb,dMq*dMa*dMb);   //Fill [Qu/Mq]*[QaQb/MaMb] with weight Mq*MaMb
 395
 396                 fHistSumOfWeightsEtaPOI[0]->Fill(dEta,dMq);        // sum of Mq
 397                 fHistSumOfWeightsEtaPOI[1]->Fill(dEta,pow(dMq,2.));// sum of Mq^2
 398                 fHistSumOfWeightsEtaPOI[2]->Fill(dEta,dMq*dMa*dMb);// sum of Mq*MaMb
 399                 fHistSumOfWeightsPtPOI[0]->Fill(dPt,dMq);          // sum of Mq
 400                 fHistSumOfWeightsPtPOI[1]->Fill(dPt,pow(dMq,2.));  // sum of Mq^2
 401                 fHistSumOfWeightsPtPOI[2]->Fill(dPt,dMq*dMa*dMb);  // sum of Mq*MaMb
 402               }
 403             }//track not in subevents
 404
 405           }//track
 406
 407         }//loop over tracks
 408
 409       fEventNumber++;
 410       //    cout<<"@@@@@ "<<fEventNumber<<" events processed"<<endl;
 411
 412     }// subevents not empty
 413     delete [] vQarray;
 414   } //event
 415 }//end of Make()
 416
 417 //--------------------------------------------------------------------
 418 void AliFlowAnalysisWithScalarProduct::GetOutputHistograms(TList *outputListHistos){
 419
 420   //get pointers to all output histograms (called before Finish())
 421
 422   if (outputListHistos) {
 423   //Get the common histograms from the output list
 424     AliFlowCommonHist *pCommonHist = dynamic_cast<AliFlowCommonHist*>
 425       (outputListHistos->FindObject("AliFlowCommonHistSP"));
 426     AliFlowCommonHistResults *pCommonHistResults = dynamic_cast<AliFlowCommonHistResults*>
 427       (outputListHistos->FindObject("AliFlowCommonHistResultsSP"));
 428     TProfile* pHistProQaQb     = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_QaQb_SP"));
 429     TProfile* pHistProUQetaRP  = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQetaRP_SP"));
 430     TProfile* pHistProUQetaPOI = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQetaPOI_SP"));
 431     TProfile* pHistProUQPtRP   = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQPtRP_SP"));
 432     TProfile* pHistProUQPtPOI  = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQPtPOI_SP"));
 433     TH1D*     pHistSumOfLinearWeights    = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_SumOfLinearWeights_SP"));
 434     TH1D*     pHistSumOfQuadraticWeights = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_SumOfQuadraticWeights_SP"));
 435     TProfile* pHistProUQQaQbPtRP    = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQQaQbPtRP_SP"));
 436     TProfile* pHistProUQQaQbEtaRP   = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQQaQbEtaRP_SP"));
 437     TProfile* pHistProUQQaQbPtPOI   = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQQaQbPtPOI_SP"));
 438     TProfile* pHistProUQQaQbEtaPOI  = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQQaQbEtaPOI_SP"));
 439
 440     TString weightFlag[3] = {"w_Qu_","w_Qu^2_","w_QuQaQb_"};
 441     TH1D* pHistSumOfWeightsPtRP[3] = {NULL};
 442     TH1D* pHistSumOfWeightsEtaRP[3] = {NULL};
 443     TH1D* pHistSumOfWeightsPtPOI[3] = {NULL};
 444     TH1D* pHistSumOfWeightsEtaPOI[3] = {NULL};
 445     for(Int_t i=0;i<3;i++)
 446     {
 447      pHistSumOfWeightsPtRP[i]   = dynamic_cast<TH1D*>(outputListHistos->FindObject(Form("Flow_SumOfWeights%sPtRP_SP",weightFlag[i].Data())));
 448      pHistSumOfWeightsEtaRP[i]  = dynamic_cast<TH1D*>(outputListHistos->FindObject(Form("Flow_SumOfWeights%sEtaRP_SP",weightFlag[i].Data())));
 449      pHistSumOfWeightsPtPOI[i]  = dynamic_cast<TH1D*>(outputListHistos->FindObject(Form("Flow_SumOfWeights%sPtPOI_SP",weightFlag[i].Data())));
 450      pHistSumOfWeightsEtaPOI[i] = dynamic_cast<TH1D*>(outputListHistos->FindObject(Form("Flow_SumOfWeights%sEtaPOI_SP",weightFlag[i].Data())));
 451     }
 452
 453     //pass the pointers to the task
 454     if (pCommonHist && pCommonHistResults && pHistProQaQb &&
 455         pHistProUQetaRP && pHistProUQetaPOI && pHistProUQPtRP && pHistProUQPtPOI
 456         && pHistSumOfLinearWeights && pHistSumOfQuadraticWeights && pHistProUQQaQbPtRP
 457         && pHistProUQQaQbEtaRP && pHistProUQQaQbPtPOI && pHistProUQQaQbEtaPOI) {
 458       this -> SetCommonHists(pCommonHist);
 459       this -> SetCommonHistsRes(pCommonHistResults);
 460       this -> SetHistProQaQb(pHistProQaQb);
 461       this -> SetHistProUQetaRP(pHistProUQetaRP);
 462       this -> SetHistProUQetaPOI(pHistProUQetaPOI);
 463       this -> SetHistProUQPtRP(pHistProUQPtRP);
 464       this -> SetHistProUQPtPOI(pHistProUQPtPOI);
 465       this -> SetHistProUQQaQbPtRP(pHistProUQQaQbPtRP);
 466       this -> SetHistProUQQaQbEtaRP(pHistProUQQaQbEtaRP);
 467       this -> SetHistProUQQaQbPtPOI(pHistProUQQaQbPtPOI);
 468       this -> SetHistProUQQaQbEtaPOI(pHistProUQQaQbEtaPOI);
 469       this -> SetHistSumOfLinearWeights(pHistSumOfLinearWeights);
 470       this -> SetHistSumOfQuadraticWeights(pHistSumOfQuadraticWeights);
 471      }
 472
 473    for(Int_t i=0;i<3;i++)
 474    {
 475     if(pHistSumOfWeightsPtRP[i]) this -> SetHistSumOfWeightsPtRP(pHistSumOfWeightsPtRP[i],i);
 476     if(pHistSumOfWeightsEtaRP[i]) this -> SetHistSumOfWeightsEtaRP(pHistSumOfWeightsEtaRP[i],i);
 477     if(pHistSumOfWeightsPtPOI[i]) this -> SetHistSumOfWeightsPtPOI(pHistSumOfWeightsPtPOI[i],i);
 478     if(pHistSumOfWeightsEtaPOI[i]) this -> SetHistSumOfWeightsEtaPOI(pHistSumOfWeightsEtaPOI[i],i);
 479    }
 480   } // end of if(outputListHistos)
 481 }
 482
 483 //--------------------------------------------------------------------
 484 void AliFlowAnalysisWithScalarProduct::Finish() {
 485
 486   //calculate flow and fill the AliFlowCommonHistResults\r
 487   if (fDebug) cout<<"AliFlowAnalysisWithScalarProduct::Finish()"<<endl;
 488
 489   cout<<"*************************************"<<endl;
 490   cout<<"*************************************"<<endl;
 491   cout<<"      Integrated flow from           "<<endl;
 492   cout<<"         Scalar product              "<<endl;
 493   cout<<endl;
 494
 495   Int_t iNbinsPt  = AliFlowCommonConstants::GetMaster()->GetNbinsPt();
 496   Int_t iNbinsEta = AliFlowCommonConstants::GetMaster()->GetNbinsEta();
 497
 498
 499   //Calculate reference flow (noname)
 500   //----------------------------------
 501   //weighted average over (QaQb/MaMb) with weight (MaMb)
 502   Double_t dQaQb  = fHistProQaQb->GetBinContent(1);
 503   Double_t dV = TMath::Sqrt(dQaQb);
 504   //statistical error of dQaQb:
 505   //  statistical error = term1 * spread * term2:
 506   //  term1 = sqrt{sum_{i=1}^{N} w^2}/(sum_{i=1}^{N} w)
 507   //  term2 = 1/sqrt(1-term1^2)
 508   Double_t dSpreadQaQb = fHistProQaQb->GetBinError(1);
 509   Double_t dSumOfLinearWeights = fHistSumOfLinearWeights->GetBinContent(1);
 510   Double_t dSumOfQuadraticWeights = fHistSumOfQuadraticWeights->GetBinContent(1);
 511   Double_t dTerm1 = 0.;
 512   Double_t dTerm2 = 0.;
 513   if(dSumOfLinearWeights) {
 514     dTerm1 = pow(dSumOfQuadraticWeights,0.5)/dSumOfLinearWeights;
 515   }
 516   if(1.-pow(dTerm1,2.) > 0.) {
 517     dTerm2 = 1./pow(1-pow(dTerm1,2.),0.5);
 518   }
 519   Double_t dStatErrorQaQb = dTerm1 * dSpreadQaQb * dTerm2;
 520   //calculate the statistical error
 521   Double_t dVerr = 0.;
 522   if(dQaQb > 0.) {
 523     dVerr = (1./(2.*pow(dQaQb,0.5)))*dStatErrorQaQb;
 524   }
 525   fCommonHistsRes->FillIntegratedFlow(dV,dVerr);
 526   cout<<"dV = "<<dV<<" +- "<<dVerr<<endl;
 527
 528   //Calculate differential flow and integrated flow (RP, POI)
 529   //---------------------------------------------------------
 530   //v as a function of eta for RP selection
 531   for(Int_t b=0;b<iNbinsEta;b++) {
 532     Double_t duQpro = fHistProUQetaRP->GetBinContent(b);
 533     Double_t dv2pro = 0.;
 534     if (dV!=0.) { dv2pro = duQpro/dV; }
 535     //calculate the statistical error
 536     Double_t dv2ProErr = CalculateStatisticalError(b, dStatErrorQaQb, fHistProUQetaRP, fHistProUQQaQbEtaRP, fHistSumOfWeightsEtaRP);
 537     //fill TH1D
 538     fCommonHistsRes->FillDifferentialFlowEtaRP(b, dv2pro, dv2ProErr);
 539   } //loop over bins b
 540
 541
 542   //v as a function of eta for POI selection
 543   for(Int_t b=0;b<iNbinsEta;b++) {
 544     Double_t duQpro = fHistProUQetaPOI->GetBinContent(b);
 545     Double_t dv2pro = 0.;
 546     if (dV!=0.) { dv2pro = duQpro/dV; }
 547     //calculate the statistical error
 548     Double_t dv2ProErr = CalculateStatisticalError(b, dStatErrorQaQb, fHistProUQetaPOI, fHistProUQQaQbEtaPOI, fHistSumOfWeightsEtaPOI);
 549
 550     //fill TH1D
 551     fCommonHistsRes->FillDifferentialFlowEtaPOI(b, dv2pro, dv2ProErr);
 552   } //loop over bins b
 553
 554
 555   //v as a function of Pt for RP selection
 556   TH1F* fHistPtRP = fCommonHists->GetHistPtRP(); //for calculating integrated flow
 557   Double_t dVRP = 0.;
 558   Double_t dSumRP = 0.;
 559   Double_t dErrVRP =0.;
 560
 561   for(Int_t b=0;b<iNbinsPt;b++) {
 562     Double_t duQpro = fHistProUQPtRP->GetBinContent(b);
 563     Double_t dv2pro = 0.;
 564     if (dV!=0.) { dv2pro = duQpro/dV; }
 565     //calculate the statistical error
 566     Double_t dv2ProErr = CalculateStatisticalError(b, dStatErrorQaQb, fHistProUQPtRP, fHistProUQQaQbPtRP, fHistSumOfWeightsPtRP);
 567
 568     //fill TH1D
 569     fCommonHistsRes->FillDifferentialFlowPtRP(b, dv2pro, dv2ProErr);
 570
 571     //calculate integrated flow for RP selection
 572     if (fHistPtRP){
 573       Double_t dYieldPt = fHistPtRP->GetBinContent(b);
 574       dVRP += dv2pro*dYieldPt;
 575       dSumRP +=dYieldPt;
 576       dErrVRP += dYieldPt*dYieldPt*dv2ProErr*dv2ProErr;
 577     } else { cout<<"fHistPtRP is NULL"<<endl; }
 578   } //loop over bins b
 579
 580   if (dSumRP != 0.) {
 581     dVRP /= dSumRP; //the pt distribution should be normalised
 582     dErrVRP /= (dSumRP*dSumRP);
 583     dErrVRP = TMath::Sqrt(dErrVRP);
 584   }
 585   fCommonHistsRes->FillIntegratedFlowRP(dVRP,dErrVRP);
 586   cout<<"dV(RP) = "<<dVRP<<" +- "<<dErrVRP<<endl;
 587
 588
 589   //v as a function of Pt for POI selection
 590   TH1F* fHistPtPOI = fCommonHists->GetHistPtPOI(); //for calculating integrated flow
 591   Double_t dVPOI = 0.;
 592   Double_t dSumPOI = 0.;
 593   Double_t dErrVPOI =0.;
 594
 595   for(Int_t b=0;b<iNbinsPt;b++) {
 596     Double_t duQpro = fHistProUQPtPOI->GetBinContent(b);
 597     Double_t dv2pro = 0.;
 598     if (dV!=0.) { dv2pro = duQpro/dV; }
 599     //calculate the statistical error
 600     Double_t dv2ProErr = CalculateStatisticalError(b, dStatErrorQaQb, fHistProUQPtPOI, fHistProUQQaQbPtPOI, fHistSumOfWeightsPtPOI);
 601
 602     //fill TH1D
 603     fCommonHistsRes->FillDifferentialFlowPtPOI(b, dv2pro, dv2ProErr);
 604
 605     //calculate integrated flow for POI selection
 606     if (fHistPtPOI){
 607       Double_t dYieldPt = fHistPtPOI->GetBinContent(b);
 608       dVPOI += dv2pro*dYieldPt;
 609       dSumPOI +=dYieldPt;
 610       dErrVPOI += dYieldPt*dYieldPt*dv2ProErr*dv2ProErr;
 611     } else { cout<<"fHistPtPOI is NULL"<<endl; }
 612   } //loop over bins b
 613
 614   if (dSumPOI != 0.) {
 615     dVPOI /= dSumPOI; //the pt distribution should be normalised\r
 616     dErrVPOI /= (dSumPOI*dSumPOI);
 617     dErrVPOI = TMath::Sqrt(dErrVPOI);
 618   }
 619   fCommonHistsRes->FillIntegratedFlowPOI(dVPOI,dErrVPOI);
 620   cout<<"dV(POI) = "<<dVPOI<<" +- "<<dErrVPOI<<endl;
 621
 622   cout<<endl;
 623   cout<<"*************************************"<<endl;
 624   cout<<"*************************************"<<endl;
 625
 626   //cout<<".....finished"<<endl;
 627 }
 628
 629
 630 //--------------------------------------------------------------------
 631 Double_t AliFlowAnalysisWithScalarProduct::CalculateStatisticalError(Int_t b, Double_t aStatErrorQaQb, TProfile* pHistProUQ, TProfile* pHistProUQQaQb, TH1D** pHistSumOfWeights) {
 632   //calculate the statistical error for differential flow for bin b
 633   Double_t duQproSpread = pHistProUQ->GetBinError(b);
 634   Double_t sumOfMq = pHistSumOfWeights[0]->GetBinContent(b);
 635   Double_t sumOfMqSquared = pHistSumOfWeights[1]->GetBinContent(b);
 636   Double_t dTerm1 = 0.;
 637   Double_t dTerm2 = 0.;
 638   if(sumOfMq) {
 639     dTerm1 = (pow(sumOfMqSquared,0.5)/sumOfMq);
 640   }
 641   if(1.-pow(dTerm1,2.)>0.) {
 642     dTerm2 = 1./pow(1.-pow(dTerm1,2.),0.5);
 643   }
 644   Double_t duQproErr = dTerm1*duQproSpread*dTerm2;
 645   // covariances:
 646   Double_t dTerm1Cov = pHistSumOfWeights[2]->GetBinContent(b);
 647   Double_t dTerm2Cov = fHistSumOfLinearWeights->GetBinContent(1);
 648   Double_t dTerm3Cov = sumOfMq;
 649   Double_t dWeightedCovariance = 0.;
 650   if(dTerm2Cov*dTerm3Cov>0.) {
 651     Double_t dDenominator = 1.-dTerm1Cov/(dTerm2Cov*dTerm3Cov);
 652     Double_t dPrefactor = dTerm1Cov/(dTerm2Cov*dTerm3Cov);
 653     if(dDenominator!=0) {
 654       Double_t dCovariance = (pHistProUQQaQb->GetBinContent(b)-fHistProQaQb->GetBinContent(1)*pHistProUQ->GetBinContent(b))/dDenominator;
 655       dWeightedCovariance = dCovariance*dPrefactor;
 656     }
 657   }
 658
 659   Double_t dv2ProErr = 0.; // final statitical error:
 660   if(fHistProQaQb->GetBinContent(1)>0.) {
 661     Double_t dv2ProErrorSquared = (1./4.)*pow(fHistProQaQb->GetBinContent(1),-3.)*
 662       (pow(pHistProUQ->GetBinContent(b),2.)*pow(aStatErrorQaQb,2.)
 663        + 4.*pow(fHistProQaQb->GetBinContent(1),2.)*pow(duQproErr,2.)
 664        - 4.*fHistProQaQb->GetBinContent(1)*pHistProUQ->GetBinContent(b)*dWeightedCovariance);
 665     if(dv2ProErrorSquared>0.) dv2ProErr = pow(dv2ProErrorSquared,0.5);
 666   }
 667
 668   return dv2ProErr;
 669 }
 670 //--------------------------------------------------------------------