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 "AliAODEvent.h"
  16 #include "AliAODForwardMult.h"
  17 #include "AliAODCentralMult.h"
  18 #include "AliGenEventHeader.h"
  19 #include "AliAnalysisManager.h"
  20 #include "AliInputEventHandler.h"
  21 #include "AliGenEventHeaderTunedPbPb.h"
  22 #include "AliCollisionGeometry.h"
  23
  24 ClassImp(AliForwardMCFlowTaskQC)
  25 #if 0
  26 ;
  27 #endif
  28 //_____________________________________________________________________
  29 AliForwardMCFlowTaskQC::AliForwardMCFlowTaskQC()
  30   : AliForwardFlowTaskQC(),
  31     fBinsForwardTR(),      // List of FMDTR analysis objects
  32     fBinsCentralTR(),      // List of SPDTR analysis objects
  33     fBinsMC(),             // List of MC particle analysis objects
  34     fAODMCHeader(),        // MC Header
  35     fHistdNdedpMC(),       // MC particle d^2N/detadphi histogram
  36     fHistFMDMCCorr(),      // FMD MC correlation
  37     fHistSPDMCCorr(),      // SPD MC correlation
  38     fWeights(0),            // Flow weights
  39     fImpactParToCent(),    // Impact parameter to centrality graph
  40     fUseImpactPar(0),      // Use impact par for centrality
  41     fUseMCVertex(0),       // Get vertex from MC header?
  42     fUseFlowWeights(0)     // Add flow to MC particles
  43 {}
  44   //
  45   // Default Constructor
  46   //
  47 //_____________________________________________________________________
  48 AliForwardMCFlowTaskQC::AliForwardMCFlowTaskQC(const char* name)
  49   : AliForwardFlowTaskQC(name),
  50     fBinsForwardTR(),      // List of FMDTR analysis objects
  51     fBinsCentralTR(),      // List of SPDTR analysis objects
  52     fBinsMC(),             // List of MC particles analysis objects
  53     fAODMCHeader(0),       // MC Header
  54     fHistdNdedpMC(),       // MC particles d^2N/detadphi histogram
  55     fHistFMDMCCorr(),      // FMD MC correlation
  56     fHistSPDMCCorr(),      // SPD MC correlation
  57     fWeights(0),            // Flow weights
  58     fImpactParToCent(),    // Impact parameter to centrality graph
  59     fUseImpactPar(0),      // Use impact par for centrality
  60     fUseMCVertex(0),       // Get vertex from MC header?
  61     fUseFlowWeights(0)     // Add flow to MC particles
  62 {
  63   //
  64   // Constructor
  65   // Parameters:
  66   //  name: Name of task
  67   //
  68
  69   //  Double_t impactParam[] = {0.,1.75,4.225,5.965,7.765,9.215,10.46,
  70   //                            11.565,12.575,13.515,16.679};
  71   //  Double_t centrality[] = {0.,2.5,7.5,15,25,35,45,55,65,75,90};
  72   Double_t impactParam[] = { 0.00,  3.72,  5.23,  7.31,  8.88, 10.20,
  73                             11.38, 12.47, 13.50, 14.51, 16.679};
  74   Double_t centrality[]  = { 0.,    5.,   10.,   20.,   30.,   40.,
  75                             50.,   60.,   70.,   80.,  100.};
  76
  77   Int_t nPoints = sizeof(impactParam)/sizeof(Double_t);
  78   fImpactParToCent = new TGraph(nPoints, impactParam, centrality);
  79
  80 }
  81 //_____________________________________________________________________
  82 AliForwardMCFlowTaskQC::AliForwardMCFlowTaskQC(const AliForwardMCFlowTaskQC& o)
  83   : AliForwardFlowTaskQC(o),
  84     fBinsForwardTR(),                      // List of FMDTR analysis objects
  85     fBinsCentralTR(),                      // List of SPDTR analysis objects
  86     fBinsMC(),                             // List of MC particles analysis objects
  87     fAODMCHeader(o.fAODMCHeader),          // MC Header
  88     fHistdNdedpMC(o.fHistdNdedpMC),        // MC particles d^2N/detadphi histogram
  89     fHistFMDMCCorr(o.fHistFMDMCCorr),      // FMD MC correlation
  90     fHistSPDMCCorr(o.fHistSPDMCCorr),      // SPD MC correlation
  91     fWeights(o.fWeights),                  // Flow weights
  92     fImpactParToCent(o.fImpactParToCent),  // Impact parameter to centrality
  93     fUseImpactPar(o.fUseImpactPar),        // Use impact par for centrality
  94     fUseMCVertex(o.fUseMCVertex),          // Get vertex from MC header?
  95     fUseFlowWeights(o.fUseFlowWeights)     // Add flow to MC particles
  96 {
  97   //
  98   // Copy Constructor
  99   //
 100 }
 101 //_____________________________________________________________________
 102 AliForwardMCFlowTaskQC&
 103 AliForwardMCFlowTaskQC::operator=(const AliForwardMCFlowTaskQC& o)
 104 {
 105   //
 106   // Assignment operator
 107   // Parameters:
 108   //    o Object to copy from
 109   //
 110   if (&o == this) return *this;
 111   fAODMCHeader     = o.fAODMCHeader;
 112   fHistdNdedpMC    = o.fHistdNdedpMC;
 113   fHistFMDMCCorr   = o.fHistFMDMCCorr;
 114   fHistSPDMCCorr   = o.fHistSPDMCCorr;
 115   fWeights         = o.fWeights;
 116   fImpactParToCent = o.fImpactParToCent;
 117   fUseImpactPar    = o.fUseImpactPar;
 118   fUseMCVertex     = o.fUseMCVertex;
 119   fUseFlowWeights  = o.fUseFlowWeights;
 120   return *this;
 121 }
 122 //_____________________________________________________________________
 123 void AliForwardMCFlowTaskQC::InitVertexBins()
 124 {
 125   //
 126   // Initiate VertexBin lists
 127   //
 128   AliForwardFlowTaskQC::InitVertexBins();
 129
 130   Bool_t isNUA = (fFlowFlags & kNUAcorr);
 131   for (Int_t v = 1; v <= fVtxAxis->GetNbins(); v++) {
 132     Int_t vL = Int_t(fVtxAxis->GetBinLowEdge(v));
 133     Int_t vH = Int_t(fVtxAxis->GetBinUpEdge(v));
 134     // FMD
 135     if ((fFlowFlags & kFMD)) {
 136       fBinsForwardTR.Add(new VertexBin(vL, vH, fMaxMoment, "FMDTR", fFlowFlags, fFMDCut, fEtaGap));
 137       if (!(fFlowFlags & k3Cor))
 138         fBinsCentralTR.Add(new VertexBin(vL, vH, fMaxMoment, "SPDTR", fFlowFlags|kSPD, fSPDCut, fEtaGap));
 139       if (isNUA) fFlowFlags ^= kNUAcorr;
 140       fBinsMC.Add(new VertexBin(vL, vH, fMaxMoment, "MC-FMD", fFlowFlags|kMC, -1, fEtaGap));
 141       if ((fFlowFlags & kStdQC))
 142         fBinsMC.Add(new VertexBin(vL, vH, fMaxMoment, "MC-SPD", fFlowFlags|kMC|kSPD, -1, fEtaGap));
 143       if (isNUA) fFlowFlags ^= kNUAcorr;
 144     }
 145     // VZERO
 146     else if ((fFlowFlags & kVZERO)) {
 147       fBinsMC.Add(new VertexBin(vL, vH, fMaxMoment, "MC-VZERO", fFlowFlags|kMC, -1, fEtaGap));
 148     }
 149   }
 150 }
 151 //_____________________________________________________________________
 152 void AliForwardMCFlowTaskQC::InitHists()
 153 {
 154   //
 155   // Initiate diagnostics hists and add to outputlist
 156   //
 157   AliForwardFlowTaskQC::InitHists();
 158
 159   TString subDetName = ((fFlowFlags & kFMD) ? "FMD" : ((fFlowFlags & kVZERO) ? "VZERO" : "none"));
 160   fHistdNdedpMC = TH2D(Form("fdNdedpMC%s%s", subDetName.Data(), GetQCType(fFlowFlags)),
 161                    Form("fdNdedpMC%s%s", subDetName.Data(), GetQCType(fFlowFlags)),
 162                    240, -6., 6., 200, 0., TMath::TwoPi());
 163
 164   fHistFMDMCCorr = new TH2D("hFMDMCCorr", "hFMDMCCorr", 200, 0., 15000., 200, 0, 20000);
 165   fHistSPDMCCorr = new TH2D("hSPDMCCorr", "hSPDMCCorr", 200, 0., 7500., 200, 0, 20000);
 166   TList* dList = (TList*)fSumList->FindObject("Diagnostics");
 167   if (!dList) {
 168     dList = new TList();
 169     dList->SetName("Diagnostics");
 170     fSumList->Add(dList);
 171   }
 172   dList->Add(fHistFMDMCCorr);
 173   dList->Add(fHistSPDMCCorr);
 174
 175   TIter nextForwardTR(&fBinsForwardTR);
 176   VertexBin* bin = 0;
 177   while ((bin = static_cast<VertexBin*>(nextForwardTR()))) {
 178     bin->AddOutput(fSumList, fCentAxis);
 179   }
 180   TIter nextCentralTR(&fBinsCentralTR);
 181   while ((bin = static_cast<VertexBin*>(nextCentralTR()))) {
 182     bin->AddOutput(fSumList, fCentAxis);
 183   }
 184   TIter nextMC(&fBinsMC);
 185   while ((bin = static_cast<VertexBin*>(nextMC()))) {
 186     bin->AddOutput(fSumList, fCentAxis);
 187   }
 188   if (fWeights) {
 189     TList* wList = new TList();
 190     wList->SetName("FlowWeights");
 191     fWeights->Init(wList);
 192     fSumList->Add(wList);
 193   }
 194
 195 }
 196 //_____________________________________________________________________
 197 Bool_t AliForwardMCFlowTaskQC::Analyze()
 198 {
 199   //
 200   // Load FMD and SPD MC objects from aod tree and call the cumulants
 201   // calculation for the correct vertexbin
 202   //
 203   if (!AliForwardFlowTaskQC::Analyze()) return kFALSE;
 204
 205   // Run analysis on trackrefs from FMD and SPD
 206   AliAODForwardMult* aodfmult = static_cast<AliAODForwardMult*>(fAOD->FindListObject("ForwardMC"));
 207   AliAODCentralMult* aodcmult = static_cast<AliAODCentralMult*>(fAOD->FindListObject("CentralClustersMC"));
 208   Int_t vtx = fVtxAxis->FindBin(fVtx)-1;
 209
 210   // if objects are present, get histograms
 211   if (aodfmult) {
 212     TH2D& fmdTRdNdetadphi = aodfmult->GetHistogram();
 213     if ((fFlowFlags & kStdQC)) {
 214       FillVtxBinList(fBinsForwardTR, fmdTRdNdetadphi, vtx);
 215     } else if ((fFlowFlags & kEtaGap)) {
 216       FillVtxBinListEtaGap(fBinsForwardTR, fmdTRdNdetadphi, fmdTRdNdetadphi, vtx/*, kDoVtxCut*/);
 217     }
 218     if (aodcmult) {
 219       TH2D& spdTRdNdetadphi = aodcmult->GetHistogram();
 220       if ((fFlowFlags & kStdQC)) {
 221         FillVtxBinList(fBinsCentralTR, spdTRdNdetadphi, vtx);
 222       } else if ((fFlowFlags & kEtaGap)) {
 223         FillVtxBinListEtaGap(fBinsCentralTR, fmdTRdNdetadphi, spdTRdNdetadphi, vtx/*, kDoVtxCut*/);
 224       } else if ((fFlowFlags & k3Cor)) {
 225         FillVtxBinList3Cor(fBinsForwardTR, spdTRdNdetadphi, fmdTRdNdetadphi, vtx);
 226       }
 227     }
 228   }
 229   // Run analysis on MC branch
 230   if (!FillMCHist()) return kFALSE;
 231
 232   if ((fFlowFlags & kStdQC)) {
 233     FillVtxBinList(fBinsMC, fHistdNdedpMC, vtx, kMC);
 234   } else if ((fFlowFlags & kEtaGap)) {
 235     FillVtxBinListEtaGap(fBinsMC, fHistdNdedpMC, fHistdNdedpMC, vtx, kMC);
 236   } else if ((fFlowFlags & k3Cor)) {
 237     FillVtxBinList3Cor(fBinsMC, fHistdNdedpMC, fHistdNdedpMC, vtx);
 238   }
 239
 240   // Mult correlation diagnostics
 241   if (aodfmult && aodcmult) {
 242     AliAODForwardMult* fmult = static_cast<AliAODForwardMult*>(fAOD->FindListObject("Forward"));
 243     AliAODCentralMult* cmult = static_cast<AliAODCentralMult*>(fAOD->FindListObject("CentralClusters"));
 244     const TH2D& fhist = fmult->GetHistogram();
 245     const TH2D& chist = cmult->GetHistogram();
 246
 247     Double_t totForward = fhist.Integral(1, fhist.GetNbinsX(), 1, fhist.GetNbinsY());
 248     Double_t totSPD = chist.Integral(1, chist.GetNbinsX(), 1, chist.GetNbinsY());
 249     Double_t totMC = fHistdNdedpMC.Integral(1, fHistdNdedpMC.GetNbinsX(), 1, fHistdNdedpMC.GetNbinsY());
 250     fHistFMDMCCorr->Fill(totForward, totMC);
 251     fHistSPDMCCorr->Fill(totSPD, totMC);
 252   }
 253
 254   return kTRUE;
 255 }
 256 //_____________________________________________________________________
 257 void AliForwardMCFlowTaskQC::Finalize()
 258 {
 259   //
 260   // Finalize command, called by Terminate()
 261   //
 262   AliForwardFlowTaskQC::Finalize();
 263
 264   EndVtxBinList(fBinsForwardTR);
 265   EndVtxBinList(fBinsCentralTR);
 266   EndVtxBinList(fBinsMC);
 267
 268   return;
 269 }
 270 //_____________________________________________________________________
 271 Bool_t AliForwardMCFlowTaskQC::CheckEvent(const AliAODForwardMult* aodfm)
 272 {
 273   //
 274   //  Function to check that an AOD event meets the cuts
 275   //
 276   //  Parameters:
 277   //   AliAODForwardMult: forward mult object with trigger and vertex info
 278   //
 279   //  Return: false if there is no trigger or if the centrality or vertex
 280   //  is out of range. Otherwise true.
 281   //
 282   fAODMCHeader = static_cast<AliAODMCHeader*>(fAOD->FindListObject(AliAODMCHeader::StdBranchName()));
 283   return AliForwardFlowTaskQC::CheckEvent(aodfm);
 284 }
 285 //_____________________________________________________________________
 286 Bool_t AliForwardMCFlowTaskQC::CheckTrigger(const AliAODForwardMult* aodfm) const
 287 {
 288   //
 289   // Function to look for a trigger string in the event.
 290   //
 291   // Parameters:
 292   //  AliAODForwardMult: forward mult object with trigger and vertex info
 293   //
 294   // Returns true if B trigger is present - for some reason this is the one we use in MC
 295   //
 296   if (aodfm) return aodfm->IsTriggerBits(AliAODForwardMult::kB);
 297   else return (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
 298                  ->IsEventSelected() & AliVEvent::kMB);
 299
 300 }
 301 // _____________________________________________________________________
 302 Bool_t AliForwardMCFlowTaskQC::GetCentrality(const AliAODForwardMult* aodfm)
 303 {
 304   //
 305   // Function to use centrality parametrization from impact parameter
 306   // if flag is not set call AliForwardFlowTaskQC::GetCentrality
 307   //
 308   // Parameters:
 309   //  AliAODForwardMult: forward mult object with trigger and vertex info
 310   //
 311   // Returns true when centrality is set.
 312   //
 313   AliGenEventHeaderTunedPbPb* header =
 314     dynamic_cast<AliGenEventHeaderTunedPbPb*>(fAODMCHeader->GetCocktailHeader(0));
 315   if (header) fCent = header->GetCentrality();
 316   else if (fUseImpactPar) fCent = GetCentFromB();
 317   else return AliForwardFlowTaskQC::GetCentrality(aodfm);
 318
 319   if (fCentAxis->GetXmin() > fCent || fCent >= fCentAxis->GetXmax()) {
 320     fHistEventSel->Fill(kInvCent);
 321     return kFALSE;
 322   }
 323   if (fCent != 0) return kTRUE;
 324   else {
 325     fHistEventSel->Fill(kNoCent);
 326     return kFALSE;
 327   }
 328 }
 329 //_____________________________________________________________________
 330 Bool_t AliForwardMCFlowTaskQC::GetVertex(const AliAODForwardMult* aodfm)
 331 {
 332   //
 333   // Function to look for vertex determination in the event using the MC header.
 334   //
 335   // Parameters:
 336   //  AliAODForwardMult: Not used
 337   //
 338   // Returns true if vertex is determined
 339   //
 340   if (fUseMCVertex) {
 341     if (fAODMCHeader) {
 342       fVtx = fAODMCHeader->GetVtxZ();
 343       if (fVtx < fVtxAxis->GetXmin() || fVtx > fVtxAxis->GetXmax()) {
 344         fHistEventSel->Fill(kInvVtx);
 345         return kFALSE;
 346       }
 347       return kTRUE;
 348     } else {
 349       fHistEventSel->Fill(kNoVtx);
 350       return kFALSE;
 351     }
 352   } else
 353     return AliForwardFlowTaskQC::GetVertex(aodfm);
 354 }
 355 //_____________________________________________________________________
 356 Bool_t AliForwardMCFlowTaskQC::FillMCHist()
 357 {
 358   //
 359   // Loop over AliAODParticle branch and fill d^2N/detadphi-histograms.
 360   // Add flow if set to do so in AddTask function
 361   //
 362   fHistdNdedpMC.Reset();
 363   Double_t minEta = -3.75;
 364   Double_t maxEta = 5.;
 365
 366   //retreive MC particles from event
 367   TClonesArray* mcArray =
 368     static_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
 369   if(!mcArray){
 370     AliWarning("No MC array found in AOD. Try making it again.");
 371     return kFALSE;
 372   }
 373
 374   if (!fAODMCHeader) AliWarning("No MC header found.");
 375
 376   Int_t ntracks = mcArray->GetEntriesFast();
 377   Double_t rp = -1, b = -1;
 378   if (fAODMCHeader) {
 379     rp = fAODMCHeader->GetReactionPlaneAngle();
 380     b = fAODMCHeader->GetImpactParameter();
 381     if (fAODMCHeader->GetNCocktailHeaders() > 1) {
 382       ntracks = fAODMCHeader->GetCocktailHeader(0)->NProduced();
 383     }
 384   }
 385
 386   for (Int_t it = 0; it < ntracks; it++) { // Track loop
 387     Double_t weight = 1;
 388     AliAODMCParticle* particle = (AliAODMCParticle*)mcArray->At(it);
 389     if (!particle) {
 390       AliError(Form("Could not receive track %d", it));
 391       continue;
 392     }
 393     if (!particle->IsPrimary()) continue;
 394     if (particle->Charge() == 0) continue;
 395     Double_t pT = particle->Pt();
 396     Double_t eta = particle->Eta();
 397     Double_t phi = particle->Phi();
 398     if (eta >= minEta && eta < maxEta) {
 399       // Add flow if it is in the argument
 400       if (fUseFlowWeights && fWeights) {
 401         weight = fWeights->CalcWeight(eta, pT, phi, particle->PdgCode(), rp, b);
 402 //      Printf("%f", weight);
 403       }
 404       fHistdNdedpMC.Fill(eta, phi, weight);
 405     }
 406   }
 407   // Set underflow bins for acceptance checks
 408   Int_t sBin = fHistdNdedpMC.GetXaxis()->FindBin(minEta);
 409   Int_t eBin = fHistdNdedpMC.GetXaxis()->FindBin(maxEta);
 410   for ( ; sBin <= eBin; sBin++) {
 411     fHistdNdedpMC.SetBinContent(sBin, 0, 1);
 412   } // End of eta bin loop
 413
 414   return kTRUE;
 415 }
 416 //_____________________________________________________________________
 417 Double_t AliForwardMCFlowTaskQC::GetCentFromB() const
 418 {
 419   //
 420   // Get centrality from MC impact parameter.
 421   //
 422   Double_t cent = -1.;
 423   if (!fAODMCHeader) return cent;
 424   Double_t b = fAODMCHeader->GetImpactParameter();
 425   cent = fImpactParToCent->Eval(b);
 426
 427   return cent;
 428 }
 429 //_____________________________________________________________________
 430 //
 431 //
 432 // EOF