Updated QA (Sylwester)
[u/mrichter/AliRoot.git] / TRD / qaAnalysis / AliTRDqaJPsi.cxx
CommitLineData
8e2f611a 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 iLs" without express or implied warranty. *
14 **************************************************************************/
15
16/* $Id$ */
17
18//
19// This class is a part of a package of high level QA monitoring for TRD.
20//
21// S. Radomski
22// radomski@physi.uni-heidelberg.de
23// March 2008
24//
25
26#include "AliTRDqaJPsi.h"
27#include "AliTRDqaAT.h"
28
29#include "TFile.h"
30#include "TH1D.h"
31#include "TH2D.h"
32#include "TChain.h"
33
34#include "AliESDEvent.h"
35#include "AliESDtrack.h"
36#include "AliKFParticle.h"
37
38#include "TLorentzVector.h"
39
40//const AliTRDqaJPsi::knSteps = 5;
41
42//______________________________________________________________________________
43
44AliTRDqaJPsi::AliTRDqaJPsi()
45 : AliAnalysisTask("",""),
46 fChain(0),
47 fESD(0),
48 fOutputContainer(0),
49 fnKFtracks(0)
50{
51 //
52 // default dummy constructor
53 //
54
55}
56//______________________________________________________________________________
57
58AliTRDqaJPsi:: AliTRDqaJPsi(AliTRDqaJPsi& /*trd*/)
59 : AliAnalysisTask("",""),
60 fChain(0),
61 fESD(0),
62 fOutputContainer(0),
63 fnKFtracks(0)
64{
65 //
66 // Dummy copy constructor
67 //
68
69 //return *this;
70}
71
72
73//______________________________________________________________________________
74AliTRDqaJPsi::AliTRDqaJPsi(const char *name)
75 : AliAnalysisTask(name,""),
76 fChain(0),
77 fESD(0),
78 fOutputContainer(0),
79 fnKFtracks(0)
80
81{
82 // Constructor.
83 // Input slot #0 works with an Ntuple
84 DefineInput(0, TChain::Class());
85 // Output slot #0 writes into a TH1 container
86 DefineOutput(0, TObjArray::Class()) ;
87}
88
89//______________________________________________________________________________
90void AliTRDqaJPsi::ConnectInputData(const Option_t *)
91{
92 // Initialisation of branch container and histograms
93
94 //AliInfo(Form("*** Initialization of %s", GetName())) ;
95
96 fChain = (TChain*)GetInputData(0);
97 fESD = new AliESDEvent();
98 fESD->ReadFromTree(fChain);
99}
100
101//________________________________________________________________________
102void AliTRDqaJPsi::CreateOutputObjects()
103{
104
105 Int_t c = 0;
106 fOutputContainer = new TObjArray(100);
107
108 const char *charge[2] = {"Neg", "Pos"};
109
110 // build histograms
111
112 for(Int_t i=0; i<knSteps; i++) {
113
114 fStatus[i] = new TH1D(Form("status_%d", i), "status", 32, -0.5, 31.5);
115 fOutputContainer->AddAt(fStatus[i], c++);
116
117 fInvMass[i] = new TH1D(Form("mass_%d", i), ";m_{inv} (GeV);", 100, 0, 5);
118 fOutputContainer->AddAt(fInvMass[i], c++);
119
120 fInvMassVec[i] = new TH1D(Form("massVec_%d", i), ";m_{inv} (GeV);", 100, 0, 5);
121 fOutputContainer->AddAt(fInvMassVec[i], c++);
122
123 fInvMassDiff[i] = new TH1D(Form("massDiff_%d", i), ";m_{inv} (GeV);", 100, -1, 1);
124 fOutputContainer->AddAt(fInvMassDiff[i], c++);
125
126 fAngleSM[i] = new TH1D(Form("angleSM_%d", i), ";#delta SM", 19, -0.5, 18.5);
127 fOutputContainer->AddAt(fAngleSM[i], c++);
128
129 fPtAngle[i] = new TH2D(Form("ptAngle_%d", i), ";p_{T} (GeV/c);#delta SM",
130 20, 0, 5, 10, -0.5, 9.5);
131 fOutputContainer->AddAt(fPtAngle[i], c++);
132
133
134 for(Int_t j=0; j<2; j++) {
135 fnTracks[j*knSteps+i] =
136 new TH1D(Form("nTracks%s_%d", charge[j],i), Form("%s;number of tracks",charge[j]), 100, -0.5, 99.5);
137 fPt[j*knSteps+i] = new TH1D(Form("pt%s_%d", charge[j], i), Form("%s;p_{T} (GeV/c)", charge[j]), 100, 0, 5);
138 fPID[j*knSteps+i] = new TH1D(Form("pid%s_%d", charge[j], i), ";electron LQ", 100, 0, 1);
139
140 fOutputContainer->AddAt(fnTracks[j*knSteps+i], c++);
141 fOutputContainer->AddAt(fPt[j*knSteps+i], c++);
142 fOutputContainer->AddAt(fPID[j*knSteps+i], c++);
143 }
144 }
145
146 //TH2D *fnGoodTracks;
147
148 printf("n hist = %d\n", c);
149}
150//______________________________________________________________________________
151void AliTRDqaJPsi::Exec(Option_t *)
152{
153 /*
154 Selection steps:
155 - Parameters In and Out
156 - TRDrefit bit
157 - TRDpid and quality
158 - pt > 0.8 GeV/c
159 - PID > 0.9
160 */
161
162
163 // Process one event
164 Long64_t entry = fChain->GetReadEntry() ;
165 if (!(entry%100)) Info("Exec", "Entry = %ld", entry);
166
167 // Processing of one event
168
169 if (!fESD) {
170 //AliError("fESD is not connected to the input!") ;
171 return ;
172 }
173
174
175 Int_t nTracks = fESD->GetNumberOfTracks();
176 Int_t cTracks[2*knSteps] = {0,0,0,0,0,0,0,0,0,0};
177 fnKFtracks = 0;
178
179 // track loop
180 for(Int_t i=0; i<nTracks; i++) {
181
182 //
183 // track selection
184 //
185 // param in and Out
186 // TRDrefit and TRDPid bit
187 //
188
189 AliESDtrack *track = fESD->GetTrack(i);
190 const AliExternalTrackParam *paramOut = track->GetOuterParam();
191 const AliExternalTrackParam *paramIn = track->GetInnerParam();
192
193 // long track ..
194 if (!paramIn) continue;
195 if (!paramOut) continue;
196
197 Int_t step = 0;
198 Int_t charge = (track->Charge() > 0) ? 1 : 0;
199 UInt_t status = track->GetStatus();
200 Double_t pt = track->Pt();
201 Double_t pid = track->GetTRDpid(AliPID::kElectron);
202
203 Double_t esdPid[5];
204 track->GetESDpid(esdPid);
205
206 // create a kalman particle
207 Int_t pdg = (charge == 0)? -11 : 11;
208 for(Int_t k=0; k<knSteps; k++) fInSample[fnKFtracks][k] = 0;
209
210 vec[fnKFtracks] = CreateVector(track);
211 tracks[fnKFtracks] = new AliKFParticle(*track, pdg);
212 fSM[fnKFtracks] = AliTRDqaAT::GetSector(paramOut->GetAlpha());
213 fnKFtracks++;
214
215 //AliTRDqaAT::PrintPID(track);
216
217 // apply the cuts
218
219 cTracks[knSteps *charge + step]++;
220 FillHist(track, step++);
221
222 if (!(status & AliESDtrack::kTRDrefit)) continue;
223
224 cTracks[knSteps *charge + step]++;
225 FillHist(track, step++);
226
227 if (!(status & AliESDtrack::kTRDpid)) continue;
228 if (track->GetTRDpidQuality() < 6) continue;
229
230 cTracks[knSteps *charge + step]++;
231 FillHist(track, step++);
232
233 if (pt < 0.8) continue;
234
235 cTracks[knSteps *charge + step]++;
236 FillHist(track, step++);
237
238 if (pid < 0.3) continue; //
239 //if (esdPid[AliPID::kElectron] < 0.5) continue;
240
241 cTracks[knSteps *charge + step]++;
242 FillHist(track, step);
243
244 for(Int_t k=0; k<2*knSteps; k++) fnTracks[k]->Fill(cTracks[k]);
245 }
246
247 // calculate invariant mass
248
249 for(Int_t k=0; k<knSteps; k++) {
250 for(Int_t i=0; i<fnKFtracks; i++) {
251 if (!fInSample[i][k]) continue;
252 for(Int_t j=i+1; j<fnKFtracks; j++) {
253 if (!fInSample[j][k]) continue;
254 AliKFParticle jpsi(*(tracks[i]), *(tracks[j]));
255 TLorentzVector jpsiVec = (*(vec[i])) + (*vec[j]);
256 fInvMass[k]->Fill(jpsi.GetMass());
257 fInvMassVec[k]->Fill(jpsiVec.M());
258 fInvMassDiff[k]->Fill(jpsiVec.M() - jpsi.GetMass());
259
260 if (jpsi.GetMass() > 2.5 && jpsi.GetMass() < 3.5) {
261 Int_t dSM = TMath::Abs(fSM[i] - fSM[j]);
262 if (dSM > 9) dSM = (18-dSM);
263 fAngleSM[k]->Fill(dSM);
264 fPtAngle[k]->Fill(TMath::Hypot(jpsi.GetPx(), jpsi.GetPy()), dSM);
265 }
266
267 }
268 }
269 }
270
271 PostData(0, fOutputContainer);
272}
273
274//______________________________________________________________________________
275void AliTRDqaJPsi::Terminate(Option_t *)
276{
277 // save histograms
278 fOutputContainer = (TObjArray*)GetOutputData(0);
279
280 TFile *file = new TFile("outJPsi.root", "RECREATE");
281 fOutputContainer->Write();
282
283 file->Flush();
284 file->Close();
285 delete file;
286
287 //for(Int_t i=0; i<fOutputContainer->GetEntries(); i++) {
288 // TObject *obj = fOu
289 // }
290}
291
292//______________________________________________________________________________
293
294void AliTRDqaJPsi::FillHist(AliESDtrack *track, Int_t step) {
295
296 Int_t charge = (track->Charge() > 0) ? 1 : 0;
297 UInt_t status = track->GetStatus();
298 Double_t pt = track->Pt();
299 Double_t pid = track->GetTRDpid(AliPID::kElectron);
300
301 Double_t esdPid[5];
302 track->GetESDpid(esdPid);
303
304 Int_t id = charge * knSteps + step;
305 AliTRDqaAT::FillStatus(fStatus[step], status);
306 fPt[id]->Fill(pt);
307 fPID[id]->Fill(pid);
308 //fPID[id]->Fill(esdPid[AliPID::kElectron]);
309
310 fInSample[fnKFtracks-1][step] = 1;
311}
312
313//______________________________________________________________________________
314
315TLorentzVector *AliTRDqaJPsi::CreateVector(AliESDtrack *track) {
316
317 TLorentzVector *vec = new TLorentzVector();
318 vec->SetXYZM(track->Px(), track->Py(), track->Pz(), 0);
319 return vec;
320}
321
322//______________________________________________________________________________