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