PWG2/FLOW/AliFlowCommon/AliFlowAnalysisWithScalarProduct.cxx

   1 /*************************************************************************
   2 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15
  16 #define AliFlowAnalysisWithScalarProduct_cxx
  17
  18 #include "Riostream.h"  //needed as include
  19 #include "TFile.h"      //needed as include
  20 #include "TList.h"
  21 #include "TMath.h"
  22 #include "TProfile.h"
  23 #include "TVector2.h"
  24
  25 class TH1F;
  26
  27 #include "AliFlowCommonConstants.h"    //needed as include
  28 #include "AliFlowEventSimple.h"
  29 #include "AliFlowTrackSimple.h"
  30 #include "AliFlowCommonHist.h"
  31 #include "AliFlowCommonHistResults.h"
  32 #include "AliFlowAnalysisWithScalarProduct.h"
  33
  34 class AliFlowVector;
  35
  36 // AliFlowAnalysisWithScalarProduct:
  37 // Description:
  38 // Maker to analyze Flow with the Scalar product method.
  39 //
  40 // author: N. van der Kolk (kolk@nikhef.nl)
  41
  42 ClassImp(AliFlowAnalysisWithScalarProduct)
  43
  44   //-----------------------------------------------------------------------
  45
  46  AliFlowAnalysisWithScalarProduct::AliFlowAnalysisWithScalarProduct():
  47    fEventNumber(0),
  48    fDebug(kFALSE),
  49    fHistList(NULL),
  50    fHistProUQetaRP(NULL),
  51    fHistProUQetaPOI(NULL),
  52    fHistProUQPtRP(NULL),
  53    fHistProUQPtPOI(NULL),
  54    fHistProQaQb(NULL),
  55    fHistProM(NULL),
  56    fCommonHists(NULL),
  57    fCommonHistsRes(NULL)
  58 {
  59   // Constructor.
  60   fHistList = new TList();
  61 }
  62  //-----------------------------------------------------------------------
  63
  64
  65  AliFlowAnalysisWithScalarProduct::~AliFlowAnalysisWithScalarProduct()
  66  {
  67    //destructor
  68    delete fHistList;
  69  }
  70
  71
  72 //-----------------------------------------------------------------------
  73
  74 void AliFlowAnalysisWithScalarProduct::WriteHistograms(TString* outputFileName)
  75 {
  76  //store the final results in output .root file
  77
  78   TFile *output = new TFile(outputFileName->Data(),"RECREATE");
  79   //output->WriteObject(fHistList, "cobjSP","SingleKey");
  80   fHistList->SetName("cobjSP");
  81   fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
  82   delete output;
  83 }
  84
  85 //-----------------------------------------------------------------------
  86
  87 void AliFlowAnalysisWithScalarProduct::WriteHistograms(TString outputFileName)
  88 {
  89  //store the final results in output .root file
  90
  91   TFile *output = new TFile(outputFileName.Data(),"RECREATE");
  92   //output->WriteObject(fHistList, "cobjSP","SingleKey");
  93   fHistList->SetName("cobjSP");
  94   fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
  95   delete output;
  96 }
  97
  98 //-----------------------------------------------------------------------
  99 void AliFlowAnalysisWithScalarProduct::Init() {
 100
 101   //Define all histograms
 102   cout<<"---Analysis with the Scalar Product Method--- Init"<<endl;
 103
 104   Int_t iNbinsPt   = AliFlowCommonConstants::GetNbinsPt();
 105   Double_t dPtMin  = AliFlowCommonConstants::GetPtMin();
 106   Double_t dPtMax  = AliFlowCommonConstants::GetPtMax();
 107   Int_t iNbinsEta  = AliFlowCommonConstants::GetNbinsEta();
 108   Double_t dEtaMin = AliFlowCommonConstants::GetEtaMin();
 109   Double_t dEtaMax = AliFlowCommonConstants::GetEtaMax();
 110
 111   fHistProUQetaRP = new TProfile("Flow_UQetaRP_SP","Flow_UQetaRP_SP",iNbinsEta,dEtaMin,dEtaMax);
 112   fHistProUQetaRP->SetXTitle("{eta}");
 113   fHistProUQetaRP->SetYTitle("<uQ>");
 114   fHistList->Add(fHistProUQetaRP);
 115
 116   fHistProUQetaPOI = new TProfile("Flow_UQetaPOI_SP","Flow_UQetaPOI_SP",iNbinsEta,dEtaMin,dEtaMax);
 117   fHistProUQetaPOI->SetXTitle("{eta}");
 118   fHistProUQetaPOI->SetYTitle("<uQ>");
 119   fHistList->Add(fHistProUQetaPOI);
 120
 121   fHistProUQPtRP = new TProfile("Flow_UQPtRP_SP","Flow_UQPtRP_SP",iNbinsPt,dPtMin,dPtMax);
 122   fHistProUQPtRP->SetXTitle("p_t (GeV)");
 123   fHistProUQPtRP->SetYTitle("<uQ>");
 124   fHistList->Add(fHistProUQPtRP);
 125
 126   fHistProUQPtPOI = new TProfile("Flow_UQPtPOI_SP","Flow_UQPtPOI_SP",iNbinsPt,dPtMin,dPtMax);
 127   fHistProUQPtPOI->SetXTitle("p_t (GeV)");
 128   fHistProUQPtPOI->SetYTitle("<uQ>");
 129   fHistList->Add(fHistProUQPtPOI);
 130
 131   fHistProQaQb = new TProfile("Flow_QaQb_SP","Flow_QaQb_SP", 1, -0.5, 0.5);
 132   fHistProQaQb->SetYTitle("<QaQb>");
 133   fHistList->Add(fHistProQaQb);
 134
 135   fHistProM = new TProfile("Flow_M_SP","Flow_M_SP",2,0.5, 2.5);
 136   fHistProM -> SetYTitle("<*>");
 137   fHistProM -> SetXTitle("<M-1>, <Ma*Mb>");
 138   fHistList->Add(fHistProM);
 139
 140   fCommonHists = new AliFlowCommonHist("AliFlowCommonHistSP");
 141   fHistList->Add(fCommonHists);
 142   fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsSP");
 143   fHistList->Add(fCommonHistsRes);
 144
 145   fEventNumber = 0;  //set number of events to zero
 146 }
 147
 148 //-----------------------------------------------------------------------
 149
 150 void AliFlowAnalysisWithScalarProduct::Make(AliFlowEventSimple* anEvent) {
 151
 152   //Fill histogram
 153   if (anEvent) {
 154
 155     //fill control histograms
 156     fCommonHists->FillControlHistograms(anEvent);
 157
 158     //get the Q vector from the FlowEvent
 159     AliFlowVector vQ = anEvent->GetQ();
 160     fHistProM -> Fill(1,vQ.GetMult()-1);
 161     //get Q vectors for the subevents
 162     AliFlowVector vQa = anEvent->GetQsub(-0.9,-0.1);
 163     AliFlowVector vQb = anEvent->GetQsub(0.1,0.9);
 164     fHistProM -> Fill(2,vQa.GetMult()*vQb.GetMult());
 165     Double_t dQaQb = vQa*vQb;
 166     fHistProQaQb -> Fill(0.,dQaQb);
 167
 168     //loop over the tracks of the event
 169     AliFlowTrackSimple*   pTrack = NULL;
 170     Int_t iNumberOfTracks = anEvent->NumberOfTracks();
 171     for (Int_t i=0;i<iNumberOfTracks;i++)
 172       {
 173         pTrack = anEvent->GetTrack(i) ;
 174         if (pTrack){
 175           Double_t dPhi = pTrack->Phi();
 176           //calculate vU
 177           TVector2 vU;
 178           Double_t dUX = TMath::Cos(2*dPhi);
 179           Double_t dUY = TMath::Sin(2*dPhi);
 180           vU.Set(dUX,dUY);
 181           Double_t dModulus = vU.Mod();
 182           if (dModulus!=0.) vU.Set(dUX/dModulus,dUY/dModulus);  // make length 1
 183           else cerr<<"dModulus is zero!"<<endl;
 184
 185           TVector2 vQm = vQ;
 186           //subtrackt particle from the flowvector if used to define it
 187           if (pTrack->InRPSelection()) {
 188             Double_t dQmX = vQm.X() - dUX;
 189             Double_t dQmY = vQm.Y() - dUY;
 190             vQm.Set(dQmX,dQmY);
 191           }
 192
 193           //dUQ = scalar product of vU and vQm
 194           Double_t dUQ = vU * vQm;
 195           Double_t dPt = pTrack->Pt();
 196           Double_t dEta = pTrack->Eta();
 197           //fill the profile histograms
 198           if (pTrack->InRPSelection()) {
 199             fHistProUQetaRP -> Fill(dEta,dUQ);
 200             fHistProUQPtRP -> Fill(dPt,dUQ);
 201           }
 202           if (pTrack->InPOISelection()) {
 203             fHistProUQetaPOI -> Fill(dEta,dUQ);
 204             fHistProUQPtPOI -> Fill(dPt,dUQ);
 205           }
 206         }//track selected
 207       }//loop over tracks
 208
 209     fEventNumber++;
 210     //    cout<<"@@@@@ "<<fEventNumber<<" events processed"<<endl;
 211   }
 212 }
 213
 214   //--------------------------------------------------------------------
 215 void AliFlowAnalysisWithScalarProduct::GetOutputHistograms(TList *outputListHistos){
 216
 217   //get pointers to all output histograms (called before Finish())
 218   if (outputListHistos) {
 219   //Get the common histograms from the output list
 220     AliFlowCommonHist *pCommonHist = dynamic_cast<AliFlowCommonHist*>
 221       (outputListHistos->FindObject("AliFlowCommonHistSP"));
 222     AliFlowCommonHistResults *pCommonHistResults = dynamic_cast<AliFlowCommonHistResults*>
 223       (outputListHistos->FindObject("AliFlowCommonHistResultsSP"));
 224     TProfile* pHistProQaQb     = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_QaQb_SP"));
 225     TProfile* pHistProM        = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_M_SP"));
 226     TProfile* pHistProUQetaRP  = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQetaRP_SP"));
 227     TProfile* pHistProUQetaPOI = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQetaPOI_SP"));
 228     TProfile* pHistProUQPtRP   = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQPtRP_SP"));
 229     TProfile* pHistProUQPtPOI  = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQPtPOI_SP"));
 230     if (pCommonHist && pCommonHistResults && pHistProQaQb && pHistProM &&
 231         pHistProUQetaRP && pHistProUQetaPOI && pHistProUQPtRP && pHistProUQPtPOI) {
 232       this -> SetCommonHists(pCommonHist);
 233       this -> SetCommonHistsRes(pCommonHistResults);
 234       this -> SetHistProQaQb(pHistProQaQb);
 235       this -> SetHistProM(pHistProM);
 236       this -> SetHistProUQetaRP(pHistProUQetaRP);
 237       this -> SetHistProUQetaPOI(pHistProUQetaPOI);
 238       this -> SetHistProUQPtRP(pHistProUQPtRP);
 239       this -> SetHistProUQPtPOI(pHistProUQPtPOI);
 240       }
 241    }
 242 }
 243
 244 //--------------------------------------------------------------------
 245 void AliFlowAnalysisWithScalarProduct::Finish() {
 246
 247   //calculate flow and fill the AliFlowCommonHistResults
 248   if (fDebug) cout<<"AliFlowAnalysisWithScalarProduct::Finish()"<<endl;
 249
 250   cout<<"*************************************"<<endl;
 251   cout<<"*************************************"<<endl;
 252   cout<<"      Integrated flow from           "<<endl;
 253   cout<<"         Scalar product              "<<endl;
 254   cout<<endl;
 255
 256   Int_t iNbinsPt  = AliFlowCommonConstants::GetNbinsPt();
 257   Int_t iNbinsEta = AliFlowCommonConstants::GetNbinsEta();
 258
 259   Double_t dMmin1    = fHistProM->GetBinContent(1);  //average over M-1
 260   Double_t dMmin1Err = fHistProM->GetBinError(1);    //error on average over M-1
 261   Double_t dMaMb    = fHistProM->GetBinContent(2);   //average over Ma*Mb
 262   Double_t dMaMbErr = fHistProM->GetBinError(2);     //error on average over Ma*Mb
 263
 264   Double_t dMcorrection = 0.;     //correction factor for Ma != Mb
 265   Double_t dMcorrectionErr = 0.;
 266   Double_t dMcorrectionErrRel = 0.;
 267   Double_t dMcorrectionErrRel2 = 0.;
 268
 269   if (dMaMb != 0. && dMmin1 != 0.) {
 270     dMcorrection    = dMmin1/(TMath::Sqrt(dMaMb));
 271     dMcorrectionErr = dMcorrection*(dMmin1Err/dMmin1 + dMaMbErr/(2*dMaMb));
 272     dMcorrectionErrRel = dMcorrectionErr/dMcorrection;
 273     dMcorrectionErrRel2 = dMcorrectionErrRel*dMcorrectionErrRel;
 274   }
 275
 276   Double_t dQaQbAv  = fHistProQaQb->GetBinContent(1); //average over events
 277   Double_t dQaQbErr = fHistProQaQb->GetBinError(1);
 278   Double_t dQaQbErrRel = 0.;
 279   if (dQaQbAv != 0.) {
 280     dQaQbErrRel = dQaQbErr/dQaQbAv; }
 281   Double_t dQaQbErrRel2 = dQaQbErrRel*dQaQbErrRel;
 282
 283   if (dQaQbAv <= 0.){
 284     //set v to -0
 285     fCommonHistsRes->FillIntegratedFlowRP(-0.,0.);
 286     fCommonHistsRes->FillIntegratedFlow(-0.,0.);
 287     cout<<"dV(RP) = -0. +- 0."<<endl;
 288     fCommonHistsRes->FillIntegratedFlowPOI(-0.,0.);
 289     cout<<"dV(POI) = -0. +- 0."<<endl;
 290   } else {
 291     Double_t dQaQbSqrt = TMath::Sqrt(dQaQbAv);
 292     if (dMaMb>0.) { dQaQbSqrt *= dMcorrection; }
 293     else { dQaQbSqrt = 0.; }
 294     Double_t dQaQbSqrtErrRel2 = dMcorrectionErrRel2 + (1/4)*dQaQbErrRel2;
 295
 296     //v as a function of eta for RP selection
 297     for(Int_t b=0;b<iNbinsEta;b++) {
 298       Double_t duQpro = fHistProUQetaRP->GetBinContent(b);
 299       Double_t duQerr = fHistProUQetaRP->GetBinError(b); //copy error for now
 300       Double_t duQerrRel = 0.;
 301       if (duQpro != 0.) {duQerrRel = duQerr/duQpro;}
 302       Double_t duQerrRel2 = duQerrRel*duQerrRel;
 303
 304       Double_t dv2pro     = 0.;
 305       if (dQaQbSqrt!=0.) { dv2pro = duQpro/dQaQbSqrt; }
 306       Double_t dv2errRel2 = duQerrRel2 + dQaQbSqrtErrRel2;
 307       Double_t dv2errRel  = 0.;
 308       if (dv2errRel2>0.) { dv2errRel  = TMath::Sqrt(dv2errRel2); }
 309       Double_t dv2err     = dv2pro*dv2errRel;
 310       //fill TH1D
 311       fCommonHistsRes->FillDifferentialFlowEtaRP(b, dv2pro, dv2err);
 312     } //loop over bins b
 313
 314     //v as a function of eta for POI selection
 315     for(Int_t b=0;b<iNbinsEta;b++) {
 316       Double_t duQpro = fHistProUQetaPOI->GetBinContent(b);
 317       Double_t duQerr = fHistProUQetaPOI->GetBinError(b); //copy error for now
 318       Double_t duQerrRel = 0.;
 319       if (duQpro != 0.) {duQerrRel = duQerr/duQpro;}
 320       Double_t duQerrRel2 = duQerrRel*duQerrRel;
 321
 322       Double_t dv2pro     = 0.;
 323       if (dQaQbSqrt!=0.) { dv2pro = duQpro/dQaQbSqrt; }
 324       Double_t dv2errRel2 = duQerrRel2 + dQaQbSqrtErrRel2;
 325       Double_t dv2errRel  = 0.;
 326       if (dv2errRel2>0.) { dv2errRel  = TMath::Sqrt(dv2errRel2); }
 327       Double_t dv2err     = dv2pro*dv2errRel;
 328       //fill TH1D
 329       fCommonHistsRes->FillDifferentialFlowEtaPOI(b, dv2pro, dv2err);
 330     } //loop over bins b
 331
 332     //v as a function of Pt for RP selection
 333     TH1F* fHistPtRP = fCommonHists->GetHistPtRP(); //for calculating integrated flow
 334     Double_t dVRP = 0.;
 335     Double_t dSum = 0.;
 336     Double_t dErrV =0.;
 337
 338     for(Int_t b=0;b<iNbinsPt;b++) {
 339       Double_t duQpro = fHistProUQPtRP->GetBinContent(b);
 340       Double_t duQerr = fHistProUQPtRP->GetBinError(b); //copy error for now
 341       Double_t duQerrRel = 0.;
 342       if (duQpro != 0.) {duQerrRel = duQerr/duQpro;}
 343       Double_t duQerrRel2 = duQerrRel*duQerrRel;
 344
 345       Double_t dv2pro     = 0.;
 346       if (dQaQbSqrt!=0.) { dv2pro = duQpro/dQaQbSqrt; }
 347       Double_t dv2errRel2 = duQerrRel2 + dQaQbSqrtErrRel2;
 348       Double_t dv2errRel  = 0.;
 349       if (dv2errRel2>0.) { dv2errRel  = TMath::Sqrt(dv2errRel2); }
 350       Double_t dv2err     = dv2pro*dv2errRel;
 351       //fill TH1D
 352       fCommonHistsRes->FillDifferentialFlowPtRP(b, dv2pro, dv2err);
 353       //calculate integrated flow for RP selection
 354       if (fHistPtRP){
 355         Double_t dYieldPt = fHistPtRP->GetBinContent(b);
 356         dVRP += dv2pro*dYieldPt;
 357         dSum +=dYieldPt;
 358         dErrV += dYieldPt*dYieldPt*dv2err*dv2err;
 359       } else { cout<<"fHistPtRP is NULL"<<endl; }
 360     } //loop over bins b
 361
 362     if (dSum != 0.) {
 363       dVRP /= dSum; //the pt distribution should be normalised
 364       dErrV /= (dSum*dSum);
 365       dErrV = TMath::Sqrt(dErrV);
 366     }
 367     fCommonHistsRes->FillIntegratedFlowRP(dVRP,dErrV);
 368     fCommonHistsRes->FillIntegratedFlow(dVRP,dErrV);
 369
 370     cout<<"dV(RP) = "<<dVRP<<" +- "<<dErrV<<endl;
 371
 372     //v as a function of Pt for POI selection
 373     TH1F* fHistPtPOI = fCommonHists->GetHistPtPOI(); //for calculating integrated flow
 374     Double_t dVPOI = 0.;
 375     dSum = 0.;
 376     dErrV =0.;
 377
 378     for(Int_t b=0;b<iNbinsPt;b++) {
 379       Double_t duQpro = fHistProUQPtPOI->GetBinContent(b);
 380       Double_t duQerr = fHistProUQPtPOI->GetBinError(b); //copy error for now
 381       Double_t duQerrRel = 0.;
 382       if (duQpro != 0.) {duQerrRel = duQerr/duQpro;}
 383       Double_t duQerrRel2 = duQerrRel*duQerrRel;
 384
 385       Double_t dv2pro     = 0.;
 386       if (dQaQbSqrt!=0.) { dv2pro = duQpro/dQaQbSqrt; }
 387       Double_t dv2errRel2 = duQerrRel2 + dQaQbSqrtErrRel2;
 388       Double_t dv2errRel  = 0.;
 389       if (dv2errRel2>0.) { dv2errRel  = TMath::Sqrt(dv2errRel2); }
 390       Double_t dv2err     = dv2pro*dv2errRel;
 391       //fill TH1D
 392       fCommonHistsRes->FillDifferentialFlowPtPOI(b, dv2pro, dv2err);
 393
 394       //calculate integrated flow for POI selection
 395       if (fHistPtPOI){
 396         Double_t dYieldPt = fHistPtPOI->GetBinContent(b);
 397         dVPOI += dv2pro*dYieldPt;
 398         dSum +=dYieldPt;
 399         dErrV += dYieldPt*dYieldPt*dv2err*dv2err;
 400       } else { cout<<"fHistPtPOI is NULL"<<endl; }
 401     } //loop over bins b
 402
 403     if (dSum != 0.) {
 404       dVPOI /= dSum; //the pt distribution should be normalised
 405       dErrV /= (dSum*dSum);
 406       dErrV = TMath::Sqrt(dErrV);
 407     }
 408     fCommonHistsRes->FillIntegratedFlowPOI(dVPOI,dErrV);
 409
 410     cout<<"dV(POI) = "<<dVPOI<<" +- "<<dErrV<<endl;
 411   }
 412   cout<<endl;
 413   cout<<"*************************************"<<endl;
 414   cout<<"*************************************"<<endl;
 415
 416   //cout<<".....finished"<<endl;
 417  }
 418
 419