[u/mrichter/AliRoot.git] / PWG1 / 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 "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
  AliInfo("Constructor AliAnaFwdDetsQA");
  // 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());
}

TH1F * AliAnaFwdDetsQA::CreateHisto(const char* name, const char* title,Int_t nBins, 
					    Double_t xMin, Double_t xMax,
					    const char* xLabel, const char* yLabel)
{
  // 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)
{
  // 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 == 0) 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
  AliInfo("AliAnaFwdDetsQA::UserCreateOutputObjects");
  // Create output container
  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*/)
{

  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 = 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);
  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 *)
{
  fListOfHistos = dynamic_cast<TList*>(GetOutputData(1));
  if (!fListOfHistos) {
    Printf("ERROR: fListOfHistos not available");
    return;
  }
	
  fT0vtxRec = dynamic_cast<TH1F*>(fListOfHistos->At(0));
  fT0time = dynamic_cast<TH1F*>(fListOfHistos->At(1));
  fT0mult = dynamic_cast<TH1F*>(fListOfHistos->At(2));
  fT0vtxRecGen = dynamic_cast<TH2F*>(fListOfHistos->At(3));
  fT0vtxRes = dynamic_cast<TH1F*>(fListOfHistos->At(4));

  fV0a = dynamic_cast<TH1F*>(fListOfHistos->At(5));
  fV0c = dynamic_cast<TH1F*>(fListOfHistos->At(6));
  fV0multA = dynamic_cast<TH1F*>(fListOfHistos->At(7));
  fV0multC = dynamic_cast<TH1F*>(fListOfHistos->At(8));
  fV0multAcorr = dynamic_cast<TH2F*>(fListOfHistos->At(9));
  fV0multCcorr = dynamic_cast<TH2F*>(fListOfHistos->At(10));
  fV0Acorr = dynamic_cast<TH2F*>(fListOfHistos->At(11));
  fV0Ccorr = dynamic_cast<TH2F*>(fListOfHistos->At(12));

  fT0ampl = dynamic_cast<TH1F*>(fListOfHistos->At(13));
  fV0ampl = dynamic_cast<TH1F*>(fListOfHistos->At(14));
  fT0time2 = dynamic_cast<TH1F*>(fListOfHistos->At(15));


  // draw the histograms if not in batch mode
  if (!gROOT->IsBatch()) {
    new TCanvas;
    fT0vtxRec->DrawCopy("E");
    new TCanvas;
    fT0time->DrawCopy("E");
    new TCanvas;
    fT0mult->DrawCopy("E");
    new TCanvas;
    fT0vtxRes->DrawCopy("E");
    new TCanvas;
    fT0vtxRecGen->DrawCopy();
    new TCanvas;
    fV0a->DrawCopy("E");
    new TCanvas;
    fV0c->DrawCopy("E");
    new TCanvas;
    fV0multA->DrawCopy("E");
    new TCanvas;
    fV0multC->DrawCopy("E");
    new TCanvas;
    fV0multAcorr->DrawCopy();
    new TCanvas;
    fV0multCcorr->DrawCopy();
    new TCanvas;
    fV0Acorr->DrawCopy();
    new TCanvas;
    fV0Ccorr->DrawCopy();
    new TCanvas;
    fT0ampl->DrawCopy("E");
    new TCanvas;
    fV0ampl->DrawCopy("E");
    new TCanvas;
    fT0time2->DrawCopy("E");
  }

  // write the output histograms to a file
  TFile* outputFile = TFile::Open("FwdDetsQA.root", "recreate");
  if (!outputFile || !outputFile->IsOpen())
    {
      Error("AliAnaFwdDetsQA", "opening output file FwdDetsQA.root failed");
      return;
    }
  fT0vtxRec->Write();
  fT0time->Write();
  fT0mult->Write();
  fT0vtxRecGen->Write();
  fT0vtxRes->Write();
  fT0ampl->Write();
  fT0time2->Write();

  fV0a->Write();
  fV0c->Write();
  fV0multA->Write();
  fV0multC->Write();
  fV0multAcorr->Write();
  fV0multCcorr->Write();
  fV0Acorr->Write();
  fV0Ccorr->Write();
  fV0ampl->Write();

  outputFile->Close();
  delete outputFile;

  //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"
	32	#include "AliESDEvent.h"
	33	#include "AliESDv0.h"
	34	#include "AliESDcascade.h"
	35	#include "AliESDMuonTrack.h"
	36	#include "AliESDCaloCluster.h"
	37	#include "AliRun.h"
	38	#include "AliMCEvent.h"
	39	#include "AliStack.h"
	40	#include "AliGenEventHeader.h"
	41	#include "AliPID.h"
	42	#include "AliESDVZERO.h"
	43
	44	#include "AliAnaFwdDetsQA.h"
	45
	46	ClassImp(AliAnaFwdDetsQA)
	47
	48	AliAnaFwdDetsQA::AliAnaFwdDetsQA():
	49	AliAnalysisTaskSE("AliAnaFwdDetsQA"),
	50	fListOfHistos(0),
	51	fT0vtxRec(0),
	52	fT0vtxRecGen(0),
	53	fT0time(0),
c543bd05	54	fT0time2(0),
684b93f0	55	fT0mult(0),
684b93f0	56	fT0vtxRes(0),
c543bd05	57	fT0ampl(0),
684b93f0	58	fV0a(0),
	59	fV0c(0),
	60	fV0multA(0),
	61	fV0multC(0),
	62	fV0multAcorr(0),
	63	fV0multCcorr(0),
	64	fV0Acorr(0),
c543bd05	65	fV0Ccorr(0),
c543bd05	66	fV0ampl(0)
684b93f0	67	{
	68	// Default constructor
	69	// Define input and output slots here
	70	// Input slot #0 works with a TChain
	71	DefineInput(0, TChain::Class());
	72	// Output slot #1 TList
	73	DefineOutput(1, TList::Class());
	74	}
	75
	76	AliAnaFwdDetsQA::AliAnaFwdDetsQA(const char* name):
	77	AliAnalysisTaskSE(name),
	78	fListOfHistos(0),
	79	fT0vtxRec(0),
	80	fT0vtxRecGen(0),
	81	fT0time(0),
c543bd05	82	fT0time2(0),
684b93f0	83	fT0mult(0),
684b93f0	84	fT0vtxRes(0),
c543bd05	85	fT0ampl(0),
684b93f0	86	fV0a(0),
	87	fV0c(0),
	88	fV0multA(0),
	89	fV0multC(0),
	90	fV0multAcorr(0),
	91	fV0multCcorr(0),
	92	fV0Acorr(0),
c543bd05	93	fV0Ccorr(0),
c543bd05	94	fV0ampl(0)
684b93f0	95	{
	96	// Constructor
	97	AliInfo("Constructor AliAnaFwdDetsQA");
	98	// Define input and output slots here
	99	// Input slot #0 works with a TChain
	100	DefineInput(0, TChain::Class());
	101	// Output slot #1 TList
	102	DefineOutput(1, TList::Class());
	103	}
	104
	105	TH1F * AliAnaFwdDetsQA::CreateHisto(const char* name, const char* title,Int_t nBins,
	106	Double_t xMin, Double_t xMax,
	107	const char* xLabel, const char* yLabel)
	108	{
	109	// create a histogram
	110	TH1F* result = new TH1F(name, title, nBins, xMin, xMax);
	111	result->SetOption("E");
	112	if (xLabel) result->GetXaxis()->SetTitle(xLabel);
	113	if (yLabel) result->GetYaxis()->SetTitle(yLabel);
	114	result->SetMarkerStyle(kFullCircle);
	115	return result;
	116	}
	117
	118	TH1F AliAnaFwdDetsQA::CreateEffHisto(const TH1F hGen, const TH1F* hRec)
	119	{
	120	// create an efficiency histogram
	121	Int_t nBins = hGen->GetNbinsX();
	122	TH1F* hEff = (TH1F*) hGen->Clone("hEff");
	123	hEff->SetTitle("");
	124	hEff->SetStats(kFALSE);
	125	hEff->SetMinimum(0.);
	126	hEff->SetMaximum(110.);
	127	hEff->GetYaxis()->SetTitle("#epsilon [%]");
	128
	129	for (Int_t iBin = 0; iBin <= nBins; iBin++) {
	130	Double_t nGenEff = hGen->GetBinContent(iBin);
	131	Double_t nRecEff = hRec->GetBinContent(iBin);
	132	if (nGenEff > 0) {
	133	Double_t eff = nRecEff/nGenEff;
	134	hEff->SetBinContent(iBin, 100. * eff);
	135	Double_t error = sqrt(eff*(1.-eff) / nGenEff);
	136	if (error == 0) error = 0.0001;
	137	hEff->SetBinError(iBin, 100. * error);
	138	}
	139	else {
	140	hEff->SetBinContent(iBin, -100.);
	141	hEff->SetBinError(iBin, 0);
	142	}
	143	}
	144	return hEff;
	145	}
	146
	147
	148	Bool_t AliAnaFwdDetsQA::FitHisto(TH1* histo, Double_t& res, Double_t& resError)
	149	{
	150	// fit a gaussian to a histogram
	151	static TF1* fitFunc = new TF1("fitFunc", "gaus");
	152	fitFunc->SetLineWidth(2);
	153	fitFunc->SetFillStyle(0);
	154	Double_t maxFitRange = 2;
	155
	156	if (histo->Integral() > 50) {
	157	Float_t mean = histo->GetMean();
	158	Float_t rms = histo->GetRMS();
159	fitFunc->SetRange(mean - maxFitRangerms, mean + maxFitRangerms);
160	fitFunc->SetParameters(mean, rms);
161	histo->Fit(fitFunc, "QRI0");
162	histo->GetFunction("fitFunc")->ResetBit(1<<9);
163	res = TMath::Abs(fitFunc->GetParameter(2));
164	resError = TMath::Abs(fitFunc->GetParError(2));
165	return kTRUE;
166	}
167	return kFALSE;
168	}
169
170	void AliAnaFwdDetsQA::UserCreateOutputObjects()
171	{
172	// Create histograms
173	AliInfo("AliAnaFwdDetsQA::UserCreateOutputObjects");
174	// Create output container
175	fListOfHistos = new TList();
176
177	fT0vtxRec = CreateHisto("hT0vtxRec", "Z vertex reconstructed with T0", 100, -25, 25, "Z_{vtx} [cm]", "");
178	fT0time = CreateHisto("hT0time", "Time0 reconstruction with T0", 5000, 10000, 20000, "t_{0} [ps]", "");
c543bd05	179	fT0time2 = CreateHisto("hT0time2", "Time0 reconstruction with T0 (mult > 10)", 5000, 10000, 20000, "t_{0} [ps]", "");
684b93f0	180	fT0mult = CreateHisto("hT0mult", "Total reconstructed multiplicity in T0", 100, -0.5, 99.5, "Multiplicity", "");
684b93f0	181	fT0vtxRes = CreateHisto("hT0vtxRes", "T0 Z vertex resolution", 100, -25, 25, "Delta(Z_{vtx}) [cm]", "");
c543bd05	182	fT0ampl = CreateHisto("hT0ampl","T0 multiplicity in single channel (all T0 channels)",400,-0.5,99.5);
684b93f0	183
	184	fT0vtxRecGen = new TH2F("hT0vtxRecGen", "T0 reconstructed vs generated Z vertex", 100, -25, 25, 100, -25, 25);
	185	fT0vtxRecGen->GetXaxis()->SetTitle("Generated Z vertex");
	186	fT0vtxRecGen->GetYaxis()->SetTitle("Reconstructed Z vertex");
	187	fT0vtxRecGen->SetMarkerStyle(kFullCircle);
	188	fT0vtxRecGen->SetMarkerSize(0.4);
	189
	190	fV0a = CreateHisto("hV0a","Number of fired PMTs (V0A)",65,-0.5,64.5);
	191	fV0c = CreateHisto("hV0c","Number of fired PMTs (V0C)",65,-0.5,64.5);
	192	fV0multA = CreateHisto("hV0multA","Total reconstructed multiplicity (V0A)",100,0.,1000.);
	193	fV0multC = CreateHisto("hV0multC","Total reconstructed multiplicity (V0C)",100,0.,1000.);
c543bd05	194	fV0ampl = CreateHisto("hV0ampl","V0 multiplicity in single channel (all V0 channels)",400,-0.5,99.5);
684b93f0	195
	196	fV0multAcorr = new TH2F("hV0multAcorr", "Reconstructed vs generated (primaries only) multiplicity (V0A)",100,0.,500.,100,0.,1000.);
	197	fV0multAcorr->GetXaxis()->SetTitle("# of primaries in V0A acceptance");
	198	fV0multAcorr->GetYaxis()->SetTitle("Reconstructed mutiplicity in V0A");
	199	fV0multCcorr = new TH2F("hV0multCcorr", "Reconstructed vs generated (primaries only) multiplicity (V0C)",100,0.,500.,100,0.,1000.);
	200	fV0multCcorr->GetXaxis()->SetTitle("# of primaries in V0C acceptance");
	201	fV0multCcorr->GetYaxis()->SetTitle("Reconstructed mutiplicity in V0C");
	202
	203	fV0Acorr = new TH2F("hV0Acorr", "Number of fired PMTs vs generated (primaries only) multiplicity (V0A)",100,0.,500.,65,-0.5,64.5);
	204	fV0Acorr->GetXaxis()->SetTitle("# of primaries in V0A acceptance");
	205	fV0Acorr->GetYaxis()->SetTitle("Number of fired PMTs in V0A");
	206	fV0Ccorr = new TH2F("hV0Ccorr", "Number of fired PMTs vs generated (primaries only) multiplicity (V0C)",100,0.,500.,65,-0.5,64.5);
	207	fV0Ccorr->GetXaxis()->SetTitle("# of primaries in V0C acceptance");
	208	fV0Ccorr->GetYaxis()->SetTitle("Number of fired PMTs in V0C");
	209
	210	fListOfHistos->Add(fT0vtxRec);
	211	fListOfHistos->Add(fT0time);
	212	fListOfHistos->Add(fT0mult);
	213	fListOfHistos->Add(fT0vtxRecGen);
	214	fListOfHistos->Add(fT0vtxRes);
	215	fListOfHistos->Add(fV0a);
	216	fListOfHistos->Add(fV0c);
	217	fListOfHistos->Add(fV0multA);
	218	fListOfHistos->Add(fV0multC);
	219	fListOfHistos->Add(fV0multAcorr);
	220	fListOfHistos->Add(fV0multCcorr);
	221	fListOfHistos->Add(fV0Acorr);
	222	fListOfHistos->Add(fV0Ccorr);
c543bd05	223	fListOfHistos->Add(fT0ampl);
	224	fListOfHistos->Add(fV0ampl);
	225	fListOfHistos->Add(fT0time2);
684b93f0	226	}
	227
	228	void AliAnaFwdDetsQA::UserExec(Option_t /option*/)
	229	{
	230
	231	AliMCEvent* mcEvent = MCEvent();
	232	if (!mcEvent) {
	233	Printf("ERROR: Could not retrieve MC event");
	234	return;
	235	}
	236
	237	//Primary vertex needed
	238	TArrayF mcvtx(3);
	239	mcEvent->GenEventHeader()->PrimaryVertex(mcvtx);
	240
	241	AliStack *stack = mcEvent->Stack();
	242	if (!stack) {
	243	Printf("ERROR: Could not retrieve MC stack");
	244	return;
	245	}
	246
	247	Int_t nV0A = 0;
	248	Int_t nV0C = 0;
	249	for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {//Check this loop again
	250	if (!stack->IsPhysicalPrimary(iTracks)) continue;
	251	AliMCParticle* track = mcEvent->GetTrack(iTracks);
	252	TParticle* particle = track->Particle();
	253	if (!particle) continue;
c543bd05	254	if (track->Charge() == 0) continue;
684b93f0	255	Double_t eta = particle->Eta();
	256	if (eta > 2.8 && eta < 5.1) {
	257	nV0A++;
	258	}
	259	if (eta > -3.7 && eta < -1.7) {
	260	nV0C++;
	261	}
	262	}
	263
	264	// ESD information
	265	AliVEvent* event = InputEvent();
	266	if (!event) {
	267	Printf("ERROR: Could not retrieve event");
	268	return;
	269	}
	270
	271	AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
	272	const AliESDTZERO* esdT0 = esd->GetESDTZERO();
	273	Double_t t0zvtx = esdT0->GetT0zVertex();
	274	Double_t t0time = esdT0->GetT0();
	275
	276	fT0vtxRec->Fill(t0zvtx);
	277	fT0time->Fill(t0time);
	278	fT0vtxRecGen->Fill(mcvtx[2],t0zvtx);
	279	t0zvtx *= -1.0;
	280	fT0vtxRes->Fill(t0zvtx - mcvtx[2]);
	281
	282	Double_t t0mult = 0;
	283	for(Int_t i = 0; i < 24; i++) {
	284	t0mult += esdT0->GetT0amplitude()[i];
c543bd05	285	fT0ampl->Fill(esdT0->GetT0amplitude()[i]);
684b93f0	286	}
	287
	288	fT0mult->Fill(t0mult);
c543bd05	289	if (t0mult > 10)
c543bd05	290	fT0time2->Fill(t0time);
684b93f0	291
	292	AliESDVZERO* esdV0 = esd->GetVZEROData();
	293	fV0a->Fill(esdV0->GetNbPMV0A());
	294	fV0c->Fill(esdV0->GetNbPMV0C());
	295	fV0multA->Fill(esdV0->GetMTotV0A());
	296	fV0multC->Fill(esdV0->GetMTotV0C());
	297
	298	fV0multAcorr->Fill((Float_t)nV0A,esdV0->GetMTotV0A());
	299	fV0multCcorr->Fill((Float_t)nV0C,esdV0->GetMTotV0C());
	300	fV0Acorr->Fill((Float_t)nV0A,(Float_t)esdV0->GetNbPMV0A());
	301	fV0Ccorr->Fill((Float_t)nV0C,(Float_t)esdV0->GetNbPMV0C());
	302
c543bd05	303	for(Int_t i = 0; i < 64; i++) {
	304	fV0ampl->Fill(esdV0->GetMultiplicity(i));
	305	}
684b93f0	306	// Post output data.
	307	PostData(1, fListOfHistos);
	308	}
	309
	310	void AliAnaFwdDetsQA::Terminate(Option_t *)
	311	{
	312	fListOfHistos = dynamic_cast<TList*>(GetOutputData(1));
	313	if (!fListOfHistos) {
	314	Printf("ERROR: fListOfHistos not available");
	315	return;
	316	}
	317
	318	fT0vtxRec = dynamic_cast<TH1F*>(fListOfHistos->At(0));
	319	fT0time = dynamic_cast<TH1F*>(fListOfHistos->At(1));
	320	fT0mult = dynamic_cast<TH1F*>(fListOfHistos->At(2));
	321	fT0vtxRecGen = dynamic_cast<TH2F*>(fListOfHistos->At(3));
	322	fT0vtxRes = dynamic_cast<TH1F*>(fListOfHistos->At(4));
	323
	324	fV0a = dynamic_cast<TH1F*>(fListOfHistos->At(5));
	325	fV0c = dynamic_cast<TH1F*>(fListOfHistos->At(6));
	326	fV0multA = dynamic_cast<TH1F*>(fListOfHistos->At(7));
	327	fV0multC = dynamic_cast<TH1F*>(fListOfHistos->At(8));
	328	fV0multAcorr = dynamic_cast<TH2F*>(fListOfHistos->At(9));
	329	fV0multCcorr = dynamic_cast<TH2F*>(fListOfHistos->At(10));
	330	fV0Acorr = dynamic_cast<TH2F*>(fListOfHistos->At(11));
	331	fV0Ccorr = dynamic_cast<TH2F*>(fListOfHistos->At(12));
	332
c543bd05	333	fT0ampl = dynamic_cast<TH1F*>(fListOfHistos->At(13));
	334	fV0ampl = dynamic_cast<TH1F*>(fListOfHistos->At(14));
	335	fT0time2 = dynamic_cast<TH1F*>(fListOfHistos->At(15));
	336
	337
684b93f0	338	// draw the histograms if not in batch mode
	339	if (!gROOT->IsBatch()) {
	340	new TCanvas;
	341	fT0vtxRec->DrawCopy("E");
	342	new TCanvas;
	343	fT0time->DrawCopy("E");
	344	new TCanvas;
	345	fT0mult->DrawCopy("E");
	346	new TCanvas;
	347	fT0vtxRes->DrawCopy("E");
	348	new TCanvas;
	349	fT0vtxRecGen->DrawCopy();
	350	new TCanvas;
	351	fV0a->DrawCopy("E");
	352	new TCanvas;
	353	fV0c->DrawCopy("E");
	354	new TCanvas;
	355	fV0multA->DrawCopy("E");
	356	new TCanvas;
	357	fV0multC->DrawCopy("E");
	358	new TCanvas;
	359	fV0multAcorr->DrawCopy();
	360	new TCanvas;
	361	fV0multCcorr->DrawCopy();
	362	new TCanvas;
	363	fV0Acorr->DrawCopy();
	364	new TCanvas;
	365	fV0Ccorr->DrawCopy();
c543bd05	366	new TCanvas;
	367	fT0ampl->DrawCopy("E");
	368	new TCanvas;
	369	fV0ampl->DrawCopy("E");
	370	new TCanvas;
	371	fT0time2->DrawCopy("E");
684b93f0	372	}
	373
	374	// write the output histograms to a file
	375	TFile* outputFile = TFile::Open("FwdDetsQA.root", "recreate");
	376	if (!outputFile \|\| !outputFile->IsOpen())
	377	{
	378	Error("AliAnaFwdDetsQA", "opening output file FwdDetsQA.root failed");
	379	return;
	380	}
	381	fT0vtxRec->Write();
	382	fT0time->Write();
	383	fT0mult->Write();
	384	fT0vtxRecGen->Write();
	385	fT0vtxRes->Write();
c543bd05	386	fT0ampl->Write();
c543bd05	387	fT0time2->Write();
684b93f0	388
	389	fV0a->Write();
	390	fV0c->Write();
	391	fV0multA->Write();
	392	fV0multC->Write();
	393	fV0multAcorr->Write();
	394	fV0multCcorr->Write();
	395	fV0Acorr->Write();
	396	fV0Ccorr->Write();
c543bd05	397	fV0ampl->Write();
684b93f0	398
	399	outputFile->Close();
	400	delete outputFile;
	401
	402	//delete esd;
	403	Info("AliAnaFwdDetsQA", "Successfully finished");
	404	}
	405
	406
	407
	408
	409
	410
	411
	412
	413
	414
	415
	416
	417
	418
	419
	420
	421
	422
	423
	424
	425
	426
	427
	428
	429
	430
	431
	432
	433
	434
	435
	436
	437
	438
	439
	440
	441
	442
	443
	444
	445
	446
	447
	448
	449
	450
	451
	452
	453
	454
	455
	456
	457
	458
	459
	460
	461
462
463
464
465
466