PWG/FLOW/Base/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 "TFile.h"
  19 #include "TList.h"
  20 #include "TMath.h"
  21 #include "TString.h"
  22 #include "TProfile.h"
  23 #include "TVector2.h"
  24 #include "TH1D.h"
  25 #include "TH1F.h"
  26 #include "TH2D.h"
  27
  28 #include "AliFlowCommonConstants.h"
  29 #include "AliFlowEventSimple.h"
  30 #include "AliFlowVector.h"
  31 #include "AliFlowTrackSimple.h"
  32 #include "AliFlowCommonHist.h"
  33 #include "AliFlowCommonHistResults.h"
  34 #include "AliFlowAnalysisWithScalarProduct.h"
  35
  36 //////////////////////////////////////////////////////////////////////////////
  37 // AliFlowAnalysisWithScalarProduct:
  38 // Description: Maker to analyze Flow from the Scalar Product method.
  39 // authors: Naomi van del Kolk (kolk@nikhef.nl)
  40 //          Ante Bilandzic (anteb@nikhef.nl)
  41 // mods:    Carlos Perez (cperez@nikhef.nl)
  42 //////////////////////////////////////////////////////////////////////////////
  43
  44 ClassImp(AliFlowAnalysisWithScalarProduct)
  45
  46 //-----------------------------------------------------------------------
  47 AliFlowAnalysisWithScalarProduct::AliFlowAnalysisWithScalarProduct():
  48 fDebug(0),
  49 fUsePhiWeights(0),
  50 fApplyCorrectionForNUA(0),
  51 fHarmonic(2),
  52 fNormalizationType(1),
  53 fTotalQvector(3),
  54 fWeightsList(NULL),
  55 fHistList(NULL),
  56 fHistProConfig(NULL),
  57 fHistProQaQbNorm(NULL),
  58 fHistSumOfWeights(NULL),
  59 fHistProNUAq(NULL),
  60 fHistProQNorm(NULL),
  61 fHistProQaQb(NULL),
  62 fHistProQaQbM(NULL),
  63 fHistMaMb(NULL),
  64 fHistQNormQaQbNorm(NULL),
  65 fHistQaNormMa(NULL),
  66 fHistQbNormMb(NULL),
  67 fResolution(NULL),
  68 fHistQaQb(NULL),
  69 fHistQaQbCos(NULL),
  70 fCommonHists(NULL),
  71 fCommonHistsuQ(NULL),
  72 fCommonHistsRes(NULL)
  73 {
  74   //ctor
  75   for(int i=0; i!=2; ++i) {
  76     fPhiWeightsSub[i] = NULL;
  77     for(int j=0; j!=2; ++j) {
  78       fHistProUQ[i][j] = NULL;
  79       fHistProUQQaQb[i][j] = NULL;
  80       fHistProNUAu[i][j][0] = NULL;
  81       fHistProNUAu[i][j][1] = NULL;
  82       for(int k=0; k!=3; ++k)
  83         fHistSumOfWeightsu[i][j][k] = NULL;
  84     }
  85   }
  86 }
  87 //-----------------------------------------------------------------------
  88 AliFlowAnalysisWithScalarProduct::~AliFlowAnalysisWithScalarProduct()
  89 {
  90   //destructor
  91   delete fWeightsList;
  92   delete fHistList;
  93 }
  94 //-----------------------------------------------------------------------
  95 void AliFlowAnalysisWithScalarProduct::Init() {
  96   //Define all histograms
  97   printf("---Analysis with the Scalar Product Method--- Init\n");
  98   printf("--- fNormalizationType %d ---\n", fNormalizationType);
  99   //save old value and prevent histograms from being added to directory
 100   //to avoid name clashes in case multiple analaysis objects are used
 101   //in an analysis
 102   Bool_t oldHistAddStatus = TH1::AddDirectoryStatus();
 103   TH1::AddDirectory(kFALSE);
 104
 105   fHistList = new TList();
 106   fHistList->SetName("cobjSP");
 107   fHistList->SetOwner();
 108
 109   TList *uQRelated = new TList();
 110   uQRelated->SetName("uQ");
 111   uQRelated->SetOwner();
 112
 113   TList *nuaRelated = new TList();
 114   nuaRelated->SetName("NUA");
 115   nuaRelated->SetOwner();
 116
 117   TList *errorRelated = new TList();
 118   errorRelated->SetName("error");
 119   errorRelated->SetOwner();
 120
 121   TList *tQARelated = new TList();
 122   tQARelated->SetName("QA");
 123   tQARelated->SetOwner();
 124
 125   fCommonHists = new AliFlowCommonHist("AliFlowCommonHist_SP");
 126   (fCommonHists->GetHarmonic())->Fill(0.5,fHarmonic); // store harmonic
 127   fHistList->Add(fCommonHists);
 128
 129   fCommonHistsuQ = new AliFlowCommonHist("AliFlowCommonHist_uQ");
 130   (fCommonHistsuQ->GetHarmonic())->Fill(0.5,fHarmonic); // store harmonic
 131   fHistList->Add(fCommonHistsuQ);
 132
 133   fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResults_SP","",fHarmonic);
 134   fHistList->Add(fCommonHistsRes);
 135
 136   fHistProConfig = new TProfile("FlowPro_Flags_SP","Flow_Flags_SP",4,0.5,4.5,"s");
 137   fHistProConfig->GetXaxis()->SetBinLabel(1,"fApplyCorrectionForNUA");
 138   fHistProConfig->GetXaxis()->SetBinLabel(2,"fNormalizationType");
 139   fHistProConfig->GetXaxis()->SetBinLabel(3,"fUsePhiWeights");
 140   fHistProConfig->GetXaxis()->SetBinLabel(4,"fHarmonic");
 141   fHistProConfig->Fill(1,fApplyCorrectionForNUA);
 142   fHistProConfig->Fill(2,fNormalizationType);
 143   fHistProConfig->Fill(3,fUsePhiWeights);
 144   fHistProConfig->Fill(4,fHarmonic);
 145   fHistList->Add(fHistProConfig);
 146
 147   fHistProQaQbNorm = new TProfile("FlowPro_QaQbNorm_SP","FlowPro_QaQbNorm_SP", 1, 0.5, 1.5,"s");
 148   fHistProQaQbNorm->SetYTitle("<QaQb/Na/Nb>");
 149   errorRelated->Add(fHistProQaQbNorm);
 150
 151   fHistProNUAq = new TProfile("FlowPro_NUAq_SP","FlowPro_NUAq_SP", 6, 0.5, 6.5,"s");
 152   fHistProNUAq->GetXaxis()->SetBinLabel( 1,"<<sin(#Phi_{a})>>");
 153   fHistProNUAq->GetXaxis()->SetBinLabel( 2,"<<cos(#Phi_{a})>>");
 154   fHistProNUAq->GetXaxis()->SetBinLabel( 3,"<<sin(#Phi_{b})>>");
 155   fHistProNUAq->GetXaxis()->SetBinLabel( 4,"<<cos(#Phi_{b})>>");
 156   fHistProNUAq->GetXaxis()->SetBinLabel( 5,"<<sin(#Phi_{t})>>");
 157   fHistProNUAq->GetXaxis()->SetBinLabel( 6,"<<cos(#Phi_{t})>>");
 158   nuaRelated->Add(fHistProNUAq);
 159
 160   fHistSumOfWeights = new TH1D("Flow_SumOfWeights_SP","Flow_SumOfWeights_SP",2,0.5, 2.5);
 161   fHistSumOfWeights->GetXaxis()->SetBinLabel( 1,"#Sigma Na*Nb");
 162   fHistSumOfWeights->GetXaxis()->SetBinLabel( 2,"#Sigma (Na*Nb)^2");
 163   errorRelated->Add(fHistSumOfWeights);
 164
 165   TString sPOI[2] = {"RP","POI"}; // backward compatibility
 166   TString sEta[2] = {"Pt","eta"}; // backward compatibility
 167   TString sTitle[2] = {"p_{T} [GeV]","#eta"};
 168   TString sWeights[3] = {"uQ","uQuQ","uQQaQb"};
 169   Int_t iNbins[2];
 170   Double_t dMin[2], dMax[2];
 171   iNbins[0] = AliFlowCommonConstants::GetMaster()->GetNbinsPt();
 172   iNbins[1] = AliFlowCommonConstants::GetMaster()->GetNbinsEta();
 173   dMin[0]   = AliFlowCommonConstants::GetMaster()->GetPtMin();
 174   dMin[1]   = AliFlowCommonConstants::GetMaster()->GetEtaMin();
 175   dMax[0]   = AliFlowCommonConstants::GetMaster()->GetPtMax();
 176   dMax[1]   = AliFlowCommonConstants::GetMaster()->GetEtaMax();
 177   for(Int_t iPOI=0; iPOI!=2; ++iPOI) for(Int_t iSpace=0; iSpace!=2; ++iSpace) {
 178     // uQ
 179     fHistProUQ[iPOI][iSpace] = new TProfile( Form( "FlowPro_UQ_%s%s_SP", sEta[iSpace].Data(), sPOI[iPOI].Data() ),
 180                                        Form( "FlowPro_UQ%s%s_SP", sEta[iSpace].Data(), sPOI[iPOI].Data() ),
 181                                        iNbins[iSpace], dMin[iSpace], dMax[iSpace], "s");
 182     fHistProUQ[iPOI][iSpace]->SetXTitle(sTitle[iSpace].Data());
 183     fHistProUQ[iPOI][iSpace]->SetYTitle("<uQ>");
 184     uQRelated->Add(fHistProUQ[iPOI][iSpace]);
 185
 186     // NUAu
 187     fHistProNUAu[iPOI][iSpace][0] = new TProfile( Form("FlowPro_NUAu_%s%s_IM_SP", sEta[iSpace].Data(), sPOI[iPOI].Data() ),
 188                                            Form("FlowPro_NUAu_%s%s_IM_SP", sEta[iSpace].Data(), sPOI[iPOI].Data() ),
 189                                            iNbins[iSpace], dMin[iSpace], dMax[iSpace]);
 190     fHistProNUAu[iPOI][iSpace][0]->SetXTitle(sTitle[iSpace].Data());
 191     nuaRelated->Add(fHistProNUAu[iPOI][iSpace][0]);
 192     fHistProNUAu[iPOI][iSpace][1] = new TProfile( Form("FlowPro_NUAu_%s%s_RE_SP", sEta[iSpace].Data(), sPOI[iPOI].Data() ),
 193                                            Form("FlowPro_NUAu_%s%s_RE_SP", sEta[iSpace].Data(), sPOI[iPOI].Data() ),
 194                                            iNbins[iSpace], dMin[iSpace], dMax[iSpace]);
 195     fHistProNUAu[iPOI][iSpace][1]->SetXTitle(sTitle[iSpace].Data());
 196     nuaRelated->Add(fHistProNUAu[iPOI][iSpace][1]);
 197
 198     // uQ QaQb
 199     fHistProUQQaQb[iPOI][iSpace] = new TProfile( Form("FlowPro_UQQaQb_%s%s_SP", sEta[iSpace].Data(), sPOI[iPOI].Data() ),
 200                                            Form("FlowPro_UQQaQb_%s%s_SP", sEta[iSpace].Data(), sPOI[iPOI].Data() ),
 201                                            iNbins[iSpace], dMin[iSpace], dMax[iSpace]);
 202     fHistProUQQaQb[iPOI][iSpace]->SetXTitle(sTitle[iSpace].Data());
 203     fHistProUQQaQb[iPOI][iSpace]-> SetYTitle("<Qu QaQb>");
 204     errorRelated->Add(fHistProUQQaQb[iPOI][iSpace]);
 205
 206     // uWeights
 207     for(Int_t i=0; i!=3; ++i) {
 208       fHistSumOfWeightsu[iPOI][iSpace][i] = new TH1D( Form("Flow_SumOfWeights_%s%s_%s_SP",sWeights[i].Data(),sPOI[iPOI].Data(),sEta[iSpace].Data()),
 209                                                       Form("Flow_SumOfWeights_%s%s_%s_SP",sWeights[i].Data(),sPOI[iPOI].Data(),sEta[iSpace].Data()),
 210                                                       iNbins[iSpace], dMin[iSpace], dMax[iSpace]);
 211       fHistSumOfWeightsu[iPOI][iSpace][i]->SetYTitle(sWeights[i].Data());
 212       fHistSumOfWeightsu[iPOI][iSpace][i]->SetXTitle(sTitle[iSpace].Data());
 213       errorRelated->Add(fHistSumOfWeightsu[iPOI][iSpace][i]);
 214     }
 215   }
 216   //weights
 217   if(fUsePhiWeights) {
 218     if(!fWeightsList) {
 219       printf( "WARNING: fWeightsList is NULL in the Scalar Product method.\n" );
 220       exit(0);
 221     }
 222     fPhiWeightsSub[0] = dynamic_cast<TH1F*>
 223                         (fWeightsList->FindObject("phi_weights_sub0"));
 224     if(!fPhiWeightsSub[0]) {
 225       printf( "WARNING: phi_weights_sub0 not found in the Scalar Product method.\n" );
 226       exit(0);
 227     }
 228     nuaRelated->Add( fPhiWeightsSub[0] );
 229     fPhiWeightsSub[1] = dynamic_cast<TH1F*>
 230                       (fWeightsList->FindObject("phi_weights_sub1"));
 231     if(!fPhiWeightsSub[1]) {
 232       printf( "WARNING: phi_weights_sub1 not found in the Scalar Product method.\n" );
 233       exit(0);
 234     }
 235     nuaRelated->Add( fPhiWeightsSub[1] );
 236   }
 237
 238
 239   fHistProQNorm = new TProfile("FlowPro_QNorm_SP","FlowPro_QNorm_SP",       1,0.5,1.5,"s");
 240   fHistProQNorm->SetYTitle("<|Qa+Qb|>");
 241   tQARelated->Add(fHistProQNorm);
 242
 243   fHistProQaQb  = new TProfile("FlowPro_QaQb_SP","FlowPro_QaQb_SP",         1,0.5,1.5,"s");
 244   fHistProQaQb->SetYTitle("<QaQb>");
 245   tQARelated->Add(fHistProQaQb);
 246
 247   fHistProQaQbM = new TProfile("FlowPro_QaQbvsM_SP","FlowPro_QaQbvsM_SP",1000,0.0,10000);
 248   fHistProQaQbM->SetYTitle("<QaQb>");
 249   fHistProQaQbM->SetXTitle("M");
 250   fHistProQaQbM->Sumw2();
 251   tQARelated->Add(fHistProQaQbM);
 252
 253   fHistMaMb = new TH2D("Flow_MavsMb_SP","Flow_MavsMb_SP",100,0.,100.,100,0.,100.);
 254   fHistMaMb->SetYTitle("Ma");
 255   fHistMaMb->SetXTitle("Mb");
 256   tQARelated->Add(fHistMaMb);
 257
 258   fHistQNormQaQbNorm = new TH2D("Flow_QNormvsQaQbNorm_SP","Flow_QNormvsQaQbNorm_SP",88,-1.1,1.1,22,0.,1.1);
 259   fHistQNormQaQbNorm->SetYTitle("|Q/Mq|");
 260   fHistQNormQaQbNorm->SetXTitle("QaQb/MaMb");
 261   tQARelated->Add(fHistQNormQaQbNorm);
 262
 263   fHistQaNormMa = new TH2D("Flow_QaNormvsMa_SP","Flow_QaNormvsMa_SP",100,0.,100.,22,0.,1.1);
 264   fHistQaNormMa->SetYTitle("|Qa/Ma|");
 265   fHistQaNormMa->SetXTitle("Ma");
 266   tQARelated->Add(fHistQaNormMa);
 267
 268   fHistQbNormMb = new TH2D("Flow_QbNormvsMb_SP","Flow_QbNormvsMb_SP",100,0.,100.,22,0.,1.1);
 269   fHistQbNormMb->SetYTitle("|Qb/Mb|");
 270   fHistQbNormMb->SetXTitle("Mb");
 271   tQARelated->Add(fHistQbNormMb);
 272
 273   fResolution = new TH1D("Flow_resolution_SP","Flow_resolution_SP",100,-1.0,1.0);
 274   fResolution->SetYTitle("dN/d(Cos2(#phi_a - #phi_b))");
 275   fResolution->SetXTitle("Cos2(#phi_a - #phi_b)");
 276   tQARelated->Add(fResolution);
 277
 278   fHistQaQb = new TH1D("Flow_QaQb_SP","Flow_QaQb_SP",200,-100.,100.);
 279   fHistQaQb->SetYTitle("dN/dQaQb");
 280   fHistQaQb->SetXTitle("dQaQb");
 281   tQARelated->Add(fHistQaQb);
 282
 283   fHistQaQbCos = new TH1D("Flow_QaQbCos_SP","Flow_QaQbCos_SP",63,0.,TMath::Pi());
 284   fHistQaQbCos->SetYTitle("dN/d#phi");
 285   fHistQaQbCos->SetXTitle("#phi");
 286   tQARelated->Add(fHistQaQbCos);
 287
 288   fHistList->Add(uQRelated);
 289   fHistList->Add(nuaRelated);
 290   fHistList->Add(errorRelated);
 291   fHistList->Add(tQARelated);
 292
 293   TH1::AddDirectory(oldHistAddStatus);
 294 }
 295
 296 //-----------------------------------------------------------------------
 297 void AliFlowAnalysisWithScalarProduct::Make(AliFlowEventSimple* anEvent) {
 298   // Scalar Product method
 299   if (!anEvent) return; // for coverity
 300
 301   // Get Q vectors for the subevents
 302   AliFlowVector* vQarray = new AliFlowVector[2];
 303   if (fUsePhiWeights)
 304     anEvent->Get2Qsub(vQarray,fHarmonic,fWeightsList,kTRUE);
 305   else
 306     anEvent->Get2Qsub(vQarray,fHarmonic);
 307   // Subevent a
 308   AliFlowVector vQa = vQarray[0];
 309   // Subevent b
 310   AliFlowVector vQb = vQarray[1];
 311   delete [] vQarray;
 312
 313   Double_t dMa = vQa.GetMult();
 314   if( dMa < 2 ) return;
 315   Double_t dMb = vQb.GetMult();
 316   if( dMb < 2 ) return;
 317   //fill control histograms
 318   if (fUsePhiWeights) {
 319     fCommonHists->FillControlHistograms(anEvent,fWeightsList,kTRUE);
 320   } else {
 321     fCommonHists->FillControlHistograms(anEvent);
 322   }
 323
 324   //Normalizing: weight the Q vectors for the subevents
 325   Double_t dNa = fNormalizationType ? dMa: vQa.Mod(); // SP corresponds to true
 326   Double_t dNb = fNormalizationType ? dMb: vQb.Mod(); // SP corresponds to true
 327   Double_t dWa = fNormalizationType ? dMa: 1; // SP corresponds to true
 328   Double_t dWb = fNormalizationType ? dMb: 1; // SP corresponds to true
 329
 330   //Scalar product of the two subevents vectors
 331   Double_t dQaQb = (vQa*vQb);
 332
 333   //printf("==============\n");
 334   //printf("vQa: { %f, %f }\n",vQa.X(),vQa.Y());
 335   //printf("QaQb/|Qa||Qb|: %f\n",dQaQb/vQa.Mod()/vQb.Mod());
 336   //printf("QaQb/|Ma||Mb|: %f\n",dQaQb/dMa/dMb);
 337
 338   //      01    10     11   <===   fTotalQVector
 339   // Q ?= Qa or Qb or QaQb
 340   AliFlowVector vQm;
 341   vQm.Set(0.0,0.0);
 342   Double_t dNq=0;
 343   if( (fTotalQvector%2)>0 ) { // 01 or 11
 344     vQm += vQa;
 345     dNq += dMa;
 346   }
 347   if( fTotalQvector>1 ) { // 10 or 11
 348     vQm += vQb;
 349     dNq += dMb;
 350   }
 351   Double_t dWq = fNormalizationType ? dNq: 1; // SP corresponds to true
 352   dNq = fNormalizationType ? dNq: vQm.Mod(); // SP corresponds to true
 353
 354   //fill some EP control histograms
 355   fHistProQaQbNorm->Fill(1.,dQaQb/dNa/dNb,dWa*dWb);  //Fill (QaQb/NaNb) with weight (WaWb).
 356   //needed for the error calculation:
 357   fHistSumOfWeights -> Fill(1.,dWa*dWb);
 358   fHistSumOfWeights -> Fill(2.,pow(dWa*dWb,2.));
 359   //needed for correcting non-uniform acceptance:
 360   fHistProNUAq->Fill(1.,vQa.Y()/dNa,dWa); // to get <<sin(phi_a)>>
 361   fHistProNUAq->Fill(2.,vQa.X()/dNa,dWa); // to get <<cos(phi_a)>>
 362   fHistProNUAq->Fill(3.,vQb.Y()/dNb,dWb); // to get <<sin(phi_b)>>
 363   fHistProNUAq->Fill(4.,vQb.X()/dNb,dWb); // to get <<cos(phi_b)>>
 364   fHistProNUAq->Fill(5.,vQm.Y()/dNq,dWq);
 365   fHistProNUAq->Fill(6.,vQm.X()/dNq,dWq);
 366
 367   //loop over the tracks of the event
 368   AliFlowTrackSimple*   pTrack = NULL;
 369   Int_t iNumberOfTracks = anEvent->NumberOfTracks();
 370   for (Int_t i=0;i<iNumberOfTracks;i++) {
 371     pTrack = anEvent->GetTrack(i) ;
 372     if (!pTrack) continue;
 373     Double_t dPhi = pTrack->Phi();
 374     Double_t dPt  = pTrack->Pt();
 375     Double_t dEta = pTrack->Eta();
 376
 377     //calculate vU
 378     TVector2 vU;
 379     Double_t dUX = TMath::Cos(fHarmonic*dPhi);
 380     Double_t dUY = TMath::Sin(fHarmonic*dPhi);
 381     vU.Set(dUX,dUY);
 382
 383     //      01    10     11   <===   fTotalQVector
 384     // Q ?= Qa or Qb or QaQb
 385     vQm.Set(0.0,0.0); //start the loop fresh
 386     Double_t dMq=0;
 387     if( (fTotalQvector%2)>0 ) { // 01 or 11
 388       vQm += vQa;
 389       dMq += dMa;
 390     }
 391     if( fTotalQvector>1 ) { // 10 or 11
 392       vQm += vQb;
 393       dMq += dMb;
 394     }
 395
 396     //remove track if in subevent
 397     for(Int_t inSubEvent=0; inSubEvent!=2; ++inSubEvent) {
 398       if( !pTrack->InSubevent( inSubEvent ) )
 399         continue;
 400       if(inSubEvent==0)
 401         if( (fTotalQvector%2)!=1 )
 402           continue;
 403       if(inSubEvent==1)
 404         if( fTotalQvector>0 )
 405           continue;
 406       Double_t dW=1;
 407       Double_t dPhiCenter = dPhi;
 408       //if phi weights are used
 409       if(fUsePhiWeights && fPhiWeightsSub[inSubEvent]) {
 410         Int_t iNBinsPhiSub = fPhiWeightsSub[inSubEvent]->GetNbinsX();
 411         Int_t phiBin = 1+(Int_t)(TMath::Floor(dPhi*iNBinsPhiSub/TMath::TwoPi()));
 412         dW = fPhiWeightsSub[inSubEvent]->GetBinContent(phiBin);
 413         dPhiCenter = fPhiWeightsSub[inSubEvent]->GetBinCenter(phiBin);
 414       }
 415       Double_t dQmX = vQm.X() - dW*(pTrack->Weight())* TMath::Cos(fHarmonic*dPhiCenter);
 416       Double_t dQmY = vQm.Y() - dW*(pTrack->Weight())* TMath::Sin(fHarmonic*dPhiCenter);
 417       vQm.Set(dQmX,dQmY);
 418       dMq = dMq-dW*pTrack->Weight();
 419     }
 420     dNq = fNormalizationType ? dMq : vQm.Mod();
 421     dWq = fNormalizationType ? dMq : 1;
 422
 423     //Filling QA (for compatibility with previous version)
 424     fHistProQaQb->Fill(1,vQa*vQb,1);
 425     fHistProQaQbM->Fill(anEvent->GetNumberOfRPs()+0.5,(vQa*vQb)/dMa/dMb,dMa*dMb);
 426     fHistQaQbCos->Fill(TMath::ACos((vQa*vQb)/vQa.Mod()/vQb.Mod()));
 427     fResolution->Fill( TMath::Cos( vQa.Phi()-vQb.Phi() ) );
 428     fHistQaQb->Fill(vQa*vQb);
 429     fHistMaMb->Fill(dMb,dMa);
 430     fHistProQNorm->Fill(1,vQm.Mod()/dMq,dMq);
 431     fHistQNormQaQbNorm->Fill((vQa*vQb)/dMa/dMb,vQm.Mod()/dMq);
 432     fHistQaNormMa->Fill(dMa,vQa.Mod()/dMa);
 433     fHistQbNormMb->Fill(dMb,vQb.Mod()/dMb);
 434
 435     Double_t dUQ = vU*vQm;
 436
 437     //fill the profile histograms
 438     for(Int_t iPOI=0; iPOI!=2; ++iPOI) {
 439       if( (iPOI==0)&&(!pTrack->InRPSelection()) )
 440         continue;
 441       if( (iPOI==1)&&(!pTrack->InPOISelection()) )
 442         continue;
 443       fHistProUQ[iPOI][0]->Fill(dPt ,dUQ/dNq,dWq); //Fill (uQ/Nq') with weight (Nq')
 444       fHistProUQ[iPOI][1]->Fill(dEta,dUQ/dNq,dWq); //Fill (uQ/Nq') with weight (Nq')
 445       //needed for the error calculation:
 446       fHistProUQQaQb[iPOI][0]-> Fill(dPt ,dUQ/dNq*dQaQb/dNa/dNb,dWq*dWa*dWb); //Fill [Qu/Nq']*[QaQb/NaNb] with weight (Nq')NaNb
 447       fHistProUQQaQb[iPOI][1]-> Fill(dEta,dUQ/dNq*dQaQb/dNa/dNb,dWq*dWa*dWb); //Fill [Qu/Nq']*[QaQb/NaNb] with weight (Nq')NaNb
 448       fHistSumOfWeightsu[iPOI][0][0]->Fill(dPt ,dWq);        // sum of Nq'
 449       fHistSumOfWeightsu[iPOI][0][1]->Fill(dPt ,pow(dWq,2.));// sum of Nq'^2
 450       fHistSumOfWeightsu[iPOI][0][2]->Fill(dPt ,dWq*dWa*dWb);// sum of Nq'*Na*Nb
 451       fHistSumOfWeightsu[iPOI][1][0]->Fill(dEta,dWq);        // sum of Nq'
 452       fHistSumOfWeightsu[iPOI][1][1]->Fill(dEta,pow(dWq,2.));// sum of Nq'^2
 453       fHistSumOfWeightsu[iPOI][1][2]->Fill(dEta,dNq*dWa*dWb);// sum of Nq'*Na*Nb
 454       //NUA:
 455       fHistProNUAu[iPOI][0][0]->Fill(dPt,dUY,1.); //sin u
 456       fHistProNUAu[iPOI][0][1]->Fill(dPt,dUX,1.); //cos u
 457       fHistProNUAu[iPOI][1][0]->Fill(dEta,dUY,1.); //sin u
 458       fHistProNUAu[iPOI][1][1]->Fill(dEta,dUX,1.); //cos u
 459     }
 460   }//loop over tracks
 461
 462 }
 463
 464 //--------------------------------------------------------------------
 465 void AliFlowAnalysisWithScalarProduct::GetOutputHistograms(TList *outputListHistos){
 466   //get pointers to all output histograms (called before Finish())
 467   fHistList = outputListHistos;
 468
 469   fCommonHists = (AliFlowCommonHist*) fHistList->FindObject("AliFlowCommonHist_SP");
 470   fCommonHistsuQ = (AliFlowCommonHist*) fHistList->FindObject("AliFlowCommonHist_uQ");
 471   fCommonHistsRes = (AliFlowCommonHistResults*) fHistList->FindObject("AliFlowCommonHistResults_SP");
 472   fHistProConfig = (TProfile*) fHistList->FindObject("FlowPro_Flags_SP");
 473   if(!fHistProConfig) printf("Error loading fHistProConfig\n");
 474   TList *uQ = (TList*) fHistList->FindObject("uQ");
 475   TList *nua = (TList*) fHistList->FindObject("NUA");
 476   TList *error = (TList*) fHistList->FindObject("error");
 477
 478   fHistProQaQbNorm = (TProfile*) error->FindObject("FlowPro_QaQbNorm_SP");
 479   if(!fHistProQaQbNorm) printf("Error loading fHistProQaQbNorm\n");
 480   fHistProNUAq = (TProfile*) nua->FindObject("FlowPro_NUAq_SP");
 481   if(!fHistProNUAq) printf("Error loading fHistProNUAq\n");
 482   fHistSumOfWeights = (TH1D*) error->FindObject("Flow_SumOfWeights_SP");
 483   if(!fHistSumOfWeights) printf("Error loading fHistSumOfWeights\n");
 484
 485   TString sPOI[2] = {"RP","POI"};
 486   TString sEta[2] = {"Pt","eta"};
 487   TString sWeights[3] = {"uQ","uQuQ","uQQaQb"};
 488   for(Int_t iPOI=0; iPOI!=2; ++iPOI) for(Int_t iSpace=0; iSpace!=2; ++iSpace) {
 489     fHistProUQ[iPOI][iSpace] = (TProfile*) uQ->FindObject( Form( "FlowPro_UQ_%s%s_SP", sEta[iSpace].Data(), sPOI[iPOI].Data() ) );
 490     if(!fHistProUQ[iPOI][iSpace]) printf("Error loading fHistProUQ[%d][%d]\n",iPOI,iSpace);
 491     fHistProNUAu[iPOI][iSpace][0] = (TProfile*) nua->FindObject( Form("FlowPro_NUAu_%s%s_IM_SP", sEta[iSpace].Data(), sPOI[iPOI].Data() ) );
 492     if(!fHistProNUAu[iPOI][iSpace][0]) printf("Error loading fHistProNUAu[%d][%d][0]\n",iPOI,iSpace);
 493     fHistProNUAu[iPOI][iSpace][1] = (TProfile*) nua->FindObject( Form("FlowPro_NUAu_%s%s_RE_SP", sEta[iSpace].Data(), sPOI[iPOI].Data() ) );
 494     if(!fHistProNUAu[iPOI][iSpace][1]) printf("Error loading fHistProNUAu[%d][%d][1]\n",iPOI,iSpace);
 495     fHistProUQQaQb[iPOI][iSpace] = (TProfile*) error->FindObject( Form("FlowPro_UQQaQb_%s%s_SP", sEta[iSpace].Data(), sPOI[iPOI].Data() ) );
 496     for(Int_t i=0; i!=3; ++i){
 497       fHistSumOfWeightsu[iPOI][iSpace][i] = (TH1D*) error->FindObject( Form("Flow_SumOfWeights_%s%s_%s_SP",sWeights[i].Data(),sPOI[iPOI].Data(),sEta[iSpace].Data()) );
 498       if(!fHistSumOfWeightsu[iPOI][iSpace][i]) printf("Error loading fHistSumOfWeightsu[%d][%d][%d]\n",iPOI,iSpace,i);
 499     }
 500   }
 501   if(fHistProConfig) {
 502     fApplyCorrectionForNUA = (Int_t) fHistProConfig->GetBinContent(1);
 503     fNormalizationType  = (Int_t) fHistProConfig->GetBinContent(2);
 504     fUsePhiWeights = (Int_t) fHistProConfig->GetBinContent(3);
 505     fHarmonic = (Int_t) fHistProConfig->GetBinContent(4);
 506   }
 507 }
 508
 509 //--------------------------------------------------------------------
 510 void AliFlowAnalysisWithScalarProduct::Finish() {
 511   //calculate flow and fill the AliFlowCommonHistResults
 512   printf("AliFlowAnalysisWithScalarProduct::Finish()\n");
 513
 514   // access harmonic:
 515   fApplyCorrectionForNUA = fHistProConfig->GetBinContent(1);
 516   fNormalizationType = fHistProConfig->GetBinContent(2);
 517   fHarmonic = fHistProConfig->GetBinContent(4);
 518
 519   printf("*************************************\n");
 520   printf("*************************************\n");
 521   printf("      Integrated flow from           \n");
 522   printf("         Scalar Product              \n\n");
 523   if(!fNormalizationType)
 524     printf("          (BehaveAsEP)               \n\n");
 525
 526   //Calculate reference flow
 527   //----------------------------------
 528   //weighted average over (QaQb/NaNb) with weight (NaNb)
 529   Double_t dEntriesQaQb = fHistProQaQbNorm->GetEntries();
 530   if( dEntriesQaQb < 1 )
 531     return;
 532   Double_t dQaQb  = fHistProQaQbNorm->GetBinContent(1);
 533   if( dQaQb < 0 )
 534     return;
 535   Double_t dSpreadQaQb = fHistProQaQbNorm->GetBinError(1);
 536
 537   //NUA qq:
 538   Double_t dImQa = fHistProNUAq->GetBinContent(1);
 539   Double_t dReQa = fHistProNUAq->GetBinContent(2);
 540   Double_t dImQb = fHistProNUAq->GetBinContent(3);
 541   Double_t dReQb = fHistProNUAq->GetBinContent(4);
 542   if(fApplyCorrectionForNUA)
 543     dQaQb = dQaQb - dImQa*dImQb - dReQa*dReQb;
 544   printf("QaQb = %f +- %f\n", dQaQb, (dSpreadQaQb/TMath::Sqrt(dEntriesQaQb)) );
 545   Double_t dV = TMath::Sqrt(dQaQb);
 546
 547   printf("ResSub = %f\n", dV );
 548   printf("fTotalQvector %d \n",fTotalQvector);
 549   if(!fNormalizationType) {
 550     if(fTotalQvector>2) {
 551       dV = ComputeResolution( TMath::Sqrt2()*FindXi(dV,1e-6) );
 552       printf("An estimate of the event plane resolution is: %f\n", dV );
 553     }
 554   }
 555   printf("ResTot = %f\n", dV );
 556   //statistical error of dQaQb:
 557   //  statistical error = term1 * spread * term2:
 558   //  term1 = sqrt{sum_{i=1}^{N} w^2}/(sum_{i=1}^{N} w)
 559   //  term2 = 1/sqrt(1-term1^2)
 560   Double_t dSumOfLinearWeights = fHistSumOfWeights->GetBinContent(1);
 561   Double_t dSumOfQuadraticWeights = fHistSumOfWeights->GetBinContent(2);
 562   Double_t dTerm1 = 0.;
 563   Double_t dTerm2 = 0.;
 564   if(dSumOfLinearWeights)
 565     dTerm1 = pow(dSumOfQuadraticWeights,0.5)/dSumOfLinearWeights;
 566   if(1.-pow(dTerm1,2.) > 0.)
 567     dTerm2 = 1./pow(1-pow(dTerm1,2.),0.5);
 568   Double_t dStatErrorQaQb = dTerm1 * dSpreadQaQb * dTerm2;
 569   Double_t dVerr = 0.;
 570   if(dQaQb > 0.)
 571     dVerr = (1./(2.*pow(dQaQb,0.5)))*dStatErrorQaQb;
 572   fCommonHistsRes->FillIntegratedFlow(dV,dVerr);
 573   printf("v%d(subevents) = %f +- %f\n",fHarmonic,dV,dVerr);
 574
 575   Int_t iNbins[2];
 576   iNbins[0] = AliFlowCommonConstants::GetMaster()->GetNbinsPt();
 577   iNbins[1] = AliFlowCommonConstants::GetMaster()->GetNbinsEta();
 578
 579   //Calculate differential flow and integrated flow (RP, POI)
 580   //---------------------------------------------------------
 581   //v as a function of eta for RP selection
 582   for(Int_t iRFPorPOI=0; iRFPorPOI != 2; ++iRFPorPOI)
 583   for(Int_t iPTorETA=0; iPTorETA != 2; ++iPTorETA)
 584   for(Int_t b=1; b != iNbins[iPTorETA]+1; ++b) {
 585     Double_t duQpro = fHistProUQ[iRFPorPOI][iPTorETA]->GetBinContent(b);
 586     if(fApplyCorrectionForNUA)
 587       duQpro = duQpro
 588         - fHistProNUAu[iRFPorPOI][iPTorETA][1]->GetBinContent(b)*fHistProNUAq->GetBinContent(6)
 589         - fHistProNUAu[iRFPorPOI][iPTorETA][0]->GetBinContent(b)*fHistProNUAq->GetBinContent(5);
 590     Double_t dv2pro = -999.;
 591     if( TMath::Abs(dV!=0.) ) { dv2pro = duQpro/dV; }
 592     //calculate the statistical error
 593     Double_t dv2ProErr = CalculateStatisticalError(iRFPorPOI, iPTorETA, b, dStatErrorQaQb);
 594     if( (iRFPorPOI==0)&&(iPTorETA==0) )
 595       fCommonHistsRes->FillDifferentialFlowPtRP(  b, dv2pro, dv2ProErr);
 596     if( (iRFPorPOI==0)&&(iPTorETA==1) )
 597       fCommonHistsRes->FillDifferentialFlowEtaRP( b, dv2pro, dv2ProErr);
 598     if( (iRFPorPOI==1)&&(iPTorETA==0) )
 599       fCommonHistsRes->FillDifferentialFlowPtPOI( b, dv2pro, dv2ProErr);
 600     if( (iRFPorPOI==1)&&(iPTorETA==1) )
 601       fCommonHistsRes->FillDifferentialFlowEtaPOI(b, dv2pro, dv2ProErr);
 602     //printf("POI %d | PT %d >>> %f +- %f\n",iRFPorPOI,iPTorETA,dv2pro,dv2ProErr);
 603   }
 604
 605   printf("\n");
 606   printf("*************************************\n");
 607   printf("*************************************\n");
 608 }
 609
 610 //-----------------------------------------------------------------------
 611 void AliFlowAnalysisWithScalarProduct::WriteHistograms(TDirectoryFile *outputFileName) const {
 612  //store the final results in output .root file
 613  outputFileName->Add(fHistList);
 614  outputFileName->Write(outputFileName->GetName(), TObject::kSingleKey);
 615 }
 616
 617 //--------------------------------------------------------------------
 618 Double_t AliFlowAnalysisWithScalarProduct::CalculateStatisticalError(Int_t iRFPorPOI, Int_t iPTorETA, Int_t b, Double_t aStatErrorQaQb) const {
 619   //calculate the statistical error for differential flow for bin b
 620   Double_t duQproSpread = fHistProUQ[iRFPorPOI][iPTorETA]->GetBinError(b);
 621   Double_t sumOfMq = fHistSumOfWeightsu[iRFPorPOI][iPTorETA][0]->GetBinContent(b);
 622   Double_t sumOfMqSquared = fHistSumOfWeightsu[iRFPorPOI][iPTorETA][1]->GetBinContent(b);
 623   Double_t dQaQb = fHistProQaQbNorm->GetBinContent(1);
 624
 625   //non-isotropic terms:
 626   if(fApplyCorrectionForNUA) {
 627     Double_t dImQa = fHistProNUAq->GetBinContent(1);  // <<sin(phi_a)>>
 628     Double_t dReQa = fHistProNUAq->GetBinContent(2);  // <<cos(phi_a)>>
 629     Double_t dImQb = fHistProNUAq->GetBinContent(3);  // <<sin(phi_b)>>
 630     Double_t dReQb = fHistProNUAq->GetBinContent(4);  // <<cos(phi_b)>>
 631     dQaQb = dQaQb - dImQa*dImQb - dReQa*dReQb;
 632   }
 633
 634   Double_t dTerm1 = 0.;
 635   Double_t dTerm2 = 0.;
 636   if(sumOfMq) {
 637     dTerm1 = (pow(sumOfMqSquared,0.5)/sumOfMq);
 638   }
 639   if(1.-pow(dTerm1,2.)>0.) {
 640     dTerm2 = 1./pow(1.-pow(dTerm1,2.),0.5);
 641   }
 642   Double_t duQproErr = dTerm1*duQproSpread*dTerm2;
 643   // covariances:
 644   Double_t dTerm1Cov = fHistSumOfWeightsu[iRFPorPOI][iPTorETA][2]->GetBinContent(b);
 645   Double_t dTerm2Cov = fHistSumOfWeights->GetBinContent(1);
 646   Double_t dTerm3Cov = sumOfMq;
 647   Double_t dWeightedCovariance = 0.;
 648   if(dTerm2Cov*dTerm3Cov>0.) {
 649     Double_t dDenominator = 1.-dTerm1Cov/(dTerm2Cov*dTerm3Cov);
 650     Double_t dPrefactor = dTerm1Cov/(dTerm2Cov*dTerm3Cov);
 651     if(dDenominator!=0) {
 652       Double_t dCovariance = ( fHistProUQQaQb[iRFPorPOI][iPTorETA]->GetBinContent(b)-dQaQb*fHistProUQ[iRFPorPOI][iPTorETA]->GetBinContent(b))/dDenominator;
 653       dWeightedCovariance = dCovariance*dPrefactor;
 654     }
 655   }
 656   Double_t dv2ProErr = 0.; // final statitical error:
 657   if(dQaQb>0.) {
 658     Double_t dv2ProErrorSquared = (1./4.)*pow(dQaQb,-3.)*
 659       (pow(fHistProUQ[iRFPorPOI][iPTorETA]->GetBinContent(b),2.)*pow(aStatErrorQaQb,2.)
 660        + 4.*pow(dQaQb,2.)*pow(duQproErr,2.)
 661        - 4.*dQaQb*fHistProUQ[iRFPorPOI][iPTorETA]->GetBinContent(b)*dWeightedCovariance);
 662     if(dv2ProErrorSquared>0.) dv2ProErr = pow(dv2ProErrorSquared,0.5);
 663   }
 664   return dv2ProErr;
 665 }
 666
 667 Double_t AliFlowAnalysisWithScalarProduct::ComputeResolution( Double_t x ) const {
 668   // Computes resolution for Event Plane method
 669   if(x > 51.3) {
 670     printf("Warning: Estimation of total resolution might be WRONG. Please check!");
 671     return 0.99981;
 672   }
 673   Double_t a = x*x/4;
 674   Double_t b = TMath::Exp(-a)*TMath::BesselI0(a)+TMath::Exp(-a)*TMath::BesselI1(a);
 675   return TMath::Sqrt(TMath::PiOver2())/2*x*b;
 676 }
 677
 678 Double_t AliFlowAnalysisWithScalarProduct::FindXi( Double_t res, Double_t prec ) const {
 679   // Computes x(res) for Event Plane method
 680   if(res > 0.99981) {
 681     printf("Warning: Resolution for subEvent is high. You reached the precision limit.");
 682     return 51.3;
 683   }
 684   int nSteps =0;
 685   Double_t xtmp=0, xmin=0, xmax=51.3, rtmp=0, delta=2*prec;
 686   while( delta > prec ) {
 687     xtmp = 0.5*(xmin+xmax);
 688     rtmp = ComputeResolution(xtmp);
 689     delta = TMath::Abs( res-rtmp );
 690     if(rtmp>res) xmax = xtmp;
 691     if(rtmp<res) xmin = xtmp;
 692     nSteps++;
 693   }
 694   return xtmp;
 695 }
 696