[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / AliForwardMCFlowTaskQC.cxx

//
// Calculate flow on MC data in the forward and central regions using the Q cumulants method.
//
// Inputs:
//  - AliAODEvent
//
// Outputs:
//  - AnalysisResults.root
//
#include "AliForwardMCFlowTaskQC.h"
#include "AliAODMCParticle.h"
#include "AliAODMCHeader.h"
#include "TGraph.h"
#include "TF1.h"
#include "AliAODEvent.h"
#include "AliAODForwardMult.h"
#include "AliAODCentralMult.h"
#include "AliGenEventHeader.h"
#include "AliAnalysisManager.h"
#include "AliInputEventHandler.h"

ClassImp(AliForwardMCFlowTaskQC)
#if 0
;
#endif
//_____________________________________________________________________
AliForwardMCFlowTaskQC::AliForwardMCFlowTaskQC() 
  : AliForwardFlowTaskQC(), 
    fBinsForwardTR(),      // List of FMDTR analysis objects
    fBinsCentralTR(),      // List of SPDTR analysis objects
    fBinsMC(),             // List of MC particle analysis objects
    fHistdNdedpMC(),       // MC particle d^2N/detadphi histogram
    fHistFMDMCCorr(),      // FMD MC correlation
    fHistSPDMCCorr(),      // SPD MC correlation
    fWeights(),            // Flow weights
    fImpactParToCent(),    // Impact parameter to centrality graph
    fUseImpactPar(0),      // Use impact par for centrality
    fUseMCVertex(0),       // Get vertex from MC header?
    fAddFlow(0),           // Add flow to MC particles
    fAddType(0),           // Add type of flow to MC particles
    fAddOrder(0)           // Add order of flow to MC particles        
{} 
  //
  // Default Constructor
  //
//_____________________________________________________________________
AliForwardMCFlowTaskQC::AliForwardMCFlowTaskQC(const char* name) 
  : AliForwardFlowTaskQC(name),
    fBinsForwardTR(),      // List of FMDTR analysis objects
    fBinsCentralTR(),      // List of SPDTR analysis objects
    fBinsMC(),             // List of MC particles analysis objects
    fHistdNdedpMC(),       // MC particles d^2N/detadphi histogram
    fHistFMDMCCorr(),      // FMD MC correlation
    fHistSPDMCCorr(),      // SPD MC correlation
    fWeights(),            // Flow weights
    fImpactParToCent(),    // Impact parameter to centrality graph
    fUseImpactPar(0),      // Use impact par for centrality
    fUseMCVertex(0),       // Get vertex from MC header?
    fAddFlow(0),           // Add flow to MC particles
    fAddType(0),           // Add type of flow to MC particles
    fAddOrder(0)           // Add order of flow to MC particles        
{ 
  // 
  // Constructor
  // Parameters:
  //  name: Name of task
  //
  
  //  Double_t impactParam[] = {0.,1.75,4.225,5.965,7.765,9.215,10.46,
  //                            11.565,12.575,13.515,16.679};
  //  Double_t centrality[] = {0.,2.5,7.5,15,25,35,45,55,65,75,90};
  Double_t impactParam[] = { 0.00,  3.72,  5.23,  7.31,  8.88, 10.20, 
			    11.38, 12.47, 13.50, 14.51, 16.679};
  Double_t centrality[]  = { 0.,    5.,   10.,   20.,   30.,   40., 
			    50.,   60.,   70.,   80.,  100.};

  Int_t nPoints = sizeof(impactParam)/sizeof(Double_t);
  fImpactParToCent = new TGraph(nPoints, impactParam, centrality);

}
//_____________________________________________________________________
AliForwardMCFlowTaskQC::AliForwardMCFlowTaskQC(const AliForwardMCFlowTaskQC& o) 
  : AliForwardFlowTaskQC(o), 
    fBinsForwardTR(),                      // List of FMDTR analysis objects
    fBinsCentralTR(),                      // List of SPDTR analysis objects
    fBinsMC(),                             // List of MC particles analysis objects
    fHistdNdedpMC(o.fHistdNdedpMC),        // MC particles d^2N/detadphi histogram
    fHistFMDMCCorr(o.fHistFMDMCCorr),      // FMD MC correlation
    fHistSPDMCCorr(o.fHistSPDMCCorr),      // SPD MC correlation
    fWeights(o.fWeights),                  // Flow weights
    fImpactParToCent(o.fImpactParToCent),  // Impact parameter to centrality
    fUseImpactPar(o.fUseImpactPar),        // Use impact par for centrality
    fUseMCVertex(o.fUseMCVertex),          // Get vertex from MC header?
    fAddFlow(o.fAddFlow),                  // Add flow to MC particles
    fAddType(o.fAddType),                  // Add type of flow to MC particles
    fAddOrder(o.fAddOrder)                 // Add order of flow to MC particles
{
  //
  // Copy Constructor
  //
} 
//_____________________________________________________________________
AliForwardMCFlowTaskQC&
AliForwardMCFlowTaskQC::operator=(const AliForwardMCFlowTaskQC& o)
{
  // 
  // Assignment operator
  // Parameters:
  //    o Object to copy from 
  //
  if (&o == this) return *this;
  fHistdNdedpMC    = o.fHistdNdedpMC;
  fHistFMDMCCorr   = o.fHistFMDMCCorr;
  fHistSPDMCCorr   = o.fHistSPDMCCorr;
  fWeights         = o.fWeights;
  fImpactParToCent = o.fImpactParToCent;
  fUseImpactPar    = o.fUseImpactPar;
  fUseMCVertex     = o.fUseMCVertex;
  fAddFlow         = o.fAddFlow;
  fAddType         = o.fAddType;
  fAddOrder        = o.fAddOrder;
  return *this;
}
//_____________________________________________________________________
void AliForwardMCFlowTaskQC::InitVertexBins()
{
  //
  // Initiate VertexBin lists
  //
  AliForwardFlowTaskQC::InitVertexBins();

  Bool_t isNUA = (fFlowFlags & kNUAcorr);
  for (Int_t v = 1; v <= fVtxAxis->GetNbins(); v++) {
    Int_t vL = Int_t(fVtxAxis->GetBinLowEdge(v));
    Int_t vH = Int_t(fVtxAxis->GetBinUpEdge(v));
    // FMD
    if ((fFlowFlags & kFMD)) {
      fBinsForwardTR.Add(new VertexBin(vL, vH, fMaxMoment, "FMDTR", fFlowFlags, fFMDCut, fEtaGap));
      if (!(fFlowFlags & k3Cor)) 
        fBinsCentralTR.Add(new VertexBin(vL, vH, fMaxMoment, "SPDTR", fFlowFlags|kSPD, fSPDCut, fEtaGap));
      if (isNUA) fFlowFlags ^= kNUAcorr;
      fBinsMC.Add(new VertexBin(vL, vH, fMaxMoment, "MC-FMD", fFlowFlags|kMC, -1, fEtaGap));
      if ((fFlowFlags & kStdQC)) 
	fBinsMC.Add(new VertexBin(vL, vH, fMaxMoment, "MC-SPD", fFlowFlags|kMC|kSPD, -1, fEtaGap));
      if (isNUA) fFlowFlags ^= kNUAcorr;
    }
    // VZERO
    else if ((fFlowFlags & kVZERO)) {
      fBinsMC.Add(new VertexBin(vL, vH, fMaxMoment, "MC-VZERO", fFlowFlags|kMC, -1, fEtaGap));
    }
  }
}
//_____________________________________________________________________
void AliForwardMCFlowTaskQC::InitHists()
{
  //
  // Initiate diagnostics hists and add to outputlist
  //
  AliForwardFlowTaskQC::InitHists();

  TString subDetName = ((fFlowFlags & kFMD) ? "FMD" : ((fFlowFlags & kVZERO) ? "VZERO" : "none"));
  fHistdNdedpMC = TH2D(Form("fdNdedpMC%s%s", subDetName.Data(), GetQCType(fFlowFlags)),
		   Form("fdNdedpMC%s%s", subDetName.Data(), GetQCType(fFlowFlags)),
		   240, -6., 6., 200, 0., TMath::TwoPi());
  
  fHistFMDMCCorr = new TH2D("hFMDMCCorr", "hFMDMCCorr", 200, 0., 15000., 200, 0, 20000);
  fHistSPDMCCorr = new TH2D("hSPDMCCorr", "hSPDMCCorr", 200, 0., 7500., 200, 0, 20000);
  TList* dList = (TList*)fSumList->FindObject("Diagnostics");
  if (!dList) {
    dList = new TList();
    dList->SetName("Diagnostics");
    fSumList->Add(dList);
  }
  dList->Add(fHistFMDMCCorr);
  dList->Add(fHistSPDMCCorr);

  TIter nextForwardTR(&fBinsForwardTR);
  VertexBin* bin = 0;
  while ((bin = static_cast<VertexBin*>(nextForwardTR()))) {
    bin->AddOutput(fSumList, fCentAxis);
  }
  TIter nextCentralTR(&fBinsCentralTR);
  while ((bin = static_cast<VertexBin*>(nextCentralTR()))) {
    bin->AddOutput(fSumList, fCentAxis);
  }
  TIter nextMC(&fBinsMC);
  while ((bin = static_cast<VertexBin*>(nextMC()))) {
    bin->AddOutput(fSumList, fCentAxis);
  }

  TList* wList = new TList();
  wList->SetName("FlowWeights");
  fWeights.Init(wList);
  fSumList->Add(wList);

}
//_____________________________________________________________________
Bool_t AliForwardMCFlowTaskQC::Analyze() 
{
  // 
  // Load FMD and SPD MC objects from aod tree and call the cumulants 
  // calculation for the correct vertexbin
  //
  if (!AliForwardFlowTaskQC::Analyze()) return kFALSE;

  // Run analysis on trackrefs from FMD and SPD
  AliAODForwardMult* aodfmult = static_cast<AliAODForwardMult*>(fAOD->FindListObject("ForwardMC"));
  AliAODCentralMult* aodcmult = static_cast<AliAODCentralMult*>(fAOD->FindListObject("CentralClustersMC"));
  Int_t vtx = fVtxAxis->FindBin(fVtx)-1;
  
  // if objects are present, get histograms
  if (aodfmult) {
    TH2D& fmdTRdNdetadphi = aodfmult->GetHistogram();
    if ((fFlowFlags & kStdQC)) {
      FillVtxBinList(fBinsForwardTR, fmdTRdNdetadphi, vtx);
    } else if ((fFlowFlags & kEtaGap)) {
      FillVtxBinListEtaGap(fBinsForwardTR, fmdTRdNdetadphi, fmdTRdNdetadphi, vtx/*, kDoVtxCut*/);
    }
    if (aodcmult) {
      TH2D& spdTRdNdetadphi = aodcmult->GetHistogram();
      if ((fFlowFlags & kStdQC)) {
	FillVtxBinList(fBinsCentralTR, spdTRdNdetadphi, vtx);
      } else if ((fFlowFlags & kEtaGap)) {
	FillVtxBinListEtaGap(fBinsCentralTR, fmdTRdNdetadphi, spdTRdNdetadphi, vtx/*, kDoVtxCut*/);
      } else if ((fFlowFlags & k3Cor)) {
	FillVtxBinList3Cor(fBinsForwardTR, spdTRdNdetadphi, fmdTRdNdetadphi, vtx);
      }
    }
  }
  // Run analysis on MC branch
  if (!FillMCHist()) return kFALSE;
  if ((fFlowFlags & kStdQC)) {
    FillVtxBinList(fBinsMC, fHistdNdedpMC, vtx);
  } else if ((fFlowFlags & kEtaGap)) {
    FillVtxBinListEtaGap(fBinsMC, fHistdNdedpMC, fHistdNdedpMC, vtx);
  } else if ((fFlowFlags & k3Cor)) {
    FillVtxBinList3Cor(fBinsMC, fHistdNdedpMC, fHistdNdedpMC, vtx);
  }

  // Mult correlation diagnostics
  if (aodfmult && aodcmult) {
    AliAODForwardMult* fmult = static_cast<AliAODForwardMult*>(fAOD->FindListObject("Forward"));
    AliAODCentralMult* cmult = static_cast<AliAODCentralMult*>(fAOD->FindListObject("CentralClusters"));
    const TH2D& fhist = fmult->GetHistogram();
    const TH2D& chist = cmult->GetHistogram();

    Double_t totForward = fhist.Integral(1, fhist.GetNbinsX(), 1, fhist.GetNbinsY());
    Double_t totSPD = chist.Integral(1, chist.GetNbinsX(), 1, chist.GetNbinsY());
    Double_t totMC = fHistdNdedpMC.Integral(1, fHistdNdedpMC.GetNbinsX(), 1, fHistdNdedpMC.GetNbinsY());
    fHistFMDMCCorr->Fill(totForward, totMC);
    fHistSPDMCCorr->Fill(totSPD, totMC);
  }

  return kTRUE;
}
//_____________________________________________________________________
void AliForwardMCFlowTaskQC::Finalize() 
{
  //
  // Finalize command, called by Terminate()
  //
  AliForwardFlowTaskQC::Finalize();

  EndVtxBinList(fBinsForwardTR);
  EndVtxBinList(fBinsCentralTR);
  EndVtxBinList(fBinsMC);

  return;
}
//_____________________________________________________________________
Bool_t AliForwardMCFlowTaskQC::CheckTrigger(const AliAODForwardMult* aodfm) const 
{
  //
  // Function to look for a trigger string in the event.
  //
  // Parameters: 
  //  AliAODForwardMult: forward mult object with trigger and vertex info
  //
  // Returns true if B trigger is present - for some reason this is the one we use in MC
  //
  if (aodfm) return aodfm->IsTriggerBits(AliAODForwardMult::kB);
  else return (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
                 ->IsEventSelected() & AliVEvent::kMB);

}
// _____________________________________________________________________
Bool_t AliForwardMCFlowTaskQC::GetCentrality(const AliAODForwardMult* aodfm)
{
  // 
  // Function to use centrality parametrization from impact parameter
  // if flag is not set call AliForwardFlowTaskQC::GetCentrality
  //
  // Parameters:
  //  AliAODForwardMult: forward mult object with trigger and vertex info
  //
  // Returns true when centrality is set.
  //
//  fCent = 2.5;
//  return kTRUE;
  if (fUseImpactPar) {
    fCent = GetCentFromB();
    if (fCentAxis->GetXmin() > fCent || fCent >= fCentAxis->GetXmax()) {
      fHistEventSel->Fill(kInvCent);
      return kFALSE;
    }
    if (fCent != 0) return kTRUE;
    else {
      fHistEventSel->Fill(kNoCent);
      return kFALSE;
    }
  } else 
    return AliForwardFlowTaskQC::GetCentrality(aodfm);
}
//_____________________________________________________________________
Bool_t AliForwardMCFlowTaskQC::GetVertex(const AliAODForwardMult* aodfm)
{
  //
  // Function to look for vertex determination in the event using the MC header.
  //
  // Parameters: 
  //  AliAODForwardMult: Not used
  //
  // Returns true if vertex is determined
  //
  if (fUseMCVertex) {
    AliAODMCHeader* header = 
      static_cast<AliAODMCHeader*>(fAOD->FindListObject(AliAODMCHeader::StdBranchName()));
    if (header) {
      fVtx = header->GetVtxZ();
      if (fVtx < fVtxAxis->GetXmin() || fVtx > fVtxAxis->GetXmax()) {
        fHistEventSel->Fill(kInvVtx);
	return kFALSE;
      }
      return kTRUE;
    } else {
      fHistEventSel->Fill(kNoVtx);
      return kFALSE;
    }
  } else 
    return AliForwardFlowTaskQC::GetVertex(aodfm);
}
//_____________________________________________________________________
Bool_t AliForwardMCFlowTaskQC::FillMCHist()  
{
  // 
  // Loop over AliAODParticle branch and fill d^2N/detadphi-histograms.
  // Add flow if set to do so in AddTask function
  //
  fHistdNdedpMC.Reset();
  Double_t minEta = -3.75;
  Double_t maxEta = 5.;

  //retreive MC particles from event
  TClonesArray* mcArray = 
    static_cast<TClonesArray*>(fAOD->FindListObject(
                               AliAODMCParticle::StdBranchName()));
  if(!mcArray){
    AliWarning("No MC array found in AOD. Try making it again.");
    return kFALSE;
  }

  AliAODMCHeader* header = 
    static_cast<AliAODMCHeader*>(fAOD->FindListObject(
                                  AliAODMCHeader::StdBranchName()));
  if (!header) 
    AliWarning("No header file found.");

  Int_t ntracks = mcArray->GetEntriesFast();
//  Double_t rp = -1, b = -1;
  if (header) {
//    rp = header->GetReactionPlaneAngle();
//    b = header->GetImpactParameter();
    if (header->GetNCocktailHeaders() > 1) {
      ntracks = header->GetCocktailHeader(0)->NProduced();
    }
  }
  
  UShort_t flowFlags = 0;
  if (fAddFlow.Length() > 1) {
    if (fAddFlow.Contains("pt"))   flowFlags |= AliForwardFlowWeights::kPt;
    if (fAddFlow.Contains("pid"))  flowFlags |= AliForwardFlowWeights::kPt;
    if (fAddFlow.Contains("eta"))  flowFlags |= AliForwardFlowWeights::kEta;
    if (fAddFlow.Contains("cent")) flowFlags |= AliForwardFlowWeights::kCent;
  }

  for (Int_t it = 0; it < ntracks; it++) { // Track loop
    Double_t weight = 1;
    AliAODMCParticle* particle = (AliAODMCParticle*)mcArray->At(it);
    if (!particle) {
      AliError(Form("Could not receive track %d", it));
      continue;
    }
    if (!particle->IsPrimary()) continue;
    if (particle->Charge() == 0) continue;
    //Double_t pT = particle->Pt();
    Double_t eta = particle->Eta();
    Double_t phi = particle->Phi();
    if (eta >= minEta && eta < maxEta) {
      // Add flow if it is in the argument
      // FLOW WEIGHTS DISABLED IN THE VERSION - COMING BACK SOON
//      if (flowFlags != 0) { 
//	weight = fWeights.CalcWeight(eta, pT, phi, particle->PdgCode(), 
//				     rp, fCent, fAddType, fAddOrder, 
//				     flowFlags) + 1;
//        weight = fWeights.CalcWeight(eta, pT, phi, particle->PdgCode(),
//                                    rp, b); 
//      Printf("%f", weight);
//      }
      fHistdNdedpMC.Fill(eta, phi, weight);
    }
  }
  // Set underflow bins for acceptance checks
  Int_t sBin = fHistdNdedpMC.GetXaxis()->FindBin(minEta);
  Int_t eBin = fHistdNdedpMC.GetXaxis()->FindBin(maxEta);
  for ( ; sBin <= eBin; sBin++) {
    fHistdNdedpMC.SetBinContent(sBin, 0, 1);
  } // End of eta bin loop

  return kTRUE;
}
//_____________________________________________________________________
Double_t AliForwardMCFlowTaskQC::GetCentFromB() const
{
  //
  // Get centrality from MC impact parameter.
  //
  Double_t cent = -1.;
  AliAODMCHeader* header = 
    static_cast<AliAODMCHeader*>(fAOD->FindListObject(AliAODMCHeader::
						      StdBranchName()));
  if (!header) return cent;
  Double_t b = header->GetImpactParameter();
  cent = fImpactParToCent->Eval(b);

  return cent;
}
//_____________________________________________________________________
//
//
// EOF
Commit	Line	Data
2b556440	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"
2b556440	15	#include "AliAODEvent.h"
	16	#include "AliAODForwardMult.h"
	17	#include "AliAODCentralMult.h"
008eda2b	18	#include "AliGenEventHeader.h"
87f694ab AH	19	#include "AliAnalysisManager.h"
87f694ab AH	20	#include "AliInputEventHandler.h"
2b556440	21
2b556440	22	ClassImp(AliForwardMCFlowTaskQC)
936b0a6c	23	#if 0
	24	;
	25	#endif
2b556440	26	//_____________________________________________________________________
936b0a6c	27	AliForwardMCFlowTaskQC::AliForwardMCFlowTaskQC()
936b0a6c	28	: AliForwardFlowTaskQC(),
87f694ab AH	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
936b0a6c	35	fWeights(), // Flow weights
	36	fImpactParToCent(), // Impact parameter to centrality graph
	37	fUseImpactPar(0), // Use impact par for centrality
87f694ab AH	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
936b0a6c	42	{}
2b556440	43	//
	44	// Default Constructor
	45	//
	46	//_____________________________________________________________________
936b0a6c	47	AliForwardMCFlowTaskQC::AliForwardMCFlowTaskQC(const char* name)
936b0a6c	48	: AliForwardFlowTaskQC(name),
87f694ab AH	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
936b0a6c	55	fWeights(), // Flow weights
	56	fImpactParToCent(), // Impact parameter to centrality graph
	57	fUseImpactPar(0), // Use impact par for centrality
87f694ab AH	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
2b556440	62	{
	63	//
	64	// Constructor
	65	// Parameters:
	66	// name: Name of task
	67	//
2b556440	68
936b0a6c	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.};
2b556440	76
	77	Int_t nPoints = sizeof(impactParam)/sizeof(Double_t);
	78	fImpactParToCent = new TGraph(nPoints, impactParam, centrality);
	79
	80	}
	81	//_____________________________________________________________________
936b0a6c	82	AliForwardMCFlowTaskQC::AliForwardMCFlowTaskQC(const AliForwardMCFlowTaskQC& o)
936b0a6c	83	: AliForwardFlowTaskQC(o),
87f694ab AH	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
936b0a6c	90	fWeights(o.fWeights), // Flow weights
	91	fImpactParToCent(o.fImpactParToCent), // Impact parameter to centrality
	92	fUseImpactPar(o.fUseImpactPar), // Use impact par for centrality
87f694ab AH	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
936b0a6c	97	{
2b556440	98	//
	99	// Copy Constructor
	100	//
936b0a6c	101	}
2b556440	102	//_____________________________________________________________________
	103	AliForwardMCFlowTaskQC&
	104	AliForwardMCFlowTaskQC::operator=(const AliForwardMCFlowTaskQC& o)
	105	{
	106	//
	107	// Assignment operator
	108	// Parameters:
	109	// o Object to copy from
	110	//
eeb6c967	111	if (&o == this) return *this;
87f694ab AH	112	fHistdNdedpMC = o.fHistdNdedpMC;
	113	fHistFMDMCCorr = o.fHistFMDMCCorr;
	114	fHistSPDMCCorr = o.fHistSPDMCCorr;
936b0a6c	115	fWeights = o.fWeights;
2b556440	116	fImpactParToCent = o.fImpactParToCent;
4b5b52b7	117	fUseImpactPar = o.fUseImpactPar;
87f694ab	118	fUseMCVertex = o.fUseMCVertex;
2b556440	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
87f694ab AH	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));
02738f97	139	if (!(fFlowFlags & k3Cor))
02738f97	140	fBinsCentralTR.Add(new VertexBin(vL, vH, fMaxMoment, "SPDTR", fFlowFlags\|kSPD, fSPDCut, fEtaGap));
87f694ab	141	if (isNUA) fFlowFlags ^= kNUAcorr;
02738f97	142	fBinsMC.Add(new VertexBin(vL, vH, fMaxMoment, "MC-FMD", fFlowFlags\|kMC, -1, fEtaGap));
	143	if ((fFlowFlags & kStdQC))
	144	fBinsMC.Add(new VertexBin(vL, vH, fMaxMoment, "MC-SPD", fFlowFlags\|kMC\|kSPD, -1, fEtaGap));
87f694ab AH	145	if (isNUA) fFlowFlags ^= kNUAcorr;
	146	}
	147	// VZERO
	148	else if ((fFlowFlags & kVZERO)) {
02738f97	149	fBinsMC.Add(new VertexBin(vL, vH, fMaxMoment, "MC-VZERO", fFlowFlags\|kMC, -1, fEtaGap));
2b556440	150	}
2b556440	151	}
2b556440	152	}
	153	//_____________________________________________________________________
	154	void AliForwardMCFlowTaskQC::InitHists()
	155	{
	156	//
	157	// Initiate diagnostics hists and add to outputlist
	158	//
	159	AliForwardFlowTaskQC::InitHists();
	160
87f694ab AH	161	TString subDetName = ((fFlowFlags & kFMD) ? "FMD" : ((fFlowFlags & kVZERO) ? "VZERO" : "none"));
	162	fHistdNdedpMC = TH2D(Form("fdNdedpMC%s%s", subDetName.Data(), GetQCType(fFlowFlags)),
	163	Form("fdNdedpMC%s%s", subDetName.Data(), GetQCType(fFlowFlags)),
	164	240, -6., 6., 200, 0., TMath::TwoPi());
	165
	166	fHistFMDMCCorr = new TH2D("hFMDMCCorr", "hFMDMCCorr", 200, 0., 15000., 200, 0, 20000);
	167	fHistSPDMCCorr = new TH2D("hSPDMCCorr", "hSPDMCCorr", 200, 0., 7500., 200, 0, 20000);
	168	TList* dList = (TList*)fSumList->FindObject("Diagnostics");
	169	if (!dList) {
	170	dList = new TList();
	171	dList->SetName("Diagnostics");
	172	fSumList->Add(dList);
	173	}
	174	dList->Add(fHistFMDMCCorr);
	175	dList->Add(fHistSPDMCCorr);
68589651	176
87f694ab	177	TIter nextForwardTR(&fBinsForwardTR);
2b556440	178	VertexBin* bin = 0;
87f694ab AH	179	while ((bin = static_cast<VertexBin*>(nextForwardTR()))) {
87f694ab AH	180	bin->AddOutput(fSumList, fCentAxis);
2b556440	181	}
87f694ab AH	182	TIter nextCentralTR(&fBinsCentralTR);
	183	while ((bin = static_cast<VertexBin*>(nextCentralTR()))) {
	184	bin->AddOutput(fSumList, fCentAxis);
2b556440	185	}
	186	TIter nextMC(&fBinsMC);
	187	while ((bin = static_cast<VertexBin*>(nextMC()))) {
87f694ab	188	bin->AddOutput(fSumList, fCentAxis);
2b556440	189	}
2b556440	190
3112bd03	191	TList* wList = new TList();
	192	wList->SetName("FlowWeights");
	193	fWeights.Init(wList);
	194	fSumList->Add(wList);
	195
2b556440	196	}
	197	//_____________________________________________________________________
	198	Bool_t AliForwardMCFlowTaskQC::Analyze()
	199	{
	200	//
	201	// Load FMD and SPD MC objects from aod tree and call the cumulants
	202	// calculation for the correct vertexbin
	203	//
	204	if (!AliForwardFlowTaskQC::Analyze()) return kFALSE;
	205
	206	// Run analysis on trackrefs from FMD and SPD
87f694ab AH	207	AliAODForwardMult* aodfmult = static_cast<AliAODForwardMult*>(fAOD->FindListObject("ForwardMC"));
87f694ab AH	208	AliAODCentralMult* aodcmult = static_cast<AliAODCentralMult*>(fAOD->FindListObject("CentralClustersMC"));
d420e249	209	Int_t vtx = fVtxAxis->FindBin(fVtx)-1;
2b556440	210
4b5b52b7	211	// if objects are present, get histograms
d420e249	212	if (aodfmult) {
87f694ab AH	213	TH2D& fmdTRdNdetadphi = aodfmult->GetHistogram();
	214	if ((fFlowFlags & kStdQC)) {
	215	FillVtxBinList(fBinsForwardTR, fmdTRdNdetadphi, vtx);
	216	} else if ((fFlowFlags & kEtaGap)) {
	217	FillVtxBinListEtaGap(fBinsForwardTR, fmdTRdNdetadphi, fmdTRdNdetadphi, vtx/, kDoVtxCut/);
68589651	218	}
68589651	219	if (aodcmult) {
87f694ab AH	220	TH2D& spdTRdNdetadphi = aodcmult->GetHistogram();
	221	if ((fFlowFlags & kStdQC)) {
	222	FillVtxBinList(fBinsCentralTR, spdTRdNdetadphi, vtx);
	223	} else if ((fFlowFlags & kEtaGap)) {
	224	FillVtxBinListEtaGap(fBinsCentralTR, fmdTRdNdetadphi, spdTRdNdetadphi, vtx/, kDoVtxCut/);
	225	} else if ((fFlowFlags & k3Cor)) {
	226	FillVtxBinList3Cor(fBinsForwardTR, spdTRdNdetadphi, fmdTRdNdetadphi, vtx);
68589651	227	}
68589651	228	}
d420e249	229	}
2b556440	230	// Run analysis on MC branch
87f694ab AH	231	if (!FillMCHist()) return kFALSE;
	232	if ((fFlowFlags & kStdQC)) {
	233	FillVtxBinList(fBinsMC, fHistdNdedpMC, vtx);
	234	} else if ((fFlowFlags & kEtaGap)) {
	235	FillVtxBinListEtaGap(fBinsMC, fHistdNdedpMC, fHistdNdedpMC, vtx);
	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);
68589651	252	}
2b556440	253
	254	return kTRUE;
	255	}
	256	//_____________________________________________________________________
	257	void AliForwardMCFlowTaskQC::Finalize()
	258	{
	259	//
	260	// Finalize command, called by Terminate()
	261	//
	262	AliForwardFlowTaskQC::Finalize();
	263
87f694ab AH	264	EndVtxBinList(fBinsForwardTR);
87f694ab AH	265	EndVtxBinList(fBinsCentralTR);
4b5b52b7	266	EndVtxBinList(fBinsMC);
2b556440	267
4b5b52b7	268	return;
	269	}
	270	//_____________________________________________________________________
d420e249	271	Bool_t AliForwardMCFlowTaskQC::CheckTrigger(const AliAODForwardMult* aodfm) const
	272	{
	273	//
	274	// Function to look for a trigger string in the event.
	275	//
	276	// Parameters:
	277	// AliAODForwardMult: forward mult object with trigger and vertex info
	278	//
	279	// Returns true if B trigger is present - for some reason this is the one we use in MC
	280	//
87f694ab AH	281	if (aodfm) return aodfm->IsTriggerBits(AliAODForwardMult::kB);
	282	else return (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
	283	->IsEventSelected() & AliVEvent::kMB);
	284
d420e249	285	}
d420e249	286	// _____________________________________________________________________
4b5b52b7	287	Bool_t AliForwardMCFlowTaskQC::GetCentrality(const AliAODForwardMult* aodfm)
	288	{
	289	//
	290	// Function to use centrality parametrization from impact parameter
	291	// if flag is not set call AliForwardFlowTaskQC::GetCentrality
	292	//
	293	// Parameters:
	294	// AliAODForwardMult: forward mult object with trigger and vertex info
	295	//
	296	// Returns true when centrality is set.
	297	//
87f694ab AH	298	// fCent = 2.5;
87f694ab AH	299	// return kTRUE;
4b5b52b7	300	if (fUseImpactPar) {
4b5b52b7	301	fCent = GetCentFromB();
87f694ab AH	302	if (fCentAxis->GetXmin() > fCent \|\| fCent >= fCentAxis->GetXmax()) {
	303	fHistEventSel->Fill(kInvCent);
	304	return kFALSE;
	305	}
	306	if (fCent != 0) return kTRUE;
	307	else {
	308	fHistEventSel->Fill(kNoCent);
	309	return kFALSE;
	310	}
	311	} else
	312	return AliForwardFlowTaskQC::GetCentrality(aodfm);
2b556440	313	}
2b556440	314	//_____________________________________________________________________
87f694ab	315	Bool_t AliForwardMCFlowTaskQC::GetVertex(const AliAODForwardMult* aodfm)
008eda2b	316	{
87f694ab AH	317	//
	318	// Function to look for vertex determination in the event using the MC header.
	319	//
	320	// Parameters:
	321	// AliAODForwardMult: Not used
	322	//
	323	// Returns true if vertex is determined
	324	//
	325	if (fUseMCVertex) {
	326	AliAODMCHeader* header =
	327	static_cast<AliAODMCHeader*>(fAOD->FindListObject(AliAODMCHeader::StdBranchName()));
	328	if (header) {
	329	fVtx = header->GetVtxZ();
	330	if (fVtx < fVtxAxis->GetXmin() \|\| fVtx > fVtxAxis->GetXmax()) {
	331	fHistEventSel->Fill(kInvVtx);
	332	return kFALSE;
	333	}
	334	return kTRUE;
	335	} else {
	336	fHistEventSel->Fill(kNoVtx);
	337	return kFALSE;
008eda2b	338	}
87f694ab AH	339	} else
87f694ab AH	340	return AliForwardFlowTaskQC::GetVertex(aodfm);
008eda2b	341	}
008eda2b	342	//_____________________________________________________________________
87f694ab	343	Bool_t AliForwardMCFlowTaskQC::FillMCHist()
2b556440	344	{
	345	//
	346	// Loop over AliAODParticle branch and fill d^2N/detadphi-histograms.
	347	// Add flow if set to do so in AddTask function
68589651	348	//
87f694ab	349	fHistdNdedpMC.Reset();
5bf83cfd	350	Double_t minEta = -3.75;
5bf83cfd	351	Double_t maxEta = 5.;
2b556440	352
2b556440	353	//retreive MC particles from event
936b0a6c	354	TClonesArray* mcArray =
d420e249	355	static_cast<TClonesArray*>(fAOD->FindListObject(
d420e249	356	AliAODMCParticle::StdBranchName()));
2b556440	357	if(!mcArray){
d420e249	358	AliWarning("No MC array found in AOD. Try making it again.");
2b556440	359	return kFALSE;
	360	}
	361
936b0a6c	362	AliAODMCHeader* header =
87f694ab	363	static_cast<AliAODMCHeader*>(fAOD->FindListObject(
d420e249	364	AliAODMCHeader::StdBranchName()));
936b0a6c	365	if (!header)
2b556440	366	AliWarning("No header file found.");
2b556440	367
2b556440	368	Int_t ntracks = mcArray->GetEntriesFast();
87f694ab AH	369	// Double_t rp = -1, b = -1;
	370	if (header) {
	371	// rp = header->GetReactionPlaneAngle();
	372	// b = header->GetImpactParameter();
	373	if (header->GetNCocktailHeaders() > 1) {
	374	ntracks = header->GetCocktailHeader(0)->NProduced();
	375	}
008eda2b	376	}
87f694ab	377
936b0a6c	378	UShort_t flowFlags = 0;
	379	if (fAddFlow.Length() > 1) {
	380	if (fAddFlow.Contains("pt")) flowFlags \|= AliForwardFlowWeights::kPt;
	381	if (fAddFlow.Contains("pid")) flowFlags \|= AliForwardFlowWeights::kPt;
	382	if (fAddFlow.Contains("eta")) flowFlags \|= AliForwardFlowWeights::kEta;
	383	if (fAddFlow.Contains("cent")) flowFlags \|= AliForwardFlowWeights::kCent;
	384	}
	385
87f694ab AH	386	for (Int_t it = 0; it < ntracks; it++) { // Track loop
87f694ab AH	387	Double_t weight = 1;
2b556440	388	AliAODMCParticle* particle = (AliAODMCParticle*)mcArray->At(it);
	389	if (!particle) {
	390	AliError(Form("Could not receive track %d", it));
	391	continue;
	392	}
008eda2b	393	if (!particle->IsPrimary()) continue;
2b556440	394	if (particle->Charge() == 0) continue;
68589651	395	//Double_t pT = particle->Pt();
2b556440	396	Double_t eta = particle->Eta();
2b556440	397	Double_t phi = particle->Phi();
87f694ab	398	if (eta >= minEta && eta < maxEta) {
2b556440	399	// Add flow if it is in the argument
87f694ab AH	400	// FLOW WEIGHTS DISABLED IN THE VERSION - COMING BACK SOON
87f694ab AH	401	// if (flowFlags != 0) {
68589651	402	// weight = fWeights.CalcWeight(eta, pT, phi, particle->PdgCode(),
	403	// rp, fCent, fAddType, fAddOrder,
	404	// flowFlags) + 1;
87f694ab AH	405	// weight = fWeights.CalcWeight(eta, pT, phi, particle->PdgCode(),
87f694ab AH	406	// rp, b);
68589651	407	// Printf("%f", weight);
87f694ab AH	408	// }
87f694ab AH	409	fHistdNdedpMC.Fill(eta, phi, weight);
2b556440	410	}
2b556440	411	}
87f694ab AH	412	// Set underflow bins for acceptance checks
	413	Int_t sBin = fHistdNdedpMC.GetXaxis()->FindBin(minEta);
	414	Int_t eBin = fHistdNdedpMC.GetXaxis()->FindBin(maxEta);
	415	for ( ; sBin <= eBin; sBin++) {
	416	fHistdNdedpMC.SetBinContent(sBin, 0, 1);
	417	} // End of eta bin loop
2b556440	418
	419	return kTRUE;
	420	}
	421	//_____________________________________________________________________
2b556440	422	Double_t AliForwardMCFlowTaskQC::GetCentFromB() const
	423	{
	424	//
	425	// Get centrality from MC impact parameter.
2b556440	426	//
2b556440	427	Double_t cent = -1.;
936b0a6c	428	AliAODMCHeader* header =
	429	static_cast<AliAODMCHeader*>(fAOD->FindListObject(AliAODMCHeader::
	430	StdBranchName()));
2b556440	431	if (!header) return cent;
9c686cb4	432	Double_t b = header->GetImpactParameter();
2b556440	433	cent = fImpactParToCent->Eval(b);
	434
	435	return cent;
	436	}
	437	//_____________________________________________________________________
	438	//
	439	//
	440	// EOF