[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, fSPDCut, fEtaGap));
      if (isNUA) fFlowFlags ^= kNUAcorr;
      fBinsMC.Add(new VertexBin(vL, vH, fMaxMoment, "MC-FMD", fFlowFlags, -1, fEtaGap));
      if (isNUA) fFlowFlags ^= kNUAcorr;
    }
    // VZERO
    else if ((fFlowFlags & kVZERO)) {
      fBinsMC.Add(new VertexBin(vL, vH, fMaxMoment, "MC-VZERO", fFlowFlags, -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();
  Int_t minEta = -3.75;
  Int_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));
	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));
2b556440	147	}
2b556440	148	}
2b556440	149	}
	150	//_____________________________________________________________________
	151	void AliForwardMCFlowTaskQC::InitHists()
	152	{
	153	//
	154	// Initiate diagnostics hists and add to outputlist
	155	//
	156	AliForwardFlowTaskQC::InitHists();
	157
87f694ab AH	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);
68589651	173
87f694ab	174	TIter nextForwardTR(&fBinsForwardTR);
2b556440	175	VertexBin* bin = 0;
87f694ab AH	176	while ((bin = static_cast<VertexBin*>(nextForwardTR()))) {
87f694ab AH	177	bin->AddOutput(fSumList, fCentAxis);
2b556440	178	}
87f694ab AH	179	TIter nextCentralTR(&fBinsCentralTR);
	180	while ((bin = static_cast<VertexBin*>(nextCentralTR()))) {
	181	bin->AddOutput(fSumList, fCentAxis);
2b556440	182	}
	183	TIter nextMC(&fBinsMC);
	184	while ((bin = static_cast<VertexBin*>(nextMC()))) {
87f694ab	185	bin->AddOutput(fSumList, fCentAxis);
2b556440	186	}
2b556440	187
3112bd03	188	TList* wList = new TList();
	189	wList->SetName("FlowWeights");
	190	fWeights.Init(wList);
	191	fSumList->Add(wList);
	192
2b556440	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
87f694ab AH	204	AliAODForwardMult* aodfmult = static_cast<AliAODForwardMult*>(fAOD->FindListObject("ForwardMC"));
87f694ab AH	205	AliAODCentralMult* aodcmult = static_cast<AliAODCentralMult*>(fAOD->FindListObject("CentralClustersMC"));
d420e249	206	Int_t vtx = fVtxAxis->FindBin(fVtx)-1;
2b556440	207
4b5b52b7	208	// if objects are present, get histograms
d420e249	209	if (aodfmult) {
87f694ab AH	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/);
68589651	215	}
68589651	216	if (aodcmult) {
87f694ab AH	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);
68589651	224	}
68589651	225	}
d420e249	226	}
2b556440	227	// Run analysis on MC branch
87f694ab AH	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);
68589651	249	}
2b556440	250
	251	return kTRUE;
	252	}
	253	//_____________________________________________________________________
	254	void AliForwardMCFlowTaskQC::Finalize()
	255	{
	256	//
	257	// Finalize command, called by Terminate()
	258	//
	259	AliForwardFlowTaskQC::Finalize();
	260
87f694ab AH	261	EndVtxBinList(fBinsForwardTR);
87f694ab AH	262	EndVtxBinList(fBinsCentralTR);
4b5b52b7	263	EndVtxBinList(fBinsMC);
2b556440	264
4b5b52b7	265	return;
	266	}
	267	//_____________________________________________________________________
d420e249	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	//
87f694ab AH	278	if (aodfm) return aodfm->IsTriggerBits(AliAODForwardMult::kB);
	279	else return (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
	280	->IsEventSelected() & AliVEvent::kMB);
	281
d420e249	282	}
d420e249	283	// _____________________________________________________________________
4b5b52b7	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	//
87f694ab AH	295	// fCent = 2.5;
87f694ab AH	296	// return kTRUE;
4b5b52b7	297	if (fUseImpactPar) {
4b5b52b7	298	fCent = GetCentFromB();
87f694ab AH	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);
2b556440	310	}
2b556440	311	//_____________________________________________________________________
87f694ab	312	Bool_t AliForwardMCFlowTaskQC::GetVertex(const AliAODForwardMult* aodfm)
008eda2b	313	{
87f694ab AH	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;
008eda2b	335	}
87f694ab AH	336	} else
87f694ab AH	337	return AliForwardFlowTaskQC::GetVertex(aodfm);
008eda2b	338	}
008eda2b	339	//_____________________________________________________________________
87f694ab	340	Bool_t AliForwardMCFlowTaskQC::FillMCHist()
2b556440	341	{
	342	//
	343	// Loop over AliAODParticle branch and fill d^2N/detadphi-histograms.
	344	// Add flow if set to do so in AddTask function
68589651	345	//
87f694ab AH	346	fHistdNdedpMC.Reset();
	347	Int_t minEta = -3.75;
	348	Int_t maxEta = 5.;
2b556440	349
2b556440	350	//retreive MC particles from event
936b0a6c	351	TClonesArray* mcArray =
d420e249	352	static_cast<TClonesArray*>(fAOD->FindListObject(
d420e249	353	AliAODMCParticle::StdBranchName()));
2b556440	354	if(!mcArray){
d420e249	355	AliWarning("No MC array found in AOD. Try making it again.");
2b556440	356	return kFALSE;
	357	}
	358
936b0a6c	359	AliAODMCHeader* header =
87f694ab	360	static_cast<AliAODMCHeader*>(fAOD->FindListObject(
d420e249	361	AliAODMCHeader::StdBranchName()));
936b0a6c	362	if (!header)
2b556440	363	AliWarning("No header file found.");
2b556440	364
2b556440	365	Int_t ntracks = mcArray->GetEntriesFast();
87f694ab AH	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	}
008eda2b	373	}
87f694ab	374
936b0a6c	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
87f694ab AH	383	for (Int_t it = 0; it < ntracks; it++) { // Track loop
87f694ab AH	384	Double_t weight = 1;
2b556440	385	AliAODMCParticle* particle = (AliAODMCParticle*)mcArray->At(it);
	386	if (!particle) {
	387	AliError(Form("Could not receive track %d", it));
	388	continue;
	389	}
008eda2b	390	if (!particle->IsPrimary()) continue;
2b556440	391	if (particle->Charge() == 0) continue;
68589651	392	//Double_t pT = particle->Pt();
2b556440	393	Double_t eta = particle->Eta();
2b556440	394	Double_t phi = particle->Phi();
87f694ab	395	if (eta >= minEta && eta < maxEta) {
2b556440	396	// Add flow if it is in the argument
87f694ab AH	397	// FLOW WEIGHTS DISABLED IN THE VERSION - COMING BACK SOON
87f694ab AH	398	// if (flowFlags != 0) {
68589651	399	// weight = fWeights.CalcWeight(eta, pT, phi, particle->PdgCode(),
	400	// rp, fCent, fAddType, fAddOrder,
	401	// flowFlags) + 1;
87f694ab AH	402	// weight = fWeights.CalcWeight(eta, pT, phi, particle->PdgCode(),
87f694ab AH	403	// rp, b);
68589651	404	// Printf("%f", weight);
87f694ab AH	405	// }
87f694ab AH	406	fHistdNdedpMC.Fill(eta, phi, weight);
2b556440	407	}
2b556440	408	}
87f694ab AH	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
2b556440	415
	416	return kTRUE;
	417	}
	418	//_____________________________________________________________________
2b556440	419	Double_t AliForwardMCFlowTaskQC::GetCentFromB() const
	420	{
	421	//
	422	// Get centrality from MC impact parameter.
2b556440	423	//
2b556440	424	Double_t cent = -1.;
936b0a6c	425	AliAODMCHeader* header =
	426	static_cast<AliAODMCHeader*>(fAOD->FindListObject(AliAODMCHeader::
	427	StdBranchName()));
2b556440	428	if (!header) return cent;
9c686cb4	429	Double_t b = header->GetImpactParameter();
2b556440	430	cent = fImpactParToCent->Eval(b);
	431
	432	return cent;
	433	}
	434	//_____________________________________________________________________
	435	//
	436	//
	437	// EOF