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");
262 for (Int_t i=0;i<ntracks;++i) {
264 AliESDtrack *track = event->GetTrack(i);
265 if (!track) continue;
268 AliExternalTrackParam * trackIn = (AliExternalTrackParam *)track->GetInnerParam();
269 AliExternalTrackParam * trackOut = (AliExternalTrackParam *)track->GetOuterParam();
270 if (!trackIn) continue;
271 if (!trackOut) continue;
273 // calculate necessary track parameters
274 //Double_t meanP = 0.5*(trackIn->GetP() + trackOut->GetP());
275 Double_t meanP = trackIn->GetP();
276 //Double_t d = trackIn->GetLinearD(0,0);
277 Int_t NclsDeDx = track->GetTPCNcls();
279 //if (meanP > 0.7 || meanP < 0.2) continue;
280 if (NclsDeDx < 60) continue;
282 // exclude tracks which do not look like primaries or are simply too short or on wrong sectors
284 //if (TMath::Abs(trackIn->GetSnp()) > 3*0.4) continue;
285 //if (TMath::Abs(trackIn->GetZ()) > 150) continue;
286 //if (seed->CookShape(1) > 1) continue;
287 //if (TMath::Abs(trackIn->GetY()) > 20) continue;
288 //if (TMath::Abs(d)>20) continue; // distance to the 0,0; select only tracks which cross chambers under proper angle
289 //if (TMath::Abs(trackIn->GetSnp()) > 0.6) continue;
290 //if (TMath::Abs(trackOut->GetSnp()) > 0.2) continue;
291 if (TMath::Abs(trackIn->GetAlpha()+0.872665)<0.01 || TMath::Abs(trackOut->GetAlpha()+0.872665)<0.01) continue; // Funny sector !
294 AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
295 if (!friendTrack) continue;
296 TObject *calibObject;
297 AliTPCseed *seed = 0;
298 for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
299 if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
303 if (meanP > 30 && TMath::Abs(trackIn->GetSnp()) < 0.2) fClusters->Fill(track->GetTPCNcls());
305 Double_t TPCsignalMax = (1/fMIP)*seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,0,159,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
308 //Double_t driftMismatch = seed->GetDriftTimeMismatch(0,0.95,0.5);
309 Double_t driftMismatch = 0;
310 // Double_t drift = 1 - (TMath::Abs(trackIn->GetZ()) + TMath::Abs(trackOut->GetZ()))/500.;
312 // particle identification
313 Double_t mass = 0.105658369;// muon
316 fPileUp->Fill(driftMismatch,TPCsignalMax);
318 fLandau->Fill(0.1,0.6);
320 //var. of interest, z,sin(phi),tan(theta),p,betaGamma,ncls
321 Double_t snpIn = TMath::Abs(trackIn->GetSnp());
322 Double_t snpOut = TMath::Abs(trackOut->GetSnp());
323 Double_t tglIn = trackIn->GetTgl();
324 Double_t tglOut = trackOut->GetTgl();
325 Double_t driftIn = trackIn->GetZ()/250.;
326 Double_t driftOut= trackIn->GetZ()/250.;
330 Float_t dEdxTot[5],dEdxTotFull[5];
331 Float_t dEdxMax[5],dEdxMaxFull[5];
333 for (Int_t itype=0; itype<5;itype++){
334 Int_t row0=0, row1 =159;
335 if (itype==1) {row0=0; row1 = 63;};
336 if (itype==2) {row0= 64; row1=63+64;}
337 if (itype==3) {row0= 63+64; row1=159;}
339 dEdxTot[itype]= (1/fMIP)*seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,0,row0,row1,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
340 dEdxMax[itype]= (1/fMIP)*seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,row0,row1,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
342 dEdxTotFull[itype]= (1/fMIP)*seed->CookdEdxNorm(0.0,0.99,0,row0,row1,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
343 dEdxMaxFull[itype]= (1/fMIP)*seed->CookdEdxNorm(0.0,0.99,1,row0,row1,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
345 dEdxTot[itype]= (1/fMIP)*seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,0,row0,row1);
346 dEdxMax[itype]= (1/fMIP)*seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,1,row0,row1);
348 dEdxTotFull[itype]= (1/fMIP)*seed->CookdEdxAnalytical(0.0,0.99,0,row0,row1);
349 dEdxMaxFull[itype]= (1/fMIP)*seed->CookdEdxAnalytical(0.0,0.99,1,row0,row1);
352 ncl[itype]=(seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,1,row0,row1,2));
357 Float_t wmeanTot=0, wmeanMax=0, sumW=0;
358 Double_t length[3] = {0.75,1,1.5};
361 for (Int_t ipad=0; ipad<3; ipad++){
362 if (ncl[1+ipad]<3) continue;
363 Double_t weight = Double_t(ncl[1+ipad])*TMath::Sqrt(length[ipad]);
364 wmeanTot+=dEdxTot[1+ipad]*weight;
365 wmeanMax+=dEdxMax[1+ipad]*weight;
369 dEdxTot[4]= wmeanTot/sumW;
370 dEdxMax[4]= wmeanMax/sumW;
372 for (Int_t itype=0;itype<5;itype++){
374 Float_t snp=(TMath::Abs(snpIn)+TMath::Abs(snpOut))*0.5;
375 Float_t tgl=(tglIn+tglOut)*0.5;
376 Float_t drift = (driftIn+driftOut)*0.5;
377 if (itype==1) {snp = TMath::Abs(snpIn); tgl = tglIn; drift= driftIn;};
378 if (itype==3) {snp = TMath::Abs(snpOut); tgl = tglOut;drift=driftOut;};
379 if (ncl[itype]<kMinClustersTracklet) continue;
380 Float_t deltaL = TMath::Sqrt(1+snp*snp+tgl*tgl);
382 Double_t vecQmax[8] = {dEdxMax[itype],drift,snp,tgl,meanP,meanP/mass,NclsDeDx, itype};
383 Double_t vecQtot[8] = {dEdxTot[itype],drift,snp,tgl,meanP,meanP/mass,NclsDeDx, itype};
387 Double_t ratioMaxTot = (dEdxTot[itype]>0) ? dEdxMax[itype]/dEdxTot[itype]:0;
388 Double_t ratioTrackletTot = (dEdxTot[0]>0) ? dEdxTot[itype]/dEdxTot[0]:0;
389 Double_t ratioTrackletMax = (dEdxMax[0]>0) ? dEdxMax[itype]/dEdxMax[0]:0;
390 Double_t ratioTrackletTruncTot = (dEdxTot[0]>0) ? dEdxTotFull[itype]/dEdxTot[itype]:0;
391 Double_t ratioTrackletTruncMax = (dEdxMax[0]>0) ? dEdxMaxFull[itype]/dEdxMax[itype]:0;
393 Double_t vecRatioMaxTot[8] = {ratioMaxTot, drift,snp,tgl,dEdxTot[0], dEdxTot[0]*deltaL,NclsDeDx,itype};
394 Double_t vecRatioTrackletTot[8] = {ratioTrackletTot, drift,snp,tgl,dEdxTot[0], dEdxTot[0]*deltaL,NclsDeDx,itype};
395 Double_t vecRatioTrackletMax[8] = {ratioTrackletMax, drift,snp,tgl,dEdxMax[0], dEdxMax[0]*deltaL,NclsDeDx,itype};
396 Double_t vecRatioTrackletTruncTot[8] = {ratioTrackletTruncTot, drift,snp,tgl,dEdxTot[0], dEdxTot[0]*deltaL,NclsDeDx,itype};
397 Double_t vecRatioTrackletTruncMax[8] = {ratioTrackletTruncMax, drift,snp,tgl,dEdxMax[0], dEdxMax[0]*deltaL,NclsDeDx,itype};
398 fDeDxQmax->Fill(vecQmax);
399 fDeDxQtot->Fill(vecQtot);
400 fDeDxRatioMaxTot->Fill(vecRatioMaxTot);
401 fDeDxRatioQmax->Fill(vecRatioTrackletTot);
402 fDeDxRatioQtot->Fill(vecRatioTrackletMax);
403 fDeDxRatioTruncQmax->Fill(vecRatioTrackletTruncTot);
404 fDeDxRatioTruncQtot->Fill(vecRatioTrackletTruncMax);
406 TTreeSRedirector * cstream = GetDebugStreamer();
408 TVectorD vQT(9,vecQtot);
409 TVectorD vQM(9,vecQmax);
410 TVectorD vQMT(9, vecRatioMaxTot);
411 TVectorD vQRT(9,vecRatioTrackletTot);
412 TVectorD vQRM(9,vecRatioTrackletMax);
413 TVectorD vQRTT(9,vecRatioTrackletTruncTot);
414 TVectorD vQRTM(9,vecRatioTrackletTruncMax);
415 TVectorF vNcl(5,ncl);
416 (*cstream) << "dEdx" <<
417 "run="<<fRun<< // run number
418 "event="<<fEvent<< // event number
419 "time="<<fTime<< // time stamp of event
420 "trigger="<<fTrigger<< // trigger
421 "mag="<<fMagF<< // magnetic field
422 "track.="<<seed<< // original seed
423 "trin.="<<trackIn<< // inner param
424 "trout.="<<trackOut<< // outer param
425 "ncl.="<<&vNcl<< // number of clusters
428 "vQT.="<<&vQT<< // trunc mean - total charge
429 "vQM.="<<&vQM<< // trunc mean - max charge
431 "vQMT.="<<&vQMT<< // ratio max/tot
432 "vQRT.="<<&vQRT<< // ratio tracklet/track - total charge
433 "vQRM.="<<&vQRM<< // ratio tracklet/track - max charge
434 "vQRTT.="<<&vQRTT<< // ratio trunc/full - total charge
435 "vQRTM.="<<&vQRTM<< // ratio trunc/full - max charge
445 void AliTPCcalibPID::Analyze() {
453 Long64_t AliTPCcalibPID::Merge(TCollection *li) {
455 TIterator* iter = li->MakeIterator();
456 AliTPCcalibPID* cal = 0;
458 while ((cal = (AliTPCcalibPID*)iter->Next())) {
459 if (!cal->InheritsFrom(AliTPCcalibPID::Class())) {
460 Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
464 if (cal->GetHistNTracks()) fHistNTracks->Add(cal->GetHistNTracks());
465 if (cal->GetHistClusters()) fClusters->Add(cal->GetHistClusters());
466 if (cal->GetHistPileUp()) fPileUp->Add(cal->GetHistPileUp());
467 if (cal->GetHistLandau()) fLandau->Add(cal->GetHistLandau());
469 if (cal->GetHistQmax()) fDeDxQmax->Add(cal->GetHistQmax());
470 if (cal->GetHistQtot()) fDeDxQtot->Add(cal->GetHistQtot());
471 if (cal->GetHistRatioMaxTot()) fDeDxRatioMaxTot->Add(cal->GetHistRatioMaxTot());
472 if (cal->GetHistRatioQmax()) fDeDxRatioQmax->Add(cal->GetHistRatioQmax());
473 if (cal->GetHistRatioQtot()) fDeDxRatioQtot->Add(cal->GetHistRatioQtot());
474 if (cal->GetHistRatioTruncQmax()) fDeDxRatioTruncQmax->Add(cal->GetHistRatioTruncQmax());
475 if (cal->GetHistRatioTruncQtot()) fDeDxRatioTruncQtot->Add(cal->GetHistRatioTruncQtot());
484 void AliTPCcalibPID::BinLogX(THnSparse *h, Int_t axisDim) {
486 // Method for the correct logarithmic binning of histograms
488 TAxis *axis = h->GetAxis(axisDim);
489 int bins = axis->GetNbins();
491 Double_t from = axis->GetXmin();
492 Double_t to = axis->GetXmax();
493 Double_t *new_bins = new Double_t[bins + 1];
496 Double_t factor = pow(to/from, 1./bins);
498 for (int i = 1; i <= bins; i++) {
499 new_bins[i] = factor * new_bins[i-1];
501 axis->Set(bins, new_bins);
510 void AliTPCcalibPID::MakeReport(const char *outputpath) {
512 // Make a standard QA plots
514 for (Int_t ipad=0;ipad<4;ipad++){
515 DrawRatioTot(ipad,outputpath);
516 DrawRatioMax(ipad,outputpath);
518 DrawRatiodEdx(0.5,3,outputpath);
519 DrawResolBGQtot(140,160,1,40,outputpath);
520 DrawResolBGQmax(140,160,1,40,outputpath);
524 void AliTPCcalibPID::DrawRatioTot(Int_t ipad, const char* outputpath){
526 // Draw default ratio histogram for given pad type
527 // ipad - 0 - short pads
531 // Int_t kmicolors[10]={1,2,3,4,6,7,8,9,10,11};
532 Int_t kmimarkers[10]={21,22,23,24,25,26,27,28,29,30};
533 AliTPCcalibPID * pid = this;
534 TCanvas *canvas= new TCanvas(Form("QtotRatio%d)",ipad),Form("QtotRatio%d)",ipad),600,800);
536 pid->GetHistRatioQtot()->GetAxis(7)->SetRange(ipad+2,ipad+2);
538 TH1 * his0 = pid->GetHistRatioQtot()->Projection(0);
539 his0->SetXTitle("dEdx_{tracklet}/dEdx_{track}");
541 his0->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
545 TH1 * hprofz = (TH1*) (pid->GetHistRatioQtot()->Projection(0,1)->ProfileX());
546 hprofz->SetXTitle("drift length");
547 hprofz->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
548 hprofz->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
549 hprofz->SetMarkerStyle(kmimarkers[0]);
550 hprofz->SetMaximum(1.1);
551 hprofz->SetMinimum(0.9);
555 TH1 * hprofphi = (TH1*) (pid->GetHistRatioQtot()->Projection(0,2)->ProfileX());
556 hprofphi->SetXTitle("sin(#phi)");
557 hprofphi->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
558 hprofphi->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
559 hprofphi->SetMarkerStyle(kmimarkers[0]);
560 hprofz->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
561 hprofphi->SetMaximum(1.1);
562 hprofphi->SetMinimum(0.9);
566 TH1 * hproftheta = (TH1*) (pid->GetHistRatioQtot()->Projection(0,3)->ProfileX());
567 hproftheta->SetXTitle("tan(#theta)");
568 hproftheta->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
569 hproftheta->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
570 hproftheta->SetMarkerStyle(kmimarkers[0]);
571 hproftheta->SetMaximum(1.1);
572 hproftheta->SetMinimum(0.9);
576 TH1 * hprofdedx = (TH1*) (pid->GetHistRatioQtot()->Projection(0,4)->ProfileX());
577 hprofdedx->SetXTitle("dEdx_{track}");
578 hprofdedx->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
579 hprofdedx->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
580 hprofdedx->SetMarkerStyle(kmimarkers[0]);
581 hprofdedx->SetMaximum(1.1);
582 hprofdedx->SetMinimum(0.9);
586 TH1 * hprofdedxL = (TH1*) (pid->GetHistRatioQtot()->Projection(0,5)->ProfileX());
587 hprofdedxL->SetXTitle("dEdx_{track}#Delta_{x}");
588 hprofdedxL->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
589 hprofdedxL->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
590 hprofdedxL->SetMarkerStyle(kmimarkers[0]);
591 hprofdedxL->SetMaximum(1.1);
592 hprofdedxL->SetMinimum(0.9);
597 canvas->SaveAs(Form("%s/QtotRatioType%d.eps",outputpath,ipad));
598 canvas->SaveAs(Form("%s/QtotRatioType%d.gif",outputpath,ipad));
601 void AliTPCcalibPID::DrawRatioMax(Int_t ipad, const char* outputpath){
603 // Draw default ration histogram for given pad type
604 // ipad - 0 - short pads
608 // Int_t kmicolors[10]={1,2,3,4,6,7,8,9,10,11};
609 Int_t kmimarkers[10]={21,22,23,24,25,26,27,28,29,30};
610 AliTPCcalibPID * pid = this;
611 TCanvas *canvas= new TCanvas(Form("QmaxRatio%d)",ipad),Form("QmaxRatio%d)",ipad),600,800);
613 pid->GetHistRatioQmax()->GetAxis(7)->SetRange(ipad+2,ipad+2);
615 TH1 * his0 = pid->GetHistRatioQmax()->Projection(0);
616 his0->SetXTitle("dEdx_{tracklet}/dEdx_{track}");
618 his0->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
622 TH1 * hprofz = (TH1*) (pid->GetHistRatioQmax()->Projection(0,1)->ProfileX());
623 hprofz->SetXTitle("drift length");
624 hprofz->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
625 hprofz->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
626 hprofz->SetMarkerStyle(kmimarkers[0]);
627 hprofz->SetMaximum(1.1);
628 hprofz->SetMinimum(0.9);
632 TH1 * hprofphi = (TH1*) (pid->GetHistRatioQmax()->Projection(0,2)->ProfileX());
633 hprofphi->SetXTitle("sin(#phi)");
634 hprofphi->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
635 hprofphi->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
636 hprofphi->SetMarkerStyle(kmimarkers[0]);
637 hprofphi->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
638 hprofphi->SetMaximum(1.1);
639 hprofphi->SetMinimum(0.9);
643 TH1 * hproftheta = (TH1*) (pid->GetHistRatioQmax()->Projection(0,3)->ProfileX());
644 hproftheta->SetXTitle("tan(#theta)");
645 hproftheta->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
646 hproftheta->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
647 hproftheta->SetMarkerStyle(kmimarkers[0]);
648 hproftheta->SetMaximum(1.1);
649 hproftheta->SetMinimum(0.9);
653 TH1 * hprofdedx = (TH1*) (pid->GetHistRatioQmax()->Projection(0,4)->ProfileX());
654 hprofdedx->SetXTitle("dEdx_{track}");
655 hprofdedx->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
656 hprofdedx->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
657 hprofdedx->SetMarkerStyle(kmimarkers[0]);
658 hprofdedx->SetMaximum(1.1);
659 hprofdedx->SetMinimum(0.9);
663 TH1 * hprofdedxL = (TH1*) (pid->GetHistRatioQmax()->Projection(0,5)->ProfileX());
664 hprofdedxL->SetXTitle("dEdx_{track}#Delta_{x}");
665 hprofdedxL->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
666 hprofdedxL->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
667 hprofdedxL->SetMarkerStyle(kmimarkers[0]);
668 hprofdedxL->SetMaximum(1.1);
669 hprofdedxL->SetMinimum(0.9);
673 canvas->SaveAs(Form("%s/QmaxRatioType%d.eps",outputpath,ipad));
674 canvas->SaveAs(Form("%s/QmaxRatioType%d.gif",outputpath,ipad));
679 void AliTPCcalibPID::DrawRatiodEdx(Float_t demin, Float_t demax, const char* outputpath){
684 // Int_t kmicolors[10]={1,2,3,4,6,7,8,9,10,11};
685 Int_t kmimarkers[10]={21,22,23,24,25,26,27,28,29,30};
686 AliTPCcalibPID * pid = this;
687 TCanvas *canvas= new TCanvas("QRatiodEdx","QRatiodEdx",600,800);
689 canvas->SetLogx(kTRUE);
691 for (Int_t ipad=0;ipad<4;ipad++){
692 pid->GetHistRatioQmax()->GetAxis(7)->SetRange(ipad+2,ipad+2);
693 pid->GetHistRatioQtot()->GetAxis(7)->SetRange(ipad+2,ipad+2);
694 pid->GetHistRatioQmax()->GetAxis(5)->SetRangeUser(demin,demax);
695 pid->GetHistRatioQtot()->GetAxis(5)->SetRangeUser(demin,demax);
697 canvas->cd(ipad*2+1)->SetLogx(kFALSE);
698 hprofP = (TH1*) (pid->GetHistRatioQmax()->Projection(0,5)->ProfileX());
699 hprofP->SetXTitle("dEdx_{track}");
700 hprofP->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
701 hprofP->SetTitle(Form("Q_{max} dEdx_{tracklet}/dEdx_{track} type %d",ipad));
702 hprofP->SetMarkerStyle(kmimarkers[0]);
703 hprofP->SetMaximum(1.1);
704 hprofP->SetMinimum(0.9);
707 canvas->cd(ipad*2+2)->SetLogx(kFALSE);
708 hprofP = (TH1*) (pid->GetHistRatioQtot()->Projection(0,5)->ProfileX());
709 hprofP->SetXTitle("dEdx_{track}");
710 hprofP->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
711 hprofP->SetTitle(Form("Q_{tot} dEdx_{tracklet}/dEdx_{track} type %d",ipad));
712 hprofP->SetMarkerStyle(kmimarkers[0]);
713 hprofP->SetMaximum(1.1);
714 hprofP->SetMinimum(0.9);
719 canvas->SaveAs(Form("%s/QratiodEdx.eps",outputpath));
720 canvas->SaveAs(Form("%s/QratiodEdx.gif",outputpath));
725 void AliTPCcalibPID::DrawResolBGQtot(Int_t minClusters, Int_t maxClusters, Float_t minp, Float_t maxp, const char *outputpath, Bool_t resol){
729 AliTPCcalibPID * pid = this;
730 Int_t kmimarkers[10]={21,22,23,24,25,26,27,28,29,30};
738 pid->GetHistQtot()->GetAxis(6)->SetRangeUser(minClusters,maxClusters);
739 pid->GetHistQtot()->GetAxis(5)->SetRangeUser(1,10000);
740 pid->GetHistQtot()->GetAxis(4)->SetRangeUser(minp,maxp);
741 TCanvas *canvas= new TCanvas("dEdxResolQ_{Tot}","dEdxResolQ_{Tot}",800,600);
746 for (Int_t ipad=0;ipad<5;ipad++){
747 canvas->cd(1+ipad)->SetLogx(kTRUE);
748 if (ipad<4) pid->GetHistQtot()->GetAxis(7)->SetRange(ipad+2,ipad+2);
749 if (ipad==4) pid->GetHistQtot()->GetAxis(7)->SetRange(1,1);
750 his = (TH2*)(pid->GetHistQtot()->Projection(0,5));
751 his->FitSlicesY(0,0,-1,0,"QNR",&arr);
753 hmean = (TH1*)arr.At(1);
754 hsigma = (TH1*)arr.At(2)->Clone();
755 hsigma->SetMaximum(0.11);
756 hsigma->SetMinimum(0.4);
757 hsigma->Divide(hmean);
759 hsigma = (TH1*)arr.At(1)->Clone();
760 hsigma->SetMaximum(2);
761 hsigma->SetMinimum(0.5);
767 hsigma->SetXTitle("#beta#gamma");
768 hsigma->SetYTitle("#sigma_{dEdx}/dEdx");
769 hsigma->SetTitle(Form("#sigma_{dEdx}/dEdx_{tot} Pad %d",ipad));
770 hsigma->SetName(Form("#sigma_{dEdx}/dEdx_{tot} Pad %d",ipad));
772 hsigma->SetTitle(Form("#sigma_{dEdx}/dEdx_{tot} Full"));
773 hsigma->SetName(Form("#sigma_{dEdx}/dEdx_{tot} Full"));
775 hsigma->SetMarkerStyle(kmimarkers[0]);
779 canvas->SaveAs(Form("%s/dEdxResolTot.eps",outputpath));
780 canvas->SaveAs(Form("%s/dEdxResolTot.gif",outputpath));
782 canvas->SaveAs(Form("%s/dEdxBGTot.eps",outputpath));
783 canvas->SaveAs(Form("%s/dEdxBGTot.gif",outputpath));
787 void AliTPCcalibPID::DrawResolBGQmax(Int_t minClusters, Int_t maxClusters, Float_t minp, Float_t maxp, const char *outputpath, Bool_t resol){
789 // Int_t minClusters=140; Int_t maxClusters=200; const char *outputpath="picPID06"
791 Int_t kmimarkers[10]={21,22,23,24,25,26,27,28,29,30};
792 AliTPCcalibPID * pid = this;
799 pid->GetHistQmax()->GetAxis(6)->SetRangeUser(minClusters,maxClusters);
800 pid->GetHistQmax()->GetAxis(4)->SetRangeUser(minp,maxp);
801 pid->GetHistQmax()->GetAxis(5)->SetRangeUser(1,10000);
802 TCanvas *canvas= new TCanvas("dEdxResolQ_{Max}","dEdxResolQ_{Max}",800,600);
807 for (Int_t ipad=0;ipad<5;ipad++){
808 canvas->cd(1+ipad)->SetLogx(kTRUE);
809 if (ipad<4) pid->GetHistQmax()->GetAxis(7)->SetRange(ipad+2,ipad+2);
810 if (ipad==4) pid->GetHistQmax()->GetAxis(7)->SetRange(1,1);
811 his = (TH2*)(pid->GetHistQmax()->Projection(0,5));
812 his->FitSlicesY(0,0,-1,0,"QNR",&arr);
814 hmean = (TH1*)arr.At(1);
815 hsigma = (TH1*)arr.At(2)->Clone();
816 hsigma->SetMaximum(0.11);
817 hsigma->SetMinimum(0.4);
818 hsigma->Divide(hmean);
820 hsigma = (TH1*)arr.At(1)->Clone();
821 hsigma->SetMaximum(2.);
822 hsigma->SetMinimum(0.5);
827 hsigma->SetXTitle("#beta#gamma");
828 hsigma->SetYTitle("#sigma_{dEdx}/dEdx");
829 hsigma->SetTitle(Form("#sigma_{dEdx}/dEdx_{max} Pad %d",ipad));
830 hsigma->SetName(Form("#sigma_{dEdx}/dEdx_{max} Pad %d",ipad));
832 hsigma->SetTitle(Form("#sigma_{dEdx}/dEdx_{max} Full"));
833 hsigma->SetName(Form("#sigma_{dEdx}/dEdx_{max} Full"));
835 hsigma->SetMarkerStyle(kmimarkers[0]);
840 canvas->SaveAs(Form("%s/dEdxResolMax.eps",outputpath));
841 canvas->SaveAs(Form("%s/dEdxResolMax.gif",outputpath));
843 canvas->SaveAs(Form("%s/dEdxBGMax.eps",outputpath));
844 canvas->SaveAs(Form("%s/dEdxBGMax.gif",outputpath));
852 void AliTPCcalibPID::DumpTree(THnSparse * hndim, const char * outname){
855 // the tree will be written to the file - outname
860 TTreeSRedirector *pcstream = new TTreeSRedirector(outname);
863 Double_t position[10];
865 Int_t *bins = new Int_t[10];
868 const Float_t rmsy0=0.5, rmsz0=1;
870 Int_t ndim = hndim->GetNdimensions();
872 AliTPCParam * param = AliTPCcalibDB::Instance()->GetParameters();
873 for (Long64_t i = 0; i < hndim->GetNbins(); ++i) {
874 value = hndim->GetBinContent(i, bins);
875 for (Int_t idim = 0; idim < ndim; idim++) {
876 position[idim] = hndim->GetAxis(idim)->GetBinCenter(bins[idim]);
880 if (TMath::Abs(position[7]-1.5)<0.1) sector=0; // inner sector
881 if (TMath::Abs(position[7]-3.5)<0.1) row=96; // long pads
883 Double_t ty = TMath::Tan(TMath::ASin(position[2]));
884 Double_t tz = position[3]*TMath::Sqrt(1+ty*ty);
885 Double_t padLength= param->GetPadPitchLength(sector,row);
886 Double_t padWidth = param->GetPadPitchWidth(sector);
887 Double_t zwidth = param->GetZWidth();
888 Double_t zlength=250-TMath::Abs(position[1])*250.;
889 // diffusion in pad, time bin units
890 Double_t diffT=TMath::Sqrt(zlength)*param->GetDiffT()/padWidth;
891 Double_t diffL=TMath::Sqrt(zlength)*param->GetDiffL()/zwidth;
893 // transform angular effect to pad units
894 Double_t pky = ty*padLength/padWidth;
895 Double_t pkz = tz*padLength/zwidth;
898 Double_t sy = TMath::Sqrt(rmsy0*rmsy0+diffT*diffT+pky*pky/12.);
899 Double_t sz = TMath::Sqrt(rmsz0*rmsz0+diffL*diffL+pkz*pkz/12.);
900 Int_t ipad = TMath::Nint(position[7]-1.5);
902 (*pcstream)<<"Dump"<<
903 "bincont="<<value<< // bin content
904 "val="<<position[0]<< // parameter difference
905 "ipad="<<ipad<< // pad type
906 "dr="<<position[1]<< //drift
909 "p2="<<position[2]<< //p2
910 "p3="<<position[3]<< //p3
911 "p="<<position[4]<< //p
912 "bg="<<position[5]<< //bg
913 "ncl="<<position[6]<< //ncl
914 "type="<<position[7]<< //tracklet
915 "tot="<<position[8]<< //is tot
916 "sy="<<sy<< // sigma y
917 "sz="<<sz<< // sigma z
918 "pky="<<pky<< // ky - angle in pad units
919 "pkz="<<pkz<< // kz - angle in pad units
928 void AliTPCcalibPID::DumpTrees(){
930 // dump the content of histogram to the tree
932 AliTPCcalibPID *pid = this;
933 DumpTree(pid->GetHistQtot(),"dumpQtot.root");
934 DumpTree(pid->GetHistQmax(),"dumpQmax.root");
935 DumpTree(pid->GetHistRatioQtot(),"dumpRatioQtot.root");
936 DumpTree(pid->GetHistRatioQmax(),"dumpRatioQmax.root");