[u/mrichter/AliRoot.git] / PWG4 / AliAnaESDSpectraQA.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.                  *
 **************************************************************************/

/* $Id: $ */
#include "AliAnaESDSpectraQA.h"

#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TList.h"
#include "TChain.h"
#include "TDirectory.h"

#include "AliAnalysisManager.h"
#include "AliLog.h"
#include "AliESDInputHandler.h"
#include "AliESDEvent.h"
#include "AliESDVertex.h"
#include "AliESDtrack.h"
#include "AliESDtrackCuts.h"


const Int_t AliAnaESDSpectraQA::fgkNPtBins=38;
const Float_t AliAnaESDSpectraQA::fgkPtMin=2;
const Float_t AliAnaESDSpectraQA::fgkPtMax=40;
const Int_t AliAnaESDSpectraQA::fgkNPhiBins=18;

ClassImp(AliAnaESDSpectraQA)

AliAnaESDSpectraQA::AliAnaESDSpectraQA(): AliAnalysisTask("AliAnaESDSpectraQA", ""), 
  fESD(0), 
  fTrackCuts(0), 
  fNEvent(0), // just to avoid warnings, inititialized in InitPointers too
  fPtAll(0),  //
  fPtSel(0),  //
  fHistList(0)
{
  InitHistPointers();
}

AliAnaESDSpectraQA::AliAnaESDSpectraQA(const char *name): 
  AliAnalysisTask(name, ""), 
  fESD(0),
  fNEvent(0), // just to avoid warnings, inititialized in InitPointers too
  fPtAll(0),  //
  fPtSel(0),  // 
  fHistList(0) {
  // Input slot #0 works with a TChain ESD
  DefineInput(0, TChain::Class());
  // Output slot #0 writes into a TList
  DefineOutput(0, TList::Class());
  InitHistPointers();
  fTrackCuts = new AliESDtrackCuts;
  fTrackCuts->SetAcceptKingDaughters(kFALSE);
  fTrackCuts->SetRequireTPCRefit(kTRUE);
  fTrackCuts->SetEtaRange(-1,1);
  // Add chisq criterium to reject 'bad' tracks that might not make it as prim
}

//________________________________________________________________________
void AliAnaESDSpectraQA::ConnectInputData(Option_t *) 
{
  // Connect ESD here
  // Called once
  TObject* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
  
  if( handler && handler->InheritsFrom("AliESDInputHandler") ) {
     fESD  =  ((AliESDInputHandler*)handler)->GetEvent();
  } 
  else {
    AliFatal("I can't get any ESD Event Handler");
  }
}

void AliAnaESDSpectraQA::InitHistPointers() {
  fNEvent = 0;
  fPtAll = 0;
  fPtSel = 0;
  for (Int_t i=0; i< 4; i++) {
    fHists[i].PhiPtNPointTPC = 0;
    fHists[i].PhiPtNPointITS = 0;
    fHists[i].PhiPtChisqC = 0;
    fHists[i].PhiPtChisqTPC = 0;
    fHists[i].PhiPtDCAR = 0;
    fHists[i].PhiPtDCAZ = 0;
  }
}

void AliAnaESDSpectraQA::CreateOutputObjects() {
  fHistList = new TList();
  //fDirectory = new TDirectory("trig_hists","Trigger histos");
  //fDirectory->cd();
  TString labels[4];
  labels[kNegA]="NegA";
  labels[kPosA]="PosA";
  labels[kNegC]="NegC";
  labels[kPosC]="PosC";

  static const Float_t kMinPhi = 0;
  static const Float_t kMaxPhi = 2*TMath::Pi();
  fNEvent = new TH1F("NEvent","NEvent",1,-0.5,0.5);
  fHistList->Add(fNEvent);
  fPtAll = new TH1F("PtAll","PtAll",fgkNPtBins, fgkPtMin, fgkPtMax);
  fHistList->Add(fPtAll);
  fPtSel = new TH1F("PtSel","PtSel",fgkNPtBins, fgkPtMin, fgkPtMax);
  fHistList->Add(fPtSel);
  for (Int_t iSide = 0; iSide < 4; iSide++) {
    TString hname="PhiPtNpointTPC";
    hname += labels[iSide];
    fHists[iSide].PhiPtNPointTPC = new TH3F(hname,hname+";#phi;p_{T} (GeV);N_{point,TPC}",fgkNPhiBins,kMinPhi,kMaxPhi,fgkNPtBins,fgkPtMin,fgkPtMax,160,0.5,160.5);
    fHistList->Add(fHists[iSide].PhiPtNPointTPC);

    hname="PhiPtNpointITS";
    hname += labels[iSide];
    fHists[iSide].PhiPtNPointITS = new TH3F(hname,hname+";#phi;p_{T} (GeV);N_{point,ITS}",fgkNPhiBins,kMinPhi,kMaxPhi,fgkNPtBins,fgkPtMin,fgkPtMax,9,-0.5,8.5);
    fHistList->Add(fHists[iSide].PhiPtNPointITS);

    hname="PhiPtChisqC";
    hname += labels[iSide];
    fHists[iSide].PhiPtChisqC = new TH3F(hname,hname+";#phi;p_{T} (GeV);#Chi^{2}/NDF",fgkNPhiBins,kMinPhi,kMaxPhi,fgkNPtBins,fgkPtMin,fgkPtMax,160,0,80);
    fHistList->Add(fHists[iSide].PhiPtChisqC);

    hname="PhiPtChisqTPC";
    hname += labels[iSide];
    fHists[iSide].PhiPtChisqTPC = new TH3F(hname,hname+";#phi;p_{T} (GeV);#Chi^{2}/NDF",fgkNPhiBins,kMinPhi,kMaxPhi,fgkNPtBins,fgkPtMin,fgkPtMax,50,0,5);
    fHistList->Add(fHists[iSide].PhiPtChisqTPC);

    hname="PhiPtDCAR";
    hname += labels[iSide];
    fHists[iSide].PhiPtDCAR = new TH3F(hname,hname+";#phi;p_{T} (GeV);DCAR",fgkNPhiBins,kMinPhi,kMaxPhi,fgkNPtBins,fgkPtMin,fgkPtMax,200,-1,1);
    fHistList->Add(fHists[iSide].PhiPtDCAR);

    hname="PhiPtDCAZ";
    hname += labels[iSide];
    fHists[iSide].PhiPtDCAZ = new TH3F(hname,hname+";#phi;p_{T} (GeV);DCAZ",fgkNPhiBins,kMinPhi,kMaxPhi,fgkNPtBins,fgkPtMin,fgkPtMax,200,-2,2);
    fHistList->Add(fHists[iSide].PhiPtDCAZ);

    hname="PhiPtSigmaToVertex";
    hname += labels[iSide];
    fHists[iSide].PhiPtSigmaToVertex = new TH3F(hname,hname+";#phi;p_{T} (GeV);n#sigma to vtx",fgkNPhiBins,kMinPhi,kMaxPhi,fgkNPtBins,fgkPtMin,fgkPtMax,50,0,5);
    fHistList->Add(fHists[iSide].PhiPtSigmaToVertex);
  }
}

void AliAnaESDSpectraQA::Exec(Option_t *) {  
  // Main loop
  // Called for each event

  if (!fESD) {
    Printf("ERROR: fESD not available");
    return;
  }

  const AliESDVertex *vtx = fESD->GetPrimaryVertex();

  // Need vertex cut
  if (vtx->GetNContributors() < 2)
    return;

  printf("Vertex title %s, status %d, nCont %d\n",vtx->GetTitle(), vtx->GetStatus(), vtx->GetNContributors());
  // Need to keep track of evts without vertex

  Int_t nTracks = fESD->GetNumberOfTracks();
  AliDebug(2,Form("nTracks %d", nTracks));
  printf("nTracks %d\n", nTracks);
  static Int_t Mult = 0;
  Mult = 0;   // Need extra init bc of static
  for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
    AliESDtrack *track = fESD->GetTrack(iTrack);
    hists *curTypeHists = 0;

    if (fTrackCuts->AcceptTrack(track)) {
      Mult++;

      Float_t dca2D, dcaZ;
      track->GetImpactParameters(dca2D,dcaZ);
      fPtAll->Fill(track->Pt());

      if (track->Eta() > 0) { // A side (crude for now)
	if (track->Charge() > 0) 
	  curTypeHists = &fHists[kPosA];
	else
	  curTypeHists = &fHists[kNegA];	
      }
      else { // C side
	if (track->Charge() > 0) 
	  curTypeHists = &fHists[kPosC];
	else
	  curTypeHists = &fHists[kNegC];	
      }

      Float_t phi = track->Phi();
      Float_t pt = track->Pt();

      UChar_t itsMap = track->GetITSClusterMap();
      Int_t nPointITS = 0;
      for (Int_t i=0; i < 6; i++) {
	if (itsMap & (1 << i))
	  nPointITS ++;
      }

      Float_t sigToVertex = fTrackCuts->GetSigmaToVertex(track);
      Float_t chisqC = 1000;
      if (track->GetConstrainedParam())
	chisqC = track->GetConstrainedChi2();

      curTypeHists->PhiPtNPointTPC->Fill(phi,pt,track->GetTPCNcls());
      curTypeHists->PhiPtNPointITS->Fill(phi,pt,nPointITS);
      curTypeHists->PhiPtChisqC->Fill(phi,pt,chisqC);
      curTypeHists->PhiPtChisqTPC->Fill(phi,pt,track->GetTPCchi2()/(track->GetTPCNclsF()-5));
      curTypeHists->PhiPtDCAR->Fill(phi,pt,dca2D);
      curTypeHists->PhiPtDCAZ->Fill(phi,pt,dcaZ);
      curTypeHists->PhiPtSigmaToVertex->Fill(phi,pt,sigToVertex);
    }
  }
  
  // Post output data
  PostData(0, fHistList);
}
Commit	Line	Data
8a907e35	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	/* $Id: $ */
	17	#include "AliAnaESDSpectraQA.h"
	18
	19	#include "TH1.h"
	20	#include "TH2.h"
	21	#include "TH3.h"
	22	#include "TList.h"
	23	#include "TChain.h"
	24	#include "TDirectory.h"
	25
	26	#include "AliAnalysisManager.h"
	27	#include "AliLog.h"
	28	#include "AliESDInputHandler.h"
	29	#include "AliESDEvent.h"
	30	#include "AliESDVertex.h"
	31	#include "AliESDtrack.h"
	32	#include "AliESDtrackCuts.h"
	33
	34
	35	const Int_t AliAnaESDSpectraQA::fgkNPtBins=38;
	36	const Float_t AliAnaESDSpectraQA::fgkPtMin=2;
	37	const Float_t AliAnaESDSpectraQA::fgkPtMax=40;
	38	const Int_t AliAnaESDSpectraQA::fgkNPhiBins=18;
	39
	40	ClassImp(AliAnaESDSpectraQA)
	41
	42	AliAnaESDSpectraQA::AliAnaESDSpectraQA(): AliAnalysisTask("AliAnaESDSpectraQA", ""),
	43	fESD(0),
	44	fTrackCuts(0),
	45	fNEvent(0), // just to avoid warnings, inititialized in InitPointers too
	46	fPtAll(0), //
	47	fPtSel(0), //
	48	fHistList(0)
	49	{
	50	InitHistPointers();
	51	}
	52
	53	AliAnaESDSpectraQA::AliAnaESDSpectraQA(const char *name):
	54	AliAnalysisTask(name, ""),
	55	fESD(0),
	56	fNEvent(0), // just to avoid warnings, inititialized in InitPointers too
	57	fPtAll(0), //
	58	fPtSel(0), //
	59	fHistList(0) {
	60	// Input slot #0 works with a TChain ESD
	61	DefineInput(0, TChain::Class());
	62	// Output slot #0 writes into a TList
	63	DefineOutput(0, TList::Class());
	64	InitHistPointers();
65	fTrackCuts = new AliESDtrackCuts;
66	fTrackCuts->SetAcceptKingDaughters(kFALSE);
67	fTrackCuts->SetRequireTPCRefit(kTRUE);
68	fTrackCuts->SetEtaRange(-1,1);
69	// Add chisq criterium to reject 'bad' tracks that might not make it as prim
70	}
71
72	//________________________________________________________________________
73	void AliAnaESDSpectraQA::ConnectInputData(Option_t *)
74	{
75	// Connect ESD here
76	// Called once
77	TObject* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
78
79	if( handler && handler->InheritsFrom("AliESDInputHandler") ) {
80	fESD = ((AliESDInputHandler*)handler)->GetEvent();
81	}
82	else {
83	AliFatal("I can't get any ESD Event Handler");
84	}
85	}
86
87	void AliAnaESDSpectraQA::InitHistPointers() {
88	fNEvent = 0;
89	fPtAll = 0;
90	fPtSel = 0;
91	for (Int_t i=0; i< 4; i++) {
92	fHists[i].PhiPtNPointTPC = 0;
93	fHists[i].PhiPtNPointITS = 0;
94	fHists[i].PhiPtChisqC = 0;
95	fHists[i].PhiPtChisqTPC = 0;
96	fHists[i].PhiPtDCAR = 0;
97	fHists[i].PhiPtDCAZ = 0;
98	}
99	}
100
101	void AliAnaESDSpectraQA::CreateOutputObjects() {
102	fHistList = new TList();
103	//fDirectory = new TDirectory("trig_hists","Trigger histos");
104	//fDirectory->cd();
105	TString labels[4];
106	labels[kNegA]="NegA";
107	labels[kPosA]="PosA";
108	labels[kNegC]="NegC";
109	labels[kPosC]="PosC";
110
111	static const Float_t kMinPhi = 0;
112	static const Float_t kMaxPhi = 2*TMath::Pi();
113	fNEvent = new TH1F("NEvent","NEvent",1,-0.5,0.5);
114	fHistList->Add(fNEvent);
115	fPtAll = new TH1F("PtAll","PtAll",fgkNPtBins, fgkPtMin, fgkPtMax);
116	fHistList->Add(fPtAll);
117	fPtSel = new TH1F("PtSel","PtSel",fgkNPtBins, fgkPtMin, fgkPtMax);
118	fHistList->Add(fPtSel);
119	for (Int_t iSide = 0; iSide < 4; iSide++) {
120	TString hname="PhiPtNpointTPC";
121	hname += labels[iSide];
122	fHists[iSide].PhiPtNPointTPC = new TH3F(hname,hname+";#phi;p_{T} (GeV);N_{point,TPC}",fgkNPhiBins,kMinPhi,kMaxPhi,fgkNPtBins,fgkPtMin,fgkPtMax,160,0.5,160.5);
123	fHistList->Add(fHists[iSide].PhiPtNPointTPC);
124
125	hname="PhiPtNpointITS";
126	hname += labels[iSide];
127	fHists[iSide].PhiPtNPointITS = new TH3F(hname,hname+";#phi;p_{T} (GeV);N_{point,ITS}",fgkNPhiBins,kMinPhi,kMaxPhi,fgkNPtBins,fgkPtMin,fgkPtMax,9,-0.5,8.5);
128	fHistList->Add(fHists[iSide].PhiPtNPointITS);
129
130	hname="PhiPtChisqC";
131	hname += labels[iSide];
132	fHists[iSide].PhiPtChisqC = new TH3F(hname,hname+";#phi;p_{T} (GeV);#Chi^{2}/NDF",fgkNPhiBins,kMinPhi,kMaxPhi,fgkNPtBins,fgkPtMin,fgkPtMax,160,0,80);
133	fHistList->Add(fHists[iSide].PhiPtChisqC);
134
135	hname="PhiPtChisqTPC";
136	hname += labels[iSide];
137	fHists[iSide].PhiPtChisqTPC = new TH3F(hname,hname+";#phi;p_{T} (GeV);#Chi^{2}/NDF",fgkNPhiBins,kMinPhi,kMaxPhi,fgkNPtBins,fgkPtMin,fgkPtMax,50,0,5);
138	fHistList->Add(fHists[iSide].PhiPtChisqTPC);
139
140	hname="PhiPtDCAR";
141	hname += labels[iSide];
142	fHists[iSide].PhiPtDCAR = new TH3F(hname,hname+";#phi;p_{T} (GeV);DCAR",fgkNPhiBins,kMinPhi,kMaxPhi,fgkNPtBins,fgkPtMin,fgkPtMax,200,-1,1);
143	fHistList->Add(fHists[iSide].PhiPtDCAR);
144
145	hname="PhiPtDCAZ";
146	hname += labels[iSide];
147	fHists[iSide].PhiPtDCAZ = new TH3F(hname,hname+";#phi;p_{T} (GeV);DCAZ",fgkNPhiBins,kMinPhi,kMaxPhi,fgkNPtBins,fgkPtMin,fgkPtMax,200,-2,2);
148	fHistList->Add(fHists[iSide].PhiPtDCAZ);
149
150	hname="PhiPtSigmaToVertex";
151	hname += labels[iSide];
152	fHists[iSide].PhiPtSigmaToVertex = new TH3F(hname,hname+";#phi;p_{T} (GeV);n#sigma to vtx",fgkNPhiBins,kMinPhi,kMaxPhi,fgkNPtBins,fgkPtMin,fgkPtMax,50,0,5);
153	fHistList->Add(fHists[iSide].PhiPtSigmaToVertex);
154	}
155	}
156
157	void AliAnaESDSpectraQA::Exec(Option_t *) {
158	// Main loop
159	// Called for each event
160
161	if (!fESD) {
162	Printf("ERROR: fESD not available");
163	return;
164	}
165
166	const AliESDVertex *vtx = fESD->GetPrimaryVertex();
167
168	// Need vertex cut
169	if (vtx->GetNContributors() < 2)
170	return;
171
172	printf("Vertex title %s, status %d, nCont %d\n",vtx->GetTitle(), vtx->GetStatus(), vtx->GetNContributors());
173	// Need to keep track of evts without vertex
174
175	Int_t nTracks = fESD->GetNumberOfTracks();
176	AliDebug(2,Form("nTracks %d", nTracks));
177	printf("nTracks %d\n", nTracks);
178	static Int_t Mult = 0;
179	Mult = 0; // Need extra init bc of static
180	for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
181	AliESDtrack *track = fESD->GetTrack(iTrack);
182	hists *curTypeHists = 0;
183
184	if (fTrackCuts->AcceptTrack(track)) {
185	Mult++;
186
187	Float_t dca2D, dcaZ;
188	track->GetImpactParameters(dca2D,dcaZ);
189	fPtAll->Fill(track->Pt());
190
191	if (track->Eta() > 0) { // A side (crude for now)
192	if (track->Charge() > 0)
193	curTypeHists = &fHists[kPosA];
194	else
195	curTypeHists = &fHists[kNegA];
196	}
197	else { // C side
198	if (track->Charge() > 0)
199	curTypeHists = &fHists[kPosC];
200	else
201	curTypeHists = &fHists[kNegC];
202	}
203
204	Float_t phi = track->Phi();
205	Float_t pt = track->Pt();
206
207	UChar_t itsMap = track->GetITSClusterMap();
208	Int_t nPointITS = 0;
209	for (Int_t i=0; i < 6; i++) {
210	if (itsMap & (1 << i))
211	nPointITS ++;
212	}
213
214	Float_t sigToVertex = fTrackCuts->GetSigmaToVertex(track);
215	Float_t chisqC = 1000;
216	if (track->GetConstrainedParam())
217	chisqC = track->GetConstrainedChi2();
218
219	curTypeHists->PhiPtNPointTPC->Fill(phi,pt,track->GetTPCNcls());
220	curTypeHists->PhiPtNPointITS->Fill(phi,pt,nPointITS);
221	curTypeHists->PhiPtChisqC->Fill(phi,pt,chisqC);
222	curTypeHists->PhiPtChisqTPC->Fill(phi,pt,track->GetTPCchi2()/(track->GetTPCNclsF()-5));
223	curTypeHists->PhiPtDCAR->Fill(phi,pt,dca2D);
224	curTypeHists->PhiPtDCAZ->Fill(phi,pt,dcaZ);
225	curTypeHists->PhiPtSigmaToVertex->Fill(phi,pt,sigToVertex);
226	}
227	}
228
229	// Post output data
230	PostData(0, fHistList);
231	}