11813a2e3360905a4732f488671acdb6f74e9f0f
[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   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 }