[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 "TProfile2D.h"
#include "AliAODEvent.h"
#include "AliAODForwardMult.h"
#include "AliAODCentralMult.h"

ClassImp(AliForwardMCFlowTaskQC)

//_____________________________________________________________________
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
  fAliceCent4th(),       // Alice QC4 vs. centrality data points
  fAlicePt2nd4050(),     // Alice QC2 vs. pT data points
  fAlicePt4th3040(),     // Alice QC4 vs. pT data points
  fAlicePt4th4050(),     // Alice QC4 vs. pT data points
  fImpactParToCent(),    // Impact parameter to centrality graph
  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
  fAliceCent4th(),       // Alice QC4 vs. centrality data points
  fAlicePt2nd4050(),     // Alice QC2 vs. pT data points
  fAlicePt4th3040(),     // Alice QC4 vs. pT data points
  fAlicePt4th4050(),     // Alice QC4 vs. pT data points
  fImpactParToCent(),    // Impact parameter to centrality graph
  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
  //
  fdNdedpMC = TH2D("fdNdedpMC", "fdNdedpMC", 48, -6., 6., 200, 0., 2.*TMath::Pi());
  fdNdedpMC.Sumw2();
  
  // Add parametrizations of ALICE data
  Double_t xCumulant4thTPCrefMultTPConlyAll[] = {2.5,7.5,15,25,35,45,55,65};
  Double_t yCumulant4thTPCrefMultTPConlyAll[] = {0.017855,0.032440,0.055818,0.073137,0.083898,0.086690,0.082040,0.077777};
  Int_t nPointsCumulant4thTPCrefMultTPConlyAll = sizeof(xCumulant4thTPCrefMultTPConlyAll)/sizeof(Double_t);
  fAliceCent4th = new TGraph(nPointsCumulant4thTPCrefMultTPConlyAll,xCumulant4thTPCrefMultTPConlyAll,
                                                        yCumulant4thTPCrefMultTPConlyAll);
 
  Double_t xCumulant2nd4050ALICE[] = {0.000000,0.250000,0.350000,0.450000,0.550000,0.650000,0.750000,0.850000,0.950000,
  1.100000,1.300000,1.500000,1.700000,1.900000,2.250000,2.750000,3.250000,3.750000,4.500000};
  Double_t yCumulant2nd4050ALICE[] = {0.000000,0.043400,0.059911,0.073516,0.089756,0.105486,0.117391,0.128199,0.138013,
  0.158271,0.177726,0.196383,0.208277,0.216648,0.242954,0.249961,0.240131,0.269006,0.207796};
  Int_t nPointsCumulant2nd4050ALICE = sizeof(xCumulant2nd4050ALICE)/sizeof(Double_t);                                      
  fAlicePt2nd4050 = new TGraph(nPointsCumulant2nd4050ALICE,xCumulant2nd4050ALICE,yCumulant2nd4050ALICE);

  Double_t xCumulant4th3040ALICE[] = {0.00000,0.250000,0.350000,0.450000,0.550000,0.650000,0.750000,0.850000,0.950000,
  1.100000,1.300000,1.500000,1.700000,1.900000,2.250000,2.750000,3.250000,3.750000,4.500000,
  5.500000,7.000000,9.000000};
  Double_t yCumulant4th3040ALICE[] = {0.000000,0.037071,0.048566,0.061083,0.070910,0.078831,0.091396,0.102026,0.109691,
  0.124449,0.139819,0.155561,0.165701,0.173678,0.191149,0.202015,0.204540,0.212560,0.195885,
  0.000000,0.000000,0.000000};
  Int_t nPointsCumulant4th3040ALICE = sizeof(xCumulant4th3040ALICE)/sizeof(Double_t);                                      
  fAlicePt4th3040 = new TGraph(nPointsCumulant4th3040ALICE,xCumulant4th3040ALICE,yCumulant4th3040ALICE);

  Double_t xCumulant4th4050ALICE[] = {0.000000,0.250000,0.350000,0.450000,0.550000,0.650000,0.750000,0.850000,0.950000,
  1.100000,1.300000,1.500000,1.700000,1.900000,2.250000,2.750000,3.250000,3.750000,4.500000};
  Double_t yCumulant4th4050ALICE[] = {0.000000,0.038646,0.049824,0.066662,0.075856,0.081583,0.099778,0.104674,0.118545,
  0.131874,0.152959,0.155348,0.169751,0.179052,0.178532,0.198851,0.185737,0.239901,0.186098};
  Int_t nPointsCumulant4th4050ALICE = sizeof(xCumulant4th4050ALICE)/sizeof(Double_t);   
  fAlicePt4th4050 = new TGraph(nPointsCumulant4th4050ALICE, xCumulant4th4050ALICE, yCumulant4th4050ALICE);

  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};

  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
  fAliceCent4th(o.fAliceCent4th),        // Alice QC4 vs. centrality data points
  fAlicePt2nd4050(o.fAlicePt2nd4050),    // Alice QC2 vs. pT data points
  fAlicePt4th3040(o.fAlicePt4th3040),    // Alice QC4 vs. pT data points
  fAlicePt4th4050(o.fAlicePt4th4050),    // Alice QC4 vs. pT data points
  fImpactParToCent(o.fImpactParToCent),  // Impact parameter to centrality graph
  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 
  //
  fdNdedpMC        = o.fdNdedpMC;
  fAliceCent4th    = o.fAliceCent4th;
  fAlicePt2nd4050  = o.fAlicePt2nd4050;
  fAlicePt4th3040  = o.fAlicePt4th3040;
  fAlicePt4th4050  = o.fAlicePt4th4050;
  fImpactParToCent = o.fImpactParToCent;
  fAddFlow         = o.fAddFlow;
  fAddType         = o.fAddType;
  fAddOrder        = o.fAddOrder;
  return *this;
}
//_____________________________________________________________________
void AliForwardMCFlowTaskQC::InitVertexBins()
{
  //
  // Initiate VertexBin lists
  //
  AliForwardFlowTaskQC::InitVertexBins();

  for(Int_t n = 1; n <= 6; n++) {
    if (!fv[n]) continue;
    for (Int_t v = -10; v < 10; v++) {
      fBinsFMDTR.Add(new VertexBin(v, v+1, n, "FMDTR"));
      fBinsSPDTR.Add(new VertexBin(v, v+1, n, "SPDTR"));
      fBinsMC.Add(new VertexBin(v, v+1, n, "MC"));
    }
  }

}
//_____________________________________________________________________
void AliForwardMCFlowTaskQC::InitHists()
{
  //
  // Initiate diagnostics hists and add to outputlist
  //
  AliForwardFlowTaskQC::InitHists();

  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);
  }

}
//_____________________________________________________________________
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"));
  if (!aodfmult) return kFALSE;
  TH2D fmdTRdNdetadphi = aodfmult->GetHistogram();

  AliAODCentralMult* aodcmult = static_cast<AliAODCentralMult*>(fAOD->FindListObject("CentralClustersMC"));
  if (!aodcmult) return kFALSE;
  TH2D spdTRdNdetadphi = aodcmult->GetHistogram();
  
  TIter nextFMDTR(&fBinsFMDTR);
  VertexBin* bin = 0;
  while ((bin = static_cast<VertexBin*>(nextFMDTR()))) {
    if (bin->CheckVertex(fZvertex)) {
      if (!bin->FillHists(fmdTRdNdetadphi)) return kFALSE;
      bin->CumulantsAccumulate(fCent);
    }
  }

  TIter nextSPDTR(&fBinsSPDTR);
  while ((bin = static_cast<VertexBin*>(nextSPDTR()))) {
    if (bin->CheckVertex(fZvertex)) {
      if (!bin->FillHists(spdTRdNdetadphi)) return kFALSE;
      bin->CumulantsAccumulate(fCent);
    }
  }

  // Run analysis on MC branch
  if (!LoopAODMC()) return kFALSE;

  TIter nextMC(&fBinsMC);
  while ((bin = static_cast<VertexBin*>(nextMC()))) {
    if (bin->CheckVertex(fZvertex)) {
      if (!bin->FillHists(fdNdedpMC)) return kFALSE;
      bin->CumulantsAccumulate(fCent);
    }
  }

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

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

  TProfile2D* fmdHist = 0;
  TProfile2D* spdHist = 0;
  TProfile2D* mcHist = 0;
  TList* correctionList = new TList();
  correctionList->SetName("CorrectionList");
//  fOutputList->Add(correctionList);

  for (Int_t i = 2; i <= 4; i += 2) {
    for (Int_t n = 1; n <= 6; n++) {
      if (!fv[n]) continue;
      fmdHist = (TProfile2D*)fOutputList->FindObject(Form("FMDQC%d_v%d_unCorr", i, n))
		  ->Clone(Form("FMDQC%d_v%d_Correction", i, n));
      spdHist = (TProfile2D*)fOutputList->FindObject(Form("SPDQC%d_v%d_unCorr", i, n))
		  ->Clone(Form("SPDQC%d_v%d_Correction", i, n));
      mcHist = (TProfile2D*)fOutputList->FindObject(Form("MCQC%d_v%d_unCorr", i, n));
     
      if (!fmdHist || !spdHist || !mcHist) {
	AliError(Form("Histogram missing, correction object not created for v%d", n));
	continue;
      }

      fmdHist->Divide(mcHist);
      spdHist->Divide(mcHist);
      fmdHist->SetTitle(Form("FMD QC{%d} v_{%d} Correction Object", i, n));
      fmdHist->SetTitle(Form("SPD QC{%d} v_{%d} Correction Object", i, n));

//      correctionList->Add(fmdHist);
//      correctionList->Add(spdHist);
    }
  }

}
//_____________________________________________________________________
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();

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

  Double_t rp = 0;
  AliAODMCHeader* header = dynamic_cast<AliAODMCHeader*>(fAOD->FindListObject(AliAODMCHeader::StdBranchName()));
  if (!header) {
    AliWarning("No header file found.");
  }
  else {
    rp = header->GetReactionPlaneAngle();
  }

  Int_t ntracks = mcArray->GetEntriesFast();

  // 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->IsPhysicalPrimary()) continue;
    if (particle->Charge() == 0) continue;
    Double_t pT = particle->Pt();
    Double_t eta = particle->Eta();
    Double_t phi = particle->Phi();
    if (TMath::Abs(eta) < 6.) {
      // Add flow if it is in the argument
      if (fAddFlow.Length() > 1) {
        if (fAddFlow.Contains("pt"))
          weight *= AddptFlow(pT);
        if (fAddFlow.Contains("pid"))
          weight *= AddpidFlow(particle->PdgCode());
        if (fAddFlow.Contains("eta"))
          weight *= AddetaFlow(eta);
        if (fAddFlow.Contains("cent")) 
          weight *= fAliceCent4th->Eval(fCent)/fAliceCent4th->Eval(45);
        
        weight *= 20*2.*TMath::Cos((Double_t)fAddOrder*(phi-rp)); 
        weight += 1;
      }
      fdNdedpMC.Fill(eta, phi, weight);
    }
  }

  return kTRUE;
}
//_____________________________________________________________________
Double_t AliForwardMCFlowTaskQC::AddptFlow(Double_t pt = 0) const 
{
  //
  // Add pt dependent flow factor
  // Parameters:
  //  pt: pT parametrization to use
  //
  Double_t weight = 0;

  switch(fAddType) 
  {
    case 1: weight = fAlicePt2nd4050->Eval(pt)*0.5+fAlicePt4th4050->Eval(pt)*0.5;
            break;
    case 2: weight = fAlicePt2nd4050->Eval(pt);
            break;
    case 3: weight = fAlicePt4th3040->Eval(pt);
            break;
    case 4: weight = fAlicePt4th4050->Eval(pt);
            break;
  } 
  
  return weight;
}
//_____________________________________________________________________
Double_t AliForwardMCFlowTaskQC::AddpidFlow(Int_t id = 0) const 
{
  //
  // Add pid dependent flow factor 
  // Parameters:
  //  id: choose PID dependent setup
  //
  Double_t weight = 0;

  switch(fAddType)
  {
    case 1: if (TMath::Abs(id) ==  211) // pion flow
              weight = 1.3;
            else if (TMath::Abs(id) ==  2212) // proton flow
              weight = 1.;
            else 
              weight = 0.7;
            break;
    case 2: weight = 1.207;
            break;
  }

  return weight;
}
//_____________________________________________________________________
Double_t AliForwardMCFlowTaskQC::AddetaFlow(Double_t eta = 0) const 
{
  //
  // Add eta dependent flow factor 
  // Parameters:
  //  eta: choose v_n(eta) shape
  //
  Double_t weight = 0;

  TF1 gaus = TF1("gaus", "gaus", -6, 6);

  switch(fAddType)
  {
     case 1: gaus.SetParameters(0.1, 0., 9);
             break;
     case 2: gaus.SetParameters(0.1, 0., 3);
             break;
     case 3: gaus.SetParameters(0.1, 0., 15);
             break;
  }

  weight = gaus.Eval(eta);

  return weight;
}
//_____________________________________________________________________
Double_t AliForwardMCFlowTaskQC::GetCentFromB() const
{
  //
  // Get centrality from MC impact parameter.
  // Values taken from: https://twiki.cern.ch/twiki/bin/viewauth/ALICE/CentStudies
  //
  Double_t cent = -1.;
  Double_t b = -1.;
  AliAODMCHeader* header = (AliAODMCHeader*)fAOD->FindListObject(AliAODMCHeader::StdBranchName());
  if (!header) return cent;
  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"
	15	#include "TProfile2D.h"
	16	#include "AliAODEvent.h"
	17	#include "AliAODForwardMult.h"
	18	#include "AliAODCentralMult.h"
	19
	20	ClassImp(AliForwardMCFlowTaskQC)
	21
	22	//_____________________________________________________________________
	23	AliForwardMCFlowTaskQC::AliForwardMCFlowTaskQC() :
	24	AliForwardFlowTaskQC(),
	25	fBinsFMDTR(), // List of FMDTR analysis objects
	26	fBinsSPDTR(), // List of SPDTR analysis objects
	27	fBinsMC(), // List of MC truth analysis objects
	28	fdNdedpMC(), // MC truth d^2N/detadphi histogram
	29	fAliceCent4th(), // Alice QC4 vs. centrality data points
	30	fAlicePt2nd4050(), // Alice QC2 vs. pT data points
	31	fAlicePt4th3040(), // Alice QC4 vs. pT data points
	32	fAlicePt4th4050(), // Alice QC4 vs. pT data points
	33	fImpactParToCent(), // Impact parameter to centrality graph
	34	fAddFlow(0), // Add flow to MC truth
	35	fAddType(0), // Add type of flow to MC truth
	36	fAddOrder(0) // Add order of flow to MC truth
	37	{}
	38	//
	39	// Default Constructor
	40	//
	41	//_____________________________________________________________________
	42	AliForwardMCFlowTaskQC::AliForwardMCFlowTaskQC(const char* name) :
	43	AliForwardFlowTaskQC(name),
	44	fBinsFMDTR(), // List of FMDTR analysis objects
	45	fBinsSPDTR(), // List of SPDTR analysis objects
	46	fBinsMC(), // List of MC truth analysis objects
	47	fdNdedpMC(), // MC truth d^2N/detadphi histogram
	48	fAliceCent4th(), // Alice QC4 vs. centrality data points
	49	fAlicePt2nd4050(), // Alice QC2 vs. pT data points
	50	fAlicePt4th3040(), // Alice QC4 vs. pT data points
	51	fAlicePt4th4050(), // Alice QC4 vs. pT data points
	52	fImpactParToCent(), // Impact parameter to centrality graph
	53	fAddFlow(0), // Add flow to MC truth
	54	fAddType(0), // Add type of flow to MC truth
	55	fAddOrder(0) // Add order of flow to MC truth
	56	{
	57	//
	58	// Constructor
	59	// Parameters:
	60	// name: Name of task
	61	//
	62	fdNdedpMC = TH2D("fdNdedpMC", "fdNdedpMC", 48, -6., 6., 200, 0., 2.*TMath::Pi());
	63	fdNdedpMC.Sumw2();
	64
65	// Add parametrizations of ALICE data
66	Double_t xCumulant4thTPCrefMultTPConlyAll[] = {2.5,7.5,15,25,35,45,55,65};
67	Double_t yCumulant4thTPCrefMultTPConlyAll[] = {0.017855,0.032440,0.055818,0.073137,0.083898,0.086690,0.082040,0.077777};
68	Int_t nPointsCumulant4thTPCrefMultTPConlyAll = sizeof(xCumulant4thTPCrefMultTPConlyAll)/sizeof(Double_t);
69	fAliceCent4th = new TGraph(nPointsCumulant4thTPCrefMultTPConlyAll,xCumulant4thTPCrefMultTPConlyAll,
70	yCumulant4thTPCrefMultTPConlyAll);
71
72	Double_t xCumulant2nd4050ALICE[] = {0.000000,0.250000,0.350000,0.450000,0.550000,0.650000,0.750000,0.850000,0.950000,
73	1.100000,1.300000,1.500000,1.700000,1.900000,2.250000,2.750000,3.250000,3.750000,4.500000};
74	Double_t yCumulant2nd4050ALICE[] = {0.000000,0.043400,0.059911,0.073516,0.089756,0.105486,0.117391,0.128199,0.138013,
75	0.158271,0.177726,0.196383,0.208277,0.216648,0.242954,0.249961,0.240131,0.269006,0.207796};
76	Int_t nPointsCumulant2nd4050ALICE = sizeof(xCumulant2nd4050ALICE)/sizeof(Double_t);
77	fAlicePt2nd4050 = new TGraph(nPointsCumulant2nd4050ALICE,xCumulant2nd4050ALICE,yCumulant2nd4050ALICE);
78
79	Double_t xCumulant4th3040ALICE[] = {0.00000,0.250000,0.350000,0.450000,0.550000,0.650000,0.750000,0.850000,0.950000,
80	1.100000,1.300000,1.500000,1.700000,1.900000,2.250000,2.750000,3.250000,3.750000,4.500000,
81	5.500000,7.000000,9.000000};
82	Double_t yCumulant4th3040ALICE[] = {0.000000,0.037071,0.048566,0.061083,0.070910,0.078831,0.091396,0.102026,0.109691,
83	0.124449,0.139819,0.155561,0.165701,0.173678,0.191149,0.202015,0.204540,0.212560,0.195885,
84	0.000000,0.000000,0.000000};
85	Int_t nPointsCumulant4th3040ALICE = sizeof(xCumulant4th3040ALICE)/sizeof(Double_t);
86	fAlicePt4th3040 = new TGraph(nPointsCumulant4th3040ALICE,xCumulant4th3040ALICE,yCumulant4th3040ALICE);
87
88	Double_t xCumulant4th4050ALICE[] = {0.000000,0.250000,0.350000,0.450000,0.550000,0.650000,0.750000,0.850000,0.950000,
89	1.100000,1.300000,1.500000,1.700000,1.900000,2.250000,2.750000,3.250000,3.750000,4.500000};
90	Double_t yCumulant4th4050ALICE[] = {0.000000,0.038646,0.049824,0.066662,0.075856,0.081583,0.099778,0.104674,0.118545,
91	0.131874,0.152959,0.155348,0.169751,0.179052,0.178532,0.198851,0.185737,0.239901,0.186098};
92	Int_t nPointsCumulant4th4050ALICE = sizeof(xCumulant4th4050ALICE)/sizeof(Double_t);
93	fAlicePt4th4050 = new TGraph(nPointsCumulant4th4050ALICE, xCumulant4th4050ALICE, yCumulant4th4050ALICE);
94
95	Double_t impactParam[] = {0.,1.75,4.225,5.965,7.765,9.215,10.46,11.565,12.575,13.515,16.679};
96	Double_t centrality[] = {0.,2.5,7.5,15,25,35,45,55,65,75,90};
97
98	Int_t nPoints = sizeof(impactParam)/sizeof(Double_t);
99	fImpactParToCent = new TGraph(nPoints, impactParam, centrality);
100
101	}
102	//_____________________________________________________________________
103	AliForwardMCFlowTaskQC::AliForwardMCFlowTaskQC(const AliForwardMCFlowTaskQC& o) :
104	AliForwardFlowTaskQC(o),
105	fBinsFMDTR(), // List of FMDTR analysis objects
106	fBinsSPDTR(), // List of SPDTR analysis objects
107	fBinsMC(), // List of MC truth analysis objects
108	fdNdedpMC(o.fdNdedpMC), // MC truth d^2N/detadphi histogram
109	fAliceCent4th(o.fAliceCent4th), // Alice QC4 vs. centrality data points
110	fAlicePt2nd4050(o.fAlicePt2nd4050), // Alice QC2 vs. pT data points
111	fAlicePt4th3040(o.fAlicePt4th3040), // Alice QC4 vs. pT data points
112	fAlicePt4th4050(o.fAlicePt4th4050), // Alice QC4 vs. pT data points
113	fImpactParToCent(o.fImpactParToCent), // Impact parameter to centrality graph
114	fAddFlow(o.fAddFlow), // Add flow to MC truth
115	fAddType(o.fAddType), // Add type of flow to MC truth
116	fAddOrder(o.fAddOrder) // Add order of flow to MC truth
117	{}
118	//
119	// Copy Constructor
120	//
121	//_____________________________________________________________________
122	AliForwardMCFlowTaskQC&
123	AliForwardMCFlowTaskQC::operator=(const AliForwardMCFlowTaskQC& o)
124	{
125	//
126	// Assignment operator
127	// Parameters:
128	// o Object to copy from
129	//
130	fdNdedpMC = o.fdNdedpMC;
131	fAliceCent4th = o.fAliceCent4th;
132	fAlicePt2nd4050 = o.fAlicePt2nd4050;
133	fAlicePt4th3040 = o.fAlicePt4th3040;
134	fAlicePt4th4050 = o.fAlicePt4th4050;
135	fImpactParToCent = o.fImpactParToCent;
136	fAddFlow = o.fAddFlow;
137	fAddType = o.fAddType;
138	fAddOrder = o.fAddOrder;
139	return *this;
140	}
141	//_____________________________________________________________________
142	void AliForwardMCFlowTaskQC::InitVertexBins()
143	{
144	//
145	// Initiate VertexBin lists
146	//
147	AliForwardFlowTaskQC::InitVertexBins();
148
149	for(Int_t n = 1; n <= 6; n++) {
150	if (!fv[n]) continue;
151	for (Int_t v = -10; v < 10; v++) {
152	fBinsFMDTR.Add(new VertexBin(v, v+1, n, "FMDTR"));
153	fBinsSPDTR.Add(new VertexBin(v, v+1, n, "SPDTR"));
154	fBinsMC.Add(new VertexBin(v, v+1, n, "MC"));
155	}
156	}
157
158	}
159	//_____________________________________________________________________
160	void AliForwardMCFlowTaskQC::InitHists()
161	{
162	//
163	// Initiate diagnostics hists and add to outputlist
164	//
165	AliForwardFlowTaskQC::InitHists();
166
167	TIter nextFMDTR(&fBinsFMDTR);
168	VertexBin* bin = 0;
169	while ((bin = static_cast<VertexBin*>(nextFMDTR()))) {
170	bin->AddOutput(fSumList);
171	}
172	TIter nextSPDTR(&fBinsSPDTR);
173	while ((bin = static_cast<VertexBin*>(nextSPDTR()))) {
174	bin->AddOutput(fSumList);
175	}
176	TIter nextMC(&fBinsMC);
177	while ((bin = static_cast<VertexBin*>(nextMC()))) {
178	bin->AddOutput(fSumList);
179	}
180
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
192	AliAODForwardMult* aodfmult = static_cast<AliAODForwardMult*>(fAOD->FindListObject("ForwardMC"));
193	if (!aodfmult) return kFALSE;
194	TH2D fmdTRdNdetadphi = aodfmult->GetHistogram();
195
196	AliAODCentralMult* aodcmult = static_cast<AliAODCentralMult*>(fAOD->FindListObject("CentralClustersMC"));
197	if (!aodcmult) return kFALSE;
198	TH2D spdTRdNdetadphi = aodcmult->GetHistogram();
199
200	TIter nextFMDTR(&fBinsFMDTR);
201	VertexBin* bin = 0;
202	while ((bin = static_cast<VertexBin*>(nextFMDTR()))) {
203	if (bin->CheckVertex(fZvertex)) {
204	if (!bin->FillHists(fmdTRdNdetadphi)) return kFALSE;
205	bin->CumulantsAccumulate(fCent);
206	}
207	}
208
209	TIter nextSPDTR(&fBinsSPDTR);
210	while ((bin = static_cast<VertexBin*>(nextSPDTR()))) {
211	if (bin->CheckVertex(fZvertex)) {
212	if (!bin->FillHists(spdTRdNdetadphi)) return kFALSE;
213	bin->CumulantsAccumulate(fCent);
214	}
215	}
216
217	// Run analysis on MC branch
218	if (!LoopAODMC()) return kFALSE;
219
220	TIter nextMC(&fBinsMC);
221	while ((bin = static_cast<VertexBin*>(nextMC()))) {
222	if (bin->CheckVertex(fZvertex)) {
223	if (!bin->FillHists(fdNdedpMC)) return kFALSE;
224	bin->CumulantsAccumulate(fCent);
225	}
226	}
227
228	return kTRUE;
229	}
230	//_____________________________________________________________________
231	void AliForwardMCFlowTaskQC::Finalize()
232	{
233	//
234	// Finalize command, called by Terminate()
235	//
236	AliForwardFlowTaskQC::Finalize();
237
238	TIter nextFMDTR(&fBinsFMDTR);
239	VertexBin* bin = 0;
240	while ((bin = static_cast<VertexBin*>(nextFMDTR()))) {
241	bin->CumulantsTerminate(fSumList, fOutputList);
242	}
243	TIter nextSPDTR(&fBinsSPDTR);
244	while ((bin = static_cast<VertexBin*>(nextSPDTR()))) {
245	bin->CumulantsTerminate(fSumList, fOutputList);
246	}
247	TIter nextMC(&fBinsMC);
248	while ((bin = static_cast<VertexBin*>(nextMC()))) {
249	bin->CumulantsTerminate(fSumList, fOutputList);
250	}
251
252	TProfile2D* fmdHist = 0;
253	TProfile2D* spdHist = 0;
254	TProfile2D* mcHist = 0;
255	TList* correctionList = new TList();
256	correctionList->SetName("CorrectionList");
257	// fOutputList->Add(correctionList);
258
259	for (Int_t i = 2; i <= 4; i += 2) {
260	for (Int_t n = 1; n <= 6; n++) {
261	if (!fv[n]) continue;
262	fmdHist = (TProfile2D*)fOutputList->FindObject(Form("FMDQC%d_v%d_unCorr", i, n))
263	->Clone(Form("FMDQC%d_v%d_Correction", i, n));
264	spdHist = (TProfile2D*)fOutputList->FindObject(Form("SPDQC%d_v%d_unCorr", i, n))
265	->Clone(Form("SPDQC%d_v%d_Correction", i, n));
266	mcHist = (TProfile2D*)fOutputList->FindObject(Form("MCQC%d_v%d_unCorr", i, n));
267
268	if (!fmdHist \|\| !spdHist \|\| !mcHist) {
269	AliError(Form("Histogram missing, correction object not created for v%d", n));
270	continue;
271	}
272
273	fmdHist->Divide(mcHist);
274	spdHist->Divide(mcHist);
275	fmdHist->SetTitle(Form("FMD QC{%d} v_{%d} Correction Object", i, n));
276	fmdHist->SetTitle(Form("SPD QC{%d} v_{%d} Correction Object", i, n));
277
278	// correctionList->Add(fmdHist);
279	// correctionList->Add(spdHist);
280	}
281	}
282
283	}
284	//_____________________________________________________________________
285	Bool_t AliForwardMCFlowTaskQC::LoopAODMC()
286	{
287	//
288	// Loop over AliAODParticle branch and fill d^2N/detadphi-histograms.
289	// Add flow if set to do so in AddTask function
290	fdNdedpMC.Reset();
291
292	//retreive MC particles from event
293	TClonesArray* mcArray = (TClonesArray*)fAOD->FindListObject(AliAODMCParticle::StdBranchName());
294	if(!mcArray){
295	// AliWarning("No MC array found in AOD. Try making it again.");
296	return kFALSE;
297	}
298
299	Double_t rp = 0;
300	AliAODMCHeader* header = dynamic_cast<AliAODMCHeader*>(fAOD->FindListObject(AliAODMCHeader::StdBranchName()));
301	if (!header) {
302	AliWarning("No header file found.");
303	}
304	else {
305	rp = header->GetReactionPlaneAngle();
306	}
307
308	Int_t ntracks = mcArray->GetEntriesFast();
309
310	// Track loop: chek how many particles will be accepted
311	Double_t weight = 0;
312	for (Int_t it = 0; it < ntracks; it++) {
313	weight = 1;
314	AliAODMCParticle* particle = (AliAODMCParticle*)mcArray->At(it);
315	if (!particle) {
316	AliError(Form("Could not receive track %d", it));
317	continue;
318	}
319	if (!particle->IsPhysicalPrimary()) continue;
320	if (particle->Charge() == 0) continue;
321	Double_t pT = particle->Pt();
322	Double_t eta = particle->Eta();
323	Double_t phi = particle->Phi();
324	if (TMath::Abs(eta) < 6.) {
325	// Add flow if it is in the argument
326	if (fAddFlow.Length() > 1) {
327	if (fAddFlow.Contains("pt"))
328	weight *= AddptFlow(pT);
329	if (fAddFlow.Contains("pid"))
330	weight *= AddpidFlow(particle->PdgCode());
331	if (fAddFlow.Contains("eta"))
332	weight *= AddetaFlow(eta);
333	if (fAddFlow.Contains("cent"))
334	weight *= fAliceCent4th->Eval(fCent)/fAliceCent4th->Eval(45);
335
336	weight = 202.TMath::Cos((Double_t)fAddOrder(phi-rp));
337	weight += 1;
338	}
339	fdNdedpMC.Fill(eta, phi, weight);
340	}
341	}
342
343	return kTRUE;
344	}
345	//_____________________________________________________________________
346	Double_t AliForwardMCFlowTaskQC::AddptFlow(Double_t pt = 0) const
347	{
348	//
349	// Add pt dependent flow factor
350	// Parameters:
351	// pt: pT parametrization to use
352	//
353	Double_t weight = 0;
354
355	switch(fAddType)
356	{
357	case 1: weight = fAlicePt2nd4050->Eval(pt)0.5+fAlicePt4th4050->Eval(pt)0.5;
358	break;
359	case 2: weight = fAlicePt2nd4050->Eval(pt);
360	break;
361	case 3: weight = fAlicePt4th3040->Eval(pt);
362	break;
363	case 4: weight = fAlicePt4th4050->Eval(pt);
364	break;
365	}
366
367	return weight;
368	}
369	//_____________________________________________________________________
370	Double_t AliForwardMCFlowTaskQC::AddpidFlow(Int_t id = 0) const
371	{
372	//
373	// Add pid dependent flow factor
374	// Parameters:
375	// id: choose PID dependent setup
376	//
377	Double_t weight = 0;
378
379	switch(fAddType)
380	{
381	case 1: if (TMath::Abs(id) == 211) // pion flow
382	weight = 1.3;
383	else if (TMath::Abs(id) == 2212) // proton flow
384	weight = 1.;
385	else
386	weight = 0.7;
387	break;
388	case 2: weight = 1.207;
389	break;
390	}
391
392	return weight;
393	}
394	//_____________________________________________________________________
395	Double_t AliForwardMCFlowTaskQC::AddetaFlow(Double_t eta = 0) const
396	{
397	//
398	// Add eta dependent flow factor
399	// Parameters:
400	// eta: choose v_n(eta) shape
401	//
402	Double_t weight = 0;
403
404	TF1 gaus = TF1("gaus", "gaus", -6, 6);
405
406	switch(fAddType)
407	{
408	case 1: gaus.SetParameters(0.1, 0., 9);
409	break;
410	case 2: gaus.SetParameters(0.1, 0., 3);
411	break;
412	case 3: gaus.SetParameters(0.1, 0., 15);
413	break;
414	}
415
416	weight = gaus.Eval(eta);
417
418	return weight;
419	}
420	//_____________________________________________________________________
421	Double_t AliForwardMCFlowTaskQC::GetCentFromB() const
422	{
423	//
424	// Get centrality from MC impact parameter.
425	// Values taken from: https://twiki.cern.ch/twiki/bin/viewauth/ALICE/CentStudies
426	//
427	Double_t cent = -1.;
428	Double_t b = -1.;
429	AliAODMCHeader* header = (AliAODMCHeader*)fAOD->FindListObject(AliAODMCHeader::StdBranchName());
430	if (!header) return cent;
431	b = header->GetImpactParameter();
432
433	cent = fImpactParToCent->Eval(b);
434
435	return cent;
436	}
437	//_____________________________________________________________________
438	//
439	//
440	// EOF
441