PWG2/FORWARD/analysis2/AliForwardFlowUtil.cxx

   1 //
   2 // Class used to handle the input from AODs and put it into histograms
   3 // the Forward Flow tasks can run on.
   4 // It can also add flow to AliAODMCParticles.
   5 //
   6 #include <iostream>
   7 #include "AliForwardFlowUtil.h"
   8 #include "AliAODCentralMult.h"
   9 #include "AliAODMCParticle.h"
  10 #include "AliAODMCHeader.h"
  11 #include "AliLog.h"
  12 #include "AliAODForwardMult.h"
  13 #include "AliAODEvent.h"
  14
  15 ClassImp(AliForwardFlowUtil)
  16
  17 //_____________________________________________________________________
  18 AliForwardFlowUtil::AliForwardFlowUtil() :
  19   fList(0),
  20   fZvertex(0)
  21 {
  22   //
  23   // Default Constructor
  24   //
  25 }
  26 //_____________________________________________________________________
  27 AliForwardFlowUtil::AliForwardFlowUtil(TList* outputList) :
  28   fList(0),
  29   fZvertex(0)
  30 {
  31   //
  32   // Constructor
  33   //
  34   // Parameters:
  35   // TList: list containing histograms for flow analysis
  36   //
  37   fList = outputList;
  38 }
  39 //_____________________________________________________________________
  40 Bool_t AliForwardFlowUtil::AODCheck(const AliAODForwardMult* aodfm) const
  41 {
  42   //
  43   // Function to check that and AOD event meets the cuts
  44   //
  45   // Parameters:
  46   //  AliAODForwardMult: forward mult object with trigger and vertex info
  47   //
  48   if (!aodfm->IsTriggerBits(AliAODForwardMult::kInel)) return kFALSE;
  49   if (!aodfm->HasIpZ()) return kFALSE;
  50   if (!aodfm->InRange(-fZvertex,fZvertex)) return kFALSE;
  51
  52   return kTRUE;
  53 }
  54 //_____________________________________________________________________
  55 Bool_t AliForwardFlowUtil::LoopAODFMD(const AliAODEvent* aodevent) const
  56 {
  57   //
  58   // Loop over AliAODFowardMult object and fill histograms provided
  59   // by flow task.
  60   //
  61   AliAODForwardMult* aodfmult = static_cast<AliAODForwardMult*>(aodevent->FindListObject("Forward"));
  62   if (!aodfmult) return kFALSE;
  63   if (!AODCheck(aodfmult)) return kFALSE;
  64
  65   TH2D hdNdetadphi = aodfmult->GetHistogram();
  66   TH1D* dNdphi = (TH1D*)fList->FindObject("hdNdphiSE");
  67   TH2D* dNdetadphi = (TH2D*)fList->FindObject("hdNdetadphiSE");
  68   dNdphi->Reset();
  69   dNdetadphi->Reset();
  70
  71   Double_t mult = 0;
  72   Double_t eta = 0;
  73   Double_t phi = 0;
  74   Double_t weight = 0;
  75   for (Int_t etaBin = 1; etaBin<=hdNdetadphi.GetNbinsX(); etaBin++) {
  76     eta = hdNdetadphi.GetXaxis()->GetBinCenter(etaBin);
  77     for (Int_t phiBin = 1; phiBin<=hdNdetadphi.GetNbinsY(); phiBin++) {
  78       phi = hdNdetadphi.GetYaxis()->GetBinCenter(phiBin);
  79       weight = hdNdetadphi.GetBinContent(etaBin, phiBin);
  80       dNdetadphi->Fill(eta, phi, weight);
  81       if (eta < 4.) dNdphi->Fill(phi, weight);
  82       mult += weight;
  83     }
  84   }
  85   dNdphi->SetBinContent(0, mult);
  86 //  if (aodfmult->HasCentrality()) dNdphi->SetBinContent(dNdphi->GetNbinsX()+1, aodfmult->GetCentrality());
  87 //  else dNdphi->SetBinContent(dNdphi->GetNbinsX()+1, -1);
  88
  89   dNdphi->SetBinContent(dNdphi->GetNbinsX()+1, GetCentFromMC(aodevent));
  90   return kTRUE;
  91 }
  92 //_____________________________________________________________________
  93 Bool_t AliForwardFlowUtil::LoopAODSPD(const AliAODEvent* aodevent) const
  94 {
  95   //
  96   // Loop over AliAODCentralMult object and fill histograms
  97   // provided by flow task
  98   //
  99   AliAODCentralMult* aodcmult = static_cast<AliAODCentralMult*>(aodevent->FindListObject("CentralClusters"));
 100   if (!aodcmult) return kFALSE;
 101
 102   TH2D hdNdetadphi = aodcmult->GetHistogram();
 103   TH2D* dNdetadphi = (TH2D*)fList->FindObject("hdNdetadphiSESPD");
 104   dNdetadphi->Reset();
 105
 106   Double_t eta = 0;
 107   Double_t phi = 0;
 108   Double_t weight = 0;
 109   for (Int_t etaBin = 1; etaBin<=hdNdetadphi.GetNbinsX(); etaBin++) {
 110     eta = hdNdetadphi.GetXaxis()->GetBinCenter(etaBin);
 111     for (Int_t phiBin = 1; phiBin<=hdNdetadphi.GetNbinsY(); phiBin++) {
 112       phi = hdNdetadphi.GetYaxis()->GetBinCenter(phiBin);
 113       weight = hdNdetadphi.GetBinContent(etaBin, phiBin);
 114       dNdetadphi->Fill(eta, phi, weight);
 115     }
 116   }
 117   return kTRUE;
 118 }
 119 //_____________________________________________________________________
 120 Bool_t AliForwardFlowUtil::LoopAODtrrefHits(const AliAODEvent* aodevent) const
 121 {
 122   //
 123   // Loop over AliAODForwardMult object, get MC track ref information
 124   // and fill flow histograms
 125   //
 126
 127   AliAODForwardMult* aodfmult = static_cast<AliAODForwardMult*>(aodevent->FindListObject("ForwardMC"));
 128   if (!aodfmult) return kFALSE;
 129  // if (!AODCheck(aodfmult)) return kFALSE;
 130
 131   TH2D hdNdetadphi = aodfmult->GetHistogram();
 132   TH1D* dNdphi = (TH1D*)fList->FindObject("hdNdphiSETrRef");
 133   TH2D* dNdetadphi = (TH2D*)fList->FindObject("hdNdetadphiSETrRef");
 134   dNdphi->Reset();
 135   dNdetadphi->Reset();
 136
 137   Double_t mult = 0;
 138   Double_t eta = 0;
 139   Double_t phi = 0;
 140   Double_t weight = 0;
 141   for(Int_t etaBin = 1; etaBin<=hdNdetadphi.GetNbinsX(); etaBin++) {
 142     eta = hdNdetadphi.GetXaxis()->GetBinCenter(etaBin);
 143     for(Int_t phiBin = 1; phiBin<=hdNdetadphi.GetNbinsY(); phiBin++) {
 144       phi = hdNdetadphi.GetYaxis()->GetBinCenter(phiBin);
 145       weight = hdNdetadphi.GetBinContent(etaBin, phiBin);
 146       dNdetadphi->Fill(eta, phi, weight);
 147       dNdphi->Fill(phi, weight);
 148       mult += weight;
 149     }
 150   }
 151   dNdphi->SetBinContent(0, mult);
 152
 153   return kTRUE;
 154 }
 155 //_____________________________________________________________________
 156 Bool_t AliForwardFlowUtil::LoopAODmc(const AliAODEvent* aodevent,
 157                                      TString addFlow = "",
 158                                      Int_t type = 0,
 159                                      Int_t order = 2) const
 160 {
 161   //
 162   // Loop over AliAODParticle branch and fill flow histograms
 163   //
 164
 165   //retreive MC particles from event
 166   TClonesArray* mcArray = (TClonesArray*)aodevent->FindListObject(AliAODMCParticle::StdBranchName());
 167   if(!mcArray){
 168     AliWarning("No MC array found in AOD. Try making it again.");
 169     return kFALSE;
 170   }
 171
 172   Double_t rp = 0;
 173   if (addFlow.Length() > 1) {
 174     AliAODMCHeader* header = (AliAODMCHeader*)aodevent->FindListObject(AliAODMCHeader::StdBranchName());
 175     rp = header->GetReactionPlaneAngle();
 176   }
 177
 178   TH2D* dNdetadphi = (TH2D*)fList->FindObject("hdNdetadphiSEMC");
 179   TH1D* dNdphi = (TH1D*)fList->FindObject("hdNdphiSEMC");
 180   dNdphi->Reset();
 181   dNdetadphi->Reset();
 182
 183   Int_t ntracks = mcArray->GetEntriesFast();
 184
 185   // Track loop: chek how many particles will be accepted
 186   Double_t mult = 0;
 187   Double_t weight = 0;
 188   for (Int_t it = 0; it < ntracks; it++) {
 189     weight = 0;
 190     AliAODMCParticle* particle = (AliAODMCParticle*)mcArray->At(it);
 191     if (!particle) {
 192       AliError(Form("Could not receive track %d", it));
 193       continue;
 194     }
 195     if (!particle->IsPhysicalPrimary()) continue;
 196     if (particle->Charge() == 0) continue;
 197     if (particle->Eta() > -3.4 && particle->Eta() < 5) {
 198       // Add flow if it is in the argument
 199       if (addFlow.Length() > 1) {
 200         if (addFlow.Contains("pt"))
 201           weight += AddptFlow(particle->Pt(), type);
 202         if (addFlow.Contains("pid"))
 203           weight += AddpidFlow(particle->PdgCode(), type);
 204         if (addFlow.Contains("eta"))
 205           weight += AddetaFlow(particle->Eta(), type);
 206         if (addFlow.Contains("flat"))
 207           weight = 0.1*type;
 208         weight *= 2.*TMath::Cos((Double_t)order*(particle->Phi()-rp));
 209       }
 210       weight += 1;
 211
 212       dNdphi->Fill(particle->Phi(), weight);
 213       dNdetadphi->Fill(particle->Eta(), particle->Phi(), weight);
 214       mult += weight;
 215     }
 216   }
 217
 218   dNdphi->SetBinContent(0, mult);
 219
 220   return kTRUE;
 221 }
 222 //_____________________________________________________________________
 223 Double_t AliForwardFlowUtil::AddptFlow(Double_t Pt = 0, Int_t type = 0) const
 224 {
 225   //
 226   // Add pt dependent flow factor
 227   //
 228   Double_t weight = 0;
 229
 230   if (type == 1) weight = 0.125*Pt;
 231
 232   if (type == 2) {
 233     weight = 0.2*(Pt/2.0);
 234     if (Pt > 2.0) weight = 0.2;
 235   }
 236
 237   if (type == 3) {
 238       weight = 0.05*(Pt/2.0);
 239       if (Pt < 2.0) weight = 0.05;
 240   }
 241
 242   if (type == 4) {
 243       weight = 0.2*(Pt/1.0);
 244       if (Pt > 1.0) weight = 0.2;
 245   }
 246
 247   if (type == 5) {
 248       weight = 0.05*(Pt/1.0);
 249       if (Pt < 1.0) weight = 0.05;
 250   }
 251
 252   if (type == 6) {
 253       weight = 0.2*(Pt/3.0);
 254       if (Pt > 3.0) weight = 0.2;
 255   }
 256
 257   return weight;
 258 }
 259 //_____________________________________________________________________
 260 Double_t AliForwardFlowUtil::AddpidFlow(Int_t ID = 0, Int_t type = 0) const
 261 {
 262   //
 263   // Add pid dependent flow factor
 264   //
 265   Double_t weight = 0;
 266
 267   if (type == 1) {
 268     weight = 0.07;
 269     if (TMath::Abs(ID) ==  211) // pion flow
 270       weight = 0.1;
 271     if (TMath::Abs(ID) ==  2212) // proton flow
 272       weight = 0.05;
 273   }
 274
 275   if (type == 2) {
 276     weight = 0.06;
 277     if (TMath::Abs(ID) ==  211) // pion flow
 278       weight = 0.1;
 279     if (TMath::Abs(ID) ==  2212) // proton flow
 280       weight = 0.08;
 281   }
 282
 283   if (type == 3) {
 284     weight = 0.05;
 285     if (TMath::Abs(ID) ==  211) // pion flow
 286       weight = 0.1;
 287     if (TMath::Abs(ID) ==  2212) // proton flow
 288       weight = 0.07;
 289   }
 290
 291   if (type == 4) {
 292     weight = 0.07;
 293     if (TMath::Abs(ID) ==  211) // pion flow
 294       weight = 0.1;
 295     if (TMath::Abs(ID) ==  2212) // proton flow
 296       weight = 0.085;
 297   }
 298
 299   return weight;
 300 }
 301 //_____________________________________________________________________
 302 Double_t AliForwardFlowUtil::AddetaFlow(Double_t eta = 0, Int_t type = 0) const
 303 {
 304   //
 305   // Add eta dependent flow factor
 306   //
 307   Double_t weight = 0;
 308
 309   if (type == 1) weight = 0.03 + 0.07 * (1 - TMath::Abs(eta) / 6.);
 310
 311   if (type == 2) weight = 0.07 * (1 - TMath::Abs(eta) / 6.);
 312
 313   if (type == 3) weight = 0.07 * (1 - TMath::Abs(eta) / 8.);
 314
 315   if (type == 4) weight = 0.07 * (1 - TMath::Abs(eta) / 10.);
 316
 317   if (type == 5) weight = 0.07 * (1 - TMath::Abs(eta) / 12.);
 318
 319   if (type == 6) weight = 0.07 * (1 - TMath::Abs(eta) / 14.);
 320
 321   if (type == 7) weight = 0.07 * (1 - TMath::Abs(eta) / 16.);
 322
 323   return weight;
 324 }
 325 //_____________________________________________________________________
 326 Double_t AliForwardFlowUtil::GetCentFromMC(const AliAODEvent* aodevent) const
 327 {
 328   //
 329   // Get centrality from MC impact parameter.
 330   // Values taken from: https://twiki.cern.ch/twiki/bin/viewauth/ALICE/CentStudies
 331   //
 332   Double_t cent = -1.;
 333   Double_t b = -1.;
 334   AliAODMCHeader* header = (AliAODMCHeader*)aodevent->FindListObject(AliAODMCHeader::StdBranchName());
 335   if (!header) return cent;
 336   b = header->GetImpactParameter();
 337
 338   if ( 0.00 <= b && b <  3.50)
 339     cent = 2.5;
 340   if ( 3.50 <= b && b <  4.95)
 341     cent = 7.5;
 342   if ( 4.95 <= b && b <  6.98)
 343     cent = 15.;
 344   if ( 6.98 <= b && b <  8.55)
 345     cent = 25.;
 346   if ( 8.55 <= b && b <  9.88)
 347     cent = 35.;
 348   if ( 9.88 <= b && b < 11.04)
 349     cent = 45.;
 350   if (11.04 <= b && b < 12.09)
 351     cent = 55.;
 352   if (12.09 <= b && b < 13.06)
 353     cent = 65.;
 354   if (13.06 <= b && b < 13.97)
 355     cent = 75.;
 356   if (13.97 <= b && b < 19.61)
 357     cent = 90.;
 358
 359   return cent;
 360 }
 361 //_____________________________________________________________________
 362 //
 363 //
 364 // EOF
 365