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
  22 ClassImp(AliForwardMCFlowTaskQC)
  23 #if 0
  24 ;
  25 #endif
  26 //_____________________________________________________________________
  27 AliForwardMCFlowTaskQC::AliForwardMCFlowTaskQC()
  28   : AliForwardFlowTaskQC(),
  29     fBinsForwardTR(),      // List of FMDTR analysis objects
  30     fBinsCentralTR(),      // List of SPDTR analysis objects
  31     fBinsMC(),             // List of MC particle analysis objects
  32     fHistdNdedpMC(),       // MC particle d^2N/detadphi histogram
  33     fHistFMDMCCorr(),      // FMD MC correlation
  34     fHistSPDMCCorr(),      // SPD MC correlation
  35     fWeights(),            // Flow weights
  36     fImpactParToCent(),    // Impact parameter to centrality graph
  37     fUseImpactPar(0),      // Use impact par for centrality
  38     fUseMCVertex(0),       // Get vertex from MC header?
  39     fAddFlow(0),           // Add flow to MC particles
  40     fAddType(0),           // Add type of flow to MC particles
  41     fAddOrder(0)           // Add order of 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     fHistdNdedpMC(),       // MC particles d^2N/detadphi histogram
  53     fHistFMDMCCorr(),      // FMD MC correlation
  54     fHistSPDMCCorr(),      // SPD MC correlation
  55     fWeights(),            // Flow weights
  56     fImpactParToCent(),    // Impact parameter to centrality graph
  57     fUseImpactPar(0),      // Use impact par for centrality
  58     fUseMCVertex(0),       // Get vertex from MC header?
  59     fAddFlow(0),           // Add flow to MC particles
  60     fAddType(0),           // Add type of flow to MC particles
  61     fAddOrder(0)           // Add order of 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     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     fAddFlow(o.fAddFlow),                  // Add flow to MC particles
  95     fAddType(o.fAddType),                  // Add type of flow to MC particles
  96     fAddOrder(o.fAddOrder)                 // Add order of flow to MC particles
  97 {
  98   //
  99   // Copy Constructor
 100   //
 101 }
 102 //_____________________________________________________________________
 103 AliForwardMCFlowTaskQC&
 104 AliForwardMCFlowTaskQC::operator=(const AliForwardMCFlowTaskQC& o)
 105 {
 106   //
 107   // Assignment operator
 108   // Parameters:
 109   //    o Object to copy from
 110   //
 111   if (&o == this) return *this;
 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   fAddFlow         = o.fAddFlow;
 120   fAddType         = o.fAddType;
 121   fAddOrder        = o.fAddOrder;
 122   return *this;
 123 }
 124 //_____________________________________________________________________
 125 void AliForwardMCFlowTaskQC::InitVertexBins()
 126 {
 127   //
 128   // Initiate VertexBin lists
 129   //
 130   AliForwardFlowTaskQC::InitVertexBins();
 131
 132   Bool_t isNUA = (fFlowFlags & kNUAcorr);
 133   for (Int_t v = 1; v <= fVtxAxis->GetNbins(); v++) {
 134     Int_t vL = Int_t(fVtxAxis->GetBinLowEdge(v));
 135     Int_t vH = Int_t(fVtxAxis->GetBinUpEdge(v));
 136     // FMD
 137     if ((fFlowFlags & kFMD)) {
 138       fBinsForwardTR.Add(new VertexBin(vL, vH, fMaxMoment, "FMDTR", fFlowFlags, fFMDCut, fEtaGap));
 139       if (!(fFlowFlags & k3Cor)) fBinsCentralTR.Add(new VertexBin(vL, vH, fMaxMoment, "SPDTR", fFlowFlags, fSPDCut, fEtaGap));
 140       if (isNUA) fFlowFlags ^= kNUAcorr;
 141       fBinsMC.Add(new VertexBin(vL, vH, fMaxMoment, "MC-FMD", fFlowFlags, -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, -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
 188   TList* wList = new TList();
 189   wList->SetName("FlowWeights");
 190   fWeights.Init(wList);
 191   fSumList->Add(wList);
 192
 193 }
 194 //_____________________________________________________________________
 195 Bool_t AliForwardMCFlowTaskQC::Analyze()
 196 {
 197   //
 198   // Load FMD and SPD MC objects from aod tree and call the cumulants
 199   // calculation for the correct vertexbin
 200   //
 201   if (!AliForwardFlowTaskQC::Analyze()) return kFALSE;
 202
 203   // Run analysis on trackrefs from FMD and SPD
 204   AliAODForwardMult* aodfmult = static_cast<AliAODForwardMult*>(fAOD->FindListObject("ForwardMC"));
 205   AliAODCentralMult* aodcmult = static_cast<AliAODCentralMult*>(fAOD->FindListObject("CentralClustersMC"));
 206   Int_t vtx = fVtxAxis->FindBin(fVtx)-1;
 207
 208   // if objects are present, get histograms
 209   if (aodfmult) {
 210     TH2D& fmdTRdNdetadphi = aodfmult->GetHistogram();
 211     if ((fFlowFlags & kStdQC)) {
 212       FillVtxBinList(fBinsForwardTR, fmdTRdNdetadphi, vtx);
 213     } else if ((fFlowFlags & kEtaGap)) {
 214       FillVtxBinListEtaGap(fBinsForwardTR, fmdTRdNdetadphi, fmdTRdNdetadphi, vtx/*, kDoVtxCut*/);
 215     }
 216     if (aodcmult) {
 217       TH2D& spdTRdNdetadphi = aodcmult->GetHistogram();
 218       if ((fFlowFlags & kStdQC)) {
 219         FillVtxBinList(fBinsCentralTR, spdTRdNdetadphi, vtx);
 220       } else if ((fFlowFlags & kEtaGap)) {
 221         FillVtxBinListEtaGap(fBinsCentralTR, fmdTRdNdetadphi, spdTRdNdetadphi, vtx/*, kDoVtxCut*/);
 222       } else if ((fFlowFlags & k3Cor)) {
 223         FillVtxBinList3Cor(fBinsForwardTR, spdTRdNdetadphi, fmdTRdNdetadphi, vtx);
 224       }
 225     }
 226   }
 227   // Run analysis on MC branch
 228   if (!FillMCHist()) return kFALSE;
 229   if ((fFlowFlags & kStdQC)) {
 230     FillVtxBinList(fBinsMC, fHistdNdedpMC, vtx);
 231   } else if ((fFlowFlags & kEtaGap)) {
 232     FillVtxBinListEtaGap(fBinsMC, fHistdNdedpMC, fHistdNdedpMC, vtx);
 233   } else if ((fFlowFlags & k3Cor)) {
 234     FillVtxBinList3Cor(fBinsMC, fHistdNdedpMC, fHistdNdedpMC, vtx);
 235   }
 236
 237   // Mult correlation diagnostics
 238   if (aodfmult && aodcmult) {
 239     AliAODForwardMult* fmult = static_cast<AliAODForwardMult*>(fAOD->FindListObject("Forward"));
 240     AliAODCentralMult* cmult = static_cast<AliAODCentralMult*>(fAOD->FindListObject("CentralClusters"));
 241     const TH2D& fhist = fmult->GetHistogram();
 242     const TH2D& chist = cmult->GetHistogram();
 243
 244     Double_t totForward = fhist.Integral(1, fhist.GetNbinsX(), 1, fhist.GetNbinsY());
 245     Double_t totSPD = chist.Integral(1, chist.GetNbinsX(), 1, chist.GetNbinsY());
 246     Double_t totMC = fHistdNdedpMC.Integral(1, fHistdNdedpMC.GetNbinsX(), 1, fHistdNdedpMC.GetNbinsY());
 247     fHistFMDMCCorr->Fill(totForward, totMC);
 248     fHistSPDMCCorr->Fill(totSPD, totMC);
 249   }
 250
 251   return kTRUE;
 252 }
 253 //_____________________________________________________________________
 254 void AliForwardMCFlowTaskQC::Finalize()
 255 {
 256   //
 257   // Finalize command, called by Terminate()
 258   //
 259   AliForwardFlowTaskQC::Finalize();
 260
 261   EndVtxBinList(fBinsForwardTR);
 262   EndVtxBinList(fBinsCentralTR);
 263   EndVtxBinList(fBinsMC);
 264
 265   return;
 266 }
 267 //_____________________________________________________________________
 268 Bool_t AliForwardMCFlowTaskQC::CheckTrigger(const AliAODForwardMult* aodfm) const
 269 {
 270   //
 271   // Function to look for a trigger string in the event.
 272   //
 273   // Parameters:
 274   //  AliAODForwardMult: forward mult object with trigger and vertex info
 275   //
 276   // Returns true if B trigger is present - for some reason this is the one we use in MC
 277   //
 278   if (aodfm) return aodfm->IsTriggerBits(AliAODForwardMult::kB);
 279   else return (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
 280                  ->IsEventSelected() & AliVEvent::kMB);
 281
 282 }
 283 // _____________________________________________________________________
 284 Bool_t AliForwardMCFlowTaskQC::GetCentrality(const AliAODForwardMult* aodfm)
 285 {
 286   //
 287   // Function to use centrality parametrization from impact parameter
 288   // if flag is not set call AliForwardFlowTaskQC::GetCentrality
 289   //
 290   // Parameters:
 291   //  AliAODForwardMult: forward mult object with trigger and vertex info
 292   //
 293   // Returns true when centrality is set.
 294   //
 295 //  fCent = 2.5;
 296 //  return kTRUE;
 297   if (fUseImpactPar) {
 298     fCent = GetCentFromB();
 299     if (fCentAxis->GetXmin() > fCent || fCent >= fCentAxis->GetXmax()) {
 300       fHistEventSel->Fill(kInvCent);
 301       return kFALSE;
 302     }
 303     if (fCent != 0) return kTRUE;
 304     else {
 305       fHistEventSel->Fill(kNoCent);
 306       return kFALSE;
 307     }
 308   } else
 309     return AliForwardFlowTaskQC::GetCentrality(aodfm);
 310 }
 311 //_____________________________________________________________________
 312 Bool_t AliForwardMCFlowTaskQC::GetVertex(const AliAODForwardMult* aodfm)
 313 {
 314   //
 315   // Function to look for vertex determination in the event using the MC header.
 316   //
 317   // Parameters:
 318   //  AliAODForwardMult: Not used
 319   //
 320   // Returns true if vertex is determined
 321   //
 322   if (fUseMCVertex) {
 323     AliAODMCHeader* header =
 324       static_cast<AliAODMCHeader*>(fAOD->FindListObject(AliAODMCHeader::StdBranchName()));
 325     if (header) {
 326       fVtx = header->GetVtxZ();
 327       if (fVtx < fVtxAxis->GetXmin() || fVtx > fVtxAxis->GetXmax()) {
 328         fHistEventSel->Fill(kInvVtx);
 329         return kFALSE;
 330       }
 331       return kTRUE;
 332     } else {
 333       fHistEventSel->Fill(kNoVtx);
 334       return kFALSE;
 335     }
 336   } else
 337     return AliForwardFlowTaskQC::GetVertex(aodfm);
 338 }
 339 //_____________________________________________________________________
 340 Bool_t AliForwardMCFlowTaskQC::FillMCHist()
 341 {
 342   //
 343   // Loop over AliAODParticle branch and fill d^2N/detadphi-histograms.
 344   // Add flow if set to do so in AddTask function
 345   //
 346   fHistdNdedpMC.Reset();
 347   Int_t minEta = -3.75;
 348   Int_t maxEta = 5.;
 349
 350   //retreive MC particles from event
 351   TClonesArray* mcArray =
 352     static_cast<TClonesArray*>(fAOD->FindListObject(
 353                                AliAODMCParticle::StdBranchName()));
 354   if(!mcArray){
 355     AliWarning("No MC array found in AOD. Try making it again.");
 356     return kFALSE;
 357   }
 358
 359   AliAODMCHeader* header =
 360     static_cast<AliAODMCHeader*>(fAOD->FindListObject(
 361                                   AliAODMCHeader::StdBranchName()));
 362   if (!header)
 363     AliWarning("No header file found.");
 364
 365   Int_t ntracks = mcArray->GetEntriesFast();
 366 //  Double_t rp = -1, b = -1;
 367   if (header) {
 368 //    rp = header->GetReactionPlaneAngle();
 369 //    b = header->GetImpactParameter();
 370     if (header->GetNCocktailHeaders() > 1) {
 371       ntracks = header->GetCocktailHeader(0)->NProduced();
 372     }
 373   }
 374
 375   UShort_t flowFlags = 0;
 376   if (fAddFlow.Length() > 1) {
 377     if (fAddFlow.Contains("pt"))   flowFlags |= AliForwardFlowWeights::kPt;
 378     if (fAddFlow.Contains("pid"))  flowFlags |= AliForwardFlowWeights::kPt;
 379     if (fAddFlow.Contains("eta"))  flowFlags |= AliForwardFlowWeights::kEta;
 380     if (fAddFlow.Contains("cent")) flowFlags |= AliForwardFlowWeights::kCent;
 381   }
 382
 383   for (Int_t it = 0; it < ntracks; it++) { // Track loop
 384     Double_t weight = 1;
 385     AliAODMCParticle* particle = (AliAODMCParticle*)mcArray->At(it);
 386     if (!particle) {
 387       AliError(Form("Could not receive track %d", it));
 388       continue;
 389     }
 390     if (!particle->IsPrimary()) continue;
 391     if (particle->Charge() == 0) continue;
 392     //Double_t pT = particle->Pt();
 393     Double_t eta = particle->Eta();
 394     Double_t phi = particle->Phi();
 395     if (eta >= minEta && eta < maxEta) {
 396       // Add flow if it is in the argument
 397       // FLOW WEIGHTS DISABLED IN THE VERSION - COMING BACK SOON
 398 //      if (flowFlags != 0) {
 399 //      weight = fWeights.CalcWeight(eta, pT, phi, particle->PdgCode(),
 400 //                                   rp, fCent, fAddType, fAddOrder,
 401 //                                   flowFlags) + 1;
 402 //        weight = fWeights.CalcWeight(eta, pT, phi, particle->PdgCode(),
 403 //                                    rp, b);
 404 //      Printf("%f", weight);
 405 //      }
 406       fHistdNdedpMC.Fill(eta, phi, weight);
 407     }
 408   }
 409   // Set underflow bins for acceptance checks
 410   Int_t sBin = fHistdNdedpMC.GetXaxis()->FindBin(minEta);
 411   Int_t eBin = fHistdNdedpMC.GetXaxis()->FindBin(maxEta);
 412   for ( ; sBin <= eBin; sBin++) {
 413     fHistdNdedpMC.SetBinContent(sBin, 0, 1);
 414   } // End of eta bin loop
 415
 416   return kTRUE;
 417 }
 418 //_____________________________________________________________________
 419 Double_t AliForwardMCFlowTaskQC::GetCentFromB() const
 420 {
 421   //
 422   // Get centrality from MC impact parameter.
 423   //
 424   Double_t cent = -1.;
 425   AliAODMCHeader* header =
 426     static_cast<AliAODMCHeader*>(fAOD->FindListObject(AliAODMCHeader::
 427                                                       StdBranchName()));
 428   if (!header) return cent;
 429   Double_t b = header->GetImpactParameter();
 430   cent = fImpactParToCent->Eval(b);
 431
 432   return cent;
 433 }
 434 //_____________________________________________________________________
 435 //
 436 //
 437 // EOF