[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
//
/**
 * @file   AliForwardMCFlowTaskQC.cxx
 * @author Christian Holm Christensen <cholm@master.hehi.nbi.dk>
 * @date   Thu Feb  7 01:09:19 2013
 * 
 * @brief  
 * 
 * 
 * @ingroup pwglf_forward_flow
 */
#include "AliForwardMCFlowTaskQC.h"
#include "AliAODMCParticle.h"
#include "AliAODMCHeader.h"
#include "TGraph.h"
#include "TF1.h"
#include "TProfile2D.h"
#include "AliAODEvent.h"
#include "AliAODForwardMult.h"
#include "AliAODCentralMult.h"
#include "AliGenEventHeader.h"

ClassImp(AliForwardMCFlowTaskQC)
#if 0
;
#endif
//_____________________________________________________________________
AliForwardMCFlowTaskQC::AliForwardMCFlowTaskQC() 
  : AliForwardFlowTaskQC(), 
    fBinsFMDTR(),          // List of FMDTR analysis objects
    fBinsSPDTR(),          // List of SPDTR analysis objects
    fBinsMC(),             // List of MC truth analysis objects
    fdNdedpMC(),           // MC truth d^2N/detadphi histogram
    fWeights(),            // Flow weights
    fImpactParToCent(),    // Impact parameter to centrality graph
    fUseImpactPar(0),      // Use impact par for centrality
    fFMDMinEta(-6),        // FMD min eta coverage for this vtx
    fFMDMaxEta(6),         // FMD max eta coverage for this vtx
    fAddFlow(0),           // Add flow to MC truth
    fAddType(0),           // Add type of flow to MC truth
    fAddOrder(0)           // Add order of flow to MC truth        
{} 
  //
  // Default Constructor
  //
//_____________________________________________________________________
AliForwardMCFlowTaskQC::AliForwardMCFlowTaskQC(const char* name) 
  : AliForwardFlowTaskQC(name),
    fBinsFMDTR(),          // List of FMDTR analysis objects
    fBinsSPDTR(),          // List of SPDTR analysis objects
    fBinsMC(),             // List of MC truth analysis objects
    fdNdedpMC(),           // MC truth d^2N/detadphi histogram
    fWeights(),            // Flow weights
    fImpactParToCent(),    // Impact parameter to centrality graph
    fUseImpactPar(0),      // Use impact par for centrality
    fFMDMinEta(-6),        // FMD min eta coverage for this vtx
    fFMDMaxEta(6),         // FMD max eta coverage for this vtx
    fAddFlow(0),           // Add flow to MC truth
    fAddType(0),           // Add type of flow to MC truth
    fAddOrder(0)           // Add order of flow to MC truth        
{ 
  // 
  // 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), 
    fBinsFMDTR(),                          // List of FMDTR analysis objects
    fBinsSPDTR(),                          // List of SPDTR analysis objects
    fBinsMC(),                             // List of MC truth analysis objects
    fdNdedpMC(o.fdNdedpMC),                // MC truth d^2N/detadphi histogram
    fWeights(o.fWeights),                  // Flow weights
    fImpactParToCent(o.fImpactParToCent),  // Impact parameter to centrality
    fUseImpactPar(o.fUseImpactPar),        // Use impact par for centrality
    fFMDMinEta(o.fFMDMinEta),              // FMD min eta coverage for this vtx
    fFMDMaxEta(o.fFMDMaxEta),              // FMD max eta coverage for this vtx
    fAddFlow(o.fAddFlow),                  // Add flow to MC truth
    fAddType(o.fAddType),                  // Add type of flow to MC truth
    fAddOrder(o.fAddOrder)                 // Add order of flow to MC truth
{
  //
  // Copy Constructor
  //
} 
//_____________________________________________________________________
AliForwardMCFlowTaskQC&
AliForwardMCFlowTaskQC::operator=(const AliForwardMCFlowTaskQC& o)
{
  // 
  // Assignment operator
  // Parameters:
  //    o Object to copy from 
  //
  if (&o == this) return *this;
  fdNdedpMC        = o.fdNdedpMC;
  fWeights         = o.fWeights;
  fImpactParToCent = o.fImpactParToCent;
  fUseImpactPar    = o.fUseImpactPar;
  fFMDMinEta       = o.fFMDMinEta;
  fFMDMaxEta       = o.fFMDMaxEta;
  fAddFlow         = o.fAddFlow;
  fAddType         = o.fAddType;
  fAddOrder        = o.fAddOrder;
  return *this;
}
//_____________________________________________________________________
void AliForwardMCFlowTaskQC::InitVertexBins()
{
  //
  // Initiate VertexBin lists
  //
  AliForwardFlowTaskQC::InitVertexBins();

  Int_t moment = 0;
  for(UShort_t n = 0; n < fV.GetSize(); n++) {
    moment = fV.At(n);
    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));
      fBinsFMDTR.Add(new VertexBin(vL, vH, moment, "FMDTR", fFlowFlags, fFMDCut, fEtaGap));
      fBinsSPDTR.Add(new VertexBin(vL, vH, moment, "SPDTR", fFlowFlags, fSPDCut, fEtaGap));
      fBinsMC.Add(new VertexBin(vL, vH, moment, "MC", fFlowFlags, -1, fEtaGap));
    }
  }
}
//_____________________________________________________________________
void AliForwardMCFlowTaskQC::InitHists()
{
  //
  // Initiate diagnostics hists and add to outputlist
  //
  AliForwardFlowTaskQC::InitHists();

  fdNdedpMC = TH2D(Form("fdNdedpMC%s", ((fFlowFlags & kEtaGap) ? "_etaGap" : "")),
		   Form("fdNdedpMC%s", ((fFlowFlags & kEtaGap) ? "_etaGap" : "")),
		   48, -6., 6., 200, 0., 2.*TMath::Pi());
  fdNdedpMC.Sumw2();

  TIter nextFMDTR(&fBinsFMDTR);
  VertexBin* bin = 0;
  while ((bin = static_cast<VertexBin*>(nextFMDTR()))) {
    bin->AddOutput(fSumList);
  }
  TIter nextSPDTR(&fBinsSPDTR);
  while ((bin = static_cast<VertexBin*>(nextSPDTR()))) {
    bin->AddOutput(fSumList);
  }
  TIter nextMC(&fBinsMC);
  while ((bin = static_cast<VertexBin*>(nextMC()))) {
    bin->AddOutput(fSumList);
  }

  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
  const AliAODForwardMult* aodfmult = 
    static_cast<AliAODForwardMult*>(fAOD->FindListObject("ForwardMC"));
  const AliAODCentralMult* aodcmult = 
    static_cast<AliAODCentralMult*>(fAOD->FindListObject("CentralClustersMC"));
  Int_t vtx = fVtxAxis->FindBin(fVtx)-1;
  
  // if objects are present, get histograms
  if (aodfmult) {
    const TH2D& fmdTRdNdetadphi = aodfmult->GetHistogram();
    if ((fFlowFlags & kEtaGap)) {
      FillVtxBinListEtaGap(fBinsFMDTR, fmdTRdNdetadphi, fmdTRdNdetadphi, vtx);
    } else {
      FillVtxBinList(fBinsFMDTR, fmdTRdNdetadphi, vtx);
    }
    
    if (aodcmult) {
      const TH2D& spdTRdNdetadphi = aodcmult->GetHistogram();
      if ((fFlowFlags & kEtaGap)) {
	FillVtxBinListEtaGap(fBinsSPDTR, fmdTRdNdetadphi, spdTRdNdetadphi, vtx);
      } else {
	FillVtxBinList(fBinsSPDTR, spdTRdNdetadphi, vtx);
      }
    }
  }
  // Run analysis on MC branch
  if (!LoopAODMC()) return kFALSE;
  if ((fFlowFlags & kEtaGap)) {
    FillVtxBinListEtaGap(fBinsMC, fdNdedpMC, fdNdedpMC, vtx);
  } else {
    FillVtxBinList(fBinsMC, fdNdedpMC, vtx);
  }

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

  EndVtxBinList(fBinsFMDTR);
  EndVtxBinList(fBinsSPDTR);
  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
  //
  return aodfm->IsTriggerBits(AliAODForwardMult::kB);
}
// _____________________________________________________________________
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.
  //
  if (fUseImpactPar) {
    fCent = GetCentFromB();
    if (fCent != -1) return kTRUE;
  }
  return AliForwardFlowTaskQC::GetCentrality(aodfm);
}
//_____________________________________________________________________
void AliForwardMCFlowTaskQC::GetFMDLimits()
{

  const AliAODForwardMult* aodfmult = 
    static_cast<AliAODForwardMult*>(fAOD->FindListObject("Forward"));
  const TH2D& h = aodfmult->GetHistogram();

  for (Int_t e = 1; ; e++) {
    if (h.GetBinContent(e, 0) != 0) { 
      fFMDMinEta = h.GetXaxis()->GetBinLowEdge(e);
      break;
    }
  }
  for (Int_t e = h.GetNbinsX(); ; e--) {
    if (h.GetBinContent(e, 0) != 0) { 
      fFMDMaxEta = h.GetXaxis()->GetBinLowEdge(e);
      break;
    }
  }

  return;
}
//_____________________________________________________________________
Bool_t AliForwardMCFlowTaskQC::LoopAODMC()  
{
  // 
  // Loop over AliAODParticle branch and fill d^2N/detadphi-histograms.
  // Add flow if set to do so in AddTask function
  //
  fdNdedpMC.Reset();
  GetFMDLimits();

  //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 = 
    dynamic_cast<AliAODMCHeader*>(fAOD->FindListObject(
                                  AliAODMCHeader::StdBranchName()));
  if (!header) 
    AliWarning("No header file found.");
  
  // Double_t rp = header->GetReactionPlaneAngle();

  Int_t ntracks = mcArray->GetEntriesFast();
  // TODO: Make this bit smarter...
  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;
  }
  // Double_t b = header->GetImpactParameter();

  // Track loop: chek how many particles will be accepted
  Double_t weight = 0;
  for (Int_t it = 0; it < ntracks; it++) {
    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 >= fFMDMinEta && eta <= fFMDMaxEta) {
      // 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);
      }*/
        fdNdedpMC.Fill(eta, phi, weight);
    }
  }
  Int_t sBin = fdNdedpMC.GetXaxis()->FindBin(fFMDMinEta);
  Int_t eBin = fdNdedpMC.GetXaxis()->FindBin(fFMDMaxEta);
  for ( ; sBin < eBin; sBin++) fdNdedpMC.SetBinContent(sBin, 0, 1.);

  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;
}
//_____________________________________________________________________
void AliForwardMCFlowTaskQC::PrintFlowSetup() const  
{
  //
  // Print the setup of the flow task
  //
  Printf("AliForwardMCFlowTaskQC::Print");
  Printf("Number of bins in vertex axis:\t%d", fVtxAxis->GetNbins());
  Printf("Range of vertex axis         :\t[%3.1f,%3.1f]", 
			  fVtxAxis->GetXmin(), fVtxAxis->GetXmax());
  printf("Doing flow analysis for      :\t");
  for (Int_t n  = 0; n < fV.GetSize(); n++) printf("v%d ", fV.At(n));
  printf("\n");
  Printf("Satellite vertex flag           :\t%s", ((fFlowFlags & kSatVtx) ? 
						   "true" : "false"));
  Printf("Symmetrize ref. flow wrt eta = 0:\t%s", ((fFlowFlags & kSymEta) ? 
						   "true" : "false"));
  Printf("Use an eta-gap for ref. flow    :\t%s", ((fFlowFlags & kEtaGap) ? 
						   "true" : "false"));
  Printf("FMD sigma cut:               :\t%f", fFMDCut);
  Printf("SPD sigma cut:               :\t%f", fSPDCut);
  if ((fFlowFlags & kEtaGap)) 
    Printf("Eta gap:                     :\t%f", fEtaGap);
}
//_____________________________________________________________________
//
//
// 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	//
33438b4c	10	/**
	11	* @file AliForwardMCFlowTaskQC.cxx
	12	* @author Christian Holm Christensen <cholm@master.hehi.nbi.dk>
	13	* @date Thu Feb 7 01:09:19 2013
	14	*
	15	* @brief
	16	*
	17	*
	18	* @ingroup pwglf_forward_flow
	19	*/
2b556440	20	#include "AliForwardMCFlowTaskQC.h"
	21	#include "AliAODMCParticle.h"
	22	#include "AliAODMCHeader.h"
	23	#include "TGraph.h"
	24	#include "TF1.h"
	25	#include "TProfile2D.h"
	26	#include "AliAODEvent.h"
	27	#include "AliAODForwardMult.h"
	28	#include "AliAODCentralMult.h"
008eda2b	29	#include "AliGenEventHeader.h"
2b556440	30
2b556440	31	ClassImp(AliForwardMCFlowTaskQC)
936b0a6c	32	#if 0
	33	;
	34	#endif
2b556440	35	//_____________________________________________________________________
936b0a6c	36	AliForwardMCFlowTaskQC::AliForwardMCFlowTaskQC()
	37	: AliForwardFlowTaskQC(),
	38	fBinsFMDTR(), // List of FMDTR analysis objects
	39	fBinsSPDTR(), // List of SPDTR analysis objects
	40	fBinsMC(), // List of MC truth analysis objects
	41	fdNdedpMC(), // MC truth d^2N/detadphi histogram
	42	fWeights(), // Flow weights
	43	fImpactParToCent(), // Impact parameter to centrality graph
	44	fUseImpactPar(0), // Use impact par for centrality
008eda2b	45	fFMDMinEta(-6), // FMD min eta coverage for this vtx
008eda2b	46	fFMDMaxEta(6), // FMD max eta coverage for this vtx
936b0a6c	47	fAddFlow(0), // Add flow to MC truth
	48	fAddType(0), // Add type of flow to MC truth
	49	fAddOrder(0) // Add order of flow to MC truth
	50	{}
2b556440	51	//
	52	// Default Constructor
	53	//
	54	//_____________________________________________________________________
936b0a6c	55	AliForwardMCFlowTaskQC::AliForwardMCFlowTaskQC(const char* name)
	56	: AliForwardFlowTaskQC(name),
	57	fBinsFMDTR(), // List of FMDTR analysis objects
	58	fBinsSPDTR(), // List of SPDTR analysis objects
	59	fBinsMC(), // List of MC truth analysis objects
	60	fdNdedpMC(), // MC truth d^2N/detadphi histogram
	61	fWeights(), // Flow weights
	62	fImpactParToCent(), // Impact parameter to centrality graph
	63	fUseImpactPar(0), // Use impact par for centrality
008eda2b	64	fFMDMinEta(-6), // FMD min eta coverage for this vtx
008eda2b	65	fFMDMaxEta(6), // FMD max eta coverage for this vtx
936b0a6c	66	fAddFlow(0), // Add flow to MC truth
	67	fAddType(0), // Add type of flow to MC truth
	68	fAddOrder(0) // Add order of flow to MC truth
2b556440	69	{
	70	//
	71	// Constructor
	72	// Parameters:
	73	// name: Name of task
	74	//
2b556440	75
936b0a6c	76	// Double_t impactParam[] = {0.,1.75,4.225,5.965,7.765,9.215,10.46,
	77	// 11.565,12.575,13.515,16.679};
	78	// Double_t centrality[] = {0.,2.5,7.5,15,25,35,45,55,65,75,90};
	79	Double_t impactParam[] = { 0.00, 3.72, 5.23, 7.31, 8.88, 10.20,
	80	11.38, 12.47, 13.50, 14.51, 16.679};
	81	Double_t centrality[] = { 0., 5., 10., 20., 30., 40.,
	82	50., 60., 70., 80., 100.};
2b556440	83
	84	Int_t nPoints = sizeof(impactParam)/sizeof(Double_t);
	85	fImpactParToCent = new TGraph(nPoints, impactParam, centrality);
	86
	87	}
	88	//_____________________________________________________________________
936b0a6c	89	AliForwardMCFlowTaskQC::AliForwardMCFlowTaskQC(const AliForwardMCFlowTaskQC& o)
	90	: AliForwardFlowTaskQC(o),
	91	fBinsFMDTR(), // List of FMDTR analysis objects
	92	fBinsSPDTR(), // List of SPDTR analysis objects
	93	fBinsMC(), // List of MC truth analysis objects
	94	fdNdedpMC(o.fdNdedpMC), // MC truth d^2N/detadphi histogram
	95	fWeights(o.fWeights), // Flow weights
	96	fImpactParToCent(o.fImpactParToCent), // Impact parameter to centrality
	97	fUseImpactPar(o.fUseImpactPar), // Use impact par for centrality
008eda2b	98	fFMDMinEta(o.fFMDMinEta), // FMD min eta coverage for this vtx
008eda2b	99	fFMDMaxEta(o.fFMDMaxEta), // FMD max eta coverage for this vtx
936b0a6c	100	fAddFlow(o.fAddFlow), // Add flow to MC truth
	101	fAddType(o.fAddType), // Add type of flow to MC truth
	102	fAddOrder(o.fAddOrder) // Add order of flow to MC truth
	103	{
2b556440	104	//
	105	// Copy Constructor
	106	//
936b0a6c	107	}
2b556440	108	//_____________________________________________________________________
	109	AliForwardMCFlowTaskQC&
	110	AliForwardMCFlowTaskQC::operator=(const AliForwardMCFlowTaskQC& o)
	111	{
	112	//
	113	// Assignment operator
	114	// Parameters:
	115	// o Object to copy from
	116	//
eeb6c967	117	if (&o == this) return *this;
2b556440	118	fdNdedpMC = o.fdNdedpMC;
936b0a6c	119	fWeights = o.fWeights;
2b556440	120	fImpactParToCent = o.fImpactParToCent;
4b5b52b7	121	fUseImpactPar = o.fUseImpactPar;
008eda2b	122	fFMDMinEta = o.fFMDMinEta;
008eda2b	123	fFMDMaxEta = o.fFMDMaxEta;
2b556440	124	fAddFlow = o.fAddFlow;
	125	fAddType = o.fAddType;
	126	fAddOrder = o.fAddOrder;
	127	return *this;
	128	}
	129	//_____________________________________________________________________
	130	void AliForwardMCFlowTaskQC::InitVertexBins()
	131	{
	132	//
	133	// Initiate VertexBin lists
	134	//
	135	AliForwardFlowTaskQC::InitVertexBins();
	136
68589651	137	Int_t moment = 0;
	138	for(UShort_t n = 0; n < fV.GetSize(); n++) {
	139	moment = fV.At(n);
	140	for (Int_t v = 1; v <= fVtxAxis->GetNbins(); v++) {
a7989168	141	Int_t vL = Int_t(fVtxAxis->GetBinLowEdge(v));
a7989168	142	Int_t vH = Int_t(fVtxAxis->GetBinUpEdge(v));
68589651	143	fBinsFMDTR.Add(new VertexBin(vL, vH, moment, "FMDTR", fFlowFlags, fFMDCut, fEtaGap));
	144	fBinsSPDTR.Add(new VertexBin(vL, vH, moment, "SPDTR", fFlowFlags, fSPDCut, fEtaGap));
	145	fBinsMC.Add(new VertexBin(vL, vH, moment, "MC", fFlowFlags, -1, fEtaGap));
2b556440	146	}
2b556440	147	}
2b556440	148	}
	149	//_____________________________________________________________________
	150	void AliForwardMCFlowTaskQC::InitHists()
	151	{
	152	//
	153	// Initiate diagnostics hists and add to outputlist
	154	//
	155	AliForwardFlowTaskQC::InitHists();
	156
68589651	157	fdNdedpMC = TH2D(Form("fdNdedpMC%s", ((fFlowFlags & kEtaGap) ? "_etaGap" : "")),
	158	Form("fdNdedpMC%s", ((fFlowFlags & kEtaGap) ? "_etaGap" : "")),
	159	48, -6., 6., 200, 0., 2.*TMath::Pi());
	160	fdNdedpMC.Sumw2();
	161
2b556440	162	TIter nextFMDTR(&fBinsFMDTR);
	163	VertexBin* bin = 0;
	164	while ((bin = static_cast<VertexBin*>(nextFMDTR()))) {
	165	bin->AddOutput(fSumList);
	166	}
	167	TIter nextSPDTR(&fBinsSPDTR);
	168	while ((bin = static_cast<VertexBin*>(nextSPDTR()))) {
	169	bin->AddOutput(fSumList);
	170	}
	171	TIter nextMC(&fBinsMC);
	172	while ((bin = static_cast<VertexBin*>(nextMC()))) {
	173	bin->AddOutput(fSumList);
	174	}
	175
3112bd03	176	TList* wList = new TList();
	177	wList->SetName("FlowWeights");
	178	fWeights.Init(wList);
	179	fSumList->Add(wList);
	180
2b556440	181	}
	182	//_____________________________________________________________________
	183	Bool_t AliForwardMCFlowTaskQC::Analyze()
	184	{
	185	//
	186	// Load FMD and SPD MC objects from aod tree and call the cumulants
	187	// calculation for the correct vertexbin
	188	//
	189	if (!AliForwardFlowTaskQC::Analyze()) return kFALSE;
	190
	191	// Run analysis on trackrefs from FMD and SPD
936b0a6c	192	const AliAODForwardMult* aodfmult =
	193	static_cast<AliAODForwardMult*>(fAOD->FindListObject("ForwardMC"));
	194	const AliAODCentralMult* aodcmult =
	195	static_cast<AliAODCentralMult*>(fAOD->FindListObject("CentralClustersMC"));
d420e249	196	Int_t vtx = fVtxAxis->FindBin(fVtx)-1;
2b556440	197
4b5b52b7	198	// if objects are present, get histograms
d420e249	199	if (aodfmult) {
d420e249	200	const TH2D& fmdTRdNdetadphi = aodfmult->GetHistogram();
68589651	201	if ((fFlowFlags & kEtaGap)) {
	202	FillVtxBinListEtaGap(fBinsFMDTR, fmdTRdNdetadphi, fmdTRdNdetadphi, vtx);
	203	} else {
	204	FillVtxBinList(fBinsFMDTR, fmdTRdNdetadphi, vtx);
	205	}
008eda2b	206
68589651	207	if (aodcmult) {
	208	const TH2D& spdTRdNdetadphi = aodcmult->GetHistogram();
	209	if ((fFlowFlags & kEtaGap)) {
	210	FillVtxBinListEtaGap(fBinsSPDTR, fmdTRdNdetadphi, spdTRdNdetadphi, vtx);
	211	} else {
	212	FillVtxBinList(fBinsSPDTR, spdTRdNdetadphi, vtx);
	213	}
	214	}
d420e249	215	}
2b556440	216	// Run analysis on MC branch
2b556440	217	if (!LoopAODMC()) return kFALSE;
68589651	218	if ((fFlowFlags & kEtaGap)) {
	219	FillVtxBinListEtaGap(fBinsMC, fdNdedpMC, fdNdedpMC, vtx);
	220	} else {
	221	FillVtxBinList(fBinsMC, fdNdedpMC, vtx);
	222	}
2b556440	223
	224	return kTRUE;
	225	}
	226	//_____________________________________________________________________
	227	void AliForwardMCFlowTaskQC::Finalize()
	228	{
	229	//
	230	// Finalize command, called by Terminate()
	231	//
	232	AliForwardFlowTaskQC::Finalize();
	233
4b5b52b7	234	EndVtxBinList(fBinsFMDTR);
	235	EndVtxBinList(fBinsSPDTR);
	236	EndVtxBinList(fBinsMC);
2b556440	237
4b5b52b7	238	return;
	239	}
	240	//_____________________________________________________________________
d420e249	241	Bool_t AliForwardMCFlowTaskQC::CheckTrigger(const AliAODForwardMult* aodfm) const
	242	{
	243	//
	244	// Function to look for a trigger string in the event.
	245	//
	246	// Parameters:
	247	// AliAODForwardMult: forward mult object with trigger and vertex info
	248	//
	249	// Returns true if B trigger is present - for some reason this is the one we use in MC
	250	//
	251	return aodfm->IsTriggerBits(AliAODForwardMult::kB);
	252	}
	253	// _____________________________________________________________________
4b5b52b7	254	Bool_t AliForwardMCFlowTaskQC::GetCentrality(const AliAODForwardMult* aodfm)
	255	{
	256	//
	257	// Function to use centrality parametrization from impact parameter
	258	// if flag is not set call AliForwardFlowTaskQC::GetCentrality
	259	//
	260	// Parameters:
	261	// AliAODForwardMult: forward mult object with trigger and vertex info
	262	//
	263	// Returns true when centrality is set.
	264	//
	265	if (fUseImpactPar) {
	266	fCent = GetCentFromB();
d420e249	267	if (fCent != -1) return kTRUE;
2b556440	268	}
d420e249	269	return AliForwardFlowTaskQC::GetCentrality(aodfm);
2b556440	270	}
2b556440	271	//_____________________________________________________________________
008eda2b	272	void AliForwardMCFlowTaskQC::GetFMDLimits()
	273	{
	274
	275	const AliAODForwardMult* aodfmult =
	276	static_cast<AliAODForwardMult*>(fAOD->FindListObject("Forward"));
	277	const TH2D& h = aodfmult->GetHistogram();
	278
	279	for (Int_t e = 1; ; e++) {
	280	if (h.GetBinContent(e, 0) != 0) {
	281	fFMDMinEta = h.GetXaxis()->GetBinLowEdge(e);
	282	break;
	283	}
	284	}
	285	for (Int_t e = h.GetNbinsX(); ; e--) {
	286	if (h.GetBinContent(e, 0) != 0) {
	287	fFMDMaxEta = h.GetXaxis()->GetBinLowEdge(e);
	288	break;
	289	}
	290	}
	291
	292	return;
	293	}
	294	//_____________________________________________________________________
2b556440	295	Bool_t AliForwardMCFlowTaskQC::LoopAODMC()
	296	{
	297	//
	298	// Loop over AliAODParticle branch and fill d^2N/detadphi-histograms.
	299	// Add flow if set to do so in AddTask function
68589651	300	//
2b556440	301	fdNdedpMC.Reset();
008eda2b	302	GetFMDLimits();
2b556440	303
2b556440	304	//retreive MC particles from event
936b0a6c	305	TClonesArray* mcArray =
d420e249	306	static_cast<TClonesArray*>(fAOD->FindListObject(
d420e249	307	AliAODMCParticle::StdBranchName()));
2b556440	308	if(!mcArray){
d420e249	309	AliWarning("No MC array found in AOD. Try making it again.");
2b556440	310	return kFALSE;
	311	}
	312
936b0a6c	313	AliAODMCHeader* header =
d420e249	314	dynamic_cast<AliAODMCHeader*>(fAOD->FindListObject(
d420e249	315	AliAODMCHeader::StdBranchName()));
936b0a6c	316	if (!header)
2b556440	317	AliWarning("No header file found.");
008eda2b	318
9c686cb4	319	// Double_t rp = header->GetReactionPlaneAngle();
2b556440	320
2b556440	321	Int_t ntracks = mcArray->GetEntriesFast();
008eda2b	322	// TODO: Make this bit smarter...
	323	if (header->GetNCocktailHeaders() > 1) {
	324	ntracks = header->GetCocktailHeader(0)->NProduced();
	325	}
2b556440	326
936b0a6c	327	UShort_t flowFlags = 0;
	328	if (fAddFlow.Length() > 1) {
	329	if (fAddFlow.Contains("pt")) flowFlags \|= AliForwardFlowWeights::kPt;
	330	if (fAddFlow.Contains("pid")) flowFlags \|= AliForwardFlowWeights::kPt;
	331	if (fAddFlow.Contains("eta")) flowFlags \|= AliForwardFlowWeights::kEta;
	332	if (fAddFlow.Contains("cent")) flowFlags \|= AliForwardFlowWeights::kCent;
	333	}
9c686cb4	334	// Double_t b = header->GetImpactParameter();
936b0a6c	335
2b556440	336	// Track loop: chek how many particles will be accepted
	337	Double_t weight = 0;
	338	for (Int_t it = 0; it < ntracks; it++) {
	339	weight = 1;
	340	AliAODMCParticle* particle = (AliAODMCParticle*)mcArray->At(it);
	341	if (!particle) {
	342	AliError(Form("Could not receive track %d", it));
	343	continue;
	344	}
008eda2b	345	if (!particle->IsPrimary()) continue;
2b556440	346	if (particle->Charge() == 0) continue;
68589651	347	//Double_t pT = particle->Pt();
2b556440	348	Double_t eta = particle->Eta();
2b556440	349	Double_t phi = particle->Phi();
008eda2b	350	if (eta >= fFMDMinEta && eta <= fFMDMaxEta) {
2b556440	351	// Add flow if it is in the argument
68589651	352	/* FLOW WEIGHTS DISABLED IN THE VERSION - COMING BACK SOON
	353	if (flowFlags != 0) {
	354	// weight = fWeights.CalcWeight(eta, pT, phi, particle->PdgCode(),
	355	// rp, fCent, fAddType, fAddOrder,
	356	// flowFlags) + 1;
03f52e9d	357	weight = fWeights.CalcWeight(eta, pT, phi, particle->PdgCode(),
03f52e9d	358	rp, b);
68589651	359	// Printf("%f", weight);
68589651	360	}*/
03f52e9d	361	fdNdedpMC.Fill(eta, phi, weight);
2b556440	362	}
2b556440	363	}
008eda2b	364	Int_t sBin = fdNdedpMC.GetXaxis()->FindBin(fFMDMinEta);
008eda2b	365	Int_t eBin = fdNdedpMC.GetXaxis()->FindBin(fFMDMaxEta);
d420e249	366	for ( ; sBin < eBin; sBin++) fdNdedpMC.SetBinContent(sBin, 0, 1.);
2b556440	367
	368	return kTRUE;
	369	}
	370	//_____________________________________________________________________
2b556440	371	Double_t AliForwardMCFlowTaskQC::GetCentFromB() const
	372	{
	373	//
	374	// Get centrality from MC impact parameter.
2b556440	375	//
2b556440	376	Double_t cent = -1.;
936b0a6c	377	AliAODMCHeader* header =
	378	static_cast<AliAODMCHeader*>(fAOD->FindListObject(AliAODMCHeader::
	379	StdBranchName()));
2b556440	380	if (!header) return cent;
9c686cb4	381	Double_t b = header->GetImpactParameter();
2b556440	382
	383	cent = fImpactParToCent->Eval(b);
	384
	385	return cent;
	386	}
	387	//_____________________________________________________________________
d420e249	388	void AliForwardMCFlowTaskQC::PrintFlowSetup() const
	389	{
	390	//
	391	// Print the setup of the flow task
	392	//
	393	Printf("AliForwardMCFlowTaskQC::Print");
	394	Printf("Number of bins in vertex axis:\t%d", fVtxAxis->GetNbins());
	395	Printf("Range of vertex axis :\t[%3.1f,%3.1f]",
	396	fVtxAxis->GetXmin(), fVtxAxis->GetXmax());
	397	printf("Doing flow analysis for :\t");
68589651	398	for (Int_t n = 0; n < fV.GetSize(); n++) printf("v%d ", fV.At(n));
d420e249	399	printf("\n");
9c686cb4	400	Printf("Satellite vertex flag :\t%s", ((fFlowFlags & kSatVtx) ?
	401	"true" : "false"));
	402	Printf("Symmetrize ref. flow wrt eta = 0:\t%s", ((fFlowFlags & kSymEta) ?
	403	"true" : "false"));
	404	Printf("Use an eta-gap for ref. flow :\t%s", ((fFlowFlags & kEtaGap) ?
	405	"true" : "false"));
d420e249	406	Printf("FMD sigma cut: :\t%f", fFMDCut);
d420e249	407	Printf("SPD sigma cut: :\t%f", fSPDCut);
68589651	408	if ((fFlowFlags & kEtaGap))
68589651	409	Printf("Eta gap: :\t%f", fEtaGap);
d420e249	410	}
d420e249	411	//_____________________________________________________________________
2b556440	412	//
	413	//
	414	// EOF
	415