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