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