PWGLF/FORWARD/analysis2/AliForwardMCFlowTaskQC.cxx

   1 //
   2 // Calculate flow on MC data in the forward and central regions using the Q cumulants method.
   3 //
   4 // Inputs:
   5 //  - AliAODEvent
   6 //
   7 // Outputs:
   8 //  - AnalysisResults.root
   9 //
  10 #include "AliForwardMCFlowTaskQC.h"
  11 #include "AliAODMCParticle.h"
  12 #include "AliAODMCHeader.h"
  13 #include "TGraph.h"
  14 #include "TF1.h"
  15 #include "TProfile2D.h"
  16 #include "AliAODEvent.h"
  17 #include "AliAODForwardMult.h"
  18 #include "AliAODCentralMult.h"
  19
  20 ClassImp(AliForwardMCFlowTaskQC)
  21
  22 //_____________________________________________________________________
  23 AliForwardMCFlowTaskQC::AliForwardMCFlowTaskQC() :
  24   AliForwardFlowTaskQC(),
  25   fBinsFMDTR(),          // List of FMDTR analysis objects
  26   fBinsSPDTR(),          // List of SPDTR analysis objects
  27   fBinsMC(),             // List of MC truth analysis objects
  28   fdNdedpMC(),           // MC truth d^2N/detadphi histogram
  29   fAliceCent4th(),       // Alice QC4 vs. centrality data points
  30   fAlicePt2nd4050(),     // Alice QC2 vs. pT data points
  31   fAlicePt4th3040(),     // Alice QC4 vs. pT data points
  32   fAlicePt4th4050(),     // Alice QC4 vs. pT data points
  33   fImpactParToCent(),    // Impact parameter to centrality graph
  34   fAddFlow(0),           // Add flow to MC truth
  35   fAddType(0),           // Add type of flow to MC truth
  36   fAddOrder(0)           // Add order of flow to MC truth
  37    {}
  38   //
  39   // Default Constructor
  40   //
  41 //_____________________________________________________________________
  42 AliForwardMCFlowTaskQC::AliForwardMCFlowTaskQC(const char* name) :
  43   AliForwardFlowTaskQC(name),
  44   fBinsFMDTR(),          // List of FMDTR analysis objects
  45   fBinsSPDTR(),          // List of SPDTR analysis objects
  46   fBinsMC(),             // List of MC truth analysis objects
  47   fdNdedpMC(),           // MC truth d^2N/detadphi histogram
  48   fAliceCent4th(),       // Alice QC4 vs. centrality data points
  49   fAlicePt2nd4050(),     // Alice QC2 vs. pT data points
  50   fAlicePt4th3040(),     // Alice QC4 vs. pT data points
  51   fAlicePt4th4050(),     // Alice QC4 vs. pT data points
  52   fImpactParToCent(),    // Impact parameter to centrality graph
  53   fAddFlow(0),           // Add flow to MC truth
  54   fAddType(0),           // Add type of flow to MC truth
  55   fAddOrder(0)           // Add order of flow to MC truth
  56 {
  57   //
  58   // Constructor
  59   // Parameters:
  60   //  name: Name of task
  61   //
  62   fdNdedpMC = TH2D("fdNdedpMC", "fdNdedpMC", 48, -6., 6., 200, 0., 2.*TMath::Pi());
  63   fdNdedpMC.Sumw2();
  64
  65   // Add parametrizations of ALICE data
  66   Double_t xCumulant4thTPCrefMultTPConlyAll[] = {2.5,7.5,15,25,35,45,55,65};
  67   Double_t yCumulant4thTPCrefMultTPConlyAll[] = {0.017855,0.032440,0.055818,0.073137,0.083898,0.086690,0.082040,0.077777};
  68   Int_t nPointsCumulant4thTPCrefMultTPConlyAll = sizeof(xCumulant4thTPCrefMultTPConlyAll)/sizeof(Double_t);
  69   fAliceCent4th = new TGraph(nPointsCumulant4thTPCrefMultTPConlyAll,xCumulant4thTPCrefMultTPConlyAll,
  70                                                         yCumulant4thTPCrefMultTPConlyAll);
  71
  72   Double_t xCumulant2nd4050ALICE[] = {0.000000,0.250000,0.350000,0.450000,0.550000,0.650000,0.750000,0.850000,0.950000,
  73   1.100000,1.300000,1.500000,1.700000,1.900000,2.250000,2.750000,3.250000,3.750000,4.500000};
  74   Double_t yCumulant2nd4050ALICE[] = {0.000000,0.043400,0.059911,0.073516,0.089756,0.105486,0.117391,0.128199,0.138013,
  75   0.158271,0.177726,0.196383,0.208277,0.216648,0.242954,0.249961,0.240131,0.269006,0.207796};
  76   Int_t nPointsCumulant2nd4050ALICE = sizeof(xCumulant2nd4050ALICE)/sizeof(Double_t);
  77   fAlicePt2nd4050 = new TGraph(nPointsCumulant2nd4050ALICE,xCumulant2nd4050ALICE,yCumulant2nd4050ALICE);
  78
  79   Double_t xCumulant4th3040ALICE[] = {0.00000,0.250000,0.350000,0.450000,0.550000,0.650000,0.750000,0.850000,0.950000,
  80   1.100000,1.300000,1.500000,1.700000,1.900000,2.250000,2.750000,3.250000,3.750000,4.500000,
  81   5.500000,7.000000,9.000000};
  82   Double_t yCumulant4th3040ALICE[] = {0.000000,0.037071,0.048566,0.061083,0.070910,0.078831,0.091396,0.102026,0.109691,
  83   0.124449,0.139819,0.155561,0.165701,0.173678,0.191149,0.202015,0.204540,0.212560,0.195885,
  84   0.000000,0.000000,0.000000};
  85   Int_t nPointsCumulant4th3040ALICE = sizeof(xCumulant4th3040ALICE)/sizeof(Double_t);
  86   fAlicePt4th3040 = new TGraph(nPointsCumulant4th3040ALICE,xCumulant4th3040ALICE,yCumulant4th3040ALICE);
  87
  88   Double_t xCumulant4th4050ALICE[] = {0.000000,0.250000,0.350000,0.450000,0.550000,0.650000,0.750000,0.850000,0.950000,
  89   1.100000,1.300000,1.500000,1.700000,1.900000,2.250000,2.750000,3.250000,3.750000,4.500000};
  90   Double_t yCumulant4th4050ALICE[] = {0.000000,0.038646,0.049824,0.066662,0.075856,0.081583,0.099778,0.104674,0.118545,
  91   0.131874,0.152959,0.155348,0.169751,0.179052,0.178532,0.198851,0.185737,0.239901,0.186098};
  92   Int_t nPointsCumulant4th4050ALICE = sizeof(xCumulant4th4050ALICE)/sizeof(Double_t);
  93   fAlicePt4th4050 = new TGraph(nPointsCumulant4th4050ALICE, xCumulant4th4050ALICE, yCumulant4th4050ALICE);
  94
  95   Double_t impactParam[] = {0.,1.75,4.225,5.965,7.765,9.215,10.46,11.565,12.575,13.515,16.679};
  96   Double_t centrality[] = {0.,2.5,7.5,15,25,35,45,55,65,75,90};
  97
  98   Int_t nPoints = sizeof(impactParam)/sizeof(Double_t);
  99   fImpactParToCent = new TGraph(nPoints, impactParam, centrality);
 100
 101 }
 102 //_____________________________________________________________________
 103 AliForwardMCFlowTaskQC::AliForwardMCFlowTaskQC(const AliForwardMCFlowTaskQC& o) :
 104   AliForwardFlowTaskQC(o),
 105   fBinsFMDTR(),                          // List of FMDTR analysis objects
 106   fBinsSPDTR(),                          // List of SPDTR analysis objects
 107   fBinsMC(),                             // List of MC truth analysis objects
 108   fdNdedpMC(o.fdNdedpMC),                // MC truth d^2N/detadphi histogram
 109   fAliceCent4th(o.fAliceCent4th),        // Alice QC4 vs. centrality data points
 110   fAlicePt2nd4050(o.fAlicePt2nd4050),    // Alice QC2 vs. pT data points
 111   fAlicePt4th3040(o.fAlicePt4th3040),    // Alice QC4 vs. pT data points
 112   fAlicePt4th4050(o.fAlicePt4th4050),    // Alice QC4 vs. pT data points
 113   fImpactParToCent(o.fImpactParToCent),  // Impact parameter to centrality graph
 114   fAddFlow(o.fAddFlow),                  // Add flow to MC truth
 115   fAddType(o.fAddType),                  // Add type of flow to MC truth
 116   fAddOrder(o.fAddOrder)                 // Add order of flow to MC truth
 117    {}
 118   //
 119   // Copy Constructor
 120   //
 121 //_____________________________________________________________________
 122 AliForwardMCFlowTaskQC&
 123 AliForwardMCFlowTaskQC::operator=(const AliForwardMCFlowTaskQC& o)
 124 {
 125   //
 126   // Assignment operator
 127   // Parameters:
 128   //    o Object to copy from
 129   //
 130   fdNdedpMC        = o.fdNdedpMC;
 131   fAliceCent4th    = o.fAliceCent4th;
 132   fAlicePt2nd4050  = o.fAlicePt2nd4050;
 133   fAlicePt4th3040  = o.fAlicePt4th3040;
 134   fAlicePt4th4050  = o.fAlicePt4th4050;
 135   fImpactParToCent = o.fImpactParToCent;
 136   fAddFlow         = o.fAddFlow;
 137   fAddType         = o.fAddType;
 138   fAddOrder        = o.fAddOrder;
 139   return *this;
 140 }
 141 //_____________________________________________________________________
 142 void AliForwardMCFlowTaskQC::InitVertexBins()
 143 {
 144   //
 145   // Initiate VertexBin lists
 146   //
 147   AliForwardFlowTaskQC::InitVertexBins();
 148
 149   for(Int_t n = 1; n <= 6; n++) {
 150     if (!fv[n]) continue;
 151     for (Int_t v = -10; v < 10; v++) {
 152       fBinsFMDTR.Add(new VertexBin(v, v+1, n, "FMDTR"));
 153       fBinsSPDTR.Add(new VertexBin(v, v+1, n, "SPDTR"));
 154       fBinsMC.Add(new VertexBin(v, v+1, n, "MC"));
 155     }
 156   }
 157
 158 }
 159 //_____________________________________________________________________
 160 void AliForwardMCFlowTaskQC::InitHists()
 161 {
 162   //
 163   // Initiate diagnostics hists and add to outputlist
 164   //
 165   AliForwardFlowTaskQC::InitHists();
 166
 167   TIter nextFMDTR(&fBinsFMDTR);
 168   VertexBin* bin = 0;
 169   while ((bin = static_cast<VertexBin*>(nextFMDTR()))) {
 170     bin->AddOutput(fSumList);
 171   }
 172   TIter nextSPDTR(&fBinsSPDTR);
 173   while ((bin = static_cast<VertexBin*>(nextSPDTR()))) {
 174     bin->AddOutput(fSumList);
 175   }
 176   TIter nextMC(&fBinsMC);
 177   while ((bin = static_cast<VertexBin*>(nextMC()))) {
 178     bin->AddOutput(fSumList);
 179   }
 180
 181 }
 182 //_____________________________________________________________________
 183 Bool_t AliForwardMCFlowTaskQC::Analyze()
 184 {
 185   //
 186   // Load FMD and SPD MC objects from aod tree and call the cumulants
 187   // calculation for the correct vertexbin
 188   //
 189   if (!AliForwardFlowTaskQC::Analyze()) return kFALSE;
 190
 191   // Run analysis on trackrefs from FMD and SPD
 192   AliAODForwardMult* aodfmult = static_cast<AliAODForwardMult*>(fAOD->FindListObject("ForwardMC"));
 193   if (!aodfmult) return kFALSE;
 194   TH2D fmdTRdNdetadphi = aodfmult->GetHistogram();
 195
 196   AliAODCentralMult* aodcmult = static_cast<AliAODCentralMult*>(fAOD->FindListObject("CentralClustersMC"));
 197   if (!aodcmult) return kFALSE;
 198   TH2D spdTRdNdetadphi = aodcmult->GetHistogram();
 199
 200   TIter nextFMDTR(&fBinsFMDTR);
 201   VertexBin* bin = 0;
 202   while ((bin = static_cast<VertexBin*>(nextFMDTR()))) {
 203     if (bin->CheckVertex(fZvertex)) {
 204       if (!bin->FillHists(fmdTRdNdetadphi)) return kFALSE;
 205       bin->CumulantsAccumulate(fCent);
 206     }
 207   }
 208
 209   TIter nextSPDTR(&fBinsSPDTR);
 210   while ((bin = static_cast<VertexBin*>(nextSPDTR()))) {
 211     if (bin->CheckVertex(fZvertex)) {
 212       if (!bin->FillHists(spdTRdNdetadphi)) return kFALSE;
 213       bin->CumulantsAccumulate(fCent);
 214     }
 215   }
 216
 217   // Run analysis on MC branch
 218   if (!LoopAODMC()) return kFALSE;
 219
 220   TIter nextMC(&fBinsMC);
 221   while ((bin = static_cast<VertexBin*>(nextMC()))) {
 222     if (bin->CheckVertex(fZvertex)) {
 223       if (!bin->FillHists(fdNdedpMC)) return kFALSE;
 224       bin->CumulantsAccumulate(fCent);
 225     }
 226   }
 227
 228   return kTRUE;
 229 }
 230 //_____________________________________________________________________
 231 void AliForwardMCFlowTaskQC::Finalize()
 232 {
 233   //
 234   // Finalize command, called by Terminate()
 235   //
 236   AliForwardFlowTaskQC::Finalize();
 237
 238   TIter nextFMDTR(&fBinsFMDTR);
 239   VertexBin* bin = 0;
 240   while ((bin = static_cast<VertexBin*>(nextFMDTR()))) {
 241     bin->CumulantsTerminate(fSumList, fOutputList);
 242   }
 243   TIter nextSPDTR(&fBinsSPDTR);
 244   while ((bin = static_cast<VertexBin*>(nextSPDTR()))) {
 245     bin->CumulantsTerminate(fSumList, fOutputList);
 246   }
 247   TIter nextMC(&fBinsMC);
 248   while ((bin = static_cast<VertexBin*>(nextMC()))) {
 249     bin->CumulantsTerminate(fSumList, fOutputList);
 250   }
 251
 252   TProfile2D* fmdHist = 0;
 253   TProfile2D* spdHist = 0;
 254   TProfile2D* mcHist = 0;
 255   TList* correctionList = new TList();
 256   correctionList->SetName("CorrectionList");
 257 //  fOutputList->Add(correctionList);
 258
 259   for (Int_t i = 2; i <= 4; i += 2) {
 260     for (Int_t n = 1; n <= 6; n++) {
 261       if (!fv[n]) continue;
 262       fmdHist = (TProfile2D*)fOutputList->FindObject(Form("FMDQC%d_v%d_unCorr", i, n))
 263                   ->Clone(Form("FMDQC%d_v%d_Correction", i, n));
 264       spdHist = (TProfile2D*)fOutputList->FindObject(Form("SPDQC%d_v%d_unCorr", i, n))
 265                   ->Clone(Form("SPDQC%d_v%d_Correction", i, n));
 266       mcHist = (TProfile2D*)fOutputList->FindObject(Form("MCQC%d_v%d_unCorr", i, n));
 267
 268       if (!fmdHist || !spdHist || !mcHist) {
 269         AliError(Form("Histogram missing, correction object not created for v%d", n));
 270         continue;
 271       }
 272
 273       fmdHist->Divide(mcHist);
 274       spdHist->Divide(mcHist);
 275       fmdHist->SetTitle(Form("FMD QC{%d} v_{%d} Correction Object", i, n));
 276       fmdHist->SetTitle(Form("SPD QC{%d} v_{%d} Correction Object", i, n));
 277
 278 //      correctionList->Add(fmdHist);
 279 //      correctionList->Add(spdHist);
 280     }
 281   }
 282
 283 }
 284 //_____________________________________________________________________
 285 Bool_t AliForwardMCFlowTaskQC::LoopAODMC()
 286 {
 287   //
 288   // Loop over AliAODParticle branch and fill d^2N/detadphi-histograms.
 289   // Add flow if set to do so in AddTask function
 290   fdNdedpMC.Reset();
 291
 292   //retreive MC particles from event
 293   TClonesArray* mcArray = (TClonesArray*)fAOD->FindListObject(AliAODMCParticle::StdBranchName());
 294   if(!mcArray){
 295 //    AliWarning("No MC array found in AOD. Try making it again.");
 296     return kFALSE;
 297   }
 298
 299   Double_t rp = 0;
 300   AliAODMCHeader* header = dynamic_cast<AliAODMCHeader*>(fAOD->FindListObject(AliAODMCHeader::StdBranchName()));
 301   if (!header) {
 302     AliWarning("No header file found.");
 303   }
 304   else {
 305     rp = header->GetReactionPlaneAngle();
 306   }
 307
 308   Int_t ntracks = mcArray->GetEntriesFast();
 309
 310   // Track loop: chek how many particles will be accepted
 311   Double_t weight = 0;
 312   for (Int_t it = 0; it < ntracks; it++) {
 313     weight = 1;
 314     AliAODMCParticle* particle = (AliAODMCParticle*)mcArray->At(it);
 315     if (!particle) {
 316       AliError(Form("Could not receive track %d", it));
 317       continue;
 318     }
 319     if (!particle->IsPhysicalPrimary()) continue;
 320     if (particle->Charge() == 0) continue;
 321     Double_t pT = particle->Pt();
 322     Double_t eta = particle->Eta();
 323     Double_t phi = particle->Phi();
 324     if (TMath::Abs(eta) < 6.) {
 325       // Add flow if it is in the argument
 326       if (fAddFlow.Length() > 1) {
 327         if (fAddFlow.Contains("pt"))
 328           weight *= AddptFlow(pT);
 329         if (fAddFlow.Contains("pid"))
 330           weight *= AddpidFlow(particle->PdgCode());
 331         if (fAddFlow.Contains("eta"))
 332           weight *= AddetaFlow(eta);
 333         if (fAddFlow.Contains("cent"))
 334           weight *= fAliceCent4th->Eval(fCent)/fAliceCent4th->Eval(45);
 335
 336         weight *= 20*2.*TMath::Cos((Double_t)fAddOrder*(phi-rp));
 337         weight += 1;
 338       }
 339       fdNdedpMC.Fill(eta, phi, weight);
 340     }
 341   }
 342
 343   return kTRUE;
 344 }
 345 //_____________________________________________________________________
 346 Double_t AliForwardMCFlowTaskQC::AddptFlow(Double_t pt = 0) const
 347 {
 348   //
 349   // Add pt dependent flow factor
 350   // Parameters:
 351   //  pt: pT parametrization to use
 352   //
 353   Double_t weight = 0;
 354
 355   switch(fAddType)
 356   {
 357     case 1: weight = fAlicePt2nd4050->Eval(pt)*0.5+fAlicePt4th4050->Eval(pt)*0.5;
 358             break;
 359     case 2: weight = fAlicePt2nd4050->Eval(pt);
 360             break;
 361     case 3: weight = fAlicePt4th3040->Eval(pt);
 362             break;
 363     case 4: weight = fAlicePt4th4050->Eval(pt);
 364             break;
 365   }
 366
 367   return weight;
 368 }
 369 //_____________________________________________________________________
 370 Double_t AliForwardMCFlowTaskQC::AddpidFlow(Int_t id = 0) const
 371 {
 372   //
 373   // Add pid dependent flow factor
 374   // Parameters:
 375   //  id: choose PID dependent setup
 376   //
 377   Double_t weight = 0;
 378
 379   switch(fAddType)
 380   {
 381     case 1: if (TMath::Abs(id) ==  211) // pion flow
 382               weight = 1.3;
 383             else if (TMath::Abs(id) ==  2212) // proton flow
 384               weight = 1.;
 385             else
 386               weight = 0.7;
 387             break;
 388     case 2: weight = 1.207;
 389             break;
 390   }
 391
 392   return weight;
 393 }
 394 //_____________________________________________________________________
 395 Double_t AliForwardMCFlowTaskQC::AddetaFlow(Double_t eta = 0) const
 396 {
 397   //
 398   // Add eta dependent flow factor
 399   // Parameters:
 400   //  eta: choose v_n(eta) shape
 401   //
 402   Double_t weight = 0;
 403
 404   TF1 gaus = TF1("gaus", "gaus", -6, 6);
 405
 406   switch(fAddType)
 407   {
 408      case 1: gaus.SetParameters(0.1, 0., 9);
 409              break;
 410      case 2: gaus.SetParameters(0.1, 0., 3);
 411              break;
 412      case 3: gaus.SetParameters(0.1, 0., 15);
 413              break;
 414   }
 415
 416   weight = gaus.Eval(eta);
 417
 418   return weight;
 419 }
 420 //_____________________________________________________________________
 421 Double_t AliForwardMCFlowTaskQC::GetCentFromB() const
 422 {
 423   //
 424   // Get centrality from MC impact parameter.
 425   // Values taken from: https://twiki.cern.ch/twiki/bin/viewauth/ALICE/CentStudies
 426   //
 427   Double_t cent = -1.;
 428   Double_t b = -1.;
 429   AliAODMCHeader* header = (AliAODMCHeader*)fAOD->FindListObject(AliAODMCHeader::StdBranchName());
 430   if (!header) return cent;
 431   b = header->GetImpactParameter();
 432
 433   cent = fImpactParToCent->Eval(b);
 434
 435   return cent;
 436 }
 437 //_____________________________________________________________________
 438 //
 439 //
 440 // EOF
 441