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   //save old value and prevent histograms from being added to directory
 132   //to avoid name clashes in case multiple analaysis objects are used
 133   //in an analysis
 134   Bool_t oldHistAddStatus = TH1::AddDirectoryStatus();
 135   TH1::AddDirectory(kFALSE);
 136
 137   Int_t iNbinsPt   = AliFlowCommonConstants::GetMaster()->GetNbinsPt();
 138   Double_t dPtMin  = AliFlowCommonConstants::GetMaster()->GetPtMin();
 139   Double_t dPtMax  = AliFlowCommonConstants::GetMaster()->GetPtMax();
 140   Int_t iNbinsEta  = AliFlowCommonConstants::GetMaster()->GetNbinsEta();
 141   Double_t dEtaMin = AliFlowCommonConstants::GetMaster()->GetEtaMin();
 142   Double_t dEtaMax = AliFlowCommonConstants::GetMaster()->GetEtaMax();
 143
 144   fHistProUQetaRP = new TProfile("Flow_UQetaRP_SP","Flow_UQetaRP_SP",iNbinsEta,dEtaMin,dEtaMax,"s");
 145   fHistProUQetaRP->SetXTitle("{eta}");
 146   fHistProUQetaRP->SetYTitle("<uQ>");
 147   fHistList->Add(fHistProUQetaRP);
 148
 149   fHistProUQetaPOI = new TProfile("Flow_UQetaPOI_SP","Flow_UQetaPOI_SP",iNbinsEta,dEtaMin,dEtaMax,"s");
 150   fHistProUQetaPOI->SetXTitle("{eta}");
 151   fHistProUQetaPOI->SetYTitle("<uQ>");
 152   fHistList->Add(fHistProUQetaPOI);
 153
 154   fHistProUQPtRP = new TProfile("Flow_UQPtRP_SP","Flow_UQPtRP_SP",iNbinsPt,dPtMin,dPtMax,"s");
 155   fHistProUQPtRP->SetXTitle("p_t (GeV)");
 156   fHistProUQPtRP->SetYTitle("<uQ>");
 157   fHistList->Add(fHistProUQPtRP);
 158
 159   fHistProUQPtPOI = new TProfile("Flow_UQPtPOI_SP","Flow_UQPtPOI_SP",iNbinsPt,dPtMin,dPtMax,"s");
 160   fHistProUQPtPOI->SetXTitle("p_t (GeV)");
 161   fHistProUQPtPOI->SetYTitle("<uQ>");
 162   fHistList->Add(fHistProUQPtPOI);
 163
 164   fHistProQaQb = new TProfile("Flow_QaQb_SP","Flow_QaQb_SP", 1, 0.5, 1.5,"s");
 165   fHistProQaQb->SetYTitle("<QaQb>");
 166   fHistList->Add(fHistProQaQb);
 167
 168   fHistSumOfLinearWeights = new TH1D("Flow_SumOfLinearWeights_SP","Flow_SumOfLinearWeights_SP",1,-0.5, 0.5);
 169   fHistSumOfLinearWeights -> SetYTitle("sum (*)");
 170   fHistSumOfLinearWeights -> SetXTitle("sum (Ma*Mb)");
 171   fHistList->Add(fHistSumOfLinearWeights);
 172
 173   fHistSumOfQuadraticWeights = new TH1D("Flow_SumOfQuadraticWeights_SP","Flow_SumOfQuadraticWeights_SP",1,-0.5, 0.5);
 174   fHistSumOfQuadraticWeights -> SetYTitle("sum (*)");
 175   fHistSumOfQuadraticWeights -> SetXTitle("sum (Ma*Mb)^2");
 176   fHistList->Add(fHistSumOfQuadraticWeights);
 177
 178   fHistProUQQaQbPtRP = new TProfile("Flow_UQQaQbPtRP_SP","Flow_UQQaQbPtRP_SP",iNbinsPt,dPtMin,dPtMax);
 179   fHistProUQQaQbPtRP -> SetYTitle("<*>");
 180   fHistProUQQaQbPtRP -> SetXTitle("<Qu QaQb>");
 181   fHistList->Add(fHistProUQQaQbPtRP);
 182
 183   fHistProUQQaQbEtaRP = new TProfile("Flow_UQQaQbEtaRP_SP","Flow_UQQaQbEtaRP_SP",iNbinsEta,dEtaMin,dEtaMax);
 184   fHistProUQQaQbEtaRP -> SetYTitle("<*>");
 185   fHistProUQQaQbEtaRP -> SetXTitle("<Qu QaQb>");
 186   fHistList->Add(fHistProUQQaQbEtaRP);
 187
 188   fHistProUQQaQbPtPOI = new TProfile("Flow_UQQaQbPtPOI_SP","Flow_UQQaQbPtPOI_SP",iNbinsPt,dPtMin,dPtMax);
 189   fHistProUQQaQbPtPOI -> SetYTitle("<*>");
 190   fHistProUQQaQbPtPOI -> SetXTitle("<Qu QaQb>");
 191   fHistList->Add(fHistProUQQaQbPtPOI);
 192
 193   fHistProUQQaQbEtaPOI = new TProfile("Flow_UQQaQbEtaPOI_SP","Flow_UQQaQbEtaPOI_SP",iNbinsEta,dEtaMin,dEtaMax);
 194   fHistProUQQaQbEtaPOI -> SetYTitle("<*>");
 195   fHistProUQQaQbEtaPOI -> SetXTitle("<Qu QaQb>");
 196   fHistList->Add(fHistProUQQaQbEtaPOI);
 197
 198   TString weightFlag[3] = {"w_Qu_","w_Qu^2_","w_QuQaQb_"};
 199   for(Int_t i=0;i<3;i++)
 200   {
 201    fHistSumOfWeightsPtRP[i] = new TH1D(Form("Flow_SumOfWeights%sPtRP_SP",weightFlag[i].Data()),
 202                               Form("Flow_SumOfWeights%sPtRP_SP",weightFlag[i].Data()),iNbinsPt,dPtMin,dPtMax);
 203    fHistSumOfWeightsPtRP[i] -> SetYTitle("sum (*)");
 204    fHistSumOfWeightsPtRP[i] -> SetXTitle("p_{T}");
 205    fHistList->Add(fHistSumOfWeightsPtRP[i]);
 206
 207    fHistSumOfWeightsEtaRP[i] = new TH1D(Form("Flow_SumOfWeights%sEtaRP_SP",weightFlag[i].Data()),
 208                                Form("Flow_SumOfWeights%sEtaRP_SP",weightFlag[i].Data()),iNbinsEta,dEtaMin,dEtaMax);
 209    fHistSumOfWeightsEtaRP[i] -> SetYTitle("sum (*)");
 210    fHistSumOfWeightsEtaRP[i] -> SetXTitle("#eta");
 211    fHistList->Add(fHistSumOfWeightsEtaRP[i]);
 212
 213    fHistSumOfWeightsPtPOI[i] = new TH1D(Form("Flow_SumOfWeights%sPtPOI_SP",weightFlag[i].Data()),
 214                                Form("Flow_SumOfWeights%sPtPOI_SP",weightFlag[i].Data()),iNbinsPt,dPtMin,dPtMax);
 215    fHistSumOfWeightsPtPOI[i] -> SetYTitle("sum (*)");
 216    fHistSumOfWeightsPtPOI[i] -> SetXTitle("p_{T}");
 217    fHistList->Add(fHistSumOfWeightsPtPOI[i]);
 218
 219    fHistSumOfWeightsEtaPOI[i] = new TH1D(Form("Flow_SumOfWeights%sEtaPOI_SP",weightFlag[i].Data()),
 220                                 Form("Flow_SumOfWeights%sEtaPOI_SP",weightFlag[i].Data()),iNbinsEta,dEtaMin,dEtaMax);
 221    fHistSumOfWeightsEtaPOI[i] -> SetYTitle("sum (*)");
 222    fHistSumOfWeightsEtaPOI[i] -> SetXTitle("#eta");
 223    fHistList->Add(fHistSumOfWeightsEtaPOI[i]);
 224   }
 225
 226   fCommonHists = new AliFlowCommonHist("AliFlowCommonHistSP");
 227   fHistList->Add(fCommonHists);
 228   fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsSP");
 229   fHistList->Add(fCommonHistsRes);
 230
 231   //weights
 232   if(fUsePhiWeights) {
 233     if(!fWeightsList) {
 234       cout<<"WARNING: fWeightsList is NULL in the Scalar Product method."<<endl;
 235       exit(0);
 236     }
 237     if(fWeightsList->FindObject("phi_weights"))  {
 238       fPhiWeights = dynamic_cast<TH1F*>(fWeightsList->FindObject("phi_weights"));
 239       fHistList->Add(fPhiWeights);
 240     } else {
 241       cout<<"WARNING: histogram with phi weights is not accessible in Scalar Product"<<endl;
 242       exit(0);
 243     }
 244   } // end of if(fUsePhiWeights)
 245
 246   fEventNumber = 0;  //set number of events to zero
 247
 248   TH1::AddDirectory(oldHistAddStatus);
 249 }
 250
 251 //-----------------------------------------------------------------------
 252 void AliFlowAnalysisWithScalarProduct::Make(AliFlowEventSimple* anEvent) {
 253
 254   //Calculate flow based on  <QaQb/MaMb> = <v^2>
 255
 256   //Fill histograms
 257   if (anEvent) {
 258
 259     //get Q vectors for the eta-subevents
 260     AliFlowVector* vQarray = new AliFlowVector[2];
 261     if (fUsePhiWeights) {
 262       anEvent->Get2Qsub(vQarray,2,fWeightsList,kTRUE);
 263     } else {
 264       anEvent->Get2Qsub(vQarray);
 265     }
 266     //subevent a
 267     AliFlowVector vQa = vQarray[0];
 268     //subevent b
 269     AliFlowVector vQb = vQarray[1];
 270
 271     //check that the subevents are not empty:
 272     Double_t dMa = vQa.GetMult();
 273     Double_t dMb = vQb.GetMult();
 274     if (dMa != 0 && dMb != 0) {
 275
 276       //fill control histograms
 277       fCommonHists->FillControlHistograms(anEvent);
 278
 279       //get total Q vector = the sum of subevent a and subevent b
 280       AliFlowVector vQ = vQa + vQb;
 281       //weight the Q vectors for the subevents by the multiplicity
 282       //Note: Weight Q only in the particle loop when it is clear if it should be (m-1) or M
 283       Double_t dQXa = vQa.X()/dMa;
 284       Double_t dQYa = vQa.Y()/dMa;
 285       vQa.Set(dQXa,dQYa);
 286
 287       Double_t dQXb = vQb.X()/dMb;
 288       Double_t dQYb = vQb.Y()/dMb;
 289       vQb.Set(dQXb,dQYb);
 290
 291       //scalar product of the two subevents
 292       Double_t dQaQb = (vQa*vQb);
 293       fHistProQaQb -> Fill(1.,dQaQb,dMa*dMb);  //Fill (QaQb/MaMb) with weight (MaMb).
 294       //needed for the error calculation:
 295       fHistSumOfLinearWeights -> Fill(0.,dMa*dMb);
 296       fHistSumOfQuadraticWeights -> Fill(0.,pow(dMa*dMb,2.));
 297
 298       //loop over the tracks of the event
 299       AliFlowTrackSimple*   pTrack = NULL;
 300       Int_t iNumberOfTracks = anEvent->NumberOfTracks();
 301       Double_t dMq =  vQ.GetMult();
 302
 303       for (Int_t i=0;i<iNumberOfTracks;i++)
 304         {
 305           pTrack = anEvent->GetTrack(i) ;
 306           if (pTrack){
 307             Double_t dPhi = pTrack->Phi();
 308             //set default phi weight to 1
 309             Double_t dW = 1.;
 310             //phi weight of pTrack
 311             if(fUsePhiWeights && fPhiWeights) {
 312               Int_t iNBinsPhi = fPhiWeights->GetNbinsX();
 313               dW = fPhiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*iNBinsPhi/TMath::TwoPi())));
 314               //bin = 1 + value*nbins/range
 315               //TMath::Floor rounds to the lower integer
 316             }
 317
 318             //calculate vU
 319             TVector2 vU;
 320             Double_t dUX = TMath::Cos(2*dPhi);
 321             Double_t dUY = TMath::Sin(2*dPhi);
 322             vU.Set(dUX,dUY);
 323             Double_t dModulus = vU.Mod();
 324             if (dModulus!=0.) vU.Set(dUX/dModulus,dUY/dModulus);  // make length 1
 325             else cerr<<"dModulus is zero!"<<endl;
 326
 327             //redefine the Q vector and devide by its multiplicity
 328             TVector2 vQm;
 329             Double_t dQmX = 0.;
 330             Double_t dQmY = 0.;
 331             //subtract particle from the flowvector if used to define it
 332             if (pTrack->InSubevent(0) || pTrack->InSubevent(1)) {
 333               dQmX = (vQ.X() - dW*dUX)/(dMq-1);
 334               dQmY = (vQ.Y() - dW*dUY)/(dMq-1);
 335
 336               vQm.Set(dQmX,dQmY);
 337
 338               //dUQ = scalar product of vU and vQm
 339               Double_t dUQ = (vU * vQm);
 340               Double_t dPt = pTrack->Pt();
 341               Double_t dEta = pTrack->Eta();
 342               //fill the profile histograms
 343               if (pTrack->InRPSelection()) {
 344                 fHistProUQetaRP -> Fill(dEta,dUQ,(dMq-1)); //Fill (Qu/(Mq-1)) with weight (Mq-1)
 345                 //needed for the error calculation:
 346                 fHistProUQQaQbEtaRP -> Fill(dEta,dUQ*dQaQb,(dMq-1)*dMa*dMb); //Fill [Qu/(Mq-1)]*[QaQb/MaMb] with weight (Mq-1)MaMb
 347                 fHistProUQPtRP -> Fill(dPt,dUQ,(dMq-1));                     //Fill (Qu/(Mq-1)) with weight (Mq-1)
 348                 fHistProUQQaQbPtRP -> Fill(dPt,dUQ*dQaQb,(dMq-1)*dMa*dMb);   //Fill [Qu/(Mq-1)]*[QaQb/MaMb] with weight (Mq-1)MaMb
 349
 350                 fHistSumOfWeightsEtaRP[0]->Fill(dEta,(dMq-1));        // sum of Mq-1
 351                 fHistSumOfWeightsEtaRP[1]->Fill(dEta,pow((dMq-1),2.));// sum of (Mq-1)^2
 352                 fHistSumOfWeightsEtaRP[2]->Fill(dEta,(dMq-1)*dMa*dMb);// sum of (Mq-1)*MaMb
 353                 fHistSumOfWeightsPtRP[0]->Fill(dPt,(dMq-1));          // sum of Mq-1
 354                 fHistSumOfWeightsPtRP[1]->Fill(dPt,pow((dMq-1),2.));  // sum of (Mq-1)^2
 355                 fHistSumOfWeightsPtRP[2]->Fill(dPt,(dMq-1)*dMa*dMb);  // sum of (Mq-1)*MaMb
 356               }
 357               if (pTrack->InPOISelection()) {
 358                 fHistProUQetaPOI -> Fill(dEta,dUQ,(dMq-1));//Fill (Qu/(Mq-1)) with weight (Mq-1)
 359                 //needed for the error calculation:
 360                 fHistProUQQaQbEtaPOI -> Fill(dEta,dUQ*dQaQb,(dMq-1)*dMa*dMb); //Fill [Qu/(Mq-1)]*[QaQb/MaMb] with weight (Mq-1)MaMb
 361                 fHistProUQPtPOI -> Fill(dPt,dUQ,(dMq-1));                     //Fill (Qu/(Mq-1)) with weight (Mq-1)
 362                 fHistProUQQaQbPtPOI -> Fill(dPt,dUQ*dQaQb,(dMq-1)*dMa*dMb);   //Fill [Qu/(Mq-1)]*[QaQb/MaMb] with weight (Mq-1)MaMb
 363
 364                 fHistSumOfWeightsEtaPOI[0]->Fill(dEta,(dMq-1));        // sum of Mq-1
 365                 fHistSumOfWeightsEtaPOI[1]->Fill(dEta,pow((dMq-1),2.));// sum of (Mq-1)^2
 366                 fHistSumOfWeightsEtaPOI[2]->Fill(dEta,(dMq-1)*dMa*dMb);// sum of (Mq-1)*MaMb
 367                 fHistSumOfWeightsPtPOI[0]->Fill(dPt,(dMq-1));          // sum of Mq-1
 368                 fHistSumOfWeightsPtPOI[1]->Fill(dPt,pow((dMq-1),2.)); // sum of (Mq-1)^2
 369                 fHistSumOfWeightsPtPOI[2]->Fill(dPt,(dMq-1)*dMa*dMb); // sum of (Mq-1)*MaMb
 370               }
 371
 372             } else { //do not subtract the particle from the flowvector
 373               dQmX = vQ.X()/dMq;
 374               dQmY = vQ.Y()/dMq;
 375
 376               vQm.Set(dQmX,dQmY);
 377
 378               //dUQ = scalar product of vU and vQm
 379               Double_t dUQ = (vU * vQm);
 380               Double_t dPt = pTrack->Pt();
 381               Double_t dEta = pTrack->Eta();
 382               //fill the profile histograms
 383               if (pTrack->InRPSelection()) {
 384                 fHistProUQetaRP -> Fill(dEta,dUQ,dMq);                   //Fill (Qu/Mq) with weight Mq
 385                 //needed for the error calculation:
 386                 fHistProUQQaQbEtaRP -> Fill(dEta,dUQ*dQaQb,dMq*dMa*dMb); //Fill [Qu/Mq]*[QaQb/MaMb] with weight Mq*MaMb
 387                 fHistProUQPtRP -> Fill(dPt,dUQ,dMq);                     //Fill (Qu/Mq) with weight Mq
 388                 fHistProUQQaQbPtRP -> Fill(dPt,dUQ*dQaQb,dMq*dMa*dMb);   //Fill [Qu/Mq]*[QaQb/MaMb] with weight Mq*MaMb
 389
 390                 fHistSumOfWeightsEtaRP[0]->Fill(dEta,dMq);        // sum of Mq
 391                 fHistSumOfWeightsEtaRP[1]->Fill(dEta,pow(dMq,2.));// sum of Mq^2
 392                 fHistSumOfWeightsEtaRP[2]->Fill(dEta,dMq*dMa*dMb);// sum of Mq*MaMb
 393                 fHistSumOfWeightsPtRP[0]->Fill(dPt,dMq);          // sum of Mq
 394                 fHistSumOfWeightsPtRP[1]->Fill(dPt,pow(dMq,2.));  // sum of Mq^2
 395                 fHistSumOfWeightsPtRP[2]->Fill(dPt,dMq*dMa*dMb);  // sum of Mq*MaMb
 396               }
 397               if (pTrack->InPOISelection()) {
 398                 fHistProUQetaPOI -> Fill(dEta,dUQ,dMq); //Fill (Qu/Mq) with weight Mq
 399                 //needed for the error calculation:
 400                 fHistProUQQaQbEtaPOI -> Fill(dEta,dUQ*dQaQb,dMq*dMa*dMb); //Fill [Qu/Mq]*[QaQb/MaMb] with weight Mq*MaMb
 401                 fHistProUQPtPOI -> Fill(dPt,dUQ,dMq);                     //Fill (Qu/Mq) with weight Mq
 402                 fHistProUQQaQbPtPOI -> Fill(dPt,dUQ*dQaQb,dMq*dMa*dMb);   //Fill [Qu/Mq]*[QaQb/MaMb] with weight Mq*MaMb
 403
 404                 fHistSumOfWeightsEtaPOI[0]->Fill(dEta,dMq);        // sum of Mq
 405                 fHistSumOfWeightsEtaPOI[1]->Fill(dEta,pow(dMq,2.));// sum of Mq^2
 406                 fHistSumOfWeightsEtaPOI[2]->Fill(dEta,dMq*dMa*dMb);// sum of Mq*MaMb
 407                 fHistSumOfWeightsPtPOI[0]->Fill(dPt,dMq);          // sum of Mq
 408                 fHistSumOfWeightsPtPOI[1]->Fill(dPt,pow(dMq,2.));  // sum of Mq^2
 409                 fHistSumOfWeightsPtPOI[2]->Fill(dPt,dMq*dMa*dMb);  // sum of Mq*MaMb
 410               }
 411             }//track not in subevents
 412
 413           }//track
 414
 415         }//loop over tracks
 416
 417       fEventNumber++;
 418       //    cout<<"@@@@@ "<<fEventNumber<<" events processed"<<endl;
 419
 420     }// subevents not empty
 421     delete [] vQarray;
 422   } //event
 423 }//end of Make()
 424
 425 //--------------------------------------------------------------------
 426 void AliFlowAnalysisWithScalarProduct::GetOutputHistograms(TList *outputListHistos){
 427
 428   //get pointers to all output histograms (called before Finish())
 429
 430   if (outputListHistos) {
 431   //Get the common histograms from the output list
 432     AliFlowCommonHist *pCommonHist = dynamic_cast<AliFlowCommonHist*>
 433       (outputListHistos->FindObject("AliFlowCommonHistSP"));
 434     AliFlowCommonHistResults *pCommonHistResults = dynamic_cast<AliFlowCommonHistResults*>
 435       (outputListHistos->FindObject("AliFlowCommonHistResultsSP"));
 436     TProfile* pHistProQaQb     = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_QaQb_SP"));
 437     TProfile* pHistProUQetaRP  = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQetaRP_SP"));
 438     TProfile* pHistProUQetaPOI = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQetaPOI_SP"));
 439     TProfile* pHistProUQPtRP   = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQPtRP_SP"));
 440     TProfile* pHistProUQPtPOI  = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQPtPOI_SP"));
 441     TH1D*     pHistSumOfLinearWeights    = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_SumOfLinearWeights_SP"));
 442     TH1D*     pHistSumOfQuadraticWeights = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_SumOfQuadraticWeights_SP"));
 443     TProfile* pHistProUQQaQbPtRP    = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQQaQbPtRP_SP"));
 444     TProfile* pHistProUQQaQbEtaRP   = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQQaQbEtaRP_SP"));
 445     TProfile* pHistProUQQaQbPtPOI   = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQQaQbPtPOI_SP"));
 446     TProfile* pHistProUQQaQbEtaPOI  = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQQaQbEtaPOI_SP"));
 447
 448     TString weightFlag[3] = {"w_Qu_","w_Qu^2_","w_QuQaQb_"};
 449     TH1D* pHistSumOfWeightsPtRP[3] = {NULL};
 450     TH1D* pHistSumOfWeightsEtaRP[3] = {NULL};
 451     TH1D* pHistSumOfWeightsPtPOI[3] = {NULL};
 452     TH1D* pHistSumOfWeightsEtaPOI[3] = {NULL};
 453     for(Int_t i=0;i<3;i++)
 454     {
 455      pHistSumOfWeightsPtRP[i]   = dynamic_cast<TH1D*>(outputListHistos->FindObject(Form("Flow_SumOfWeights%sPtRP_SP",weightFlag[i].Data())));
 456      pHistSumOfWeightsEtaRP[i]  = dynamic_cast<TH1D*>(outputListHistos->FindObject(Form("Flow_SumOfWeights%sEtaRP_SP",weightFlag[i].Data())));
 457      pHistSumOfWeightsPtPOI[i]  = dynamic_cast<TH1D*>(outputListHistos->FindObject(Form("Flow_SumOfWeights%sPtPOI_SP",weightFlag[i].Data())));
 458      pHistSumOfWeightsEtaPOI[i] = dynamic_cast<TH1D*>(outputListHistos->FindObject(Form("Flow_SumOfWeights%sEtaPOI_SP",weightFlag[i].Data())));
 459     }
 460
 461     //pass the pointers to the task
 462     if (pCommonHist && pCommonHistResults && pHistProQaQb &&
 463         pHistProUQetaRP && pHistProUQetaPOI && pHistProUQPtRP && pHistProUQPtPOI
 464         && pHistSumOfLinearWeights && pHistSumOfQuadraticWeights && pHistProUQQaQbPtRP
 465         && pHistProUQQaQbEtaRP && pHistProUQQaQbPtPOI && pHistProUQQaQbEtaPOI) {
 466       this -> SetCommonHists(pCommonHist);
 467       this -> SetCommonHistsRes(pCommonHistResults);
 468       this -> SetHistProQaQb(pHistProQaQb);
 469       this -> SetHistProUQetaRP(pHistProUQetaRP);
 470       this -> SetHistProUQetaPOI(pHistProUQetaPOI);
 471       this -> SetHistProUQPtRP(pHistProUQPtRP);
 472       this -> SetHistProUQPtPOI(pHistProUQPtPOI);
 473       this -> SetHistProUQQaQbPtRP(pHistProUQQaQbPtRP);
 474       this -> SetHistProUQQaQbEtaRP(pHistProUQQaQbEtaRP);
 475       this -> SetHistProUQQaQbPtPOI(pHistProUQQaQbPtPOI);
 476       this -> SetHistProUQQaQbEtaPOI(pHistProUQQaQbEtaPOI);
 477       this -> SetHistSumOfLinearWeights(pHistSumOfLinearWeights);
 478       this -> SetHistSumOfQuadraticWeights(pHistSumOfQuadraticWeights);
 479      }
 480
 481    for(Int_t i=0;i<3;i++)
 482    {
 483     if(pHistSumOfWeightsPtRP[i]) this -> SetHistSumOfWeightsPtRP(pHistSumOfWeightsPtRP[i],i);
 484     if(pHistSumOfWeightsEtaRP[i]) this -> SetHistSumOfWeightsEtaRP(pHistSumOfWeightsEtaRP[i],i);
 485     if(pHistSumOfWeightsPtPOI[i]) this -> SetHistSumOfWeightsPtPOI(pHistSumOfWeightsPtPOI[i],i);
 486     if(pHistSumOfWeightsEtaPOI[i]) this -> SetHistSumOfWeightsEtaPOI(pHistSumOfWeightsEtaPOI[i],i);
 487    }
 488   } // end of if(outputListHistos)
 489 }
 490
 491 //--------------------------------------------------------------------
 492 void AliFlowAnalysisWithScalarProduct::Finish() {
 493
 494   //calculate flow and fill the AliFlowCommonHistResults
 495   if (fDebug) cout<<"AliFlowAnalysisWithScalarProduct::Finish()"<<endl;
 496
 497   cout<<"*************************************"<<endl;
 498   cout<<"*************************************"<<endl;
 499   cout<<"      Integrated flow from           "<<endl;
 500   cout<<"         Scalar product              "<<endl;
 501   cout<<endl;
 502
 503   Int_t iNbinsPt  = AliFlowCommonConstants::GetMaster()->GetNbinsPt();
 504   Int_t iNbinsEta = AliFlowCommonConstants::GetMaster()->GetNbinsEta();
 505
 506   //Calculate reference flow (noname)
 507   //----------------------------------
 508   //weighted average over (QaQb/MaMb) with weight (MaMb)
 509   Double_t dQaQb  = fHistProQaQb->GetBinContent(1);
 510   Double_t dV = 0.;
 511   if(dQaQb>=0.)
 512   {
 513    dV = TMath::Sqrt(dQaQb);
 514   }
 515   //statistical error of dQaQb:
 516   //  statistical error = term1 * spread * term2:
 517   //  term1 = sqrt{sum_{i=1}^{N} w^2}/(sum_{i=1}^{N} w)
 518   //  term2 = 1/sqrt(1-term1^2)
 519   Double_t dSpreadQaQb = fHistProQaQb->GetBinError(1);
 520   Double_t dSumOfLinearWeights = fHistSumOfLinearWeights->GetBinContent(1);
 521   Double_t dSumOfQuadraticWeights = fHistSumOfQuadraticWeights->GetBinContent(1);
 522   Double_t dTerm1 = 0.;
 523   Double_t dTerm2 = 0.;
 524   if(dSumOfLinearWeights) {
 525     dTerm1 = pow(dSumOfQuadraticWeights,0.5)/dSumOfLinearWeights;
 526   }
 527   if(1.-pow(dTerm1,2.) > 0.) {
 528     dTerm2 = 1./pow(1-pow(dTerm1,2.),0.5);
 529   }
 530   Double_t dStatErrorQaQb = dTerm1 * dSpreadQaQb * dTerm2;
 531   //calculate the statistical error
 532   Double_t dVerr = 0.;
 533   if(dQaQb > 0.) {
 534     dVerr = (1./(2.*pow(dQaQb,0.5)))*dStatErrorQaQb;
 535   }
 536   fCommonHistsRes->FillIntegratedFlow(dV,dVerr);
 537   cout<<"dV = "<<dV<<" +- "<<dVerr<<endl;
 538
 539   //Calculate differential flow and integrated flow (RP, POI)
 540   //---------------------------------------------------------
 541   //v as a function of eta for RP selection
 542   for(Int_t b=1;b<iNbinsEta+1;b++) {
 543     Double_t duQpro = fHistProUQetaRP->GetBinContent(b);
 544     Double_t dv2pro = 0.;
 545     if (dV!=0.) { dv2pro = duQpro/dV; }
 546     //calculate the statistical error
 547     Double_t dv2ProErr = CalculateStatisticalError(b, dStatErrorQaQb, fHistProUQetaRP, fHistProUQQaQbEtaRP, fHistSumOfWeightsEtaRP);
 548     //fill TH1D
 549     fCommonHistsRes->FillDifferentialFlowEtaRP(b, dv2pro, dv2ProErr);
 550   } //loop over bins b
 551
 552
 553   //v as a function of eta for POI selection
 554   for(Int_t b=1;b<iNbinsEta+1;b++) {
 555     Double_t duQpro = fHistProUQetaPOI->GetBinContent(b);
 556     Double_t dv2pro = 0.;
 557     if (dV!=0.) { dv2pro = duQpro/dV; }
 558     //calculate the statistical error
 559     Double_t dv2ProErr = CalculateStatisticalError(b, dStatErrorQaQb, fHistProUQetaPOI, fHistProUQQaQbEtaPOI, fHistSumOfWeightsEtaPOI);
 560
 561     //fill TH1D
 562     fCommonHistsRes->FillDifferentialFlowEtaPOI(b, dv2pro, dv2ProErr);
 563   } //loop over bins b
 564
 565
 566   //v as a function of Pt for RP selection
 567   TH1F* fHistPtRP = fCommonHists->GetHistPtRP(); //for calculating integrated flow
 568   Double_t dVRP = 0.;
 569   Double_t dSumRP = 0.;
 570   Double_t dErrVRP =0.;
 571
 572   for(Int_t b=1;b<iNbinsPt+1;b++) {
 573     Double_t duQpro = fHistProUQPtRP->GetBinContent(b);
 574     Double_t dv2pro = 0.;
 575     if (dV!=0.) { dv2pro = duQpro/dV; }
 576     //calculate the statistical error
 577     Double_t dv2ProErr = CalculateStatisticalError(b, dStatErrorQaQb, fHistProUQPtRP, fHistProUQQaQbPtRP, fHistSumOfWeightsPtRP);
 578
 579     //fill TH1D
 580     fCommonHistsRes->FillDifferentialFlowPtRP(b, dv2pro, dv2ProErr);
 581
 582     //calculate integrated flow for RP selection
 583     if (fHistPtRP){
 584       Double_t dYieldPt = fHistPtRP->GetBinContent(b);
 585       dVRP += dv2pro*dYieldPt;
 586       dSumRP +=dYieldPt;
 587       dErrVRP += dYieldPt*dYieldPt*dv2ProErr*dv2ProErr;
 588     } else { cout<<"fHistPtRP is NULL"<<endl; }
 589   } //loop over bins b
 590
 591   if (dSumRP != 0.) {
 592     dVRP /= dSumRP; //the pt distribution should be normalised
 593     dErrVRP /= (dSumRP*dSumRP);
 594     dErrVRP = TMath::Sqrt(dErrVRP);
 595   }
 596   fCommonHistsRes->FillIntegratedFlowRP(dVRP,dErrVRP);
 597   cout<<"dV(RP) = "<<dVRP<<" +- "<<dErrVRP<<endl;
 598
 599
 600   //v as a function of Pt for POI selection
 601   TH1F* fHistPtPOI = fCommonHists->GetHistPtPOI(); //for calculating integrated flow
 602   Double_t dVPOI = 0.;
 603   Double_t dSumPOI = 0.;
 604   Double_t dErrVPOI =0.;
 605
 606   for(Int_t b=1;b<iNbinsPt+1;b++) {
 607     Double_t duQpro = fHistProUQPtPOI->GetBinContent(b);
 608     Double_t dv2pro = 0.;
 609     if (dV!=0.) { dv2pro = duQpro/dV; }
 610     //calculate the statistical error
 611     Double_t dv2ProErr = CalculateStatisticalError(b, dStatErrorQaQb, fHistProUQPtPOI, fHistProUQQaQbPtPOI, fHistSumOfWeightsPtPOI);
 612
 613     //fill TH1D
 614     fCommonHistsRes->FillDifferentialFlowPtPOI(b, dv2pro, dv2ProErr);
 615
 616     //calculate integrated flow for POI selection
 617     if (fHistPtPOI){
 618       Double_t dYieldPt = fHistPtPOI->GetBinContent(b);
 619       dVPOI += dv2pro*dYieldPt;
 620       dSumPOI +=dYieldPt;
 621       dErrVPOI += dYieldPt*dYieldPt*dv2ProErr*dv2ProErr;
 622     } else { cout<<"fHistPtPOI is NULL"<<endl; }
 623   } //loop over bins b
 624
 625   if (dSumPOI != 0.) {
 626     dVPOI /= dSumPOI; //the pt distribution should be normalised
 627     dErrVPOI /= (dSumPOI*dSumPOI);
 628     dErrVPOI = TMath::Sqrt(dErrVPOI);
 629   }
 630   fCommonHistsRes->FillIntegratedFlowPOI(dVPOI,dErrVPOI);
 631   cout<<"dV(POI) = "<<dVPOI<<" +- "<<dErrVPOI<<endl;
 632
 633   cout<<endl;
 634   cout<<"*************************************"<<endl;
 635   cout<<"*************************************"<<endl;
 636
 637   //cout<<".....finished"<<endl;
 638 }
 639
 640
 641 //--------------------------------------------------------------------
 642 Double_t AliFlowAnalysisWithScalarProduct::CalculateStatisticalError(Int_t b, Double_t aStatErrorQaQb, TProfile* pHistProUQ, TProfile* pHistProUQQaQb, TH1D** pHistSumOfWeights) {
 643   //calculate the statistical error for differential flow for bin b
 644   Double_t duQproSpread = pHistProUQ->GetBinError(b);
 645   Double_t sumOfMq = pHistSumOfWeights[0]->GetBinContent(b);
 646   Double_t sumOfMqSquared = pHistSumOfWeights[1]->GetBinContent(b);
 647   Double_t dTerm1 = 0.;
 648   Double_t dTerm2 = 0.;
 649   if(sumOfMq) {
 650     dTerm1 = (pow(sumOfMqSquared,0.5)/sumOfMq);
 651   }
 652   if(1.-pow(dTerm1,2.)>0.) {
 653     dTerm2 = 1./pow(1.-pow(dTerm1,2.),0.5);
 654   }
 655   Double_t duQproErr = dTerm1*duQproSpread*dTerm2;
 656   // covariances:
 657   Double_t dTerm1Cov = pHistSumOfWeights[2]->GetBinContent(b);
 658   Double_t dTerm2Cov = fHistSumOfLinearWeights->GetBinContent(1);
 659   Double_t dTerm3Cov = sumOfMq;
 660   Double_t dWeightedCovariance = 0.;
 661   if(dTerm2Cov*dTerm3Cov>0.) {
 662     Double_t dDenominator = 1.-dTerm1Cov/(dTerm2Cov*dTerm3Cov);
 663     Double_t dPrefactor = dTerm1Cov/(dTerm2Cov*dTerm3Cov);
 664     if(dDenominator!=0) {
 665       Double_t dCovariance = (pHistProUQQaQb->GetBinContent(b)-fHistProQaQb->GetBinContent(1)*pHistProUQ->GetBinContent(b))/dDenominator;
 666       dWeightedCovariance = dCovariance*dPrefactor;
 667     }
 668   }
 669
 670   Double_t dv2ProErr = 0.; // final statitical error:
 671   if(fHistProQaQb->GetBinContent(1)>0.) {
 672     Double_t dv2ProErrorSquared = (1./4.)*pow(fHistProQaQb->GetBinContent(1),-3.)*
 673       (pow(pHistProUQ->GetBinContent(b),2.)*pow(aStatErrorQaQb,2.)
 674        + 4.*pow(fHistProQaQb->GetBinContent(1),2.)*pow(duQproErr,2.)
 675        - 4.*fHistProQaQb->GetBinContent(1)*pHistProUQ->GetBinContent(b)*dWeightedCovariance);
 676     if(dv2ProErrorSquared>0.) dv2ProErr = pow(dv2ProErrorSquared,0.5);
 677   }
 678
 679   return dv2ProErr;
 680 }
 681 //--------------------------------------------------------------------