[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 "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);
  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;
  }
	
  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"
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);
	278	const AliESDTZERO* esdT0 = esd->GetESDTZERO();
	279	Double_t t0zvtx = esdT0->GetT0zVertex();
	280	Double_t t0time = esdT0->GetT0();
	281
	282	fT0vtxRec->Fill(t0zvtx);
	283	fT0time->Fill(t0time);
	284	fT0vtxRecGen->Fill(mcvtx[2],t0zvtx);
	285	t0zvtx *= -1.0;
	286	fT0vtxRes->Fill(t0zvtx - mcvtx[2]);
	287
	288	Double_t t0mult = 0;
	289	for(Int_t i = 0; i < 24; i++) {
	290	t0mult += esdT0->GetT0amplitude()[i];
c543bd05	291	fT0ampl->Fill(esdT0->GetT0amplitude()[i]);
684b93f0	292	}
	293
	294	fT0mult->Fill(t0mult);
c543bd05	295	if (t0mult > 10)
c543bd05	296	fT0time2->Fill(t0time);
684b93f0	297
	298	AliESDVZERO* esdV0 = esd->GetVZEROData();
	299	fV0a->Fill(esdV0->GetNbPMV0A());
	300	fV0c->Fill(esdV0->GetNbPMV0C());
	301	fV0multA->Fill(esdV0->GetMTotV0A());
	302	fV0multC->Fill(esdV0->GetMTotV0C());
	303
	304	fV0multAcorr->Fill((Float_t)nV0A,esdV0->GetMTotV0A());
	305	fV0multCcorr->Fill((Float_t)nV0C,esdV0->GetMTotV0C());
	306	fV0Acorr->Fill((Float_t)nV0A,(Float_t)esdV0->GetNbPMV0A());
	307	fV0Ccorr->Fill((Float_t)nV0C,(Float_t)esdV0->GetNbPMV0C());
	308
c543bd05	309	for(Int_t i = 0; i < 64; i++) {
	310	fV0ampl->Fill(esdV0->GetMultiplicity(i));
	311	}
684b93f0	312	// Post output data.
	313	PostData(1, fListOfHistos);
	314	}
	315
	316	void AliAnaFwdDetsQA::Terminate(Option_t *)
	317	{
fc6c8a4a	318	// Terminate
fc6c8a4a	319	// Store the output histos
684b93f0	320	fListOfHistos = dynamic_cast<TList*>(GetOutputData(1));
	321	if (!fListOfHistos) {
	322	Printf("ERROR: fListOfHistos not available");
	323	return;
	324	}
	325
	326	fT0vtxRec = dynamic_cast<TH1F*>(fListOfHistos->At(0));
	327	fT0time = dynamic_cast<TH1F*>(fListOfHistos->At(1));
	328	fT0mult = dynamic_cast<TH1F*>(fListOfHistos->At(2));
	329	fT0vtxRecGen = dynamic_cast<TH2F*>(fListOfHistos->At(3));
	330	fT0vtxRes = dynamic_cast<TH1F*>(fListOfHistos->At(4));
	331
	332	fV0a = dynamic_cast<TH1F*>(fListOfHistos->At(5));
	333	fV0c = dynamic_cast<TH1F*>(fListOfHistos->At(6));
	334	fV0multA = dynamic_cast<TH1F*>(fListOfHistos->At(7));
	335	fV0multC = dynamic_cast<TH1F*>(fListOfHistos->At(8));
	336	fV0multAcorr = dynamic_cast<TH2F*>(fListOfHistos->At(9));
	337	fV0multCcorr = dynamic_cast<TH2F*>(fListOfHistos->At(10));
	338	fV0Acorr = dynamic_cast<TH2F*>(fListOfHistos->At(11));
	339	fV0Ccorr = dynamic_cast<TH2F*>(fListOfHistos->At(12));
	340
c543bd05	341	fT0ampl = dynamic_cast<TH1F*>(fListOfHistos->At(13));
	342	fV0ampl = dynamic_cast<TH1F*>(fListOfHistos->At(14));
	343	fT0time2 = dynamic_cast<TH1F*>(fListOfHistos->At(15));
	344
	345
684b93f0	346	// draw the histograms if not in batch mode
	347	if (!gROOT->IsBatch()) {
	348	new TCanvas;
	349	fT0vtxRec->DrawCopy("E");
	350	new TCanvas;
	351	fT0time->DrawCopy("E");
	352	new TCanvas;
	353	fT0mult->DrawCopy("E");
	354	new TCanvas;
	355	fT0vtxRes->DrawCopy("E");
	356	new TCanvas;
	357	fT0vtxRecGen->DrawCopy();
	358	new TCanvas;
	359	fV0a->DrawCopy("E");
	360	new TCanvas;
	361	fV0c->DrawCopy("E");
	362	new TCanvas;
	363	fV0multA->DrawCopy("E");
	364	new TCanvas;
	365	fV0multC->DrawCopy("E");
	366	new TCanvas;
	367	fV0multAcorr->DrawCopy();
	368	new TCanvas;
	369	fV0multCcorr->DrawCopy();
	370	new TCanvas;
	371	fV0Acorr->DrawCopy();
	372	new TCanvas;
	373	fV0Ccorr->DrawCopy();
c543bd05	374	new TCanvas;
	375	fT0ampl->DrawCopy("E");
	376	new TCanvas;
	377	fV0ampl->DrawCopy("E");
	378	new TCanvas;
	379	fT0time2->DrawCopy("E");
684b93f0	380	}
	381
	382	// write the output histograms to a file
	383	TFile* outputFile = TFile::Open("FwdDetsQA.root", "recreate");
	384	if (!outputFile \|\| !outputFile->IsOpen())
	385	{
	386	Error("AliAnaFwdDetsQA", "opening output file FwdDetsQA.root failed");
	387	return;
	388	}
	389	fT0vtxRec->Write();
	390	fT0time->Write();
	391	fT0mult->Write();
	392	fT0vtxRecGen->Write();
	393	fT0vtxRes->Write();
c543bd05	394	fT0ampl->Write();
c543bd05	395	fT0time2->Write();
684b93f0	396
	397	fV0a->Write();
	398	fV0c->Write();
	399	fV0multA->Write();
	400	fV0multC->Write();
	401	fV0multAcorr->Write();
	402	fV0multCcorr->Write();
	403	fV0Acorr->Write();
	404	fV0Ccorr->Write();
c543bd05	405	fV0ampl->Write();
684b93f0	406
	407	outputFile->Close();
	408	delete outputFile;
	409
	410	//delete esd;
	411	Info("AliAnaFwdDetsQA", "Successfully finished");
	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
	467
	468
	469
470
471
472
473
474