]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPCcalibPID.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
some histos added for TPC clusters
[u/mrichter/AliRoot.git] / TPC / AliTPCcalibPID.cxx
CommitLineData
401ad2f1 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/*
17
18
19Basic calibration and QA class for the PID of the TPC based on dEdx.
20
21The corrections which are done on the cluster level (angle,drift time etc.) are here checked by plotting the dE/dx for selected tracks as a function of these variables. In addition the Bethe-Bloch-Parameterisation is fitted to the distribution and the corresponding parameters are being extracted.
22
231.a) Fast equalization for different gain in pad regions:
24
25TFile f("CalibObjects.root")
26AliTPCcalibPID *cal = f->Get("TPCCalib")->FindObject("calibPID")
27
37fb53c1 28cal->GetHistRatioMaxTot()->Projection(0)->Fit("gaus")
29cal->GetHistRatioTracklet()->Projection(0)->Fit("gaus")
401ad2f1 30
311.b) Update OCDB:
32.x $ALICE_ROOT/TPC/macros/ConfigOCDB.C
33AliTPCClusterParam * parcl = AliTPCcalibDB::Instance()->GetClusterParam();
37fb53c1 34*parcl->fQpadMnorm)[ipad] = oldvalue*corrFactor
401ad2f1 35
36Int_t runNumber = 0;
37AliCDBMetaData *metaData= new AliCDBMetaData();
38metaData->SetObjectClassName("AliTPCClusterParam");
39metaData->SetResponsible("Alexander Kalweit");
40metaData->SetBeamPeriod(1);
41metaData->SetAliRootVersion("05-23-02");
42metaData->SetComment("October runs calibration");
43AliCDBId id1("TPC/Calib/ClusterParam", runNumber, AliCDBRunRange::Infinity());
44gStorage = AliCDBManager::Instance()->GetStorage("local://$ALICE_ROOT/OCDB");
45gStorage->Put(parcl, id1, metaData);
46
47
482.) Checks should be done with particles on the Fermi-Plateau and whoch leave a long track:
49
50cal->GetHistQmax()->GetAxis(6)->SetRangeUser(120,160)
51cal->GetHistQmax()->GetAxis(4)->SetRangeUser(20,50)
52
53variables to check:
54
55cal->GetHistQmax()->Projection(0,1)->Draw() z-dep.
56cal->GetHistQmax()->Projection(0,2)->Draw() snp-dep.
57cal->GetHistQmax()->Projection(0,3)->Draw() tgl-dep.
58
59
60
61 Comments to be written here:
62 1. What do we calibrate.
63 2. How to interpret results
64 3. Simple example
65 4. Analysis using debug streamers.
66
67Double_t initialParam[] = {3.81470,15.2608,3.25306e-07,2.51791,2.71012}
68
69
70Send comments etc. to: A.Kalweit@gsi.de, marian.ivanov@cern.ch
71*/
72
73
74#include "Riostream.h"
75#include "TChain.h"
76#include "TTree.h"
77#include "TH1F.h"
78#include "TH2F.h"
79#include "TList.h"
80#include "TMath.h"
81#include "TCanvas.h"
82#include "TFile.h"
83#include "TF1.h"
84#include "TVectorD.h"
85#include "TProfile.h"
86
87#include "AliTPCcalibDB.h"
88#include "AliTPCclusterMI.h"
89#include "AliTPCClusterParam.h"
90#include "AliTPCseed.h"
91#include "AliESDVertex.h"
92#include "AliESDEvent.h"
93#include "AliESDfriend.h"
94#include "AliESDInputHandler.h"
95#include "AliAnalysisManager.h"
96#include "AliTPCParam.h"
97
98#include "AliLog.h"
99
100#include "AliTPCcalibPID.h"
101
102#include "TTreeStream.h"
103
104
105ClassImp(AliTPCcalibPID)
106
107
108AliTPCcalibPID::AliTPCcalibPID()
109 :AliTPCcalibBase(),
110 fMIP(0),
111 fLowerTrunc(0),
112 fUpperTrunc(0),
113 fUseShapeNorm(0),
114 fUsePosNorm(0),
115 fUsePadNorm(0),
116 fIsCosmic(0),
117 fHistNTracks(0),
118 fClusters(0),
119 fPileUp(0),
120 fLandau(0),
121 fDeDxQmax(0),
122 fDeDxQtot(0),
37fb53c1 123 fDeDxRatioMaxTot(0),
124 fDeDxRatioQmax(0),
125 fDeDxRatioQtot(0),
126 fDeDxRatioTruncQtot(0),
127 fDeDxRatioTruncQmax(0)
401ad2f1 128{
129 //
130 // Empty default cosntructor
131 //
132 AliInfo("Default Constructor");
133}
134
135
136AliTPCcalibPID::AliTPCcalibPID(const Text_t *name, const Text_t *title)
137 :AliTPCcalibBase(),
138 fMIP(0),
139 fLowerTrunc(0),
140 fUpperTrunc(0),
141 fUseShapeNorm(0),
142 fUsePosNorm(0),
143 fUsePadNorm(0),
144 fIsCosmic(0),
145 fHistNTracks(0),
146 fClusters(0),
147 fPileUp(0),
148 fLandau(0),
149 fDeDxQmax(0),
150 fDeDxQtot(0),
37fb53c1 151 fDeDxRatioMaxTot(0),
152 fDeDxRatioQmax(0),
153 fDeDxRatioQtot(0) ,
154 fDeDxRatioTruncQtot(0),
155 fDeDxRatioTruncQmax(0)
401ad2f1 156{
157 //
158 //
159 //
160 SetName(name);
161 SetTitle(title);
162 //
163 fMIP = 40.;
164 fLowerTrunc = 0;
165 fUpperTrunc = 0.6;
166 //
167 fUseShapeNorm = kTRUE;
37fb53c1 168 fUsePosNorm = 0;
401ad2f1 169 fUsePadNorm = 1;
170 //
171 fIsCosmic = kTRUE;
172 //
37fb53c1 173 // dE/dx, z, phi, theta, p, bg, ncls, tracklet type
174 Int_t binsQA[8] = {150, 10, 10, 10, 50, 50, 8, 5};
175 Double_t xminQA[8] = {0., 0, 0, 0, 0.01, 0.1, 60, 0};
176 Double_t xmaxQA[8] = {10., 1, 0.6, 1.5, 50, 500, 160, 5};
401ad2f1 177 //
37fb53c1 178 //
179 //
180 // dE/dx, z, phi, theta, dEdx, dEdx*dl, ncls, tracklet type
181 Int_t binsRA[9] = {100, 10, 10, 10, 25, 25, 4, 5};
182 Double_t xminRa[9] = {0.5, 0, 0, 0, 0.2, 0.2, 60, 0};
183 Double_t xmaxRa[9] = {4.0, 1, 0.6, 1.5, 5, 5, 160, 5};
401ad2f1 184
185 // z,sin(phi),tan(theta),p,betaGamma,ncls
37fb53c1 186 fDeDxQmax = new THnSparseS("fDeDxQmax","Qmax;(z,sin(phi),tan(theta),p,betaGamma,ncls,type); TPC signal Qmax (a.u.)",8,binsQA,xminQA,xmaxQA);
187 fDeDxQtot = new THnSparseS("fDeDxQtot","Q_{tot};(z,sin(phi),tan(theta),p,betaGamma,ncls,type); TPC signal Qmax (a.u.)",8,binsQA,xminQA,xmaxQA);
188 //
189 // ratio histograms
190 //
191 fDeDxRatioMaxTot = new THnSparseS("fDeDxRatioMaxTot","Q_{max}/Q_{tot};(z,sin(phi),tan(theta),dEdx,dEdx*dl,ncls,type); TPC signal Qmax/Qtot (a.u.)",8,binsRA,xminRa,xmaxRa);
192 fDeDxRatioQmax = new THnSparseS("fDeDxRatioQmax","Q_{mtracklet}/Q_{mtrack};(z,sin(phi),tan(theta),dEdx,dEdx*dl,ncls,type,qtupe); TPC signal Tracklet/Track (a.u.)",8,binsRA,xminRa,xmaxRa);
193 fDeDxRatioQtot = new THnSparseS("fDeDxRatioQtot","Q_{ttracklet}/Q_{ttrack};(z,sin(phi),tan(theta),dEdx,dEdx*dl,ncls,type,qtupe); TPC signal Tracklet/Track (a.u.)",8,binsRA,xminRa,xmaxRa);
194 fDeDxRatioTruncQmax = new THnSparseS("fDeDxRatioTrunQmax","Q_{max}/Q_{maxtrunc};(z,sin(phi),tan(theta),dEdx,dEdx*dl,ncls,type,qtupe); TPC signal Full/Trunc. (a.u.)",8,binsRA,xminRa,xmaxRa);
195 fDeDxRatioTruncQtot = new THnSparseS("fDeDxRatioTruncQtot","Q_{tot}/Q_{tottrunc};(z,sin(phi),tan(theta),dEdx,dEdx*dl,ncls,type,qtupe); TPC signal Full/Trunc (a.u.)",8,binsRA,xminRa,xmaxRa);
196
197
401ad2f1 198 BinLogX(fDeDxQmax,4); BinLogX(fDeDxQmax,5);
199 BinLogX(fDeDxQtot,4); BinLogX(fDeDxQtot,5);
37fb53c1 200 //
201 BinLogX(fDeDxRatioMaxTot,4); BinLogX(fDeDxRatioMaxTot,5);
202 BinLogX(fDeDxRatioQmax,4); BinLogX(fDeDxRatioQmax,5);
203 BinLogX(fDeDxRatioQtot,4); BinLogX(fDeDxRatioQtot,5);
204 BinLogX(fDeDxRatioTruncQmax,4); BinLogX(fDeDxRatioTruncQmax,5);
205 BinLogX(fDeDxRatioTruncQtot,4); BinLogX(fDeDxRatioTruncQtot,5);
401ad2f1 206 //
207 fHistNTracks = new TH1F("ntracks","Number of Tracks per Event; number of tracks per event; number of tracks",1001,-0.5,1000.5);
208 fClusters = new TH1F("signal","Number of Clusters per track; number of clusters per track n_{cl}; counts",40,0,160);
209 fPileUp = new TH2F("PileUp","timing plots.; #Delta L ; dEdx signal ",400,-1,1,400,0,200);
210 fLandau = new TH2F("Landau","Landau.; pad type ; cluster charge ",3,-0.5,2.5,400,0,1000);
211 //
212 AliInfo("Non Default Constructor");
213 //
214}
215
216
217AliTPCcalibPID::~AliTPCcalibPID(){
218 //
219 //
220 //
37fb53c1 221
401ad2f1 222}
223
224
225
226void AliTPCcalibPID::Process(AliESDEvent *event) {
227 //
228 //
229 //
230 if (!event) {
231 Printf("ERROR: ESD not available");
232 return;
233 }
37fb53c1 234 const Int_t kMinClustersTracklet = 25;
401ad2f1 235 Int_t ntracks=event->GetNumberOfTracks();
236 fHistNTracks->Fill(ntracks);
237
238 AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
239 if (!ESDfriend) {
240 Printf("ERROR: ESDfriend not available");
241 return;
242 }
243 //
244 // track loop
245 //
246 for (Int_t i=0;i<ntracks;++i) {
247
248 AliESDtrack *track = event->GetTrack(i);
249 if (!track) continue;
250
251
252 const AliExternalTrackParam * trackIn = track->GetInnerParam();
253 const AliExternalTrackParam * trackOut = track->GetOuterParam();
254 if (!trackIn) continue;
255 if (!trackOut) continue;
256
257 // calculate necessary track parameters
258 //Double_t meanP = 0.5*(trackIn->GetP() + trackOut->GetP());
259 Double_t meanP = trackIn->GetP();
37fb53c1 260 //Double_t d = trackIn->GetLinearD(0,0);
401ad2f1 261 Int_t NclsDeDx = track->GetTPCNcls();
262
263 //if (meanP > 0.7 || meanP < 0.2) continue;
264 if (NclsDeDx < 60) continue;
265
266 // exclude tracks which do not look like primaries or are simply too short or on wrong sectors
267
268 //if (TMath::Abs(trackIn->GetSnp()) > 3*0.4) continue;
269 //if (TMath::Abs(trackIn->GetZ()) > 150) continue;
270 //if (seed->CookShape(1) > 1) continue;
271 //if (TMath::Abs(trackIn->GetY()) > 20) continue;
272 //if (TMath::Abs(d)>20) continue; // distance to the 0,0; select only tracks which cross chambers under proper angle
273 //if (TMath::Abs(trackIn->GetSnp()) > 0.6) continue;
274 //if (TMath::Abs(trackOut->GetSnp()) > 0.2) continue;
275 if (TMath::Abs(trackIn->GetAlpha()+0.872665)<0.01 || TMath::Abs(trackOut->GetAlpha()+0.872665)<0.01) continue; // Funny sector !
276
277 // Get seeds
278 AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
279 if (!friendTrack) continue;
280 TObject *calibObject;
281 AliTPCseed *seed = 0;
282 for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
283 if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
284 }
285
286 if (seed) {
287 if (meanP > 30 && TMath::Abs(trackIn->GetSnp()) < 0.2) fClusters->Fill(track->GetTPCNcls());
288 // calculate dEdx
289 // (Double_t low, Double_t up, Int_t type, Int_t i1, Int_t i2, Bool_t shapeNorm,Bool_t posNorm, Int_t padNorm, Int_t returnVal)
37fb53c1 290 //Double_t TPCsignalTot = (1/fMIP)*seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,0,0,159,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
401ad2f1 291 Double_t TPCsignalMax = (1/fMIP)*seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,0,159,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
292 //
37fb53c1 293 //Double_t TPCsignalShort = seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,0,63,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
294 //Double_t TPCsignalMedium = seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,63,63+64,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
295 //Double_t TPCsignalLong = seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,63+64,159,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
401ad2f1 296 //
297 //Double_t driftMismatch = seed->GetDriftTimeMismatch(0,0.95,0.5);
298 Double_t driftMismatch = 0;
37fb53c1 299 // Double_t drift = 1 - (TMath::Abs(trackIn->GetZ()) + TMath::Abs(trackOut->GetZ()))/500.;
401ad2f1 300
301 // particle identification
302 Double_t mass = 0.105658369;// muon
303
304 if (meanP > 30) {
305 fPileUp->Fill(driftMismatch,TPCsignalMax);
306 //
307 fLandau->Fill(0.1,0.6);
308 }
309 //var. of interest, z,sin(phi),tan(theta),p,betaGamma,ncls
37fb53c1 310 Double_t snpIn = TMath::Abs(trackIn->GetSnp());
311 Double_t snpOut = TMath::Abs(trackOut->GetSnp());
312 Double_t tglIn = TMath::Abs(trackIn->GetTgl());
313 Double_t tglOut = TMath::Abs(trackOut->GetTgl());
314 Double_t driftIn = 1 - (TMath::Abs(trackIn->GetZ()))/250.;
315 Double_t driftOut = 1 - (TMath::Abs(trackIn->GetZ()))/250.;
401ad2f1 316 //
37fb53c1 317 // dEdx in region
318 //
319 Float_t dEdxTot[5],dEdxTotFull[5];
320 Float_t dEdxMax[5],dEdxMaxFull[5];
321 Int_t ncl[5];
322 for (Int_t itype=0; itype<5;itype++){
323 Int_t row0=0, row1 =159;
324 if (itype==1) {row0=0; row1 = 63;};
325 if (itype==2) {row0= 64; row1=63+64;}
326 if (itype==3) {row0= 63+64; row1=159;}
327 dEdxTot[itype]= (1/fMIP)*seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,0,row0,row1,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
328 dEdxMax[itype]= (1/fMIP)*seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,row0,row1,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
329 // non trucated dEdx
330 dEdxTotFull[itype]= (1/fMIP)*seed->CookdEdxNorm(0.0,0.99,0,row0,row1,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
331 dEdxMaxFull[itype]= (1/fMIP)*seed->CookdEdxNorm(0.0,0.99,1,row0,row1,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
332
333
334 ncl[itype]=TMath::Nint(seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,row0,row1,fUseShapeNorm,fUsePosNorm,fUsePadNorm,2));
335 }
336 //
337 //
338 //
339 Float_t wmeanTot=0, wmeanMax=0, sumW=0;
340 Double_t length[3] = {0.75,1,1.5};
341
342 for (Int_t ipad=0; ipad<3; ipad++){
343 if (ncl[1+ipad]<3) continue;
344 Double_t weight = Double_t(ncl[1+ipad])*TMath::Sqrt(length[ipad]);
345 wmeanTot+=dEdxTot[1+ipad]*weight;
346 wmeanMax+=dEdxMax[1+ipad]*weight;
347 sumW+=weight;
348 }
349 if (sumW>0){
350 dEdxTot[4]= wmeanTot/sumW;
351 dEdxMax[4]= wmeanMax/sumW;
352 }
353 for (Int_t itype=0;itype<5;itype++){
354 //
355 Float_t snp=(TMath::Abs(snpIn)+TMath::Abs(snpOut))*0.5;
356 Float_t tgl=(TMath::Abs(tglIn)+TMath::Abs(tglOut))*0.5;
357 Float_t drift = (driftIn+driftOut)*0.5;
358 if (itype==1) {snp = TMath::Abs(snpIn); tgl = TMath::Abs(tglIn); drift= driftIn;};
359 if (itype==3) {snp = TMath::Abs(snpOut); tgl = TMath::Abs(tglOut);drift=driftOut;};
360 if (ncl[itype]<kMinClustersTracklet) continue;
361 Float_t deltaL = TMath::Sqrt(1+snp*snp+tgl*tgl);
362 //
363 Double_t vecQmax[8] = {dEdxMax[itype],drift,snp,tgl,meanP,meanP/mass,NclsDeDx, itype};
364 Double_t vecQtot[8] = {dEdxTot[itype],drift,snp,tgl,meanP,meanP/mass,NclsDeDx, itype};
365 //
366 //
367 //
368 Double_t ratioMaxTot = (dEdxTot[itype]>0) ? dEdxMax[itype]/dEdxTot[itype]:0;
369 Double_t ratioTrackletTot = (dEdxTot[0]>0) ? dEdxTot[itype]/dEdxTot[0]:0;
370 Double_t ratioTrackletMax = (dEdxMax[0]>0) ? dEdxMax[itype]/dEdxMax[0]:0;
371 Double_t ratioTrackletTruncTot = (dEdxTot[0]>0) ? dEdxTotFull[itype]/dEdxTot[itype]:0;
372 Double_t ratioTrackletTruncMax = (dEdxMax[0]>0) ? dEdxMaxFull[itype]/dEdxMax[itype]:0;
373
374 Double_t vecRatioMaxTot[8] = {ratioMaxTot, drift,snp,tgl,dEdxTot[0], dEdxTot[0]*deltaL,NclsDeDx,itype};
375 Double_t vecRatioTrackletTot[8] = {ratioTrackletTot, drift,snp,tgl,dEdxTot[0], dEdxTot[0]*deltaL,NclsDeDx,itype};
376 Double_t vecRatioTrackletMax[8] = {ratioTrackletMax, drift,snp,tgl,dEdxMax[0], dEdxMax[0]*deltaL,NclsDeDx,itype};
377 Double_t vecRatioTrackletTruncTot[8] = {ratioTrackletTruncTot, drift,snp,tgl,dEdxTot[0], dEdxTot[0]*deltaL,NclsDeDx,itype};
378 Double_t vecRatioTrackletTruncMax[8] = {ratioTrackletTruncMax, drift,snp,tgl,dEdxMax[0], dEdxMax[0]*deltaL,NclsDeDx,itype};
379 fDeDxQmax->Fill(vecQmax);
380 fDeDxQtot->Fill(vecQtot);
381 fDeDxRatioMaxTot->Fill(vecRatioMaxTot);
382 fDeDxRatioQmax->Fill(vecRatioTrackletTot);
383 fDeDxRatioQtot->Fill(vecRatioTrackletMax);
384 fDeDxRatioTruncQmax->Fill(vecRatioTrackletTruncTot);
385 fDeDxRatioTruncQtot->Fill(vecRatioTrackletTruncMax);
386 //
387 TTreeSRedirector * cstream = GetDebugStreamer();
388 if (cstream){
389 TVectorD vQT(9,vecQtot);
390 TVectorD vQM(9,vecQmax);
391 TVectorD vQMT(9, vecRatioMaxTot);
392 TVectorD vQRT(9,vecRatioTrackletTot);
393 TVectorD vQRM(9,vecRatioTrackletMax);
394 TVectorD vQRTT(9,vecRatioTrackletTruncTot);
395 TVectorD vQRTM(9,vecRatioTrackletTruncMax);
396 (*cstream) << "dEdx" <<
397 "run="<<fRun<< // run number
398 "event="<<fEvent<< // event number
399 "time="<<fTime<< // time stamp of event
400 "trigger="<<fTrigger<< // trigger
401 "mag="<<fMagF<< // magnetic field
402 "track.="<<seed<< // original seed
403 //
404 "vQT.="<<&vQT<< // trunc mean - total charge
405 "vQM.="<<&vQM<< // trunc mean - max charge
406 //
407 "vQMT.="<<&vQMT<< // ratio max/tot
408 "vQRT.="<<&vQRT<< // ratio tracklet/track - total charge
409 "vQRM.="<<&vQRM<< // ratio tracklet/track - max charge
410 "vQRTT.="<<&vQRTT<< // ratio trunc/full - total charge
411 "vQRTM.="<<&vQRTM<< // ratio trunc/full - max charge
412 "\n";
413 }
414 }
415 }
416 }
401ad2f1 417}
418
419
420
421void AliTPCcalibPID::Analyze() {
422
423
424 return;
425
426}
427
428
429Long64_t AliTPCcalibPID::Merge(TCollection *li) {
430
431 TIterator* iter = li->MakeIterator();
432 AliTPCcalibPID* cal = 0;
433
434 while ((cal = (AliTPCcalibPID*)iter->Next())) {
435 if (!cal->InheritsFrom(AliTPCcalibPID::Class())) {
436 Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
437 return -1;
438 }
439
440 if (cal->GetHistNTracks()) fHistNTracks->Add(cal->GetHistNTracks());
441 if (cal->GetHistClusters()) fClusters->Add(cal->GetHistClusters());
442 if (cal->GetHistPileUp()) fPileUp->Add(cal->GetHistPileUp());
443 if (cal->GetHistLandau()) fLandau->Add(cal->GetHistLandau());
444 //
445 if (cal->GetHistQmax()) fDeDxQmax->Add(cal->GetHistQmax());
446 if (cal->GetHistQtot()) fDeDxQtot->Add(cal->GetHistQtot());
37fb53c1 447 if (cal->GetHistRatioMaxTot()) fDeDxRatioMaxTot->Add(cal->GetHistRatioMaxTot());
448 if (cal->GetHistRatioQmax()) fDeDxRatioQmax->Add(cal->GetHistRatioQmax());
449 if (cal->GetHistRatioQtot()) fDeDxRatioQtot->Add(cal->GetHistRatioQtot());
450 if (cal->GetHistRatioTruncQmax()) fDeDxRatioTruncQmax->Add(cal->GetHistRatioTruncQmax());
451 if (cal->GetHistRatioTruncQtot()) fDeDxRatioTruncQtot->Add(cal->GetHistRatioTruncQtot());
401ad2f1 452 }
453
454 return 0;
455
456}
457
458
459
401ad2f1 460void AliTPCcalibPID::BinLogX(THnSparse *h, Int_t axisDim) {
461
462 // Method for the correct logarithmic binning of histograms
463
464 TAxis *axis = h->GetAxis(axisDim);
465 int bins = axis->GetNbins();
466
467 Double_t from = axis->GetXmin();
468 Double_t to = axis->GetXmax();
469 Double_t *new_bins = new Double_t[bins + 1];
470
471 new_bins[0] = from;
472 Double_t factor = pow(to/from, 1./bins);
473
474 for (int i = 1; i <= bins; i++) {
475 new_bins[i] = factor * new_bins[i-1];
476 }
477 axis->Set(bins, new_bins);
478 delete new_bins;
479
480}
481
37fb53c1 482
483
484
485
486void AliTPCcalibPID::MakeReport(const char *outputpath) {
487 //
488 // Make a standard QA plots
489 //
490 for (Int_t ipad=0;ipad<4;ipad++){
491 DrawRatioTot(ipad,outputpath);
492 DrawRatioMax(ipad,outputpath);
493 }
494 DrawRatiodEdx(0.5,3,outputpath);
495 DrawResolBGQtot(140,160,outputpath);
496 DrawResolBGQmax(140,160,outputpath);
497 return;
498}
499
500void AliTPCcalibPID::DrawRatioTot(Int_t ipad, const char* outputpath){
501 //
502 // Draw default ratio histogram for given pad type
503 // ipad - 0 - short pads
504 // 1 - medium pads
505 // 2 - long pads
506 //
507 // Int_t kmicolors[10]={1,2,3,4,6,7,8,9,10,11};
508 Int_t kmimarkers[10]={21,22,23,24,25,26,27,28,29,30};
509 AliTPCcalibPID * pid = this;
510 TCanvas *canvas= new TCanvas(Form("QtotRatio%d)",ipad),Form("QtotRatio%d)",ipad));
511 canvas->Divide(3,2);
512 pid->GetHistRatioQtot()->GetAxis(7)->SetRange(ipad+2,ipad+2);
513 canvas->cd(1);
514 TH1 * his0 = pid->GetHistRatioQtot()->Projection(0);
515 his0->SetXTitle("dEdx_{tracklet}/dEdx_{track}");
516 his0->SetYTitle("");
517 his0->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
518 his0->Draw();
519 //
520 canvas->cd(2);
521 TH1 * hprofz = (TH1*) (pid->GetHistRatioQtot()->Projection(0,1)->ProfileX());
522 hprofz->SetXTitle("drift length");
523 hprofz->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
524 hprofz->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
525 hprofz->SetMarkerStyle(kmimarkers[0]);
526 hprofz->SetMaximum(1.1);
527 hprofz->SetMinimum(0.9);
528 hprofz->Draw();
529 //
530 canvas->cd(3);
531 TH1 * hprofphi = (TH1*) (pid->GetHistRatioQtot()->Projection(0,2)->ProfileX());
532 hprofphi->SetXTitle("sin(#phi)");
533 hprofphi->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
534 hprofphi->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
535 hprofphi->SetMarkerStyle(kmimarkers[0]);
536 hprofz->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
537 hprofphi->SetMaximum(1.1);
538 hprofphi->SetMinimum(0.9);
539 hprofphi->Draw();
540 //
541 canvas->cd(4);
542 TH1 * hproftheta = (TH1*) (pid->GetHistRatioQtot()->Projection(0,3)->ProfileX());
543 hproftheta->SetXTitle("tan(#theta)");
544 hproftheta->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
545 hproftheta->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
546 hproftheta->SetMarkerStyle(kmimarkers[0]);
547 hproftheta->SetMaximum(1.1);
548 hproftheta->SetMinimum(0.9);
549 hproftheta->Draw();
550 //
551 canvas->cd(5);
552 TH1 * hprofdedx = (TH1*) (pid->GetHistRatioQtot()->Projection(0,4)->ProfileX());
553 hprofdedx->SetXTitle("dEdx_{track}");
554 hprofdedx->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
555 hprofdedx->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
556 hprofdedx->SetMarkerStyle(kmimarkers[0]);
557 hprofdedx->SetMaximum(1.1);
558 hprofdedx->SetMinimum(0.9);
559 hprofdedx->Draw();
560
561 canvas->cd(6);
562 TH1 * hprofdedxL = (TH1*) (pid->GetHistRatioQtot()->Projection(0,5)->ProfileX());
563 hprofdedxL->SetXTitle("dEdx_{track}#Delta_{x}");
564 hprofdedxL->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
565 hprofdedxL->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
566 hprofdedxL->SetMarkerStyle(kmimarkers[0]);
567 hprofdedxL->SetMaximum(1.1);
568 hprofdedxL->SetMinimum(0.9);
569 hprofdedxL->Draw();
570
571
572
573 canvas->SaveAs(Form("%s/QtotRatioType%d.eps",outputpath,ipad));
574 canvas->SaveAs(Form("%s/QtotRatioType%d.gif",outputpath,ipad));
575}
576
577void AliTPCcalibPID::DrawRatioMax(Int_t ipad, const char* outputpath){
578 //
579 // Draw default ration histogram for given pad type
580 // ipad - 0 - short pads
581 // 1 - medium pads
582 // 2 - long pads
583 //
584 // Int_t kmicolors[10]={1,2,3,4,6,7,8,9,10,11};
585 Int_t kmimarkers[10]={21,22,23,24,25,26,27,28,29,30};
586 AliTPCcalibPID * pid = this;
587 TCanvas *canvas= new TCanvas(Form("QmaxRatio%d)",ipad),Form("QmaxRatio%d)",ipad));
588 canvas->Divide(3,2);
589 pid->GetHistRatioQmax()->GetAxis(7)->SetRange(ipad+2,ipad+2);
590 canvas->cd(1);
591 TH1 * his0 = pid->GetHistRatioQmax()->Projection(0);
592 his0->SetXTitle("dEdx_{tracklet}/dEdx_{track}");
593 his0->SetYTitle("");
594 his0->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
595 his0->Draw();
596 //
597 canvas->cd(2);
598 TH1 * hprofz = (TH1*) (pid->GetHistRatioQmax()->Projection(0,1)->ProfileX());
599 hprofz->SetXTitle("drift length");
600 hprofz->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
601 hprofz->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
602 hprofz->SetMarkerStyle(kmimarkers[0]);
603 hprofz->SetMaximum(1.1);
604 hprofz->SetMinimum(0.9);
605 hprofz->Draw();
606 //
607 canvas->cd(3);
608 TH1 * hprofphi = (TH1*) (pid->GetHistRatioQmax()->Projection(0,2)->ProfileX());
609 hprofphi->SetXTitle("sin(#phi)");
610 hprofphi->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
611 hprofphi->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
612 hprofphi->SetMarkerStyle(kmimarkers[0]);
613 hprofphi->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
614 hprofphi->SetMaximum(1.1);
615 hprofphi->SetMinimum(0.9);
616 hprofphi->Draw();
617 //
618 canvas->cd(4);
619 TH1 * hproftheta = (TH1*) (pid->GetHistRatioQmax()->Projection(0,3)->ProfileX());
620 hproftheta->SetXTitle("tan(#theta)");
621 hproftheta->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
622 hproftheta->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
623 hproftheta->SetMarkerStyle(kmimarkers[0]);
624 hproftheta->SetMaximum(1.1);
625 hproftheta->SetMinimum(0.9);
626 hproftheta->Draw();
627
628 canvas->cd(5);
629 TH1 * hprofdedx = (TH1*) (pid->GetHistRatioQmax()->Projection(0,4)->ProfileX());
630 hprofdedx->SetXTitle("dEdx_{track}");
631 hprofdedx->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
632 hprofdedx->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
633 hprofdedx->SetMarkerStyle(kmimarkers[0]);
634 hprofdedx->SetMaximum(1.1);
635 hprofdedx->SetMinimum(0.9);
636 hprofdedx->Draw();
637
638 canvas->cd(6);
639 TH1 * hprofdedxL = (TH1*) (pid->GetHistRatioQmax()->Projection(0,5)->ProfileX());
640 hprofdedxL->SetXTitle("dEdx_{track}#Delta_{x}");
641 hprofdedxL->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
642 hprofdedxL->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
643 hprofdedxL->SetMarkerStyle(kmimarkers[0]);
644 hprofdedxL->SetMaximum(1.1);
645 hprofdedxL->SetMinimum(0.9);
646 hprofdedxL->Draw();
647
648
649 canvas->SaveAs(Form("%s/QmaxRatioType%d.eps",outputpath,ipad));
650 canvas->SaveAs(Form("%s/QmaxRatioType%d.gif",outputpath,ipad));
651}
652
653
654
655void AliTPCcalibPID::DrawRatiodEdx(Float_t demin, Float_t demax, const char* outputpath){
656 //
657 //
658 //
659 //
660 // Int_t kmicolors[10]={1,2,3,4,6,7,8,9,10,11};
661 Int_t kmimarkers[10]={21,22,23,24,25,26,27,28,29,30};
662 AliTPCcalibPID * pid = this;
663 TCanvas *canvas= new TCanvas("QRatiodEdx","QRatiodEdx");
664 canvas->Divide(2,4);
665 canvas->SetLogx(kTRUE);
666 TH1 * hprofP=0;
667 for (Int_t ipad=0;ipad<4;ipad++){
668 pid->GetHistRatioQmax()->GetAxis(7)->SetRange(ipad+2,ipad+2);
669 pid->GetHistRatioQtot()->GetAxis(7)->SetRange(ipad+2,ipad+2);
670 pid->GetHistRatioQmax()->GetAxis(5)->SetRangeUser(demin,demax);
671 pid->GetHistRatioQtot()->GetAxis(5)->SetRangeUser(demin,demax);
672
673 canvas->cd(ipad*2+1)->SetLogx(kFALSE);
674 hprofP = (TH1*) (pid->GetHistRatioQmax()->Projection(0,5)->ProfileX());
675 hprofP->SetXTitle("dEdx_{track}");
676 hprofP->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
677 hprofP->SetTitle(Form("Q_{max} dEdx_{tracklet}/dEdx_{track} type %d",0));
678 hprofP->SetMarkerStyle(kmimarkers[0]);
679 hprofP->SetMaximum(1.1);
680 hprofP->SetMinimum(0.9);
681 hprofP->Draw();
682 // pad Tot
683 canvas->cd(ipad*2+2)->SetLogx(kFALSE);
684 hprofP = (TH1*) (pid->GetHistRatioQtot()->Projection(0,5)->ProfileX());
685 hprofP->SetXTitle("dEdx_{track}");
686 hprofP->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
687 hprofP->SetTitle(Form("Q_{tot} dEdx_{tracklet}/dEdx_{track} type %d",0));
688 hprofP->SetMarkerStyle(kmimarkers[0]);
689 hprofP->SetMaximum(1.1);
690 hprofP->SetMinimum(0.9);
691 hprofP->Draw();
692 }
693 //
694 //
695 canvas->SaveAs(Form("%s/QratiodEdx%d.eps",outputpath));
696 canvas->SaveAs(Form("%s/QratiodEdx%d.gif",outputpath));
697}
698
699
700
701void AliTPCcalibPID::DrawResolBGQtot(Int_t minClusters, Int_t maxClusters, const char *outputpath){
702 //
703 //
704 //
705 AliTPCcalibPID * pid = this;
706 Int_t kmimarkers[10]={21,22,23,24,25,26,27,28,29,30};
707
708 TObjArray arr;
709 TH2 * his =0;
710 TH1 * hmean =0;
711 TH1 * hsigma=0;
712 //
713 // set cut
714 pid->GetHistQtot()->GetAxis(6)->SetRangeUser(minClusters,maxClusters);
715 pid->GetHistQtot()->GetAxis(5)->SetRangeUser(1,10000);
716 TCanvas *canvas= new TCanvas("dEdxResolQ_{Tot}","dEdxResolQ_{Tot}");
717 canvas->Divide(2,3);
718 //
719 //
720 //
721 for (Int_t ipad=0;ipad<5;ipad++){
722 canvas->cd(1+ipad)->SetLogx(kTRUE);
723 if (ipad<4) pid->GetHistQtot()->GetAxis(7)->SetRange(ipad+2,ipad+2);
724 if (ipad==4) pid->GetHistQtot()->GetAxis(7)->SetRange(1,1);
725 his = (TH2*)(pid->GetHistQtot()->Projection(0,5));
726 his->FitSlicesY(0,0,-1,0,"QNR",&arr);
727 hmean = (TH1*)arr.At(1);
728 hsigma = (TH1*)arr.At(2)->Clone();
729 hsigma->Divide(hmean);
730 arr.SetOwner(kTRUE);
731 arr.Delete();
732 delete his;
733
734 hsigma->SetXTitle("#beta#gamma");
735 hsigma->SetYTitle("#sigma_{dEdx}/dEdx");
736 hsigma->SetTitle(Form("#sigma_{dEdx}/dEdx_{tot} Pad %d",ipad));
737 hsigma->SetName(Form("#sigma_{dEdx}/dEdx_{tot} Pad %d",ipad));
738 if (ipad==4) {
739 hsigma->SetTitle(Form("#sigma_{dEdx}/dEdx_{tot} Full"));
740 hsigma->SetName(Form("#sigma_{dEdx}/dEdx_{tot} Full"));
741 }
742 hsigma->SetMarkerStyle(kmimarkers[0]);
743 hsigma->SetMaximum(0.15);
744 hsigma->SetMinimum(0.0);
745 hsigma->Draw();
746 }
747 canvas->SaveAs(Form("%s/dEdxResolMax.eps",outputpath));
748 canvas->SaveAs(Form("%s/dEdxResolMax.gif",outputpath));
749}
750
751void AliTPCcalibPID::DrawResolBGQmax(Int_t minClusters, Int_t maxClusters, const char *outputpath){
752 //
753 // Int_t minClusters=140; Int_t maxClusters=200; const char *outputpath="picPID06"
754 //
755 Int_t kmimarkers[10]={21,22,23,24,25,26,27,28,29,30};
756 AliTPCcalibPID * pid = this;
757 TObjArray arr;
758 TH2 * his =0;
759 TH1 * hmean =0;
760 TH1 * hsigma=0;
761 //
762 // set cut
763 pid->GetHistQmax()->GetAxis(6)->SetRangeUser(minClusters,maxClusters);
764 pid->GetHistQmax()->GetAxis(5)->SetRangeUser(1,10000);
765 TCanvas *canvas= new TCanvas("dEdxResolQ_{Max}","dEdxResolQ_{Max}");
766 canvas->Divide(2,3);
767 //
768 //
769 //
770 for (Int_t ipad=0;ipad<5;ipad++){
771 canvas->cd(1+ipad)->SetLogx(kTRUE);
772 if (ipad<4) pid->GetHistQmax()->GetAxis(7)->SetRange(ipad+2,ipad+2);
773 if (ipad==4) pid->GetHistQmax()->GetAxis(7)->SetRange(1,1);
774 his = (TH2*)(pid->GetHistQmax()->Projection(0,5));
775 his->FitSlicesY(0,0,-1,0,"QNR",&arr);
776 hmean = (TH1*)arr.At(1);
777 hsigma = (TH1*)arr.At(2)->Clone();
778 hsigma->Divide(hmean);
779 arr.SetOwner(kTRUE);
780 arr.Delete();
781 delete his;
782 hsigma->SetXTitle("#beta#gamma");
783 hsigma->SetYTitle("#sigma_{dEdx}/dEdx");
784 hsigma->SetTitle(Form("#sigma_{dEdx}/dEdx_{max} Pad %d",ipad));
785 hsigma->SetName(Form("#sigma_{dEdx}/dEdx_{max} Pad %d",ipad));
786 if (ipad==4){
787 hsigma->SetTitle(Form("#sigma_{dEdx}/dEdx_{max} Full"));
788 hsigma->SetName(Form("#sigma_{dEdx}/dEdx_{max} Full"));
789 }
790 hsigma->SetMarkerStyle(kmimarkers[0]);
791 hsigma->SetMaximum(0.15);
792 hsigma->SetMinimum(0.0);
793 hsigma->Draw();
794 }
795 canvas->SaveAs(Form("%s/dEdxResolMax.eps",outputpath));
796 canvas->SaveAs(Form("%s/dEdxResolMax.gif",outputpath));
797}
798
799
800
801void AliTPCcalibPID::DumpTree(THnSparse * hndim, const char * outname){
802 //
803 // make a tree
804 // the tree will be written to the file - outname
805 //
806 TTreeSRedirector *pcstream = new TTreeSRedirector(outname);
807 //
808 //
809 Double_t position[10];
810 Double_t value;
811 Int_t *bins = new Int_t[10];
812 //
813 //
814 Int_t ndim = hndim->GetNdimensions();
815 //
816 for (Long64_t i = 0; i < hndim->GetNbins(); ++i) {
817 value = hndim->GetBinContent(i, bins);
818 for (Int_t idim = 0; idim < ndim; idim++) {
819 position[idim] = hndim->GetAxis(idim)->GetBinCenter(bins[idim]);
820 }
821 Double_t ty = TMath::Tan(TMath::ASin(position[2]));
822 Double_t tz = position[3]*TMath::Sqrt(1+ty*ty);
823 //
824 (*pcstream)<<"Dump"<<
825 "bincont="<<value<< // bin content
826 "val="<<position[0]<< // parameter difference
827 "dr="<<position[1]<< //drift
828 "ty="<<ty<< //phi
829 "tz="<<tz<< //theta
830 "p2="<<position[2]<< //p2
831 "p3="<<position[3]<< //p3
832 "p="<<position[4]<< //p
833 "bg="<<position[5]<< //bg
834 "ncl="<<position[6]<< //ncl
835 "type="<<position[7]<< //tracklet
836 "tot="<<position[8]<< //is tot
837 "\n";
838 }
839 delete pcstream;
840}
841
842
843void AliTPCcalibPID::DumpTrees(){
844 //
845 // dump the content of histogram to the tree
846 //
847 AliTPCcalibPID *pid = this;
848 DumpTree(pid->GetHistQtot(),"dumpQtot.root");
849 DumpTree(pid->GetHistQmax(),"dumpQmax.root");
850 DumpTree(pid->GetHistRatioQtot(),"dumpRatioQtot.root");
851 DumpTree(pid->GetHistRatioQmax(),"dumpRatioQmax.root");
852}
853
854
855