1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
19 Basic calibration and QA class for the PID of the TPC based on dEdx.
21 The 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.
23 1.a) Fast equalization for different gain in pad regions:
25 TFile f("CalibObjects.root")
26 AliTPCcalibPID *cal = f->Get("TPCCalib")->FindObject("calibPID")
28 cal->GetHistRatioMaxTot()->Projection(0)->Fit("gaus")
29 cal->GetHistRatioTracklet()->Projection(0)->Fit("gaus")
32 .x $ALICE_ROOT/TPC/macros/ConfigOCDB.C
33 AliTPCClusterParam * parcl = AliTPCcalibDB::Instance()->GetClusterParam();
34 *parcl->fQpadMnorm)[ipad] = oldvalue*corrFactor
37 AliCDBMetaData *metaData= new AliCDBMetaData();
38 metaData->SetObjectClassName("AliTPCClusterParam");
39 metaData->SetResponsible("Alexander Kalweit");
40 metaData->SetBeamPeriod(1);
41 metaData->SetAliRootVersion("05-23-02");
42 metaData->SetComment("October runs calibration");
43 AliCDBId id1("TPC/Calib/ClusterParam", runNumber, AliCDBRunRange::Infinity());
44 gStorage = AliCDBManager::Instance()->GetStorage("local://$ALICE_ROOT/OCDB");
45 gStorage->Put(parcl, id1, metaData);
48 2.) Checks should be done with particles on the Fermi-Plateau and whoch leave a long track:
50 cal->GetHistQmax()->GetAxis(6)->SetRangeUser(120,160)
51 cal->GetHistQmax()->GetAxis(4)->SetRangeUser(20,50)
55 cal->GetHistQmax()->Projection(0,1)->Draw() z-dep.
56 cal->GetHistQmax()->Projection(0,2)->Draw() snp-dep.
57 cal->GetHistQmax()->Projection(0,3)->Draw() tgl-dep.
61 Comments to be written here:
62 1. What do we calibrate.
63 2. How to interpret results
65 4. Analysis using debug streamers.
67 Double_t initialParam[] = {3.81470,15.2608,3.25306e-07,2.51791,2.71012}
70 Send comments etc. to: A.Kalweit@gsi.de, marian.ivanov@cern.ch
74 #include "Riostream.h"
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"
100 #include "AliTPCcalibPID.h"
102 #include "TTreeStream.h"
105 ClassImp(AliTPCcalibPID)
108 AliTPCcalibPID::AliTPCcalibPID()
126 fDeDxRatioTruncQtot(0),
127 fDeDxRatioTruncQmax(0)
130 // Empty default cosntructor
132 AliInfo("Default Constructor");
136 AliTPCcalibPID::AliTPCcalibPID(const Text_t *name, const Text_t *title)
154 fDeDxRatioTruncQtot(0),
155 fDeDxRatioTruncQmax(0)
167 fUseShapeNorm = kTRUE;
173 // dE/dx, z, phi, theta, p, bg, ncls, tracklet type
174 Int_t binsQA[8] = {300, 20, 10, 20, 50, 50, 8, 5};
175 Double_t xminQA[8] = {0.2, -1, 0, -1.5, 0.01, 0.1, 60, 0};
176 Double_t xmaxQA[8] = {10., 1, 0.9, 1.5, 50, 500, 160, 5};
180 // dE/dx, z, phi, theta, dEdx, dEdx*dl, ncls, tracklet type
181 Int_t binsRA[9] = {100, 20, 10, 20, 25, 25, 4, 5};
182 Double_t xminRa[9] = {0.5, -1, 0, -1.5, 0.2, 0.2, 60, 0};
183 Double_t xmaxRa[9] = {2, 1, 0.9, 1.5, 5, 5, 160, 5};
185 // z,sin(phi),tan(theta),p,betaGamma,ncls
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);
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);
197 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);
200 BinLogX(fDeDxQmax,4); BinLogX(fDeDxQmax,5);BinLogX(fDeDxQmax,0);
201 BinLogX(fDeDxQtot,4); BinLogX(fDeDxQtot,5);BinLogX(fDeDxQmax,0);
203 BinLogX(fDeDxRatioMaxTot,4); BinLogX(fDeDxRatioMaxTot,5);
204 BinLogX(fDeDxRatioQmax,4); BinLogX(fDeDxRatioQmax,5);
205 BinLogX(fDeDxRatioQtot,4); BinLogX(fDeDxRatioQtot,5);
206 BinLogX(fDeDxRatioTruncQmax,4); BinLogX(fDeDxRatioTruncQmax,5);
207 BinLogX(fDeDxRatioTruncQtot,4); BinLogX(fDeDxRatioTruncQtot,5);
209 fHistNTracks = new TH1F("ntracks","Number of Tracks per Event; number of tracks per event; number of tracks",1001,-0.5,1000.5);
210 fClusters = new TH1F("signal","Number of Clusters per track; number of clusters per track n_{cl}; counts",40,0,160);
211 fPileUp = new TH2F("PileUp","timing plots.; #Delta L ; dEdx signal ",400,-1,1,400,0,200);
212 fLandau = new TH2F("Landau","Landau.; pad type ; cluster charge ",3,-0.5,2.5,400,0,1000);
214 AliInfo("Non Default Constructor");
219 AliTPCcalibPID::~AliTPCcalibPID(){
223 delete fHistNTracks; // histogram showing number of ESD tracks per event
224 delete fClusters; // histogram showing the number of clusters per track
225 delete fPileUp; // histogram which shows correlation between time mismatch and dEdx signal
226 delete fLandau; // histogran which shows Landau distribution for the three pad geometries
228 delete fDeDxQmax; // histogram which shows dEdx (Qmax) as a function of z,sin(phi),tan(theta),p,betaGamma
229 delete fDeDxQtot; // histogram which shows dEdx (Qtot) as a function of z,sin(phi),tan(theta),p,betaGamma
233 delete fDeDxRatioMaxTot; // histogram which shows dEdx ratio (Qmax/Qtot) as a function of z,sin(phi),tan(theta),dEdx,dEdx*dl
234 delete fDeDxRatioQmax; // dEdx ratio (tracklet/track) as a function of z,sin(phi),tan(theta),dEdx,dEdx*dl
235 delete fDeDxRatioQtot; // dEdx ratio (tracklet/track) as a function of z,sin(phi),tan(theta),dEdx,dEdx*dl
236 delete fDeDxRatioTruncQtot; // dEdx ratio (tracklet/track) as a function of z,sin(phi),tan(theta),dEdx,dEdx*dl
237 delete fDeDxRatioTruncQmax; // dEdx ratio (tracklet/track) as a function of z,sin(phi),tan(theta),dEdx,dEdx*d
242 void AliTPCcalibPID::Process(AliESDEvent *event) {
247 Printf("ERROR: ESD not available");
250 const Int_t kMinClustersTracklet = 25;
251 Int_t ntracks=event->GetNumberOfTracks();
252 fHistNTracks->Fill(ntracks);
254 AliESDfriend *esdFriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
256 Printf("ERROR: esdFriend not available");
259 if (esdFriend->TestSkipBit()) return;
263 for (Int_t i=0;i<ntracks;++i) {
265 AliESDtrack *track = event->GetTrack(i);
266 if (!track) continue;
269 AliExternalTrackParam * trackIn = (AliExternalTrackParam *)track->GetInnerParam();
270 AliExternalTrackParam * trackOut = (AliExternalTrackParam *)track->GetOuterParam();
271 if (!trackIn) continue;
272 if (!trackOut) continue;
274 // calculate necessary track parameters
275 //Double_t meanP = 0.5*(trackIn->GetP() + trackOut->GetP());
276 Double_t meanP = trackIn->GetP();
277 //Double_t d = trackIn->GetLinearD(0,0);
278 Int_t nclsDeDx = track->GetTPCNcls();
280 //if (meanP > 0.7 || meanP < 0.2) continue;
281 if (nclsDeDx < 60) continue;
283 // exclude tracks which do not look like primaries or are simply too short or on wrong sectors
285 //if (TMath::Abs(trackIn->GetSnp()) > 3*0.4) continue;
286 //if (TMath::Abs(trackIn->GetZ()) > 150) continue;
287 //if (seed->CookShape(1) > 1) continue;
288 //if (TMath::Abs(trackIn->GetY()) > 20) continue;
289 //if (TMath::Abs(d)>20) continue; // distance to the 0,0; select only tracks which cross chambers under proper angle
290 //if (TMath::Abs(trackIn->GetSnp()) > 0.6) continue;
291 //if (TMath::Abs(trackOut->GetSnp()) > 0.2) continue;
292 if (TMath::Abs(trackIn->GetAlpha()+0.872665)<0.01 || TMath::Abs(trackOut->GetAlpha()+0.872665)<0.01) continue; // Funny sector !
295 AliESDfriendTrack *friendTrack = esdFriend->GetTrack(i);
296 if (!friendTrack) continue;
297 TObject *calibObject;
298 AliTPCseed *seed = 0;
299 for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
300 if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
304 if (meanP > 30 && TMath::Abs(trackIn->GetSnp()) < 0.2) fClusters->Fill(track->GetTPCNcls());
306 Double_t tpcSignalMax = (1/fMIP)*seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,0,159,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
309 //Double_t driftMismatch = seed->GetDriftTimeMismatch(0,0.95,0.5);
310 Double_t driftMismatch = 0;
311 // Double_t drift = 1 - (TMath::Abs(trackIn->GetZ()) + TMath::Abs(trackOut->GetZ()))/500.;
313 // particle identification
314 Double_t mass = 0.105658369;// muon
317 fPileUp->Fill(driftMismatch,tpcSignalMax);
319 fLandau->Fill(0.1,0.6);
321 //var. of interest, z,sin(phi),tan(theta),p,betaGamma,ncls
322 Double_t snpIn = TMath::Abs(trackIn->GetSnp());
323 Double_t snpOut = TMath::Abs(trackOut->GetSnp());
324 Double_t tglIn = trackIn->GetTgl();
325 Double_t tglOut = trackOut->GetTgl();
326 Double_t driftIn = trackIn->GetZ()/250.;
327 Double_t driftOut= trackIn->GetZ()/250.;
331 Float_t dEdxTot[5],dEdxTotFull[5];
332 Float_t dEdxMax[5],dEdxMaxFull[5];
334 for (Int_t itype=0; itype<5;itype++){
335 Int_t row0=0, row1 =159;
336 if (itype==1) {row0=0; row1 = 63;};
337 if (itype==2) {row0= 64; row1=63+64;}
338 if (itype==3) {row0= 63+64; row1=159;}
340 dEdxTot[itype]= (1/fMIP)*seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,0,row0,row1,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
341 dEdxMax[itype]= (1/fMIP)*seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,row0,row1,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
343 dEdxTotFull[itype]= (1/fMIP)*seed->CookdEdxNorm(0.0,0.99,0,row0,row1,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
344 dEdxMaxFull[itype]= (1/fMIP)*seed->CookdEdxNorm(0.0,0.99,1,row0,row1,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
346 dEdxTot[itype]= (1/fMIP)*seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,0,row0,row1);
347 dEdxMax[itype]= (1/fMIP)*seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,1,row0,row1);
349 dEdxTotFull[itype]= (1/fMIP)*seed->CookdEdxAnalytical(0.0,0.99,0,row0,row1);
350 dEdxMaxFull[itype]= (1/fMIP)*seed->CookdEdxAnalytical(0.0,0.99,1,row0,row1);
353 ncl[itype]=(seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,1,row0,row1,2));
358 Float_t wmeanTot=0, wmeanMax=0, sumW=0;
359 Double_t length[3] = {0.75,1,1.5};
362 for (Int_t ipad=0; ipad<3; ipad++){
363 if (ncl[1+ipad]<3) continue;
364 Double_t weight = Double_t(ncl[1+ipad])*TMath::Sqrt(length[ipad]);
365 wmeanTot+=dEdxTot[1+ipad]*weight;
366 wmeanMax+=dEdxMax[1+ipad]*weight;
370 dEdxTot[4]= wmeanTot/sumW;
371 dEdxMax[4]= wmeanMax/sumW;
373 for (Int_t itype=0;itype<5;itype++){
375 Float_t snp=(TMath::Abs(snpIn)+TMath::Abs(snpOut))*0.5;
376 Float_t tgl=(tglIn+tglOut)*0.5;
377 Float_t drift = (driftIn+driftOut)*0.5;
378 if (itype==1) {snp = TMath::Abs(snpIn); tgl = tglIn; drift= driftIn;};
379 if (itype==3) {snp = TMath::Abs(snpOut); tgl = tglOut;drift=driftOut;};
380 if (ncl[itype]<kMinClustersTracklet) continue;
381 Float_t deltaL = TMath::Sqrt(1+snp*snp+tgl*tgl);
383 Double_t vecQmax[8] = {dEdxMax[itype],drift,snp,tgl,meanP,meanP/mass,nclsDeDx, itype};
384 Double_t vecQtot[8] = {dEdxTot[itype],drift,snp,tgl,meanP,meanP/mass,nclsDeDx, itype};
388 Double_t ratioMaxTot = (dEdxTot[itype]>0) ? dEdxMax[itype]/dEdxTot[itype]:0;
389 Double_t ratioTrackletTot = (dEdxTot[0]>0) ? dEdxTot[itype]/dEdxTot[0]:0;
390 Double_t ratioTrackletMax = (dEdxMax[0]>0) ? dEdxMax[itype]/dEdxMax[0]:0;
391 Double_t ratioTrackletTruncTot = (dEdxTot[0]>0) ? dEdxTotFull[itype]/dEdxTot[itype]:0;
392 Double_t ratioTrackletTruncMax = (dEdxMax[0]>0) ? dEdxMaxFull[itype]/dEdxMax[itype]:0;
394 Double_t vecRatioMaxTot[8] = {ratioMaxTot, drift,snp,tgl,dEdxTot[0], dEdxTot[0]*deltaL,nclsDeDx,itype};
395 Double_t vecRatioTrackletTot[8] = {ratioTrackletTot, drift,snp,tgl,dEdxTot[0], dEdxTot[0]*deltaL,nclsDeDx,itype};
396 Double_t vecRatioTrackletMax[8] = {ratioTrackletMax, drift,snp,tgl,dEdxMax[0], dEdxMax[0]*deltaL,nclsDeDx,itype};
397 Double_t vecRatioTrackletTruncTot[8] = {ratioTrackletTruncTot, drift,snp,tgl,dEdxTot[0], dEdxTot[0]*deltaL,nclsDeDx,itype};
398 Double_t vecRatioTrackletTruncMax[8] = {ratioTrackletTruncMax, drift,snp,tgl,dEdxMax[0], dEdxMax[0]*deltaL,nclsDeDx,itype};
399 fDeDxQmax->Fill(vecQmax);
400 fDeDxQtot->Fill(vecQtot);
401 fDeDxRatioMaxTot->Fill(vecRatioMaxTot);
402 fDeDxRatioQmax->Fill(vecRatioTrackletMax);
403 fDeDxRatioQtot->Fill(vecRatioTrackletTot);
404 fDeDxRatioTruncQmax->Fill(vecRatioTrackletTruncMax);
405 fDeDxRatioTruncQtot->Fill(vecRatioTrackletTruncTot);
407 TTreeSRedirector * cstream = GetDebugStreamer();
409 TVectorD vQT(9,vecQtot);
410 TVectorD vQM(9,vecQmax);
411 TVectorD vQMT(9, vecRatioMaxTot);
412 TVectorD vQRT(9,vecRatioTrackletTot);
413 TVectorD vQRM(9,vecRatioTrackletMax);
414 TVectorD vQRTT(9,vecRatioTrackletTruncTot);
415 TVectorD vQRTM(9,vecRatioTrackletTruncMax);
416 TVectorF vNcl(5,ncl);
417 (*cstream) << "dEdx" <<
418 "run="<<fRun<< // run number
419 "event="<<fEvent<< // event number
420 "time="<<fTime<< // time stamp of event
421 "trigger="<<fTrigger<< // trigger
422 "mag="<<fMagF<< // magnetic field
423 "track.="<<seed<< // original seed
424 "trin.="<<trackIn<< // inner param
425 "trout.="<<trackOut<< // outer param
426 "ncl.="<<&vNcl<< // number of clusters
429 "vQT.="<<&vQT<< // trunc mean - total charge
430 "vQM.="<<&vQM<< // trunc mean - max charge
432 "vQMT.="<<&vQMT<< // ratio max/tot
433 "vQRT.="<<&vQRT<< // ratio tracklet/track - total charge
434 "vQRM.="<<&vQRM<< // ratio tracklet/track - max charge
435 "vQRTT.="<<&vQRTT<< // ratio trunc/full - total charge
436 "vQRTM.="<<&vQRTM<< // ratio trunc/full - max charge
446 void AliTPCcalibPID::Analyze() {
454 Long64_t AliTPCcalibPID::Merge(TCollection *li) {
456 TIterator* iter = li->MakeIterator();
457 AliTPCcalibPID* cal = 0;
459 while ((cal = (AliTPCcalibPID*)iter->Next())) {
460 if (!cal->InheritsFrom(AliTPCcalibPID::Class())) {
461 Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
465 if (cal->GetHistNTracks()) fHistNTracks->Add(cal->GetHistNTracks());
466 if (cal->GetHistClusters()) fClusters->Add(cal->GetHistClusters());
467 if (cal->GetHistPileUp()) fPileUp->Add(cal->GetHistPileUp());
468 if (cal->GetHistLandau()) fLandau->Add(cal->GetHistLandau());
470 if (cal->GetHistQmax()) fDeDxQmax->Add(cal->GetHistQmax());
471 if (cal->GetHistQtot()) fDeDxQtot->Add(cal->GetHistQtot());
472 if (cal->GetHistRatioMaxTot()) fDeDxRatioMaxTot->Add(cal->GetHistRatioMaxTot());
473 if (cal->GetHistRatioQmax()) fDeDxRatioQmax->Add(cal->GetHistRatioQmax());
474 if (cal->GetHistRatioQtot()) fDeDxRatioQtot->Add(cal->GetHistRatioQtot());
475 if (cal->GetHistRatioTruncQmax()) fDeDxRatioTruncQmax->Add(cal->GetHistRatioTruncQmax());
476 if (cal->GetHistRatioTruncQtot()) fDeDxRatioTruncQtot->Add(cal->GetHistRatioTruncQtot());
485 void AliTPCcalibPID::BinLogX(THnSparse *h, Int_t axisDim) {
487 // Method for the correct logarithmic binning of histograms
489 TAxis *axis = h->GetAxis(axisDim);
490 int bins = axis->GetNbins();
492 Double_t from = axis->GetXmin();
493 Double_t to = axis->GetXmax();
494 Double_t *newBins = new Double_t[bins + 1];
497 Double_t factor = pow(to/from, 1./bins);
499 for (int i = 1; i <= bins; i++) {
500 newBins[i] = factor * newBins[i-1];
502 axis->Set(bins, newBins);
511 void AliTPCcalibPID::MakeReport(const char *outputpath) {
513 // Make a standard QA plots
515 for (Int_t ipad=0;ipad<4;ipad++){
516 DrawRatioTot(ipad,outputpath);
517 DrawRatioMax(ipad,outputpath);
519 DrawRatiodEdx(0.5,3,outputpath);
520 DrawResolBGQtot(140,160,1,40,outputpath);
521 DrawResolBGQmax(140,160,1,40,outputpath);
525 void AliTPCcalibPID::DrawRatioTot(Int_t ipad, const char* outputpath){
527 // Draw default ratio histogram for given pad type
528 // ipad - 0 - short pads
532 // Int_t kmicolors[10]={1,2,3,4,6,7,8,9,10,11};
533 Int_t kmimarkers[10]={21,22,23,24,25,26,27,28,29,30};
534 AliTPCcalibPID * pid = this;
535 TCanvas *canvas= new TCanvas(Form("QtotRatio%d)",ipad),Form("QtotRatio%d)",ipad),600,800);
537 pid->GetHistRatioQtot()->GetAxis(7)->SetRange(ipad+2,ipad+2);
539 TH1 * his0 = pid->GetHistRatioQtot()->Projection(0);
540 his0->SetXTitle("dEdx_{tracklet}/dEdx_{track}");
542 his0->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
546 TH1 * hprofz = (TH1*) (pid->GetHistRatioQtot()->Projection(0,1)->ProfileX());
547 hprofz->SetXTitle("drift length");
548 hprofz->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
549 hprofz->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
550 hprofz->SetMarkerStyle(kmimarkers[0]);
551 hprofz->SetMaximum(1.1);
552 hprofz->SetMinimum(0.9);
556 TH1 * hprofphi = (TH1*) (pid->GetHistRatioQtot()->Projection(0,2)->ProfileX());
557 hprofphi->SetXTitle("sin(#phi)");
558 hprofphi->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
559 hprofphi->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
560 hprofphi->SetMarkerStyle(kmimarkers[0]);
561 hprofz->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
562 hprofphi->SetMaximum(1.1);
563 hprofphi->SetMinimum(0.9);
567 TH1 * hproftheta = (TH1*) (pid->GetHistRatioQtot()->Projection(0,3)->ProfileX());
568 hproftheta->SetXTitle("tan(#theta)");
569 hproftheta->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
570 hproftheta->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
571 hproftheta->SetMarkerStyle(kmimarkers[0]);
572 hproftheta->SetMaximum(1.1);
573 hproftheta->SetMinimum(0.9);
577 TH1 * hprofdedx = (TH1*) (pid->GetHistRatioQtot()->Projection(0,4)->ProfileX());
578 hprofdedx->SetXTitle("dEdx_{track}");
579 hprofdedx->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
580 hprofdedx->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
581 hprofdedx->SetMarkerStyle(kmimarkers[0]);
582 hprofdedx->SetMaximum(1.1);
583 hprofdedx->SetMinimum(0.9);
587 TH1 * hprofdedxL = (TH1*) (pid->GetHistRatioQtot()->Projection(0,5)->ProfileX());
588 hprofdedxL->SetXTitle("dEdx_{track}#Delta_{x}");
589 hprofdedxL->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
590 hprofdedxL->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
591 hprofdedxL->SetMarkerStyle(kmimarkers[0]);
592 hprofdedxL->SetMaximum(1.1);
593 hprofdedxL->SetMinimum(0.9);
598 canvas->SaveAs(Form("%s/QtotRatioType%d.eps",outputpath,ipad));
599 canvas->SaveAs(Form("%s/QtotRatioType%d.gif",outputpath,ipad));
602 void AliTPCcalibPID::DrawRatioMax(Int_t ipad, const char* outputpath){
604 // Draw default ration histogram for given pad type
605 // ipad - 0 - short pads
609 // Int_t kmicolors[10]={1,2,3,4,6,7,8,9,10,11};
610 Int_t kmimarkers[10]={21,22,23,24,25,26,27,28,29,30};
611 AliTPCcalibPID * pid = this;
612 TCanvas *canvas= new TCanvas(Form("QmaxRatio%d)",ipad),Form("QmaxRatio%d)",ipad),600,800);
614 pid->GetHistRatioQmax()->GetAxis(7)->SetRange(ipad+2,ipad+2);
616 TH1 * his0 = pid->GetHistRatioQmax()->Projection(0);
617 his0->SetXTitle("dEdx_{tracklet}/dEdx_{track}");
619 his0->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
623 TH1 * hprofz = (TH1*) (pid->GetHistRatioQmax()->Projection(0,1)->ProfileX());
624 hprofz->SetXTitle("drift length");
625 hprofz->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
626 hprofz->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
627 hprofz->SetMarkerStyle(kmimarkers[0]);
628 hprofz->SetMaximum(1.1);
629 hprofz->SetMinimum(0.9);
633 TH1 * hprofphi = (TH1*) (pid->GetHistRatioQmax()->Projection(0,2)->ProfileX());
634 hprofphi->SetXTitle("sin(#phi)");
635 hprofphi->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
636 hprofphi->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
637 hprofphi->SetMarkerStyle(kmimarkers[0]);
638 hprofphi->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
639 hprofphi->SetMaximum(1.1);
640 hprofphi->SetMinimum(0.9);
644 TH1 * hproftheta = (TH1*) (pid->GetHistRatioQmax()->Projection(0,3)->ProfileX());
645 hproftheta->SetXTitle("tan(#theta)");
646 hproftheta->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
647 hproftheta->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
648 hproftheta->SetMarkerStyle(kmimarkers[0]);
649 hproftheta->SetMaximum(1.1);
650 hproftheta->SetMinimum(0.9);
654 TH1 * hprofdedx = (TH1*) (pid->GetHistRatioQmax()->Projection(0,4)->ProfileX());
655 hprofdedx->SetXTitle("dEdx_{track}");
656 hprofdedx->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
657 hprofdedx->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
658 hprofdedx->SetMarkerStyle(kmimarkers[0]);
659 hprofdedx->SetMaximum(1.1);
660 hprofdedx->SetMinimum(0.9);
664 TH1 * hprofdedxL = (TH1*) (pid->GetHistRatioQmax()->Projection(0,5)->ProfileX());
665 hprofdedxL->SetXTitle("dEdx_{track}#Delta_{x}");
666 hprofdedxL->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
667 hprofdedxL->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
668 hprofdedxL->SetMarkerStyle(kmimarkers[0]);
669 hprofdedxL->SetMaximum(1.1);
670 hprofdedxL->SetMinimum(0.9);
674 canvas->SaveAs(Form("%s/QmaxRatioType%d.eps",outputpath,ipad));
675 canvas->SaveAs(Form("%s/QmaxRatioType%d.gif",outputpath,ipad));
680 void AliTPCcalibPID::DrawRatiodEdx(Float_t demin, Float_t demax, const char* outputpath){
685 // Int_t kmicolors[10]={1,2,3,4,6,7,8,9,10,11};
686 Int_t kmimarkers[10]={21,22,23,24,25,26,27,28,29,30};
687 AliTPCcalibPID * pid = this;
688 TCanvas *canvas= new TCanvas("QRatiodEdx","QRatiodEdx",600,800);
690 canvas->SetLogx(kTRUE);
692 for (Int_t ipad=0;ipad<4;ipad++){
693 pid->GetHistRatioQmax()->GetAxis(7)->SetRange(ipad+2,ipad+2);
694 pid->GetHistRatioQtot()->GetAxis(7)->SetRange(ipad+2,ipad+2);
695 pid->GetHistRatioQmax()->GetAxis(5)->SetRangeUser(demin,demax);
696 pid->GetHistRatioQtot()->GetAxis(5)->SetRangeUser(demin,demax);
698 canvas->cd(ipad*2+1)->SetLogx(kFALSE);
699 hprofP = (TH1*) (pid->GetHistRatioQmax()->Projection(0,5)->ProfileX());
700 hprofP->SetXTitle("dEdx_{track}");
701 hprofP->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
702 hprofP->SetTitle(Form("Q_{max} dEdx_{tracklet}/dEdx_{track} type %d",ipad));
703 hprofP->SetMarkerStyle(kmimarkers[0]);
704 hprofP->SetMaximum(1.1);
705 hprofP->SetMinimum(0.9);
708 canvas->cd(ipad*2+2)->SetLogx(kFALSE);
709 hprofP = (TH1*) (pid->GetHistRatioQtot()->Projection(0,5)->ProfileX());
710 hprofP->SetXTitle("dEdx_{track}");
711 hprofP->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
712 hprofP->SetTitle(Form("Q_{tot} dEdx_{tracklet}/dEdx_{track} type %d",ipad));
713 hprofP->SetMarkerStyle(kmimarkers[0]);
714 hprofP->SetMaximum(1.1);
715 hprofP->SetMinimum(0.9);
720 canvas->SaveAs(Form("%s/QratiodEdx.eps",outputpath));
721 canvas->SaveAs(Form("%s/QratiodEdx.gif",outputpath));
726 void AliTPCcalibPID::DrawResolBGQtot(Int_t minClusters, Int_t maxClusters, Float_t minp, Float_t maxp, const char *outputpath, Bool_t resol){
730 AliTPCcalibPID * pid = this;
731 Int_t kmimarkers[10]={21,22,23,24,25,26,27,28,29,30};
739 pid->GetHistQtot()->GetAxis(6)->SetRangeUser(minClusters,maxClusters);
740 pid->GetHistQtot()->GetAxis(5)->SetRangeUser(1,10000);
741 pid->GetHistQtot()->GetAxis(4)->SetRangeUser(minp,maxp);
742 TCanvas *canvas= new TCanvas("dEdxResolQ_{Tot}","dEdxResolQ_{Tot}",800,600);
747 for (Int_t ipad=0;ipad<5;ipad++){
748 canvas->cd(1+ipad)->SetLogx(kTRUE);
749 if (ipad<4) pid->GetHistQtot()->GetAxis(7)->SetRange(ipad+2,ipad+2);
750 if (ipad==4) pid->GetHistQtot()->GetAxis(7)->SetRange(1,1);
751 his = (TH2*)(pid->GetHistQtot()->Projection(0,5));
752 his->FitSlicesY(0,0,-1,0,"QNR",&arr);
754 hmean = (TH1*)arr.At(1);
755 hsigma = (TH1*)arr.At(2)->Clone();
756 hsigma->SetMaximum(0.11);
757 hsigma->SetMinimum(0.4);
758 hsigma->Divide(hmean);
760 hsigma = (TH1*)arr.At(1)->Clone();
761 hsigma->SetMaximum(2);
762 hsigma->SetMinimum(0.5);
768 hsigma->SetXTitle("#beta#gamma");
769 hsigma->SetYTitle("#sigma_{dEdx}/dEdx");
770 hsigma->SetTitle(Form("#sigma_{dEdx}/dEdx_{tot} Pad %d",ipad));
771 hsigma->SetName(Form("#sigma_{dEdx}/dEdx_{tot} Pad %d",ipad));
773 hsigma->SetTitle(Form("#sigma_{dEdx}/dEdx_{tot} Full"));
774 hsigma->SetName(Form("#sigma_{dEdx}/dEdx_{tot} Full"));
776 hsigma->SetMarkerStyle(kmimarkers[0]);
780 canvas->SaveAs(Form("%s/dEdxResolTot.eps",outputpath));
781 canvas->SaveAs(Form("%s/dEdxResolTot.gif",outputpath));
783 canvas->SaveAs(Form("%s/dEdxBGTot.eps",outputpath));
784 canvas->SaveAs(Form("%s/dEdxBGTot.gif",outputpath));
788 void AliTPCcalibPID::DrawResolBGQmax(Int_t minClusters, Int_t maxClusters, Float_t minp, Float_t maxp, const char *outputpath, Bool_t resol){
790 // Int_t minClusters=140; Int_t maxClusters=200; const char *outputpath="picPID06"
792 Int_t kmimarkers[10]={21,22,23,24,25,26,27,28,29,30};
793 AliTPCcalibPID * pid = this;
800 pid->GetHistQmax()->GetAxis(6)->SetRangeUser(minClusters,maxClusters);
801 pid->GetHistQmax()->GetAxis(4)->SetRangeUser(minp,maxp);
802 pid->GetHistQmax()->GetAxis(5)->SetRangeUser(1,10000);
803 TCanvas *canvas= new TCanvas("dEdxResolQ_{Max}","dEdxResolQ_{Max}",800,600);
808 for (Int_t ipad=0;ipad<5;ipad++){
809 canvas->cd(1+ipad)->SetLogx(kTRUE);
810 if (ipad<4) pid->GetHistQmax()->GetAxis(7)->SetRange(ipad+2,ipad+2);
811 if (ipad==4) pid->GetHistQmax()->GetAxis(7)->SetRange(1,1);
812 his = (TH2*)(pid->GetHistQmax()->Projection(0,5));
813 his->FitSlicesY(0,0,-1,0,"QNR",&arr);
815 hmean = (TH1*)arr.At(1);
816 hsigma = (TH1*)arr.At(2)->Clone();
817 hsigma->SetMaximum(0.11);
818 hsigma->SetMinimum(0.4);
819 hsigma->Divide(hmean);
821 hsigma = (TH1*)arr.At(1)->Clone();
822 hsigma->SetMaximum(2.);
823 hsigma->SetMinimum(0.5);
828 hsigma->SetXTitle("#beta#gamma");
829 hsigma->SetYTitle("#sigma_{dEdx}/dEdx");
830 hsigma->SetTitle(Form("#sigma_{dEdx}/dEdx_{max} Pad %d",ipad));
831 hsigma->SetName(Form("#sigma_{dEdx}/dEdx_{max} Pad %d",ipad));
833 hsigma->SetTitle(Form("#sigma_{dEdx}/dEdx_{max} Full"));
834 hsigma->SetName(Form("#sigma_{dEdx}/dEdx_{max} Full"));
836 hsigma->SetMarkerStyle(kmimarkers[0]);
841 canvas->SaveAs(Form("%s/dEdxResolMax.eps",outputpath));
842 canvas->SaveAs(Form("%s/dEdxResolMax.gif",outputpath));
844 canvas->SaveAs(Form("%s/dEdxBGMax.eps",outputpath));
845 canvas->SaveAs(Form("%s/dEdxBGMax.gif",outputpath));
853 void AliTPCcalibPID::DumpTree(THnSparse * hndim, const char * outname){
856 // the tree will be written to the file - outname
861 TTreeSRedirector *pcstream = new TTreeSRedirector(outname);
864 Double_t position[10];
869 const Float_t rmsy0=0.5, rmsz0=1;
871 Int_t ndim = hndim->GetNdimensions();
873 AliTPCParam * param = AliTPCcalibDB::Instance()->GetParameters();
874 for (Long64_t i = 0; i < hndim->GetNbins(); ++i) {
875 value = hndim->GetBinContent(i, bins);
876 for (Int_t idim = 0; idim < ndim; idim++) {
877 position[idim] = hndim->GetAxis(idim)->GetBinCenter(bins[idim]);
881 if (TMath::Abs(position[7]-1.5)<0.1) sector=0; // inner sector
882 if (TMath::Abs(position[7]-3.5)<0.1) row=96; // long pads
884 Double_t ty = TMath::Tan(TMath::ASin(position[2]));
885 Double_t tz = position[3]*TMath::Sqrt(1+ty*ty);
886 Double_t padLength= param->GetPadPitchLength(sector,row);
887 Double_t padWidth = param->GetPadPitchWidth(sector);
888 Double_t zwidth = param->GetZWidth();
889 Double_t zlength=250-TMath::Abs(position[1])*250.;
890 // diffusion in pad, time bin units
891 Double_t diffT=TMath::Sqrt(zlength)*param->GetDiffT()/padWidth;
892 Double_t diffL=TMath::Sqrt(zlength)*param->GetDiffL()/zwidth;
894 // transform angular effect to pad units
895 Double_t pky = ty*padLength/padWidth;
896 Double_t pkz = tz*padLength/zwidth;
899 Double_t sy = TMath::Sqrt(rmsy0*rmsy0+diffT*diffT+pky*pky/12.);
900 Double_t sz = TMath::Sqrt(rmsz0*rmsz0+diffL*diffL+pkz*pkz/12.);
901 Int_t ipad = TMath::Nint(position[7]-1.5);
903 (*pcstream)<<"Dump"<<
904 "bincont="<<value<< // bin content
905 "val="<<position[0]<< // parameter difference
906 "ipad="<<ipad<< // pad type
907 "dr="<<position[1]<< //drift
910 "p2="<<position[2]<< //p2
911 "p3="<<position[3]<< //p3
912 "p="<<position[4]<< //p
913 "bg="<<position[5]<< //bg
914 "ncl="<<position[6]<< //ncl
915 "type="<<position[7]<< //tracklet
916 "tot="<<position[8]<< //is tot
917 "sy="<<sy<< // sigma y
918 "sz="<<sz<< // sigma z
919 "pky="<<pky<< // ky - angle in pad units
920 "pkz="<<pkz<< // kz - angle in pad units
929 void AliTPCcalibPID::DumpTrees(){
931 // dump the content of histogram to the tree
933 AliTPCcalibPID *pid = this;
934 DumpTree(pid->GetHistQtot(),"dumpQtot.root");
935 DumpTree(pid->GetHistQmax(),"dumpQmax.root");
936 DumpTree(pid->GetHistRatioQtot(),"dumpRatioQtot.root");
937 DumpTree(pid->GetHistRatioQmax(),"dumpRatioQmax.root");