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     Double_t dMass = pTrack->Mass();
 384
 385     //calculate vU
 386     TVector2 vU;
 387     Double_t dUX = TMath::Cos(fHarmonic*dPhi);
 388     Double_t dUY = TMath::Sin(fHarmonic*dPhi);
 389     vU.Set(dUX,dUY);
 390
 391     //      01    10     11   <===   fTotalQVector
 392     // Q ?= Qa or Qb or QaQb
 393     vQm.Set(0.0,0.0); //start the loop fresh
 394     Double_t dMq=0;
 395     if( (fTotalQvector%2)>0 ) { // 01 or 11
 396       vQm += vQa;
 397       dMq += dMa;
 398     }
 399     if( fTotalQvector>1 ) { // 10 or 11
 400       vQm += vQb;
 401       dMq += dMb;
 402     }
 403
 404     //remove track if in subevent
 405     for(Int_t inSubEvent=0; inSubEvent<2; ++inSubEvent) {
 406       if( !pTrack->InSubevent( inSubEvent ) )
 407         continue;
 408       if(inSubEvent==0)
 409         if( (fTotalQvector%2)!=1 )
 410           continue;
 411       if(inSubEvent==1)
 412         if( fTotalQvector<2 )
 413           continue;
 414       Double_t dW=1;
 415       //Double_t dPhiCenter = dPhi;
 416       //if phi weights are used
 417       if(fUsePhiWeights && fPhiWeightsSub[inSubEvent]) {
 418         Int_t iNBinsPhiSub = fPhiWeightsSub[inSubEvent]->GetNbinsX();
 419         Int_t phiBin = 1+(Int_t)(TMath::Floor(dPhi*iNBinsPhiSub/TMath::TwoPi()));
 420         dW = fPhiWeightsSub[inSubEvent]->GetBinContent(phiBin);
 421         //dPhiCenter = fPhiWeightsSub[inSubEvent]->GetBinCenter(phiBin);
 422       }
 423       //Double_t dQmX = vQm.X() - dW*(pTrack->Weight())* TMath::Cos(fHarmonic*dPhiCenter);
 424       //Double_t dQmY = vQm.Y() - dW*(pTrack->Weight())* TMath::Sin(fHarmonic*dPhiCenter);
 425       //vQm.Set(dQmX,dQmY);
 426
 427       //subtrack the track from the Q vector, but only if it was used to construct this
 428       //Q vector: i.e. check wether it has the same tags and is in the same subevent
 429       //this is especially important for the daughters (as for the mother it is already checked)
 430       Int_t numberOfsubtractedDaughters=vQm.SubtractTrackWithDaughters(pTrack,dW);
 431
 432       if(!fMinimalBook) {
 433         fHistNumberOfSubtractedDaughters->Fill(numberOfsubtractedDaughters);
 434       }
 435
 436       dMq = dMq-dW*pTrack->Weight();
 437     }
 438     dNq = fNormalizationType ? dMq : vQm.Mod();
 439     dWq = fNormalizationType ? dMq : 1;
 440
 441     //Filling QA (for compatibility with previous version)
 442     if(!fMinimalBook) {
 443       fHistProQaQb->Fill(1,vQa*vQb,1);
 444       fHistProQaQbM->Fill(anEvent->GetNumberOfRPs()+0.5,(vQa*vQb)/dMa/dMb,dMa*dMb);
 445       fHistQaQbCos->Fill(TMath::ACos((vQa*vQb)/vQa.Mod()/vQb.Mod()));
 446       fResolution->Fill( TMath::Cos( vQa.Phi()-vQb.Phi() ) );
 447       fHistQaQb->Fill(vQa*vQb);
 448       fHistMaMb->Fill(dMb,dMa);
 449       fHistProQNorm->Fill(1,vQm.Mod()/dMq,dMq);
 450       fHistQNormQaQbNorm->Fill((vQa*vQb)/dMa/dMb,vQm.Mod()/dMq);
 451       fHistQaNormMa->Fill(dMa,vQa.Mod()/dMa);
 452       fHistQbNormMb->Fill(dMb,vQb.Mod()/dMb);
 453     }
 454
 455     Double_t dUQ = vU*vQm;
 456
 457     //fill the profile histograms
 458     for(Int_t iPOI=0; iPOI!=2; ++iPOI) {
 459       if( (iPOI==0)&&(!pTrack->InRPSelection()) )
 460         continue;
 461       if( (iPOI==1)&&(!pTrack->InPOISelection(fPOItype)) )
 462         continue;
 463       fHistProUQ[iPOI][0]->Fill(dPt ,dUQ/dNq,dWq); //Fill (uQ/Nq') with weight (Nq')
 464       fHistProUQ[iPOI][1]->Fill(dEta,dUQ/dNq,dWq); //Fill (uQ/Nq') with weight (Nq')
 465       //needed for the error calculation:
 466       fHistProUQQaQb[iPOI][0]-> Fill(dPt ,dUQ/dNq*dQaQb/dNa/dNb,dWq*dWa*dWb); //Fill [Qu/Nq']*[QaQb/NaNb] with weight (Nq')NaNb
 467       fHistProUQQaQb[iPOI][1]-> Fill(dEta,dUQ/dNq*dQaQb/dNa/dNb,dWq*dWa*dWb); //Fill [Qu/Nq']*[QaQb/NaNb] with weight (Nq')NaNb
 468       fHistSumOfWeightsu[iPOI][0][0]->Fill(dPt ,dWq);        // sum of Nq'
 469       fHistSumOfWeightsu[iPOI][0][1]->Fill(dPt ,pow(dWq,2.));// sum of Nq'^2
 470       fHistSumOfWeightsu[iPOI][0][2]->Fill(dPt ,dWq*dWa*dWb);// sum of Nq'*Na*Nb
 471       fHistSumOfWeightsu[iPOI][1][0]->Fill(dEta,dWq);        // sum of Nq'
 472       fHistSumOfWeightsu[iPOI][1][1]->Fill(dEta,pow(dWq,2.));// sum of Nq'^2
 473       fHistSumOfWeightsu[iPOI][1][2]->Fill(dEta,dNq*dWa*dWb);// sum of Nq'*Na*Nb
 474       //NUA:
 475       fHistProNUAu[iPOI][0][0]->Fill(dPt,dUY,1.); //sin u
 476       fHistProNUAu[iPOI][0][1]->Fill(dPt,dUX,1.); //cos u
 477       fHistProNUAu[iPOI][1][0]->Fill(dEta,dUY,1.); //sin u
 478       fHistProNUAu[iPOI][1][1]->Fill(dEta,dUX,1.); //cos u
 479     }
 480   }//loop over tracks
 481
 482 }
 483
 484 //--------------------------------------------------------------------
 485 void AliFlowAnalysisWithScalarProduct::GetOutputHistograms(TList *outputListHistos){
 486   //get pointers to all output histograms (called before Finish())
 487   fHistList = outputListHistos;
 488
 489   fCommonHists = (AliFlowCommonHist*) fHistList->FindObject("AliFlowCommonHist_SP");
 490   //  fCommonHistsuQ = (AliFlowCommonHist*) fHistList->FindObject("AliFlowCommonHist_uQ");
 491   fCommonHistsRes = (AliFlowCommonHistResults*) fHistList->FindObject("AliFlowCommonHistResults_SP");
 492   fHistProConfig = (TProfile*) fHistList->FindObject("FlowPro_Flags_SP");
 493   if(!fHistProConfig) printf("Error loading fHistProConfig\n");
 494   TList *uQ = (TList*) fHistList->FindObject("uQ");
 495   TList *nua = (TList*) fHistList->FindObject("NUA");
 496   TList *error = (TList*) fHistList->FindObject("error");
 497
 498   fHistProQaQbNorm = (TProfile*) error->FindObject("FlowPro_QaQbNorm_SP");
 499   if(!fHistProQaQbNorm) printf("Error loading fHistProQaQbNorm\n");
 500   fHistProNUAq = (TProfile*) nua->FindObject("FlowPro_NUAq_SP");
 501   if(!fHistProNUAq) printf("Error loading fHistProNUAq\n");
 502   fHistSumOfWeights = (TH1D*) error->FindObject("Flow_SumOfWeights_SP");
 503   if(!fHistSumOfWeights) printf("Error loading fHistSumOfWeights\n");
 504
 505   TString sPOI[2] = {"RP","POI"};
 506   TString sEta[2] = {"Pt","eta"};
 507   TString sWeights[3] = {"uQ","uQuQ","uQQaQb"};
 508   for(Int_t iPOI=0; iPOI!=2; ++iPOI) for(Int_t iSpace=0; iSpace!=2; ++iSpace) {
 509     fHistProUQ[iPOI][iSpace] = (TProfile*) uQ->FindObject( Form( "FlowPro_UQ_%s%s_SP", sEta[iSpace].Data(), sPOI[iPOI].Data() ) );
 510     if(!fHistProUQ[iPOI][iSpace]) printf("Error loading fHistProUQ[%d][%d]\n",iPOI,iSpace);
 511     fHistProNUAu[iPOI][iSpace][0] = (TProfile*) nua->FindObject( Form("FlowPro_NUAu_%s%s_IM_SP", sEta[iSpace].Data(), sPOI[iPOI].Data() ) );
 512     if(!fHistProNUAu[iPOI][iSpace][0]) printf("Error loading fHistProNUAu[%d][%d][0]\n",iPOI,iSpace);
 513     fHistProNUAu[iPOI][iSpace][1] = (TProfile*) nua->FindObject( Form("FlowPro_NUAu_%s%s_RE_SP", sEta[iSpace].Data(), sPOI[iPOI].Data() ) );
 514     if(!fHistProNUAu[iPOI][iSpace][1]) printf("Error loading fHistProNUAu[%d][%d][1]\n",iPOI,iSpace);
 515     fHistProUQQaQb[iPOI][iSpace] = (TProfile*) error->FindObject( Form("FlowPro_UQQaQb_%s%s_SP", sEta[iSpace].Data(), sPOI[iPOI].Data() ) );
 516     for(Int_t i=0; i!=3; ++i){
 517       fHistSumOfWeightsu[iPOI][iSpace][i] = (TH1D*) error->FindObject( Form("Flow_SumOfWeights_%s%s_%s_SP",sWeights[i].Data(),sPOI[iPOI].Data(),sEta[iSpace].Data()) );
 518       if(!fHistSumOfWeightsu[iPOI][iSpace][i]) printf("Error loading fHistSumOfWeightsu[%d][%d][%d]\n",iPOI,iSpace,i);
 519     }
 520   }
 521   if(fHistProConfig) {
 522     fApplyCorrectionForNUA = (Int_t) fHistProConfig->GetBinContent(1);
 523     fNormalizationType  = (Int_t) fHistProConfig->GetBinContent(2);
 524     fUsePhiWeights = (Int_t) fHistProConfig->GetBinContent(3);
 525     fHarmonic = (Int_t) fHistProConfig->GetBinContent(4);
 526   }
 527 }
 528
 529 //--------------------------------------------------------------------
 530 void AliFlowAnalysisWithScalarProduct::Finish() {
 531   //calculate flow and fill the AliFlowCommonHistResults
 532   printf("AliFlowAnalysisWithScalarProduct::Finish()\n");
 533
 534   // access harmonic:
 535   fApplyCorrectionForNUA = (Int_t)(fHistProConfig->GetBinContent(1));
 536   fNormalizationType = (Int_t)(fHistProConfig->GetBinContent(2));
 537   fHarmonic = (Int_t)(fHistProConfig->GetBinContent(4));
 538
 539   printf("*************************************\n");
 540   printf("*************************************\n");
 541   printf("      Integrated flow from           \n");
 542   printf("         Scalar Product              \n\n");
 543   if(!fNormalizationType)
 544     printf("          (BehaveAsEP)               \n\n");
 545
 546   //Calculate reference flow
 547   //----------------------------------
 548   //weighted average over (QaQb/NaNb) with weight (NaNb)
 549   Double_t dEntriesQaQb = fHistProQaQbNorm->GetEntries();
 550   if( dEntriesQaQb < 1 )
 551     return;
 552   Double_t dQaQb  = fHistProQaQbNorm->GetBinContent(1);
 553   if( dQaQb < 0 )
 554     return;
 555   Double_t dSpreadQaQb = fHistProQaQbNorm->GetBinError(1);
 556
 557   //NUA qq:
 558   Double_t dImQa = fHistProNUAq->GetBinContent(1);
 559   Double_t dReQa = fHistProNUAq->GetBinContent(2);
 560   Double_t dImQb = fHistProNUAq->GetBinContent(3);
 561   Double_t dReQb = fHistProNUAq->GetBinContent(4);
 562   if(fApplyCorrectionForNUA)
 563     dQaQb = dQaQb - dImQa*dImQb - dReQa*dReQb;
 564   printf("QaQb = %f +- %f\n", dQaQb, (dSpreadQaQb/TMath::Sqrt(dEntriesQaQb)) );
 565   Double_t dV = TMath::Sqrt(dQaQb);
 566
 567   printf("ResSub = %f\n", dV );
 568   printf("fTotalQvector %d \n",fTotalQvector);
 569   if(!fNormalizationType) {
 570     if(fTotalQvector>2) {
 571       dV = ComputeResolution( TMath::Sqrt2()*FindXi(dV,1e-6) );
 572       printf("An estimate of the event plane resolution is: %f\n", dV );
 573     }
 574   }
 575   printf("ResTot = %f\n", dV );
 576   //statistical error of dQaQb:
 577   //  statistical error = term1 * spread * term2:
 578   //  term1 = sqrt{sum_{i=1}^{N} w^2}/(sum_{i=1}^{N} w)
 579   //  term2 = 1/sqrt(1-term1^2)
 580   Double_t dSumOfLinearWeights = fHistSumOfWeights->GetBinContent(1);
 581   Double_t dSumOfQuadraticWeights = fHistSumOfWeights->GetBinContent(2);
 582   Double_t dTerm1 = 0.;
 583   Double_t dTerm2 = 0.;
 584   if(dSumOfLinearWeights)
 585     dTerm1 = pow(dSumOfQuadraticWeights,0.5)/dSumOfLinearWeights;
 586   if(1.-pow(dTerm1,2.) > 0.)
 587     dTerm2 = 1./pow(1-pow(dTerm1,2.),0.5);
 588   Double_t dStatErrorQaQb = dTerm1 * dSpreadQaQb * dTerm2;
 589   Double_t dVerr = 0.;
 590   if(dQaQb > 0.)
 591     dVerr = (1./(2.*pow(dQaQb,0.5)))*dStatErrorQaQb;
 592   fCommonHistsRes->FillIntegratedFlow(dV,dVerr);
 593   printf("v%d(subevents) = %f +- %f\n",fHarmonic,dV,dVerr);
 594
 595   Int_t iNbins[2];
 596   iNbins[0] = AliFlowCommonConstants::GetMaster()->GetNbinsPt();
 597   iNbins[1] = AliFlowCommonConstants::GetMaster()->GetNbinsEta();
 598
 599   //Calculate differential flow and integrated flow (RP, POI)
 600   //---------------------------------------------------------
 601   //v as a function of eta for RP selection
 602   for(Int_t iRFPorPOI=0; iRFPorPOI != 2; ++iRFPorPOI)
 603   for(Int_t iPTorETA=0; iPTorETA != 2; ++iPTorETA)
 604   for(Int_t b=1; b != iNbins[iPTorETA]+1; ++b) {
 605     Double_t duQpro = fHistProUQ[iRFPorPOI][iPTorETA]->GetBinContent(b);
 606     if(fApplyCorrectionForNUA)
 607       duQpro = duQpro
 608         - fHistProNUAu[iRFPorPOI][iPTorETA][1]->GetBinContent(b)*fHistProNUAq->GetBinContent(6)
 609         - fHistProNUAu[iRFPorPOI][iPTorETA][0]->GetBinContent(b)*fHistProNUAq->GetBinContent(5);
 610     Double_t dv2pro = -999.;
 611     if( TMath::Abs(dV!=0.) ) { dv2pro = duQpro/dV; }
 612     //calculate the statistical error
 613     Double_t dv2ProErr = CalculateStatisticalError(iRFPorPOI, iPTorETA, b, dStatErrorQaQb);
 614     if( (iRFPorPOI==0)&&(iPTorETA==0) )
 615       fCommonHistsRes->FillDifferentialFlowPtRP(  b, dv2pro, dv2ProErr);
 616     if( (iRFPorPOI==0)&&(iPTorETA==1) )
 617       fCommonHistsRes->FillDifferentialFlowEtaRP( b, dv2pro, dv2ProErr);
 618     if( (iRFPorPOI==1)&&(iPTorETA==0) )
 619       fCommonHistsRes->FillDifferentialFlowPtPOI( b, dv2pro, dv2ProErr);
 620     if( (iRFPorPOI==1)&&(iPTorETA==1) )
 621       fCommonHistsRes->FillDifferentialFlowEtaPOI(b, dv2pro, dv2ProErr);
 622     //printf("POI %d | PT %d >>> %f +- %f\n",iRFPorPOI,iPTorETA,dv2pro,dv2ProErr);
 623   }
 624
 625   printf("\n");
 626   printf("*************************************\n");
 627   printf("*************************************\n");
 628 }
 629
 630 //-----------------------------------------------------------------------
 631 void AliFlowAnalysisWithScalarProduct::WriteHistograms(TDirectoryFile *outputFileName) const {
 632  //store the final results in output .root file
 633  outputFileName->Add(fHistList);
 634  outputFileName->Write(outputFileName->GetName(), TObject::kSingleKey);
 635 }
 636
 637 //--------------------------------------------------------------------
 638 Double_t AliFlowAnalysisWithScalarProduct::CalculateStatisticalError(Int_t iRFPorPOI, Int_t iPTorETA, Int_t b, Double_t aStatErrorQaQb) const {
 639   //calculate the statistical error for differential flow for bin b
 640   Double_t duQproSpread = fHistProUQ[iRFPorPOI][iPTorETA]->GetBinError(b);
 641   Double_t sumOfMq = fHistSumOfWeightsu[iRFPorPOI][iPTorETA][0]->GetBinContent(b);
 642   Double_t sumOfMqSquared = fHistSumOfWeightsu[iRFPorPOI][iPTorETA][1]->GetBinContent(b);
 643   Double_t dQaQb = fHistProQaQbNorm->GetBinContent(1);
 644
 645   //non-isotropic terms:
 646   if(fApplyCorrectionForNUA) {
 647     Double_t dImQa = fHistProNUAq->GetBinContent(1);  // <<sin(phi_a)>>
 648     Double_t dReQa = fHistProNUAq->GetBinContent(2);  // <<cos(phi_a)>>
 649     Double_t dImQb = fHistProNUAq->GetBinContent(3);  // <<sin(phi_b)>>
 650     Double_t dReQb = fHistProNUAq->GetBinContent(4);  // <<cos(phi_b)>>
 651     dQaQb = dQaQb - dImQa*dImQb - dReQa*dReQb;
 652   }
 653
 654   Double_t dTerm1 = 0.;
 655   Double_t dTerm2 = 0.;
 656   if(sumOfMq) {
 657     dTerm1 = (pow(sumOfMqSquared,0.5)/sumOfMq);
 658   }
 659   if(1.-pow(dTerm1,2.)>0.) {
 660     dTerm2 = 1./pow(1.-pow(dTerm1,2.),0.5);
 661   }
 662   Double_t duQproErr = dTerm1*duQproSpread*dTerm2;
 663   // covariances:
 664   Double_t dTerm1Cov = fHistSumOfWeightsu[iRFPorPOI][iPTorETA][2]->GetBinContent(b);
 665   Double_t dTerm2Cov = fHistSumOfWeights->GetBinContent(1);
 666   Double_t dTerm3Cov = sumOfMq;
 667   Double_t dWeightedCovariance = 0.;
 668   if(dTerm2Cov*dTerm3Cov>0.) {
 669     Double_t dDenominator = 1.-dTerm1Cov/(dTerm2Cov*dTerm3Cov);
 670     Double_t dPrefactor = dTerm1Cov/(dTerm2Cov*dTerm3Cov);
 671     if(dDenominator!=0) {
 672       Double_t dCovariance = ( fHistProUQQaQb[iRFPorPOI][iPTorETA]->GetBinContent(b)-dQaQb*fHistProUQ[iRFPorPOI][iPTorETA]->GetBinContent(b))/dDenominator;
 673       dWeightedCovariance = dCovariance*dPrefactor;
 674     }
 675   }
 676   Double_t dv2ProErr = 0.; // final statitical error:
 677   if(dQaQb>0.) {
 678     Double_t dv2ProErrorSquared = (1./4.)*pow(dQaQb,-3.)*
 679       (pow(fHistProUQ[iRFPorPOI][iPTorETA]->GetBinContent(b),2.)*pow(aStatErrorQaQb,2.)
 680        + 4.*pow(dQaQb,2.)*pow(duQproErr,2.)
 681        - 4.*dQaQb*fHistProUQ[iRFPorPOI][iPTorETA]->GetBinContent(b)*dWeightedCovariance);
 682     if(dv2ProErrorSquared>0.) dv2ProErr = pow(dv2ProErrorSquared,0.5);
 683   }
 684   return dv2ProErr;
 685 }
 686
 687 Double_t AliFlowAnalysisWithScalarProduct::ComputeResolution( Double_t x ) const {
 688   // Computes resolution for Event Plane method
 689   if(x > 51.3) {
 690     printf("Warning: Estimation of total resolution might be WRONG. Please check!");
 691     return 0.99981;
 692   }
 693   Double_t a = x*x/4;
 694   Double_t b = TMath::Exp(-a)*TMath::BesselI0(a)+TMath::Exp(-a)*TMath::BesselI1(a);
 695   return TMath::Sqrt(TMath::PiOver2())/2*x*b;
 696 }
 697
 698 Double_t AliFlowAnalysisWithScalarProduct::FindXi( Double_t res, Double_t prec ) const {
 699   // Computes x(res) for Event Plane method
 700   if(res > 0.99981) {
 701     printf("Warning: Resolution for subEvent is high. You reached the precision limit.");
 702     return 51.3;
 703   }
 704   int nSteps =0;
 705   Double_t xtmp=0, xmin=0, xmax=51.3, rtmp=0, delta=2*prec;
 706   while( delta > prec ) {
 707     xtmp = 0.5*(xmin+xmax);
 708     rtmp = ComputeResolution(xtmp);
 709     delta = TMath::Abs( res-rtmp );
 710     if(rtmp>res) xmax = xtmp;
 711     if(rtmp<res) xmin = xtmp;
 712     nSteps++;
 713   }
 714   return xtmp;
 715 }
 716