[u/mrichter/AliRoot.git] / PWGPP / AliAnaFwdDetsQA.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

//------------------------------
// Analysis task for quality-assurance
// of forward detectors ESD
//
// 12/06/2009 cvetan.cheshkov@cern.ch
//------------------------------

#include "TChain.h"
#include "TROOT.h"
#include "TFile.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TF1.h"
#include "TCanvas.h"
#include "TVector3.h"
#include "TParticle.h"
#include "AliVParticle.h"
#include "AliMCParticle.h"
#include "AliESDEvent.h"
#include "AliESDv0.h"
#include "AliESDcascade.h"
#include "AliESDMuonTrack.h"
#include "AliESDCaloCluster.h"
#include "AliRun.h"
#include "AliMCEvent.h"
#include "AliStack.h"
#include "AliGenEventHeader.h"
#include "AliPID.h"
#include "AliESDVZERO.h"

#include "AliAnaFwdDetsQA.h"

ClassImp(AliAnaFwdDetsQA)

AliAnaFwdDetsQA::AliAnaFwdDetsQA():
AliAnalysisTaskSE("AliAnaFwdDetsQA"),
  fListOfHistos(0),
  fT0vtxRec(0),
  fT0vtxRecGen(0),
  fT0time(0),
  fT0time2(0),
  fT0mult(0),
  fT0vtxRes(0),
  fT0ampl(0),
  fV0a(0),
  fV0c(0),
  fV0multA(0),
  fV0multC(0),
  fV0multAcorr(0),
  fV0multCcorr(0),
  fV0Acorr(0),
  fV0Ccorr(0),
  fV0ampl(0)
{
  // Default constructor
  // Define input and output slots here
  // Input slot #0 works with a TChain
  DefineInput(0, TChain::Class());
  // Output slot #1 TList
  DefineOutput(1, TList::Class());
}

AliAnaFwdDetsQA::AliAnaFwdDetsQA(const char* name):
AliAnalysisTaskSE(name),
  fListOfHistos(0),
  fT0vtxRec(0),
  fT0vtxRecGen(0),
  fT0time(0),
  fT0time2(0),
  fT0mult(0),
  fT0vtxRes(0),
  fT0ampl(0),
  fV0a(0),
  fV0c(0),
  fV0multA(0),
  fV0multC(0),
  fV0multAcorr(0),
  fV0multCcorr(0),
  fV0Acorr(0),
  fV0Ccorr(0),
  fV0ampl(0)
{
  // Constructor
  // Define input and output slots here
  // Input slot #0 works with a TChain
  AliInfo("Constructor AliAnaFwdDetsQA");
  DefineInput(0, TChain::Class());
  // Output slot #1 TList
  DefineOutput(1, TList::Class());
}

TH1F * AliAnaFwdDetsQA::CreateHisto(const char* name, const char* title,Int_t nBins, 
					    Double_t xMin, Double_t xMax,
					    const char* xLabel, const char* yLabel)
{
  // helper method which can be used
  // in order to create a histogram
  TH1F* result = new TH1F(name, title, nBins, xMin, xMax);
  result->SetOption("E");
  if (xLabel) result->GetXaxis()->SetTitle(xLabel);
  if (yLabel) result->GetYaxis()->SetTitle(yLabel);
  result->SetMarkerStyle(kFullCircle);
  return result;
}

TH1F *AliAnaFwdDetsQA::CreateEffHisto(const TH1F* hGen, const TH1F* hRec)
{
  // helper method which can be used
  // in order create an efficiency histogram
  Int_t nBins = hGen->GetNbinsX();
  TH1F* hEff = (TH1F*) hGen->Clone("hEff");
  hEff->SetTitle("");
  hEff->SetStats(kFALSE);
  hEff->SetMinimum(0.);
  hEff->SetMaximum(110.);
  hEff->GetYaxis()->SetTitle("#epsilon [%]");
  
  for (Int_t iBin = 0; iBin <= nBins; iBin++) {
    Double_t nGenEff = hGen->GetBinContent(iBin);
    Double_t nRecEff = hRec->GetBinContent(iBin);
    if (nGenEff > 0) {
      Double_t eff = nRecEff/nGenEff;
      hEff->SetBinContent(iBin, 100. * eff);
      Double_t error = sqrt(eff*(1.-eff) / nGenEff);
      if (error < 1e-12) error = 0.0001;
      hEff->SetBinError(iBin, 100. * error);			
    }
    else {
      hEff->SetBinContent(iBin, -100.);
      hEff->SetBinError(iBin, 0);
    }
  }
  return hEff;
}


Bool_t AliAnaFwdDetsQA::FitHisto(TH1* histo, Double_t& res, Double_t& resError)
{
  // fit a gaussian to
  // a histogram
  static TF1* fitFunc = new TF1("fitFunc", "gaus");
  fitFunc->SetLineWidth(2);
  fitFunc->SetFillStyle(0);
  Double_t maxFitRange = 2;
  
  if (histo->Integral() > 50) {
    Float_t mean = histo->GetMean();
    Float_t rms = histo->GetRMS();
    fitFunc->SetRange(mean - maxFitRange*rms, mean + maxFitRange*rms);
    fitFunc->SetParameters(mean, rms);
    histo->Fit(fitFunc, "QRI0");
    histo->GetFunction("fitFunc")->ResetBit(1<<9);
    res = TMath::Abs(fitFunc->GetParameter(2));
    resError = TMath::Abs(fitFunc->GetParError(2));
    return kTRUE;
  }
  return kFALSE;
}

void AliAnaFwdDetsQA::UserCreateOutputObjects()
{
  // Create histograms
  // Create output container
  AliInfo("AliAnaFwdDetsQA::UserCreateOutputObjects");
  fListOfHistos = new TList();
  
  fT0vtxRec = CreateHisto("hT0vtxRec", "Z vertex reconstructed with T0", 100, -25, 25, "Z_{vtx} [cm]", "");
  fT0time   = CreateHisto("hT0time", "Time0 reconstruction with T0", 5000, 10000, 20000, "t_{0} [ps]", "");
  fT0time2  = CreateHisto("hT0time2", "Time0 reconstruction with T0 (mult > 10)", 5000, 10000, 20000, "t_{0} [ps]", "");
  fT0mult   = CreateHisto("hT0mult", "Total reconstructed multiplicity in T0", 100, -0.5, 99.5, "Multiplicity", "");
  fT0vtxRes = CreateHisto("hT0vtxRes", "T0 Z vertex resolution", 100, -25, 25, "Delta(Z_{vtx}) [cm]", "");
  fT0ampl   = CreateHisto("hT0ampl","T0 multiplicity in single channel (all T0 channels)",400,-0.5,99.5);
  
  fT0vtxRecGen = new TH2F("hT0vtxRecGen", "T0 reconstructed vs generated Z vertex", 100, -25, 25, 100, -25, 25);
  fT0vtxRecGen->GetXaxis()->SetTitle("Generated Z vertex");
  fT0vtxRecGen->GetYaxis()->SetTitle("Reconstructed Z vertex");
  fT0vtxRecGen->SetMarkerStyle(kFullCircle);
  fT0vtxRecGen->SetMarkerSize(0.4);

  fV0a = CreateHisto("hV0a","Number of fired PMTs (V0A)",65,-0.5,64.5);
  fV0c = CreateHisto("hV0c","Number of fired PMTs (V0C)",65,-0.5,64.5);
  fV0multA = CreateHisto("hV0multA","Total reconstructed multiplicity (V0A)",100,0.,1000.);
  fV0multC = CreateHisto("hV0multC","Total reconstructed multiplicity (V0C)",100,0.,1000.);
  fV0ampl  = CreateHisto("hV0ampl","V0 multiplicity in single channel (all V0 channels)",400,-0.5,99.5);

  fV0multAcorr = new TH2F("hV0multAcorr", "Reconstructed vs generated (primaries only) multiplicity (V0A)",100,0.,500.,100,0.,1000.);
  fV0multAcorr->GetXaxis()->SetTitle("# of primaries in V0A acceptance");
  fV0multAcorr->GetYaxis()->SetTitle("Reconstructed mutiplicity in V0A");
  fV0multCcorr = new TH2F("hV0multCcorr", "Reconstructed vs generated (primaries only) multiplicity (V0C)",100,0.,500.,100,0.,1000.);
  fV0multCcorr->GetXaxis()->SetTitle("# of primaries in V0C acceptance");
  fV0multCcorr->GetYaxis()->SetTitle("Reconstructed mutiplicity in V0C");

  fV0Acorr = new TH2F("hV0Acorr", "Number of fired PMTs vs generated (primaries only) multiplicity (V0A)",100,0.,500.,65,-0.5,64.5);
  fV0Acorr->GetXaxis()->SetTitle("# of primaries in V0A acceptance");
  fV0Acorr->GetYaxis()->SetTitle("Number of fired PMTs in V0A");
  fV0Ccorr = new TH2F("hV0Ccorr", "Number of fired PMTs vs generated (primaries only) multiplicity (V0C)",100,0.,500.,65,-0.5,64.5);
  fV0Ccorr->GetXaxis()->SetTitle("# of primaries in V0C acceptance");
  fV0Ccorr->GetYaxis()->SetTitle("Number of fired PMTs in V0C");

  fListOfHistos->Add(fT0vtxRec);
  fListOfHistos->Add(fT0time);
  fListOfHistos->Add(fT0mult);
  fListOfHistos->Add(fT0vtxRecGen);
  fListOfHistos->Add(fT0vtxRes);
  fListOfHistos->Add(fV0a);
  fListOfHistos->Add(fV0c);
  fListOfHistos->Add(fV0multA);
  fListOfHistos->Add(fV0multC);
  fListOfHistos->Add(fV0multAcorr);
  fListOfHistos->Add(fV0multCcorr);
  fListOfHistos->Add(fV0Acorr);
  fListOfHistos->Add(fV0Ccorr);
  fListOfHistos->Add(fT0ampl);
  fListOfHistos->Add(fV0ampl);
  fListOfHistos->Add(fT0time2);
}

void AliAnaFwdDetsQA::UserExec(Option_t */*option*/)
{
  // The analysis code
  // goes here
  AliMCEvent* mcEvent = MCEvent();
  if (!mcEvent) {
    Printf("ERROR: Could not retrieve MC event");
    return;
  }	
	
  //Primary vertex needed
  TArrayF mcvtx(3);  
  mcEvent->GenEventHeader()->PrimaryVertex(mcvtx);

  AliStack *stack = mcEvent->Stack();
  if (!stack) {
    Printf("ERROR: Could not retrieve MC stack");
    return;
  }

  Int_t nV0A = 0;
  Int_t nV0C = 0;
  for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {//Check this loop again
    if (!stack->IsPhysicalPrimary(iTracks)) continue;
    AliMCParticle* track = (AliMCParticle*)mcEvent->GetTrack(iTracks);
    TParticle* particle = track->Particle();
    if (!particle) continue;
    if (track->Charge() == 0) continue;
    Double_t eta = particle->Eta();
    if (eta > 2.8 && eta < 5.1) {
      nV0A++;
    }
    if (eta > -3.7 && eta < -1.7) {
      nV0C++;
    }
  }

  // ESD information
  AliVEvent* event = InputEvent();
  if (!event) {
    Printf("ERROR: Could not retrieve event");
    return;
  }

  AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
  if (!esd) {
    Printf("ERROR: Could not retrieve esd");
    return;
  }
  const AliESDTZERO* esdT0 = esd->GetESDTZERO();
  Double_t t0zvtx = esdT0->GetT0zVertex();
  Double_t t0time = esdT0->GetT0();

  fT0vtxRec->Fill(t0zvtx);
  fT0time->Fill(t0time);
  fT0vtxRecGen->Fill(mcvtx[2],t0zvtx);
  t0zvtx *= -1.0;
  fT0vtxRes->Fill(t0zvtx - mcvtx[2]);

  Double_t t0mult = 0;
  for(Int_t i = 0; i < 24; i++) {
    t0mult += esdT0->GetT0amplitude()[i];
    fT0ampl->Fill(esdT0->GetT0amplitude()[i]);
  }

  fT0mult->Fill(t0mult);
  if (t0mult > 10)
    fT0time2->Fill(t0time);

  AliESDVZERO* esdV0 = esd->GetVZEROData();
  fV0a->Fill(esdV0->GetNbPMV0A());
  fV0c->Fill(esdV0->GetNbPMV0C());
  fV0multA->Fill(esdV0->GetMTotV0A());
  fV0multC->Fill(esdV0->GetMTotV0C());

  fV0multAcorr->Fill((Float_t)nV0A,esdV0->GetMTotV0A());
  fV0multCcorr->Fill((Float_t)nV0C,esdV0->GetMTotV0C());
  fV0Acorr->Fill((Float_t)nV0A,(Float_t)esdV0->GetNbPMV0A());
  fV0Ccorr->Fill((Float_t)nV0C,(Float_t)esdV0->GetNbPMV0C());

  for(Int_t i = 0; i < 64; i++) {
    fV0ampl->Fill(esdV0->GetMultiplicity(i));
  }
  // Post output data.
  PostData(1, fListOfHistos);
}

void AliAnaFwdDetsQA::Terminate(Option_t *)
{
  // Terminate
  // Store the output histos
  fListOfHistos = dynamic_cast<TList*>(GetOutputData(1));
  if (!fListOfHistos) {
    Printf("ERROR: fListOfHistos not available");
    return;
  }
	
  //delete esd;
  Info("AliAnaFwdDetsQA", "Successfully finished");
}
Commit	Line	Data
684b93f0	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	//------------------------------
	17	// Analysis task for quality-assurance
	18	// of forward detectors ESD
	19	//
	20	// 12/06/2009 cvetan.cheshkov@cern.ch
	21	//------------------------------
	22
	23	#include "TChain.h"
	24	#include "TROOT.h"
	25	#include "TFile.h"
	26	#include "TH1F.h"
	27	#include "TH2F.h"
	28	#include "TF1.h"
	29	#include "TCanvas.h"
	30	#include "TVector3.h"
	31	#include "TParticle.h"
8b04dae1	32	#include "AliVParticle.h"
8b04dae1	33	#include "AliMCParticle.h"
684b93f0	34	#include "AliESDEvent.h"
	35	#include "AliESDv0.h"
	36	#include "AliESDcascade.h"
	37	#include "AliESDMuonTrack.h"
	38	#include "AliESDCaloCluster.h"
	39	#include "AliRun.h"
	40	#include "AliMCEvent.h"
	41	#include "AliStack.h"
	42	#include "AliGenEventHeader.h"
	43	#include "AliPID.h"
	44	#include "AliESDVZERO.h"
	45
	46	#include "AliAnaFwdDetsQA.h"
	47
	48	ClassImp(AliAnaFwdDetsQA)
	49
	50	AliAnaFwdDetsQA::AliAnaFwdDetsQA():
	51	AliAnalysisTaskSE("AliAnaFwdDetsQA"),
	52	fListOfHistos(0),
	53	fT0vtxRec(0),
	54	fT0vtxRecGen(0),
	55	fT0time(0),
c543bd05	56	fT0time2(0),
684b93f0	57	fT0mult(0),
684b93f0	58	fT0vtxRes(0),
c543bd05	59	fT0ampl(0),
684b93f0	60	fV0a(0),
	61	fV0c(0),
	62	fV0multA(0),
	63	fV0multC(0),
	64	fV0multAcorr(0),
	65	fV0multCcorr(0),
	66	fV0Acorr(0),
c543bd05	67	fV0Ccorr(0),
c543bd05	68	fV0ampl(0)
684b93f0	69	{
	70	// Default constructor
	71	// Define input and output slots here
	72	// Input slot #0 works with a TChain
	73	DefineInput(0, TChain::Class());
	74	// Output slot #1 TList
	75	DefineOutput(1, TList::Class());
	76	}
	77
	78	AliAnaFwdDetsQA::AliAnaFwdDetsQA(const char* name):
	79	AliAnalysisTaskSE(name),
	80	fListOfHistos(0),
	81	fT0vtxRec(0),
	82	fT0vtxRecGen(0),
	83	fT0time(0),
c543bd05	84	fT0time2(0),
684b93f0	85	fT0mult(0),
684b93f0	86	fT0vtxRes(0),
c543bd05	87	fT0ampl(0),
684b93f0	88	fV0a(0),
	89	fV0c(0),
	90	fV0multA(0),
	91	fV0multC(0),
	92	fV0multAcorr(0),
	93	fV0multCcorr(0),
	94	fV0Acorr(0),
c543bd05	95	fV0Ccorr(0),
c543bd05	96	fV0ampl(0)
684b93f0	97	{
684b93f0	98	// Constructor
684b93f0	99	// Define input and output slots here
684b93f0	100	// Input slot #0 works with a TChain
fc6c8a4a	101	AliInfo("Constructor AliAnaFwdDetsQA");
684b93f0	102	DefineInput(0, TChain::Class());
	103	// Output slot #1 TList
	104	DefineOutput(1, TList::Class());
	105	}
	106
	107	TH1F * AliAnaFwdDetsQA::CreateHisto(const char* name, const char* title,Int_t nBins,
	108	Double_t xMin, Double_t xMax,
	109	const char* xLabel, const char* yLabel)
	110	{
fc6c8a4a	111	// helper method which can be used
fc6c8a4a	112	// in order to create a histogram
684b93f0	113	TH1F* result = new TH1F(name, title, nBins, xMin, xMax);
	114	result->SetOption("E");
	115	if (xLabel) result->GetXaxis()->SetTitle(xLabel);
	116	if (yLabel) result->GetYaxis()->SetTitle(yLabel);
	117	result->SetMarkerStyle(kFullCircle);
	118	return result;
	119	}
	120
	121	TH1F AliAnaFwdDetsQA::CreateEffHisto(const TH1F hGen, const TH1F* hRec)
	122	{
fc6c8a4a	123	// helper method which can be used
fc6c8a4a	124	// in order create an efficiency histogram
684b93f0	125	Int_t nBins = hGen->GetNbinsX();
	126	TH1F* hEff = (TH1F*) hGen->Clone("hEff");
	127	hEff->SetTitle("");
	128	hEff->SetStats(kFALSE);
	129	hEff->SetMinimum(0.);
	130	hEff->SetMaximum(110.);
	131	hEff->GetYaxis()->SetTitle("#epsilon [%]");
	132
	133	for (Int_t iBin = 0; iBin <= nBins; iBin++) {
	134	Double_t nGenEff = hGen->GetBinContent(iBin);
	135	Double_t nRecEff = hRec->GetBinContent(iBin);
	136	if (nGenEff > 0) {
	137	Double_t eff = nRecEff/nGenEff;
	138	hEff->SetBinContent(iBin, 100. * eff);
	139	Double_t error = sqrt(eff*(1.-eff) / nGenEff);
fc6c8a4a	140	if (error < 1e-12) error = 0.0001;
684b93f0	141	hEff->SetBinError(iBin, 100. * error);
	142	}
	143	else {
	144	hEff->SetBinContent(iBin, -100.);
	145	hEff->SetBinError(iBin, 0);
	146	}
	147	}
	148	return hEff;
	149	}
	150
	151
	152	Bool_t AliAnaFwdDetsQA::FitHisto(TH1* histo, Double_t& res, Double_t& resError)
	153	{
fc6c8a4a	154	// fit a gaussian to
fc6c8a4a	155	// a histogram
684b93f0	156	static TF1* fitFunc = new TF1("fitFunc", "gaus");
	157	fitFunc->SetLineWidth(2);
	158	fitFunc->SetFillStyle(0);
	159	Double_t maxFitRange = 2;
	160
	161	if (histo->Integral() > 50) {
	162	Float_t mean = histo->GetMean();
	163	Float_t rms = histo->GetRMS();
	164	fitFunc->SetRange(mean - maxFitRangerms, mean + maxFitRangerms);
	165	fitFunc->SetParameters(mean, rms);
	166	histo->Fit(fitFunc, "QRI0");
	167	histo->GetFunction("fitFunc")->ResetBit(1<<9);
	168	res = TMath::Abs(fitFunc->GetParameter(2));
	169	resError = TMath::Abs(fitFunc->GetParError(2));
	170	return kTRUE;
	171	}
	172	return kFALSE;
	173	}
	174
	175	void AliAnaFwdDetsQA::UserCreateOutputObjects()
	176	{
	177	// Create histograms
684b93f0	178	// Create output container
fc6c8a4a	179	AliInfo("AliAnaFwdDetsQA::UserCreateOutputObjects");
684b93f0	180	fListOfHistos = new TList();
	181
	182	fT0vtxRec = CreateHisto("hT0vtxRec", "Z vertex reconstructed with T0", 100, -25, 25, "Z_{vtx} [cm]", "");
	183	fT0time = CreateHisto("hT0time", "Time0 reconstruction with T0", 5000, 10000, 20000, "t_{0} [ps]", "");
c543bd05	184	fT0time2 = CreateHisto("hT0time2", "Time0 reconstruction with T0 (mult > 10)", 5000, 10000, 20000, "t_{0} [ps]", "");
684b93f0	185	fT0mult = CreateHisto("hT0mult", "Total reconstructed multiplicity in T0", 100, -0.5, 99.5, "Multiplicity", "");
684b93f0	186	fT0vtxRes = CreateHisto("hT0vtxRes", "T0 Z vertex resolution", 100, -25, 25, "Delta(Z_{vtx}) [cm]", "");
c543bd05	187	fT0ampl = CreateHisto("hT0ampl","T0 multiplicity in single channel (all T0 channels)",400,-0.5,99.5);
684b93f0	188
	189	fT0vtxRecGen = new TH2F("hT0vtxRecGen", "T0 reconstructed vs generated Z vertex", 100, -25, 25, 100, -25, 25);
	190	fT0vtxRecGen->GetXaxis()->SetTitle("Generated Z vertex");
	191	fT0vtxRecGen->GetYaxis()->SetTitle("Reconstructed Z vertex");
	192	fT0vtxRecGen->SetMarkerStyle(kFullCircle);
	193	fT0vtxRecGen->SetMarkerSize(0.4);
	194
	195	fV0a = CreateHisto("hV0a","Number of fired PMTs (V0A)",65,-0.5,64.5);
	196	fV0c = CreateHisto("hV0c","Number of fired PMTs (V0C)",65,-0.5,64.5);
	197	fV0multA = CreateHisto("hV0multA","Total reconstructed multiplicity (V0A)",100,0.,1000.);
	198	fV0multC = CreateHisto("hV0multC","Total reconstructed multiplicity (V0C)",100,0.,1000.);
c543bd05	199	fV0ampl = CreateHisto("hV0ampl","V0 multiplicity in single channel (all V0 channels)",400,-0.5,99.5);
684b93f0	200
	201	fV0multAcorr = new TH2F("hV0multAcorr", "Reconstructed vs generated (primaries only) multiplicity (V0A)",100,0.,500.,100,0.,1000.);
	202	fV0multAcorr->GetXaxis()->SetTitle("# of primaries in V0A acceptance");
	203	fV0multAcorr->GetYaxis()->SetTitle("Reconstructed mutiplicity in V0A");
	204	fV0multCcorr = new TH2F("hV0multCcorr", "Reconstructed vs generated (primaries only) multiplicity (V0C)",100,0.,500.,100,0.,1000.);
	205	fV0multCcorr->GetXaxis()->SetTitle("# of primaries in V0C acceptance");
	206	fV0multCcorr->GetYaxis()->SetTitle("Reconstructed mutiplicity in V0C");
	207
	208	fV0Acorr = new TH2F("hV0Acorr", "Number of fired PMTs vs generated (primaries only) multiplicity (V0A)",100,0.,500.,65,-0.5,64.5);
	209	fV0Acorr->GetXaxis()->SetTitle("# of primaries in V0A acceptance");
	210	fV0Acorr->GetYaxis()->SetTitle("Number of fired PMTs in V0A");
	211	fV0Ccorr = new TH2F("hV0Ccorr", "Number of fired PMTs vs generated (primaries only) multiplicity (V0C)",100,0.,500.,65,-0.5,64.5);
	212	fV0Ccorr->GetXaxis()->SetTitle("# of primaries in V0C acceptance");
	213	fV0Ccorr->GetYaxis()->SetTitle("Number of fired PMTs in V0C");
	214
	215	fListOfHistos->Add(fT0vtxRec);
	216	fListOfHistos->Add(fT0time);
	217	fListOfHistos->Add(fT0mult);
	218	fListOfHistos->Add(fT0vtxRecGen);
	219	fListOfHistos->Add(fT0vtxRes);
	220	fListOfHistos->Add(fV0a);
	221	fListOfHistos->Add(fV0c);
	222	fListOfHistos->Add(fV0multA);
	223	fListOfHistos->Add(fV0multC);
	224	fListOfHistos->Add(fV0multAcorr);
	225	fListOfHistos->Add(fV0multCcorr);
	226	fListOfHistos->Add(fV0Acorr);
	227	fListOfHistos->Add(fV0Ccorr);
c543bd05	228	fListOfHistos->Add(fT0ampl);
	229	fListOfHistos->Add(fV0ampl);
	230	fListOfHistos->Add(fT0time2);
684b93f0	231	}
	232
	233	void AliAnaFwdDetsQA::UserExec(Option_t /option*/)
	234	{
fc6c8a4a	235	// The analysis code
fc6c8a4a	236	// goes here
684b93f0	237	AliMCEvent* mcEvent = MCEvent();
	238	if (!mcEvent) {
	239	Printf("ERROR: Could not retrieve MC event");
	240	return;
	241	}
	242
	243	//Primary vertex needed
	244	TArrayF mcvtx(3);
	245	mcEvent->GenEventHeader()->PrimaryVertex(mcvtx);
	246
	247	AliStack *stack = mcEvent->Stack();
	248	if (!stack) {
	249	Printf("ERROR: Could not retrieve MC stack");
	250	return;
	251	}
	252
	253	Int_t nV0A = 0;
	254	Int_t nV0C = 0;
	255	for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {//Check this loop again
	256	if (!stack->IsPhysicalPrimary(iTracks)) continue;
8b04dae1	257	AliMCParticle* track = (AliMCParticle*)mcEvent->GetTrack(iTracks);
684b93f0	258	TParticle* particle = track->Particle();
684b93f0	259	if (!particle) continue;
c543bd05	260	if (track->Charge() == 0) continue;
684b93f0	261	Double_t eta = particle->Eta();
	262	if (eta > 2.8 && eta < 5.1) {
	263	nV0A++;
	264	}
	265	if (eta > -3.7 && eta < -1.7) {
	266	nV0C++;
	267	}
	268	}
	269
	270	// ESD information
	271	AliVEvent* event = InputEvent();
	272	if (!event) {
	273	Printf("ERROR: Could not retrieve event");
	274	return;
	275	}
	276
	277	AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
ae02804c	278	if (!esd) {
	279	Printf("ERROR: Could not retrieve esd");
	280	return;
	281	}
684b93f0	282	const AliESDTZERO* esdT0 = esd->GetESDTZERO();
	283	Double_t t0zvtx = esdT0->GetT0zVertex();
	284	Double_t t0time = esdT0->GetT0();
	285
	286	fT0vtxRec->Fill(t0zvtx);
	287	fT0time->Fill(t0time);
	288	fT0vtxRecGen->Fill(mcvtx[2],t0zvtx);
	289	t0zvtx *= -1.0;
	290	fT0vtxRes->Fill(t0zvtx - mcvtx[2]);
	291
	292	Double_t t0mult = 0;
	293	for(Int_t i = 0; i < 24; i++) {
	294	t0mult += esdT0->GetT0amplitude()[i];
c543bd05	295	fT0ampl->Fill(esdT0->GetT0amplitude()[i]);
684b93f0	296	}
	297
	298	fT0mult->Fill(t0mult);
c543bd05	299	if (t0mult > 10)
c543bd05	300	fT0time2->Fill(t0time);
684b93f0	301
	302	AliESDVZERO* esdV0 = esd->GetVZEROData();
	303	fV0a->Fill(esdV0->GetNbPMV0A());
	304	fV0c->Fill(esdV0->GetNbPMV0C());
	305	fV0multA->Fill(esdV0->GetMTotV0A());
	306	fV0multC->Fill(esdV0->GetMTotV0C());
	307
	308	fV0multAcorr->Fill((Float_t)nV0A,esdV0->GetMTotV0A());
	309	fV0multCcorr->Fill((Float_t)nV0C,esdV0->GetMTotV0C());
	310	fV0Acorr->Fill((Float_t)nV0A,(Float_t)esdV0->GetNbPMV0A());
	311	fV0Ccorr->Fill((Float_t)nV0C,(Float_t)esdV0->GetNbPMV0C());
	312
c543bd05	313	for(Int_t i = 0; i < 64; i++) {
	314	fV0ampl->Fill(esdV0->GetMultiplicity(i));
	315	}
684b93f0	316	// Post output data.
	317	PostData(1, fListOfHistos);
	318	}
	319
	320	void AliAnaFwdDetsQA::Terminate(Option_t *)
	321	{
fc6c8a4a	322	// Terminate
fc6c8a4a	323	// Store the output histos
684b93f0	324	fListOfHistos = dynamic_cast<TList*>(GetOutputData(1));
	325	if (!fListOfHistos) {
	326	Printf("ERROR: fListOfHistos not available");
	327	return;
	328	}
	329
684b93f0	330	//delete esd;
	331	Info("AliAnaFwdDetsQA", "Successfully finished");
	332	}
	333
	334
	335
	336
	337
	338
	339
	340
	341
	342
	343
	344
	345
	346
	347
	348
	349
	350
	351
	352
	353
	354
	355
	356
	357
	358
	359
	360
	361
	362
	363
	364
	365
	366
	367
	368
	369
	370
	371
	372
	373
	374
	375
	376
	377
	378
	379
	380
	381
	382
	383
	384
	385
	386
	387
	388
	389
	390
	391
	392
	393
394