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    fEtaMin(0.01),
  49    fEtaMax(1.),
  50    fDebug(kFALSE),
  51    fHistList(NULL),
  52    fHistProUQetaRP(NULL),
  53    fHistProUQetaPOI(NULL),
  54    fHistProUQPtRP(NULL),
  55    fHistProUQPtPOI(NULL),
  56    fHistProQaQb(NULL),
  57    fHistProM(NULL),
  58    fCommonHists(NULL),
  59    fCommonHistsRes(NULL)
  60 {
  61   // Constructor.
  62   fHistList = new TList();
  63 }
  64  //-----------------------------------------------------------------------
  65
  66
  67  AliFlowAnalysisWithScalarProduct::~AliFlowAnalysisWithScalarProduct()
  68  {
  69    //destructor
  70    delete fHistList;
  71  }
  72
  73
  74 //-----------------------------------------------------------------------
  75
  76 void AliFlowAnalysisWithScalarProduct::WriteHistograms(TString* outputFileName)
  77 {
  78  //store the final results in output .root file
  79
  80   TFile *output = new TFile(outputFileName->Data(),"RECREATE");
  81   //output->WriteObject(fHistList, "cobjSP","SingleKey");
  82   fHistList->SetName("cobjSP");
  83   fHistList->SetOwner(kTRUE);
  84   fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
  85   delete output;
  86 }
  87
  88 //-----------------------------------------------------------------------
  89
  90 void AliFlowAnalysisWithScalarProduct::WriteHistograms(TString outputFileName)
  91 {
  92  //store the final results in output .root file
  93
  94   TFile *output = new TFile(outputFileName.Data(),"RECREATE");
  95   //output->WriteObject(fHistList, "cobjSP","SingleKey");
  96   fHistList->SetName("cobjSP");
  97   fHistList->SetOwner(kTRUE);
  98   fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
  99   delete output;
 100 }
 101
 102 //-----------------------------------------------------------------------
 103 void AliFlowAnalysisWithScalarProduct::Init() {
 104
 105   //Define all histograms
 106   cout<<"---Analysis with the Scalar Product Method--- Init"<<endl;
 107
 108   Int_t iNbinsPt   = AliFlowCommonConstants::GetNbinsPt();
 109   Double_t dPtMin  = AliFlowCommonConstants::GetPtMin();
 110   Double_t dPtMax  = AliFlowCommonConstants::GetPtMax();
 111   Int_t iNbinsEta  = AliFlowCommonConstants::GetNbinsEta();
 112   Double_t dEtaMin = AliFlowCommonConstants::GetEtaMin();
 113   Double_t dEtaMax = AliFlowCommonConstants::GetEtaMax();
 114
 115   fHistProUQetaRP = new TProfile("Flow_UQetaRP_SP","Flow_UQetaRP_SP",iNbinsEta,dEtaMin,dEtaMax);
 116   fHistProUQetaRP->SetXTitle("{eta}");
 117   fHistProUQetaRP->SetYTitle("<uQ>");
 118   fHistList->Add(fHistProUQetaRP);
 119
 120   fHistProUQetaPOI = new TProfile("Flow_UQetaPOI_SP","Flow_UQetaPOI_SP",iNbinsEta,dEtaMin,dEtaMax);
 121   fHistProUQetaPOI->SetXTitle("{eta}");
 122   fHistProUQetaPOI->SetYTitle("<uQ>");
 123   fHistList->Add(fHistProUQetaPOI);
 124
 125   fHistProUQPtRP = new TProfile("Flow_UQPtRP_SP","Flow_UQPtRP_SP",iNbinsPt,dPtMin,dPtMax);
 126   fHistProUQPtRP->SetXTitle("p_t (GeV)");
 127   fHistProUQPtRP->SetYTitle("<uQ>");
 128   fHistList->Add(fHistProUQPtRP);
 129
 130   fHistProUQPtPOI = new TProfile("Flow_UQPtPOI_SP","Flow_UQPtPOI_SP",iNbinsPt,dPtMin,dPtMax);
 131   fHistProUQPtPOI->SetXTitle("p_t (GeV)");
 132   fHistProUQPtPOI->SetYTitle("<uQ>");
 133   fHistList->Add(fHistProUQPtPOI);
 134
 135   fHistProQaQb = new TProfile("Flow_QaQb_SP","Flow_QaQb_SP", 1, -0.5, 0.5);
 136   fHistProQaQb->SetYTitle("<QaQb>");
 137   fHistList->Add(fHistProQaQb);
 138
 139   fHistProM = new TProfile("Flow_M_SP","Flow_M_SP",2,0.5, 2.5);
 140   fHistProM -> SetYTitle("<*>");
 141   fHistProM -> SetXTitle("<M-1>, <Ma*Mb>");
 142   fHistList->Add(fHistProM);
 143
 144   fCommonHists = new AliFlowCommonHist("AliFlowCommonHistSP");
 145   fHistList->Add(fCommonHists);
 146   fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsSP");
 147   fHistList->Add(fCommonHistsRes);
 148
 149   fEventNumber = 0;  //set number of events to zero
 150 }
 151
 152 //-----------------------------------------------------------------------
 153
 154 void AliFlowAnalysisWithScalarProduct::Make(AliFlowEventSimple* anEvent) {
 155
 156   //Fill histogram
 157   if (anEvent) {
 158
 159     //fill control histograms
 160     fCommonHists->FillControlHistograms(anEvent);
 161
 162     //get Q vectors for the eta-subevents
 163     AliFlowVector* vQarray = new AliFlowVector[2];
 164     anEvent->GetQsub(vQarray);
 165     AliFlowVector vQa = vQarray[0];
 166     AliFlowVector vQb = vQarray[1];
 167     //get total Q vector
 168     AliFlowVector vQ = vQa + vQb;
 169
 170     //fill the multiplicity histograms for the prefactor
 171     fHistProM -> Fill(1,vQ.GetMult()-1);                //<M-1>
 172     fHistProM -> Fill(2,vQa.GetMult()*vQb.GetMult());   //<Ma*Mb>
 173     //scalar product of the two subevents
 174     Double_t dQaQb = vQa*vQb;
 175     fHistProQaQb -> Fill(0.,dQaQb);
 176
 177     //loop over the tracks of the event
 178     AliFlowTrackSimple*   pTrack = NULL;
 179     Int_t iNumberOfTracks = anEvent->NumberOfTracks();
 180     for (Int_t i=0;i<iNumberOfTracks;i++)
 181       {
 182         pTrack = anEvent->GetTrack(i) ;
 183         if (pTrack){
 184           Double_t dPhi = pTrack->Phi();
 185           //calculate vU
 186           TVector2 vU;
 187           Double_t dUX = TMath::Cos(2*dPhi);
 188           Double_t dUY = TMath::Sin(2*dPhi);
 189           vU.Set(dUX,dUY);
 190           Double_t dModulus = vU.Mod();
 191           if (dModulus!=0.) vU.Set(dUX/dModulus,dUY/dModulus);  // make length 1
 192           else cerr<<"dModulus is zero!"<<endl;
 193
 194           TVector2 vQm = vQ;
 195           //subtract particle from the flowvector if used to define it
 196           if (pTrack->InRPSelection()) {
 197             if (pTrack->InSubevent(0) || pTrack->InSubevent(1)) {
 198               Double_t dQmX = vQm.X() - dUX;
 199               Double_t dQmY = vQm.Y() - dUY;
 200               vQm.Set(dQmX,dQmY);
 201             }
 202           }
 203
 204           //dUQ = scalar product of vU and vQm
 205           Double_t dUQ = vU * vQm;
 206           Double_t dPt = pTrack->Pt();
 207           Double_t dEta = pTrack->Eta();
 208           //fill the profile histograms
 209           if (pTrack->InRPSelection()) {
 210             fHistProUQetaRP -> Fill(dEta,dUQ);
 211             fHistProUQPtRP -> Fill(dPt,dUQ);
 212           }
 213           if (pTrack->InPOISelection()) {
 214             fHistProUQetaPOI -> Fill(dEta,dUQ);
 215             fHistProUQPtPOI -> Fill(dPt,dUQ);
 216           }
 217         }//track selected
 218       }//loop over tracks
 219
 220     fEventNumber++;
 221     //    cout<<"@@@@@ "<<fEventNumber<<" events processed"<<endl;
 222     delete [] vQarray;
 223   }
 224 }
 225
 226   //--------------------------------------------------------------------
 227 void AliFlowAnalysisWithScalarProduct::GetOutputHistograms(TList *outputListHistos){
 228
 229   //get pointers to all output histograms (called before Finish())
 230   if (outputListHistos) {
 231   //Get the common histograms from the output list
 232     AliFlowCommonHist *pCommonHist = dynamic_cast<AliFlowCommonHist*>
 233       (outputListHistos->FindObject("AliFlowCommonHistSP"));
 234     AliFlowCommonHistResults *pCommonHistResults = dynamic_cast<AliFlowCommonHistResults*>
 235       (outputListHistos->FindObject("AliFlowCommonHistResultsSP"));
 236     TProfile* pHistProQaQb     = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_QaQb_SP"));
 237     TProfile* pHistProM        = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_M_SP"));
 238     TProfile* pHistProUQetaRP  = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQetaRP_SP"));
 239     TProfile* pHistProUQetaPOI = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQetaPOI_SP"));
 240     TProfile* pHistProUQPtRP   = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQPtRP_SP"));
 241     TProfile* pHistProUQPtPOI  = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQPtPOI_SP"));
 242     if (pCommonHist && pCommonHistResults && pHistProQaQb && pHistProM &&
 243         pHistProUQetaRP && pHistProUQetaPOI && pHistProUQPtRP && pHistProUQPtPOI) {
 244       this -> SetCommonHists(pCommonHist);
 245       this -> SetCommonHistsRes(pCommonHistResults);
 246       this -> SetHistProQaQb(pHistProQaQb);
 247       this -> SetHistProM(pHistProM);
 248       this -> SetHistProUQetaRP(pHistProUQetaRP);
 249       this -> SetHistProUQetaPOI(pHistProUQetaPOI);
 250       this -> SetHistProUQPtRP(pHistProUQPtRP);
 251       this -> SetHistProUQPtPOI(pHistProUQPtPOI);
 252       }
 253    }
 254 }
 255
 256 //--------------------------------------------------------------------
 257 void AliFlowAnalysisWithScalarProduct::Finish() {
 258
 259   //calculate flow and fill the AliFlowCommonHistResults
 260   if (fDebug) cout<<"AliFlowAnalysisWithScalarProduct::Finish()"<<endl;
 261
 262   cout<<"*************************************"<<endl;
 263   cout<<"*************************************"<<endl;
 264   cout<<"      Integrated flow from           "<<endl;
 265   cout<<"         Scalar product              "<<endl;
 266   cout<<endl;
 267
 268   Int_t iNbinsPt  = AliFlowCommonConstants::GetNbinsPt();
 269   Int_t iNbinsEta = AliFlowCommonConstants::GetNbinsEta();
 270
 271   Double_t dMmin1    = fHistProM->GetBinContent(1);  //average over M-1
 272   Double_t dMmin1Err = fHistProM->GetBinError(1);    //error on average over M-1
 273   Double_t dMaMb     = fHistProM->GetBinContent(2);  //average over Ma*Mb
 274   Double_t dMaMbErr  = fHistProM->GetBinError(2);    //error on average over Ma*Mb
 275
 276   Double_t dMcorrection = 0.;     //correction factor for Ma != Mb
 277   Double_t dMcorrectionErr = 0.;
 278   Double_t dMcorrectionErrRel = 0.;
 279   Double_t dMcorrectionErrRel2 = 0.;
 280
 281   if (dMaMb != 0. && dMmin1 != 0.) {
 282     dMcorrection    = dMmin1/(TMath::Sqrt(dMaMb));
 283     dMcorrectionErr = dMcorrection*(dMmin1Err/dMmin1 + dMaMbErr/(2*dMaMb));
 284     dMcorrectionErrRel = dMcorrectionErr/dMcorrection;
 285     dMcorrectionErrRel2 = dMcorrectionErrRel*dMcorrectionErrRel;
 286   }
 287
 288   Double_t dQaQbAv  = TMath::Abs(fHistProQaQb->GetBinContent(1)); //average over events //TEST TAKE ABS
 289   Double_t dQaQbErr = fHistProQaQb->GetBinError(1);
 290   Double_t dQaQbErrRel = 0.;
 291   if (dQaQbAv != 0.) {
 292     dQaQbErrRel = dQaQbErr/dQaQbAv; }
 293   Double_t dQaQbErrRel2 = dQaQbErrRel*dQaQbErrRel;
 294
 295   if (dQaQbAv <= 0.){
 296     //set v to -0
 297     fCommonHistsRes->FillIntegratedFlowRP(-0.,0.);
 298     fCommonHistsRes->FillIntegratedFlow(-0.,0.);
 299     cout<<"dV(RP) = -0. +- 0."<<endl;
 300     fCommonHistsRes->FillIntegratedFlowPOI(-0.,0.);
 301     cout<<"dV(POI) = -0. +- 0."<<endl;
 302   } else {
 303   Double_t dQaQbSqrt = TMath::Sqrt(dQaQbAv);  //DOES NOT WORK IF dQaQbAv IS NEGATIVE
 304     if (dMaMb>0.) { dQaQbSqrt *= dMcorrection; }
 305     else { dQaQbSqrt = 0.; }
 306     Double_t dQaQbSqrtErrRel2 = dMcorrectionErrRel2 + (1/4)*dQaQbErrRel2;
 307
 308     //v as a function of eta for RP selection
 309     for(Int_t b=0;b<iNbinsEta;b++) {
 310       Double_t duQpro = fHistProUQetaRP->GetBinContent(b);
 311       Double_t duQerr = fHistProUQetaRP->GetBinError(b); //copy error for now
 312       Double_t duQerrRel = 0.;
 313       if (duQpro != 0.) {duQerrRel = duQerr/duQpro;}
 314       Double_t duQerrRel2 = duQerrRel*duQerrRel;
 315
 316       Double_t dv2pro     = 0.;
 317       if (dQaQbSqrt!=0.) { dv2pro = duQpro/dQaQbSqrt; }
 318       Double_t dv2errRel2 = duQerrRel2 + dQaQbSqrtErrRel2;
 319       Double_t dv2errRel  = 0.;
 320       if (dv2errRel2>0.) { dv2errRel  = TMath::Sqrt(dv2errRel2); }
 321       Double_t dv2err     = dv2pro*dv2errRel;
 322       //fill TH1D
 323       fCommonHistsRes->FillDifferentialFlowEtaRP(b, dv2pro, dv2err);
 324     } //loop over bins b
 325
 326     //v as a function of eta for POI selection
 327     for(Int_t b=0;b<iNbinsEta;b++) {
 328       Double_t duQpro = fHistProUQetaPOI->GetBinContent(b);
 329       Double_t duQerr = fHistProUQetaPOI->GetBinError(b); //copy error for now
 330       Double_t duQerrRel = 0.;
 331       if (duQpro != 0.) {duQerrRel = duQerr/duQpro;}
 332       Double_t duQerrRel2 = duQerrRel*duQerrRel;
 333
 334       Double_t dv2pro     = 0.;
 335       if (dQaQbSqrt!=0.) { dv2pro = duQpro/dQaQbSqrt; }
 336       Double_t dv2errRel2 = duQerrRel2 + dQaQbSqrtErrRel2;
 337       Double_t dv2errRel  = 0.;
 338       if (dv2errRel2>0.) { dv2errRel  = TMath::Sqrt(dv2errRel2); }
 339       Double_t dv2err     = dv2pro*dv2errRel;
 340       //fill TH1D
 341       fCommonHistsRes->FillDifferentialFlowEtaPOI(b, dv2pro, dv2err);
 342     } //loop over bins b
 343
 344     //v as a function of Pt for RP selection
 345     TH1F* fHistPtRP = fCommonHists->GetHistPtRP(); //for calculating integrated flow
 346     Double_t dVRP = 0.;
 347     Double_t dSum = 0.;
 348     Double_t dErrV =0.;
 349
 350     for(Int_t b=0;b<iNbinsPt;b++) {
 351       Double_t duQpro = fHistProUQPtRP->GetBinContent(b);
 352       Double_t duQerr = fHistProUQPtRP->GetBinError(b); //copy error for now
 353       Double_t duQerrRel = 0.;
 354       if (duQpro != 0.) {duQerrRel = duQerr/duQpro;}
 355       Double_t duQerrRel2 = duQerrRel*duQerrRel;
 356
 357       Double_t dv2pro     = 0.;
 358       if (dQaQbSqrt!=0.) { dv2pro = duQpro/dQaQbSqrt; }
 359       Double_t dv2errRel2 = duQerrRel2 + dQaQbSqrtErrRel2;
 360       Double_t dv2errRel  = 0.;
 361       if (dv2errRel2>0.) { dv2errRel  = TMath::Sqrt(dv2errRel2); }
 362       Double_t dv2err     = dv2pro*dv2errRel;
 363       //fill TH1D
 364       fCommonHistsRes->FillDifferentialFlowPtRP(b, dv2pro, dv2err);
 365       //calculate integrated flow for RP selection
 366       if (fHistPtRP){
 367         Double_t dYieldPt = fHistPtRP->GetBinContent(b);
 368         dVRP += dv2pro*dYieldPt;
 369         dSum +=dYieldPt;
 370         dErrV += dYieldPt*dYieldPt*dv2err*dv2err;
 371       } else { cout<<"fHistPtRP is NULL"<<endl; }
 372     } //loop over bins b
 373
 374     if (dSum != 0.) {
 375       dVRP /= dSum; //the pt distribution should be normalised
 376       dErrV /= (dSum*dSum);
 377       dErrV = TMath::Sqrt(dErrV);
 378     }
 379     fCommonHistsRes->FillIntegratedFlowRP(dVRP,dErrV);
 380     fCommonHistsRes->FillIntegratedFlow(dVRP,dErrV);
 381
 382     cout<<"dV(RP) = "<<dVRP<<" +- "<<dErrV<<endl;
 383
 384     //v as a function of Pt for POI selection
 385     TH1F* fHistPtPOI = fCommonHists->GetHistPtPOI(); //for calculating integrated flow
 386     Double_t dVPOI = 0.;
 387     dSum = 0.;
 388     dErrV =0.;
 389
 390     for(Int_t b=0;b<iNbinsPt;b++) {
 391       Double_t duQpro = fHistProUQPtPOI->GetBinContent(b);
 392       Double_t duQerr = fHistProUQPtPOI->GetBinError(b); //copy error for now
 393       Double_t duQerrRel = 0.;
 394       if (duQpro != 0.) {duQerrRel = duQerr/duQpro;}
 395       Double_t duQerrRel2 = duQerrRel*duQerrRel;
 396
 397       Double_t dv2pro     = 0.;
 398       if (dQaQbSqrt!=0.) { dv2pro = duQpro/dQaQbSqrt; }
 399       Double_t dv2errRel2 = duQerrRel2 + dQaQbSqrtErrRel2;
 400       Double_t dv2errRel  = 0.;
 401       if (dv2errRel2>0.) { dv2errRel  = TMath::Sqrt(dv2errRel2); }
 402       Double_t dv2err     = dv2pro*dv2errRel;
 403       //fill TH1D
 404       fCommonHistsRes->FillDifferentialFlowPtPOI(b, dv2pro, dv2err);
 405
 406       //calculate integrated flow for POI selection
 407       if (fHistPtPOI){
 408         Double_t dYieldPt = fHistPtPOI->GetBinContent(b);
 409         dVPOI += dv2pro*dYieldPt;
 410         dSum +=dYieldPt;
 411         dErrV += dYieldPt*dYieldPt*dv2err*dv2err;
 412       } else { cout<<"fHistPtPOI is NULL"<<endl; }
 413     } //loop over bins b
 414
 415     if (dSum != 0.) {
 416       dVPOI /= dSum; //the pt distribution should be normalised
 417       dErrV /= (dSum*dSum);
 418       dErrV = TMath::Sqrt(dErrV);
 419     }
 420     fCommonHistsRes->FillIntegratedFlowPOI(dVPOI,dErrV);
 421
 422     cout<<"dV(POI) = "<<dVPOI<<" +- "<<dErrV<<endl;
 423     }
 424   cout<<endl;
 425   cout<<"*************************************"<<endl;
 426   cout<<"*************************************"<<endl;
 427
 428   //cout<<".....finished"<<endl;
 429  }
 430
 431