ANALYSIS/AliFlowAnalysis.cxx

   1 #include "AliFlowAnalysis.h"
   2 //________________________________
   3 ///////////////////////////////////////////////////////////
   4 //
   5 // class AliFlowAnalysis
   6 //
   7 // Flow Analysis
   8 //
   9 //
  10 // S.Radomski@gsi.de
  11 // Piotr.Skowronski@cern.ch
  12 //
  13 ///////////////////////////////////////////////////////////
  14 /*********************************************************/
  15
  16 #include <TString.h>
  17 #include <TParticle.h>
  18
  19 #include <AliStack.h>
  20 #include <AliAOD.h>
  21 #include <AliVAODParticle.h>
  22 #include <AliAODParticleCut.h>
  23
  24 #include <AliESDtrack.h>
  25 #include <AliESD.h>
  26
  27 ClassImp(AliFlowAnalysis)
  28
  29 AliFlowAnalysis::AliFlowAnalysis():
  30  fPartCut(0x0)
  31 {
  32  //ctor
  33 }
  34 /*********************************************************/
  35
  36 AliFlowAnalysis::~AliFlowAnalysis()
  37 {
  38  //dtor
  39   delete fPartCut;
  40 }
  41 /*********************************************************/
  42
  43 Int_t AliFlowAnalysis::Init()
  44 {
  45   //Initilizes anaysis
  46   Info("Init","");
  47   return 0;
  48 }
  49 /*********************************************************/
  50
  51 Int_t AliFlowAnalysis::ProcessEvent(AliAOD* aodrec, AliAOD* aodsim)
  52 {
  53
  54   Info("ProcessEvent","Sim AOD address %#x",aodsim);
  55   Double_t psi = 0, v2 = 0;
  56   if (aodrec)
  57    {
  58      GetFlow(aodrec,v2,psi);
  59      Info("ProcessEvent","Reconstructed Event: Event plane is %f, V2 is %f",psi,v2);
  60    }
  61
  62   if (aodsim)
  63    {
  64      GetFlow(aodsim,v2,psi);
  65      Info("ProcessEvent","Simulated Event: Event plane is %f, V2 is %f",psi,v2);
  66    }
  67
  68   return 0;
  69
  70 }
  71 /*********************************************************/
  72
  73 Int_t AliFlowAnalysis::Finish()
  74 {
  75   //Finish analysis and writes results
  76   Info("Init","Finish");
  77   return 0;
  78 }
  79 /*********************************************************/
  80
  81 Double_t AliFlowAnalysis::GetEventPlane(AliAOD* aod)
  82 {
  83   //returns event plane in degrees
  84   if (aod == 0x0)
  85    {
  86      Error("AliFlowAnalysis::GetFlow","Pointer to AOD is NULL");
  87      return -1.0;
  88    }
  89
  90   Double_t psi;
  91   Int_t mult = aod->GetNumberOfParticles();
  92
  93   Double_t ssin = 0, scos = 0;
  94
  95   for (Int_t i=0; i<mult; i++)
  96    {
  97      AliVAODParticle* aodtrack = aod->GetParticle(i);
  98      if (aodtrack == 0x0)
  99       {
 100         Error("AliFlowAnalysis::GetEventPlane","Can not get track %d", i);
 101         continue;
 102       }
 103
 104      if (fPartCut)
 105       if (fPartCut->Rejected(aodtrack))
 106         continue;
 107
 108      Double_t phi = TMath::Pi()+TMath::ATan2(-aodtrack->Py(),-aodtrack->Px());
 109
 110      ssin += TMath::Sin( 2.0 * phi );
 111      scos += TMath::Cos( 2.0 * phi );
 112    }
 113
 114   psi = atan2 (ssin, scos) / 2.0;
 115   psi = psi * 180. / TMath::Pi();
 116
 117   return psi;
 118
 119 }
 120 /*********************************************************/
 121
 122 void AliFlowAnalysis::GetFlow(AliAOD* aod,Double_t& v2,Double_t& psi)
 123 {
 124 //returns flow parameters: v2 and event plane
 125   if (aod == 0x0)
 126    {
 127      Error("AliFlowAnalysis::GetFlow","Pointer to AOD is NULL");
 128      return;
 129    }
 130
 131   psi = GetEventPlane(aod);
 132   Int_t mult = aod->GetNumberOfParticles();
 133
 134   Double_t ssin = 0, scos = 0;
 135
 136   for (Int_t i=0; i<mult; i++)
 137    {
 138      AliVAODParticle* aodtrack = aod->GetParticle(i);
 139      if (aodtrack == 0x0)
 140       {
 141         Error("AliFlowAnalysis::GetEventPlane","Can not get track %d", i);
 142         continue;
 143       }
 144      if (fPartCut)
 145       if (fPartCut->Rejected(aodtrack))
 146         continue;
 147
 148      Double_t phi = TMath::Pi()+TMath::ATan2(-aodtrack->Py(),-aodtrack->Px());
 149      ssin += TMath::Sin( 2.0 * (phi - psi));
 150      scos += TMath::Cos( 2.0 * (phi - psi));
 151    }
 152
 153   v2 = TMath::Hypot(ssin,scos);
 154 }
 155
 156
 157 /*********************************************************/
 158
 159 Double_t AliFlowAnalysis::GetEventPlane(AliESD* esd)
 160 {
 161   //returns event plane
 162   if (esd == 0x0)
 163    {
 164      ::Error("AliFlowAnalysis::GetFlow","Pointer to ESD is NULL");
 165      return -1.0;
 166    }
 167
 168   Double_t psi;
 169   Int_t mult = esd->GetNumberOfTracks();
 170
 171   Double_t ssin = 0, scos = 0;
 172
 173   for (Int_t i=0; i<mult; i++)
 174    {
 175      AliESDtrack* esdtrack = esd->GetTrack(i);
 176      if (esdtrack == 0x0)
 177       {
 178         ::Error("AliFlowAnalysis::GetEventPlane","Can not get track %d", i);
 179         continue;
 180       }
 181
 182      Double_t mom[3];//momentum
 183      esdtrack->GetPxPyPz(mom);
 184      Double_t phi = TMath::Pi()+TMath::ATan2(-mom[1],-mom[0]);
 185
 186      ssin += TMath::Sin( 2.0 * phi );
 187      scos += TMath::Cos( 2.0 * phi );
 188    }
 189
 190   psi = atan2 (ssin, scos) / 2.0;
 191   psi = psi * 180. / TMath::Pi();
 192
 193   return psi;
 194
 195 }
 196 /*********************************************************/
 197
 198 void AliFlowAnalysis::GetFlow(AliESD* esd,Double_t& v2,Double_t& psi)
 199 {
 200 //returns flow parameters: v2 and event plane
 201   if (esd == 0x0)
 202    {
 203      ::Error("AliFlowAnalysis::GetFlow","Pointer to ESD is NULL");
 204      return;
 205    }
 206
 207   psi = GetEventPlane(esd);
 208   Int_t mult = esd->GetNumberOfTracks();
 209
 210   Double_t ssin = 0, scos = 0;
 211
 212   for (Int_t i=0; i<mult; i++)
 213    {
 214      AliESDtrack* esdtrack = esd->GetTrack(i);
 215      if (esdtrack == 0x0)
 216       {
 217         ::Error("AliFlowAnalysis::GetEventPlane","Can not get track %d", i);
 218         continue;
 219       }
 220
 221      Double_t mom[3];//momentum
 222      esdtrack->GetPxPyPz(mom);
 223      Double_t phi = TMath::Pi()+TMath::ATan2(-mom[1],-mom[0]);
 224
 225      ssin += TMath::Sin( 2.0 * (phi - psi));
 226      scos += TMath::Cos( 2.0 * (phi - psi));
 227    }
 228
 229   v2 = TMath::Hypot(ssin,scos);
 230 }
 231 /*********************************************************/