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 #include "TGraph.h"
  16 #include "TF1.h"
  17
  18 ClassImp(AliForwardFlowUtil)
  19
  20 //_____________________________________________________________________
  21 AliForwardFlowUtil::AliForwardFlowUtil() :
  22   fList(0),
  23   fCent(-1),
  24   fVertex(1111),
  25   fAliceCent4th(),
  26   fAlicePt2nd4050(),
  27   fAlicePt4th3040(),
  28   fAlicePt4th4050(),
  29   fImpactParToCent()
  30    {}
  31   //
  32   // Default Constructor
  33   //
  34 //_____________________________________________________________________
  35 AliForwardFlowUtil::AliForwardFlowUtil(TList* outputList) :
  36   fList(0),
  37   fCent(-1),
  38   fVertex(1111),
  39   fAliceCent4th(),
  40   fAlicePt2nd4050(),
  41   fAlicePt4th3040(),
  42   fAlicePt4th4050(),
  43   fImpactParToCent()
  44 {
  45   //
  46   // Constructor
  47   //
  48   // Parameters:
  49   // TList: list containing histograms for flow analysis
  50   //
  51   fList = outputList;
  52
  53   Double_t xCumulant4thTPCrefMultTPConlyAll[] = {2.5,7.5,15,25,35,45,55,65};
  54   Double_t yCumulant4thTPCrefMultTPConlyAll[] = {0.017855,0.032440,0.055818,0.073137,0.083898,0.086690,0.082040,0.077777};
  55   Int_t nPointsCumulant4thTPCrefMultTPConlyAll = sizeof(xCumulant4thTPCrefMultTPConlyAll)/sizeof(Double_t);
  56   fAliceCent4th = new TGraph(nPointsCumulant4thTPCrefMultTPConlyAll,xCumulant4thTPCrefMultTPConlyAll,
  57                                                         yCumulant4thTPCrefMultTPConlyAll);
  58
  59   Double_t xCumulant2nd4050ALICE[] = {0.000000,0.250000,0.350000,0.450000,0.550000,0.650000,0.750000,0.850000,0.950000,
  60   1.100000,1.300000,1.500000,1.700000,1.900000,2.250000,2.750000,3.250000,3.750000,4.500000};
  61   Double_t yCumulant2nd4050ALICE[] = {0.000000,0.043400,0.059911,0.073516,0.089756,0.105486,0.117391,0.128199,0.138013,
  62   0.158271,0.177726,0.196383,0.208277,0.216648,0.242954,0.249961,0.240131,0.269006,0.207796};
  63   Int_t nPointsCumulant2nd4050ALICE = sizeof(xCumulant2nd4050ALICE)/sizeof(Double_t);
  64   fAlicePt2nd4050 = new TGraph(nPointsCumulant2nd4050ALICE,xCumulant2nd4050ALICE,yCumulant2nd4050ALICE);
  65
  66   Double_t xCumulant4th3040ALICE[] = {0.00000,0.250000,0.350000,0.450000,0.550000,0.650000,0.750000,0.850000,0.950000,
  67   1.100000,1.300000,1.500000,1.700000,1.900000,2.250000,2.750000,3.250000,3.750000,4.500000,
  68   5.500000,7.000000,9.000000};
  69   Double_t yCumulant4th3040ALICE[] = {0.000000,0.037071,0.048566,0.061083,0.070910,0.078831,0.091396,0.102026,0.109691,
  70   0.124449,0.139819,0.155561,0.165701,0.173678,0.191149,0.202015,0.204540,0.212560,0.195885,
  71   0.000000,0.000000,0.000000};
  72   Int_t nPointsCumulant4th3040ALICE = sizeof(xCumulant4th3040ALICE)/sizeof(Double_t);
  73   fAlicePt4th3040 = new TGraph(nPointsCumulant4th3040ALICE,xCumulant4th3040ALICE,yCumulant4th3040ALICE);
  74
  75   Double_t xCumulant4th4050ALICE[] = {0.000000,0.250000,0.350000,0.450000,0.550000,0.650000,0.750000,0.850000,0.950000,
  76   1.100000,1.300000,1.500000,1.700000,1.900000,2.250000,2.750000,3.250000,3.750000,4.500000};
  77   Double_t yCumulant4th4050ALICE[] = {0.000000,0.038646,0.049824,0.066662,0.075856,0.081583,0.099778,0.104674,0.118545,
  78   0.131874,0.152959,0.155348,0.169751,0.179052,0.178532,0.198851,0.185737,0.239901,0.186098};
  79   Int_t nPointsCumulant4th4050ALICE = sizeof(xCumulant4th4050ALICE)/sizeof(Double_t);
  80   fAlicePt4th4050 = new TGraph(nPointsCumulant4th4050ALICE, xCumulant4th4050ALICE, yCumulant4th4050ALICE);
  81
  82   Double_t impactParam[] = {0.,1.75,4.225,5.965,7.765,9.215,10.46,11.565,12.575,13.515,16.679};
  83   Double_t centrality[] = {0.,2.5,7.5,15,25,35,45,55,65,75,90};
  84
  85   Int_t nPoints = sizeof(impactParam)/sizeof(Double_t);
  86   fImpactParToCent = new TGraph(nPoints, impactParam, centrality);
  87
  88 }
  89 //_____________________________________________________________________
  90 Bool_t AliForwardFlowUtil::AODCheck(const AliAODForwardMult* aodfm)
  91 {
  92   //
  93   // Function to check that and AOD event meets the cuts
  94   //
  95   // Parameters:
  96   //  AliAODForwardMult: forward mult object with trigger and vertex info
  97   //
  98   fCent = -1;
  99   fVertex = 1111;
 100
 101   if (!aodfm->IsTriggerBits(AliAODForwardMult::kOffline)) return kFALSE;
 102   fCent = (Double_t)aodfm->GetCentrality();
 103   if (0. >= fCent || fCent >= 100.) return kFALSE;
 104   TH1D* vertex = (TH1D*)fList->FindObject("VertexAll");
 105   fVertex = aodfm->GetIpZ();
 106   vertex->Fill(fVertex);
 107   if (TMath::Abs(fVertex) >= 5.) return kFALSE;
 108   return kTRUE;
 109 }
 110 //_____________________________________________________________________
 111 Bool_t AliForwardFlowUtil::LoopAODFMD(const AliAODEvent* aodevent)
 112 {
 113   //
 114   // Loop over AliAODFowardMult object and fill histograms provided
 115   // by flow task.
 116   //
 117
 118   AliAODForwardMult* aodfmult = static_cast<AliAODForwardMult*>(aodevent->FindListObject("Forward"));
 119   if (!aodfmult) return kFALSE;
 120   if (!AODCheck(aodfmult)) return kFALSE;
 121
 122   TH2D hdNdetadphi = aodfmult->GetHistogram();
 123   TH2D* vertex = (TH2D*)fList->FindObject("CoverageVsVertex");
 124   TH2D* dNdphi = (TH2D*)fList->FindObject("hdNdphiSEFMD");
 125   TH2D* dNdetadphi = (TH2D*)fList->FindObject("hdNdetadphiSEFMD");
 126   dNdphi->Reset();
 127   dNdetadphi->Reset();
 128
 129   Double_t mult = 0;
 130   Double_t eta = 0;
 131   Double_t phi = 0;
 132   Double_t weight = 0;
 133   for (Int_t etaBin = 1; etaBin<=hdNdetadphi.GetNbinsX(); etaBin++) {
 134     eta = hdNdetadphi.GetXaxis()->GetBinCenter(etaBin);
 135     if (TMath::Abs(eta) < 1.75) continue;
 136     for (Int_t phiBin = 0; phiBin<=hdNdetadphi.GetNbinsY(); phiBin++) {
 137       phi = hdNdetadphi.GetYaxis()->GetBinCenter(phiBin);
 138       weight = hdNdetadphi.GetBinContent(etaBin, phiBin);
 139       if (phiBin == 0) {
 140         vertex->Fill(eta, fVertex, weight);
 141         continue;
 142       }
 143       dNdetadphi->Fill(eta, phi, weight);
 144       dNdphi->Fill(eta, phi, weight);
 145       dNdphi->Fill(-1.*eta, phi, weight);
 146       mult += weight;
 147     }
 148   }
 149 //  fCent = GetCentFromMC(aodevent);
 150 //  fCent = 0.5;
 151
 152   return kTRUE;
 153 }
 154 //_____________________________________________________________________
 155 Bool_t AliForwardFlowUtil::LoopAODSPD(const AliAODEvent* aodevent) const
 156 {
 157   //
 158   // Loop over AliAODCentralMult object and fill histograms
 159   // provided by flow task
 160   //
 161   AliAODCentralMult* aodcmult = static_cast<AliAODCentralMult*>(aodevent->FindListObject("CentralClusters"));
 162   if (!aodcmult) return kFALSE;
 163
 164   TH2D hdNdetadphi = aodcmult->GetHistogram();
 165   TH2D* vertex = (TH2D*)fList->FindObject("CoverageVsVertex");
 166   TH2D* dNdetadphi = (TH2D*)fList->FindObject("hdNdetadphiSESPD");
 167   TH2D* dNdphi = (TH2D*)fList->FindObject("hdNdphiSESPD");
 168   dNdphi->Reset();
 169   dNdetadphi->Reset();
 170
 171   Double_t eta = 0;
 172   Double_t phi = 0;
 173   Double_t weight = 0;
 174   Double_t mult = 0;
 175   for (Int_t etaBin = 1; etaBin<=hdNdetadphi.GetNbinsX(); etaBin++) {
 176     eta = hdNdetadphi.GetXaxis()->GetBinCenter(etaBin);
 177     if (TMath::Abs(eta) > 1.75) continue;
 178     for (Int_t phiBin = 0; phiBin<=hdNdetadphi.GetNbinsY(); phiBin++) {
 179       phi = hdNdetadphi.GetYaxis()->GetBinCenter(phiBin);
 180       weight = hdNdetadphi.GetBinContent(etaBin, phiBin);
 181       if (phiBin == 0) {
 182         vertex->Fill(eta, fVertex, weight);
 183         continue;
 184       }
 185       dNdetadphi->Fill(eta, phi, weight);
 186       if (TMath::Abs(eta) < 0.5) continue;
 187       dNdphi->Fill(eta, phi, weight);
 188 //      dNdphi->Fill(-1.*eta, phi, weight);
 189       mult += weight;
 190     }
 191   }
 192
 193   return kTRUE;
 194 }
 195 //_____________________________________________________________________
 196 Bool_t AliForwardFlowUtil::LoopAODFMDtrrefHits(const AliAODEvent* aodevent) const
 197 {
 198   //
 199   // Loop over AliAODForwardMult object, get MC track ref information
 200   // and fill flow histograms
 201   //
 202
 203   AliAODForwardMult* aodfmult = static_cast<AliAODForwardMult*>(aodevent->FindListObject("ForwardMC"));
 204   if (!aodfmult) return kFALSE;
 205
 206   TH2D hdNdetadphi = aodfmult->GetHistogram();
 207   TH2D* dNdphi = (TH2D*)fList->FindObject("hdNdphiSEFMDTR");
 208   TH2D* dNdetadphi = (TH2D*)fList->FindObject("hdNdetadphiSEFMDTR");
 209   dNdphi->Reset();
 210   dNdetadphi->Reset();
 211
 212   Double_t mult = 0;
 213   Double_t eta = 0;
 214   Double_t phi = 0;
 215   Double_t weight = 0;
 216   for(Int_t etaBin = 1; etaBin<=hdNdetadphi.GetNbinsX(); etaBin++) {
 217     eta = hdNdetadphi.GetXaxis()->GetBinCenter(etaBin);
 218     if (TMath::Abs(eta) < 1.75) continue;
 219     for(Int_t phiBin = 1; phiBin<=hdNdetadphi.GetNbinsY(); phiBin++) {
 220       phi = hdNdetadphi.GetYaxis()->GetBinCenter(phiBin);
 221       weight = hdNdetadphi.GetBinContent(etaBin, phiBin);
 222       dNdetadphi->Fill(eta, phi, weight);
 223       dNdphi->Fill(eta, phi, weight);
 224       dNdphi->Fill(-1.*eta, phi, weight);
 225       mult += weight;
 226     }
 227   }
 228
 229   return kTRUE;
 230 }
 231 //_____________________________________________________________________
 232 Bool_t AliForwardFlowUtil::LoopAODSPDtrrefHits(const AliAODEvent* aodevent) const
 233 {
 234   //
 235   // Loop over AliAODCentralMult object and fill histograms
 236   // provided by flow task
 237   //
 238   AliAODCentralMult* aodcmult = static_cast<AliAODCentralMult*>(aodevent->FindListObject("CentralClustersMC"));
 239   if (!aodcmult) return kFALSE;
 240
 241   TH2D hdNdetadphi = aodcmult->GetHistogram();
 242   TH2D* dNdphi = (TH2D*)fList->FindObject("hdNdphiSESPDTR");
 243   TH2D* dNdetadphi = (TH2D*)fList->FindObject("hdNdetadphiSESPDTR");
 244   dNdphi->Reset();
 245   dNdetadphi->Reset();
 246
 247   Double_t eta = 0;
 248   Double_t phi = 0;
 249   Double_t weight = 0;
 250   Double_t mult = 0;
 251   for (Int_t etaBin = 1; etaBin<=hdNdetadphi.GetNbinsX(); etaBin++) {
 252     eta = hdNdetadphi.GetXaxis()->GetBinCenter(etaBin);
 253     if (TMath::Abs(eta) > 1.75) continue;
 254     for (Int_t phiBin = 1; phiBin<=hdNdetadphi.GetNbinsY(); phiBin++) {
 255       phi = hdNdetadphi.GetYaxis()->GetBinCenter(phiBin);
 256       weight = hdNdetadphi.GetBinContent(etaBin, phiBin);
 257       dNdetadphi->Fill(eta, phi, weight);
 258       dNdphi->Fill(eta, phi, weight);
 259       dNdphi->Fill(-1.*eta, phi, weight);
 260       mult += weight;
 261     }
 262   }
 263
 264   return kTRUE;
 265 }
 266 //_____________________________________________________________________
 267 Bool_t AliForwardFlowUtil::LoopAODmc(const AliAODEvent* aodevent,
 268                                      TString addFlow = "",
 269                                      Int_t type = 0,
 270                                      Int_t order = 2) const
 271 {
 272   //
 273   // Loop over AliAODParticle branch and fill flow histograms
 274   //
 275
 276   //retreive MC particles from event
 277   TClonesArray* mcArray = (TClonesArray*)aodevent->FindListObject(AliAODMCParticle::StdBranchName());
 278   if(!mcArray){
 279 //    AliWarning("No MC array found in AOD. Try making it again.");
 280     return kFALSE;
 281   }
 282
 283   Double_t rp = 0;
 284   AliAODMCHeader* header = dynamic_cast<AliAODMCHeader*>(aodevent->FindListObject(AliAODMCHeader::StdBranchName()));
 285   if (!header) {
 286     AliWarning("No header file found.");
 287   }
 288   else {
 289     rp = header->GetReactionPlaneAngle();
 290   }
 291
 292   TH2D* dNdetadphi = (TH2D*)fList->FindObject("hdNdetadphiSEMC");
 293   TH2D* dNdphi = (TH2D*)fList->FindObject("hdNdphiSEMC");
 294   TProfile2D* mcTruth = (TProfile2D*)fList->FindObject("pMCTruth");
 295   dNdphi->Reset();
 296   dNdetadphi->Reset();
 297
 298   Int_t ntracks = mcArray->GetEntriesFast();
 299
 300   // Track loop: chek how many particles will be accepted
 301   Double_t mult = 0;
 302   Double_t weight = 0;
 303   for (Int_t it = 0; it < ntracks; it++) {
 304     weight = 1;
 305     AliAODMCParticle* particle = (AliAODMCParticle*)mcArray->At(it);
 306     if (!particle) {
 307       AliError(Form("Could not receive track %d", it));
 308       continue;
 309     }
 310     if (!particle->IsPhysicalPrimary()) continue;
 311     if (particle->Charge() == 0) continue;
 312     Double_t pT = particle->Pt();
 313 //    if (pT <= 0.2 || pT >= 5.) continue;
 314     Double_t eta = particle->Eta();
 315     Double_t phi = particle->Phi();
 316     if (TMath::Abs(eta) < 6.) {
 317       // Add flow if it is in the argument
 318       if (addFlow.Length() > 1) {
 319         if (addFlow.Contains("pt"))
 320           weight *= AddptFlow(pT, type);
 321 //          weight *= AddptFlow(pT, 1);
 322         if (addFlow.Contains("pid"))
 323           weight *= AddpidFlow(particle->PdgCode(), 1);
 324         if (addFlow.Contains("eta"))
 325           weight *= AddetaFlow(eta, type);
 326 //          weight *= AddetaFlow(eta, 2);
 327         if (addFlow.Contains("cent"))
 328           weight *= fAliceCent4th->Eval(fCent)/fAliceCent4th->Eval(45);
 329
 330         weight *= 20*2.*TMath::Cos((Double_t)order*(phi-rp));
 331         weight += 1;
 332       }
 333       dNdphi->Fill(eta, phi, weight);
 334       dNdphi->Fill(-1.*eta, phi, weight);
 335       dNdetadphi->Fill(eta, phi, weight);
 336       mcTruth->Fill(eta, fCent, weight*TMath::Cos(2.*(phi-rp)));
 337       mult += weight;
 338     }
 339   }
 340
 341   return kTRUE;
 342 }
 343 //_____________________________________________________________________
 344 Double_t AliForwardFlowUtil::AddptFlow(Double_t pt = 0, Int_t type = 0) const
 345 {
 346   //
 347   // Add pt dependent flow factor
 348   //
 349   Double_t weight = 0;
 350
 351   switch(type)
 352   {
 353     case 1: weight = fAlicePt2nd4050->Eval(pt)*0.5+fAlicePt4th4050->Eval(pt)*0.5;
 354             break;
 355     case 2: weight = fAlicePt2nd4050->Eval(pt);
 356             break;
 357     case 3: weight = fAlicePt4th3040->Eval(pt);
 358             break;
 359     case 4: weight = fAlicePt4th4050->Eval(pt);
 360             break;
 361   }
 362
 363   return weight;
 364 }
 365 //_____________________________________________________________________
 366 Double_t AliForwardFlowUtil::AddpidFlow(Int_t id = 0, Int_t type = 0) const
 367 {
 368   //
 369   // Add pid dependent flow factor
 370   //
 371   Double_t weight = 0;
 372
 373   switch(type)
 374   {
 375     case 1: if (TMath::Abs(id) ==  211) // pion flow
 376               weight = 1.;
 377             else if (TMath::Abs(id) ==  2212) // proton flow
 378               weight = 1.3;
 379             else
 380               weight = 0.7;
 381             break;
 382     case 2: weight = 1.207;
 383             break;
 384   }
 385
 386   return weight;
 387 }
 388 //_____________________________________________________________________
 389 Double_t AliForwardFlowUtil::AddetaFlow(Double_t eta = 0, Int_t type = 0) const
 390 {
 391   //
 392   // Add eta dependent flow factor
 393   //
 394   Double_t weight = 0;
 395
 396   TF1 gaus = TF1("gaus", "gaus", -6, 6);
 397
 398   switch(type)
 399   {
 400      case 1: gaus.SetParameters(0.1, 0., 9);
 401              break;
 402      case 2: gaus.SetParameters(0.1, 0., 3);
 403              break;
 404      case 3: gaus.SetParameters(0.1, 0., 15);
 405              break;
 406   }
 407
 408   weight = gaus.Eval(eta);
 409
 410   return weight;
 411 }
 412 //_____________________________________________________________________
 413 Double_t AliForwardFlowUtil::GetCentFromMC(const AliAODEvent* aodevent) const
 414 {
 415   //
 416   // Get centrality from MC impact parameter.
 417   // Values taken from: https://twiki.cern.ch/twiki/bin/viewauth/ALICE/CentStudies
 418   //
 419   Double_t cent = -1.;
 420   Double_t b = -1.;
 421   AliAODMCHeader* header = (AliAODMCHeader*)aodevent->FindListObject(AliAODMCHeader::StdBranchName());
 422   if (!header) return cent;
 423   b = header->GetImpactParameter();
 424
 425   cent = fImpactParToCent->Eval(b);
 426
 427   return cent;
 428 }
 429 //_____________________________________________________________________
 430 //
 431 //
 432 // EOF
 433