Fixes for coding violations
[u/mrichter/AliRoot.git] / TPC / AliTPCcalibPID.cxx
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
19 Basic calibration and QA class for the PID of the TPC based on dEdx.
20
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.
22
23 1.a) Fast equalization for different gain in pad regions:
24
25 TFile f("CalibObjects.root")
26 AliTPCcalibPID *cal = f->Get("TPCCalib")->FindObject("calibPID")
27
28 cal->GetHistRatioMaxTot()->Projection(0)->Fit("gaus")
29 cal->GetHistRatioTracklet()->Projection(0)->Fit("gaus")
30
31 1.b) Update OCDB:
32 .x $ALICE_ROOT/TPC/macros/ConfigOCDB.C
33 AliTPCClusterParam * parcl = AliTPCcalibDB::Instance()->GetClusterParam();
34 *parcl->fQpadMnorm)[ipad] = oldvalue*corrFactor
35
36 Int_t runNumber = 0;
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);
46   
47
48 2.) Checks should be done with particles on the Fermi-Plateau and whoch leave a long track:
49
50 cal->GetHistQmax()->GetAxis(6)->SetRangeUser(120,160)
51 cal->GetHistQmax()->GetAxis(4)->SetRangeUser(20,50)
52
53 variables to check:
54
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.
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
67 Double_t initialParam[] = {3.81470,15.2608,3.25306e-07,2.51791,2.71012}
68
69
70 Send 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
105 ClassImp(AliTPCcalibPID)
106
107
108 AliTPCcalibPID::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),
123    fDeDxRatioMaxTot(0),
124    fDeDxRatioQmax(0),
125    fDeDxRatioQtot(0),
126    fDeDxRatioTruncQtot(0),
127    fDeDxRatioTruncQmax(0)
128 {  
129   //
130   // Empty default cosntructor
131   //
132   AliInfo("Default Constructor");  
133 }
134
135
136 AliTPCcalibPID::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),
151    fDeDxRatioMaxTot(0),
152    fDeDxRatioQmax(0), 
153    fDeDxRatioQtot(0) ,
154    fDeDxRatioTruncQtot(0),
155    fDeDxRatioTruncQmax(0)
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;
168   fUsePosNorm = 0;
169   fUsePadNorm = 1;
170   //
171   fIsCosmic  = kTRUE;
172   //
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};
177   //
178   //
179   //
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};
184   
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);
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
196
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);
198
199
200   BinLogX(fDeDxQmax,4); BinLogX(fDeDxQmax,5);BinLogX(fDeDxQmax,0);
201   BinLogX(fDeDxQtot,4); BinLogX(fDeDxQtot,5);BinLogX(fDeDxQmax,0);
202   //
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);
208   //
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);
213   //
214   AliInfo("Non Default Constructor");  
215   //
216 }
217
218
219 AliTPCcalibPID::~AliTPCcalibPID(){
220   //
221   //
222   //
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
227   //
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
230   //
231   // ratio histograms
232   //
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
238 }
239
240
241
242 void AliTPCcalibPID::Process(AliESDEvent *event) {
243   //
244   //
245   //
246   if (!event) {
247     Printf("ERROR: ESD not available");
248     return;
249   }  
250   const Int_t kMinClustersTracklet = 25;
251   Int_t ntracks=event->GetNumberOfTracks(); 
252   fHistNTracks->Fill(ntracks);
253   
254   AliESDfriend *esdFriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
255   if (!esdFriend) {
256    Printf("ERROR: esdFriend not available");
257    return;
258   }  
259   //
260   // track loop
261   //
262   for (Int_t i=0;i<ntracks;++i) {
263
264     AliESDtrack *track = event->GetTrack(i);
265     if (!track) continue;
266     
267     
268     AliExternalTrackParam * trackIn  = (AliExternalTrackParam *)track->GetInnerParam();
269     AliExternalTrackParam * trackOut = (AliExternalTrackParam *)track->GetOuterParam();
270     if (!trackIn) continue;
271     if (!trackOut) continue;
272
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();
278
279     //if (meanP > 0.7 || meanP < 0.2) continue;
280      if (nclsDeDx < 60) continue;     
281
282     // exclude tracks which do not look like primaries or are simply too short or on wrong sectors
283
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 !
292     
293     // Get seeds
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;
300     }    
301
302     if (seed) {
303       if (meanP > 30 && TMath::Abs(trackIn->GetSnp()) < 0.2) fClusters->Fill(track->GetTPCNcls());
304       // calculate dEdx
305       Double_t tpcSignalMax    = (1/fMIP)*seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,0,159,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
306       //
307       //
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.;        
311       
312       // particle identification
313       Double_t mass = 0.105658369;// muon
314       
315       if (meanP > 30) {
316         fPileUp->Fill(driftMismatch,tpcSignalMax);
317         //
318         fLandau->Fill(0.1,0.6);
319       }
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.;
327       //
328       // dEdx in region
329       //
330       Float_t dEdxTot[5],dEdxTotFull[5];
331       Float_t dEdxMax[5],dEdxMaxFull[5];
332       Float_t   ncl[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;}
338         if (fUsePosNorm==0){
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);
341           // non trucated dEdx
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);
344         }else{
345           dEdxTot[itype]= (1/fMIP)*seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,0,row0,row1);
346           dEdxMax[itype]= (1/fMIP)*seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,1,row0,row1);
347           // non trucated dEdx
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);
350         }
351
352         ncl[itype]=(seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,1,row0,row1,2));
353       }
354       //
355       //
356       //
357       Float_t wmeanTot=0, wmeanMax=0, sumW=0;
358       Double_t length[3] = {0.75,1,1.5};
359       //       //
360       
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;
366         sumW+=weight;
367       }
368       if (sumW>0){
369         dEdxTot[4]= wmeanTot/sumW;
370         dEdxMax[4]= wmeanMax/sumW;      
371       }
372       for (Int_t itype=0;itype<5;itype++){
373         //
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);
381         //
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};
384         //
385         //
386         //
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;
392
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(vecRatioTrackletMax); 
402         fDeDxRatioQtot->Fill(vecRatioTrackletTot); 
403         fDeDxRatioTruncQmax->Fill(vecRatioTrackletTruncMax); 
404         fDeDxRatioTruncQtot->Fill(vecRatioTrackletTruncTot); 
405         //
406         TTreeSRedirector * cstream =  GetDebugStreamer();
407         if (cstream){
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
426             "itype="<<itype<<
427             //
428             "vQT.="<<&vQT<<             // trunc mean - total charge
429             "vQM.="<<&vQM<<             // trunc mean - max charge 
430             //
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
436             "\n";
437         }
438       }   
439     }    
440   }  
441 }
442
443
444
445 void AliTPCcalibPID::Analyze() {
446
447
448   return;
449
450 }
451
452
453 Long64_t AliTPCcalibPID::Merge(TCollection *li) {
454
455   TIterator* iter = li->MakeIterator();
456   AliTPCcalibPID* cal = 0;
457
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());
461       return -1;
462     }
463     
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());
468     //
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());
476   }
477   
478   return 0;
479   
480 }
481
482
483
484 void AliTPCcalibPID::BinLogX(THnSparse *h, Int_t axisDim) {
485
486   // Method for the correct logarithmic binning of histograms
487
488   TAxis *axis = h->GetAxis(axisDim);
489   int bins = axis->GetNbins();
490
491   Double_t from = axis->GetXmin();
492   Double_t to = axis->GetXmax();
493   Double_t *newBins = new Double_t[bins + 1];
494    
495   newBins[0] = from;
496   Double_t factor = pow(to/from, 1./bins);
497   
498   for (int i = 1; i <= bins; i++) {
499    newBins[i] = factor * newBins[i-1];
500   }
501   axis->Set(bins, newBins);
502   delete [] newBins;
503   
504 }
505
506
507
508
509
510 void AliTPCcalibPID::MakeReport(const char *outputpath) {
511   //
512   // Make a standard QA plots
513   //
514   for (Int_t ipad=0;ipad<4;ipad++){
515     DrawRatioTot(ipad,outputpath);
516     DrawRatioMax(ipad,outputpath);
517   }
518   DrawRatiodEdx(0.5,3,outputpath);
519   DrawResolBGQtot(140,160,1,40,outputpath);
520   DrawResolBGQmax(140,160,1,40,outputpath);
521   return;
522 }
523
524 void  AliTPCcalibPID::DrawRatioTot(Int_t ipad, const char* outputpath){
525   //
526   // Draw default ratio histogram for given pad type
527   // ipad - 0 - short pads
528   //        1 - medium pads
529   //        2 - long pads
530   //
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);
535   canvas->Divide(3,2);
536   pid->GetHistRatioQtot()->GetAxis(7)->SetRange(ipad+2,ipad+2);
537   canvas->cd(1);
538   TH1 * his0 = pid->GetHistRatioQtot()->Projection(0);
539   his0->SetXTitle("dEdx_{tracklet}/dEdx_{track}");
540   his0->SetYTitle("");  
541   his0->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
542   his0->Draw();
543   //
544   canvas->cd(2);
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);
552   hprofz->Draw();
553   //
554   canvas->cd(3);
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);
563   hprofphi->Draw();
564   //
565   canvas->cd(4);
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);
573   hproftheta->Draw();
574   //
575   canvas->cd(5);
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);
583   hprofdedx->Draw();
584
585   canvas->cd(6);
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);
593   hprofdedxL->Draw();
594
595
596
597   canvas->SaveAs(Form("%s/QtotRatioType%d.eps",outputpath,ipad));
598   canvas->SaveAs(Form("%s/QtotRatioType%d.gif",outputpath,ipad));
599 }
600
601 void  AliTPCcalibPID::DrawRatioMax(Int_t ipad, const char* outputpath){
602   //
603   // Draw default ration histogram for given pad type
604   // ipad - 0 - short pads
605   //        1 - medium pads
606   //        2 - long pads
607   //
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);
612   canvas->Divide(3,2);
613   pid->GetHistRatioQmax()->GetAxis(7)->SetRange(ipad+2,ipad+2);
614   canvas->cd(1);
615   TH1 * his0 = pid->GetHistRatioQmax()->Projection(0);
616   his0->SetXTitle("dEdx_{tracklet}/dEdx_{track}");
617   his0->SetYTitle("");  
618   his0->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
619   his0->Draw();
620   //
621   canvas->cd(2);
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);
629   hprofz->Draw();
630   //
631   canvas->cd(3);
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);
640   hprofphi->Draw();
641   //
642   canvas->cd(4);
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);
650   hproftheta->Draw();
651
652   canvas->cd(5);
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);
660   hprofdedx->Draw();
661
662   canvas->cd(6);
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);
670   hprofdedxL->Draw();
671
672
673   canvas->SaveAs(Form("%s/QmaxRatioType%d.eps",outputpath,ipad));
674   canvas->SaveAs(Form("%s/QmaxRatioType%d.gif",outputpath,ipad));
675 }
676
677
678
679 void AliTPCcalibPID::DrawRatiodEdx(Float_t demin, Float_t demax, const char* outputpath){
680   //
681   //
682   //
683   //
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);
688   canvas->Divide(2,4);
689   canvas->SetLogx(kTRUE);
690   TH1 * hprofP=0;
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);
696
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);
705     hprofP->Draw();
706     // pad  Tot
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);
715     hprofP->Draw();
716   }
717   //
718   //
719   canvas->SaveAs(Form("%s/QratiodEdx.eps",outputpath));
720   canvas->SaveAs(Form("%s/QratiodEdx.gif",outputpath));
721 }
722
723
724
725 void AliTPCcalibPID::DrawResolBGQtot(Int_t minClusters, Int_t maxClusters, Float_t minp, Float_t maxp, const char *outputpath, Bool_t resol){
726   //
727   // 
728   //
729   AliTPCcalibPID * pid = this;
730   Int_t kmimarkers[10]={21,22,23,24,25,26,27,28,29,30};
731
732   TObjArray arr;
733   TH2 * his   =0; 
734   TH1 * hmean =0;
735   TH1 * hsigma=0;
736   //
737   // set cut
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);
742   canvas->Divide(2,3);
743   //
744   //
745   //
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);
752     if (resol){
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);
758     }else{
759       hsigma = (TH1*)arr.At(1)->Clone();
760       hsigma->SetMaximum(2);
761       hsigma->SetMinimum(0.5);
762     }
763     arr.SetOwner(kTRUE);
764     arr.Delete();
765     delete his;
766
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));
771     if (ipad==4) {
772       hsigma->SetTitle(Form("#sigma_{dEdx}/dEdx_{tot} Full"));
773       hsigma->SetName(Form("#sigma_{dEdx}/dEdx_{tot} Full"));
774     }
775     hsigma->SetMarkerStyle(kmimarkers[0]);
776     hsigma->Draw();
777   }
778   if (resol){
779     canvas->SaveAs(Form("%s/dEdxResolTot.eps",outputpath));
780     canvas->SaveAs(Form("%s/dEdxResolTot.gif",outputpath));
781   }else {
782     canvas->SaveAs(Form("%s/dEdxBGTot.eps",outputpath));
783     canvas->SaveAs(Form("%s/dEdxBGTot.gif",outputpath));
784   }
785 }
786
787 void AliTPCcalibPID::DrawResolBGQmax(Int_t minClusters, Int_t maxClusters, Float_t minp, Float_t maxp, const char *outputpath, Bool_t resol){
788   //
789   // Int_t minClusters=140;  Int_t maxClusters=200; const char *outputpath="picPID06"
790   //
791   Int_t kmimarkers[10]={21,22,23,24,25,26,27,28,29,30};
792   AliTPCcalibPID * pid = this;
793   TObjArray arr;
794   TH2 * his   =0; 
795   TH1 * hmean =0;
796   TH1 * hsigma=0;
797   //
798   // set cut
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);
803   canvas->Divide(2,3);
804   //
805   //
806   //
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);
813     if (resol){
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);
819     }else{
820       hsigma = (TH1*)arr.At(1)->Clone();      
821       hsigma->SetMaximum(2.);
822       hsigma->SetMinimum(0.5);
823     }
824     arr.SetOwner(kTRUE);
825     arr.Delete();
826     delete his;
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));
831     if (ipad==4){
832       hsigma->SetTitle(Form("#sigma_{dEdx}/dEdx_{max} Full"));
833       hsigma->SetName(Form("#sigma_{dEdx}/dEdx_{max} Full"));
834     }
835     hsigma->SetMarkerStyle(kmimarkers[0]);
836     hsigma->Draw();
837   }
838
839   if (resol){
840     canvas->SaveAs(Form("%s/dEdxResolMax.eps",outputpath));
841     canvas->SaveAs(Form("%s/dEdxResolMax.gif",outputpath));
842   }else {
843     canvas->SaveAs(Form("%s/dEdxBGMax.eps",outputpath));
844     canvas->SaveAs(Form("%s/dEdxBGMax.gif",outputpath));
845   }
846
847
848 }
849
850
851
852 void AliTPCcalibPID::DumpTree(THnSparse * hndim, const char * outname){
853   //
854   // make a tree
855   // the tree will be written to the file - outname
856   //
857   //
858   
859
860   TTreeSRedirector *pcstream = new TTreeSRedirector(outname);
861   //
862   //
863   Double_t position[10];
864   Double_t value; 
865   Int_t *bins = new Int_t[10];
866   //
867   //
868   const Float_t rmsy0=0.5, rmsz0=1;
869   //
870   Int_t ndim = hndim->GetNdimensions();
871   //
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]);
877     }      
878     Int_t sector=36;
879     Int_t row=0;
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
882
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;
892     //
893     // transform angular effect to pad units
894     Double_t pky   = ty*padLength/padWidth;
895     Double_t pkz   = tz*padLength/zwidth;
896     //
897     //
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);
901      //
902     (*pcstream)<<"Dump"<<
903       "bincont="<<value<<      // bin content
904       "val="<<position[0]<<    // parameter difference
905       "ipad="<<ipad<<          // pad type
906       "dr="<<position[1]<<     //drift
907       "ty="<<ty<<              //phi
908       "tz="<<tz<<              //theta      
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
920       "diy="<<diffT<<
921       "diz="<<diffL<<
922       "\n";
923   }
924   delete pcstream;
925 }
926
927
928 void AliTPCcalibPID::DumpTrees(){
929   //
930   // dump the content of histogram to the tree
931   //
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");
937 }
938
939
940