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