PWG/FLOW/Base/AliFlowLYZEventPlane.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 AliFlowAnalysisWithLYZEventPlane_cxx
  17
  18 // AliFlowLYZEventPlane:
  19 //
  20 // Class to calculate the event plane and event weight from the LYZ method
  21 // It needs input from the standard LYZ first and second run
  22 // author: N. van der Kolk (kolk@nikhef.nl)
  23
  24 #include "Riostream.h"
  25 #include "TProfile.h"
  26 #include "TFile.h"
  27 #include "TComplex.h"
  28 #include "TList.h"
  29
  30 #include "AliFlowVector.h"
  31 #include "AliFlowLYZConstants.h"
  32 #include "AliFlowEventSimple.h"
  33 #include "AliFlowLYZEventPlane.h"
  34
  35 class AliFlowTrackSimple;
  36
  37 using std::cout;
  38 using std::endl;
  39 ClassImp(AliFlowLYZEventPlane)
  40
  41   //-----------------------------------------------------------------------
  42
  43 AliFlowLYZEventPlane::AliFlowLYZEventPlane():
  44   fSecondRunList(0),
  45   fWR(0),
  46   fPsi(0),
  47   fSecondReDtheta(0),
  48   fSecondImDtheta(0),
  49   fFirstr0theta(0)
  50 {
  51   // Constructor.
  52
  53 }
  54 //-----------------------------------------------------------------------
  55
  56
  57 AliFlowLYZEventPlane::~AliFlowLYZEventPlane()
  58 {
  59   //destructor
  60   delete fSecondRunList;
  61 }
  62
  63
  64 //-----------------------------------------------------------------------
  65 void AliFlowLYZEventPlane::Init()
  66 {
  67   //Declare histograms & get input files
  68   cout<<"---Lee Yang Zeros Event Plane Method---"<<endl;
  69
  70   //input histograms
  71   if (fSecondRunList) {
  72     fSecondReDtheta = (TProfile*)fSecondRunList->FindObject("Second_FlowPro_ReDtheta_LYZSUM");
  73     fSecondImDtheta = (TProfile*)fSecondRunList->FindObject("Second_FlowPro_ImDtheta_LYZSUM");
  74     fFirstr0theta   = (TProfile*)fSecondRunList->FindObject("First_FlowPro_r0theta_LYZSUM");
  75
  76     //warnings
  77     if (!fSecondReDtheta) {cout<<"fSecondReDtheta is NULL!"<<endl; }
  78     if (!fSecondImDtheta) {cout<<"fSecondImDtheta is NULL!"<<endl; }
  79     if (!fFirstr0theta)   {cout<<"fFirstr0theta is NULL!"<<endl; }
  80   }
  81
  82
  83 }
  84
  85 //-----------------------------------------------------------------------
  86 void AliFlowLYZEventPlane::CalculateRPandW(AliFlowVector aQ)
  87 {
  88   //declare variables
  89   Int_t iNtheta = AliFlowLYZConstants::GetMaster()->GetNtheta();
  90   Double_t dCosTerm = 0;
  91   Double_t dSinTerm = 0;
  92   TComplex cDtheta;
  93   TComplex cRatio;
  94
  95   if (aQ.Mod()==0.) { cout<<"Q vector is NULL"<<endl; }
  96   else {
  97     for (Int_t theta=0;theta<iNtheta;theta++)   {
  98       Double_t dTheta = ((float)theta/iNtheta)*TMath::Pi()/2;
  99       //Calculate Qtheta
 100       Double_t dQtheta = aQ.X()*cos(2*dTheta)+aQ.Y()*sin(2*dTheta);  //get Qtheta from Q vector
 101
 102       //get R0
 103       Double_t dR0 = fFirstr0theta->GetBinContent(theta+1);
 104       //cerr<<"dR0 = "<<dR0<<endl;
 105
 106       //get Dtheta
 107       Double_t dReDtheta = fSecondReDtheta->GetBinContent(theta+1);
 108       //cerr<<"dReDtheta = "<<dReDtheta<<endl;
 109       Double_t dImDtheta = fSecondImDtheta->GetBinContent(theta+1);
 110       //cerr<<"dImDtheta = "<<dImDtheta<<endl;
 111       cDtheta(dReDtheta,dImDtheta);
 112
 113       TComplex cExpo(0.,dR0*dQtheta);                  //Complex number: 0 +(i r0 Qtheta)
 114       if (cDtheta.Rho()!=0.) { cRatio =(TComplex::Exp(cExpo))/cDtheta; } //(e^(i r0 Qtheta))/Dtheta
 115       else { cRatio(0.,0.); }
 116
 117       //sum over theta
 118       dCosTerm += cRatio.Re() * TMath::Cos(2*dTheta);    //Re{(e^(i r0 Qtheta))/Dtheta } cos(2 theta)
 119       dSinTerm += cRatio.Re() * TMath::Sin(2*dTheta);    //Re{(e^(i r0 Qtheta))/Dtheta } sin(2 theta)
 120
 121     }//loop over theta
 122
 123     //average over theta
 124     dCosTerm /= iNtheta;
 125     dSinTerm /= iNtheta;
 126     //cerr<<"cosTerm & sinTerm are: "<<dCosTerm<<" & "<<dSinTerm<<endl;
 127
 128     //calculate fWR
 129     fWR = TMath::Sqrt(dCosTerm*dCosTerm + dSinTerm*dSinTerm);
 130
 131     //calculate fRP
 132     fPsi = 0.5*TMath::ATan2(dSinTerm,dCosTerm);   //takes care of the signs correctly!
 133     if (fPsi < 0.) { fPsi += TMath::Pi(); }       //to shift distribution from (-pi/2 to pi/2) to (0 to pi)
 134   }
 135
 136 }
 137