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