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   if (esdFriend->TestSkipBit()) return;
 260   //
 261   // track loop
 262   //
 263   for (Int_t i=0;i<ntracks;++i) {
 264
 265     AliESDtrack *track = event->GetTrack(i);
 266     if (!track) continue;
 267
 268
 269     AliExternalTrackParam * trackIn  = (AliExternalTrackParam *)track->GetInnerParam();
 270     AliExternalTrackParam * trackOut = (AliExternalTrackParam *)track->GetOuterParam();
 271     if (!trackIn) continue;
 272     if (!trackOut) continue;
 273
 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();
 279
 280     //if (meanP > 0.7 || meanP < 0.2) continue;
 281      if (nclsDeDx < 60) continue;
 282
 283     // exclude tracks which do not look like primaries or are simply too short or on wrong sectors
 284
 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 !
 293
 294     // Get seeds
 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;
 301     }
 302
 303     if (seed) {
 304       if (meanP > 30 && TMath::Abs(trackIn->GetSnp()) < 0.2) fClusters->Fill(track->GetTPCNcls());
 305       // calculate dEdx
 306       Double_t tpcSignalMax    = (1/fMIP)*seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,0,159,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
 307       //
 308       //
 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.;
 312
 313       // particle identification
 314       Double_t mass = 0.105658369;// muon
 315
 316       if (meanP > 30) {
 317         fPileUp->Fill(driftMismatch,tpcSignalMax);
 318         //
 319         fLandau->Fill(0.1,0.6);
 320       }
 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.;
 328       //
 329       // dEdx in region
 330       //
 331       Float_t dEdxTot[5],dEdxTotFull[5];
 332       Float_t dEdxMax[5],dEdxMaxFull[5];
 333       Float_t   ncl[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;}
 339         if (fUsePosNorm==0){
 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);
 342           // non trucated dEdx
 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);
 345         }else{
 346           dEdxTot[itype]= (1/fMIP)*seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,0,row0,row1);
 347           dEdxMax[itype]= (1/fMIP)*seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,1,row0,row1);
 348           // non trucated dEdx
 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);
 351         }
 352
 353         ncl[itype]=(seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,1,row0,row1,2));
 354       }
 355       //
 356       //
 357       //
 358       Float_t wmeanTot=0, wmeanMax=0, sumW=0;
 359       Double_t length[3] = {0.75,1,1.5};
 360       //       //
 361
 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;
 367         sumW+=weight;
 368       }
 369       if (sumW>0){
 370         dEdxTot[4]= wmeanTot/sumW;
 371         dEdxMax[4]= wmeanMax/sumW;
 372       }
 373       for (Int_t itype=0;itype<5;itype++){
 374         //
 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);
 382         //
 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};
 385         //
 386         //
 387         //
 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;
 393
 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);
 406         //
 407         TTreeSRedirector * cstream =  GetDebugStreamer();
 408         if (cstream){
 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
 427             "itype="<<itype<<
 428             //
 429             "vQT.="<<&vQT<<             // trunc mean - total charge
 430             "vQM.="<<&vQM<<             // trunc mean - max charge
 431             //
 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
 437             "\n";
 438         }
 439       }
 440     }
 441   }
 442 }
 443
 444
 445
 446 void AliTPCcalibPID::Analyze() {
 447
 448
 449   return;
 450
 451 }
 452
 453
 454 Long64_t AliTPCcalibPID::Merge(TCollection *li) {
 455
 456   TIterator* iter = li->MakeIterator();
 457   AliTPCcalibPID* cal = 0;
 458
 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());
 462       return -1;
 463     }
 464
 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());
 469     //
 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());
 477   }
 478
 479   return 0;
 480
 481 }
 482
 483
 484
 485 void AliTPCcalibPID::BinLogX(THnSparse *h, Int_t axisDim) {
 486
 487   // Method for the correct logarithmic binning of histograms
 488
 489   TAxis *axis = h->GetAxis(axisDim);
 490   int bins = axis->GetNbins();
 491
 492   Double_t from = axis->GetXmin();
 493   Double_t to = axis->GetXmax();
 494   Double_t *newBins = new Double_t[bins + 1];
 495
 496   newBins[0] = from;
 497   Double_t factor = pow(to/from, 1./bins);
 498
 499   for (int i = 1; i <= bins; i++) {
 500    newBins[i] = factor * newBins[i-1];
 501   }
 502   axis->Set(bins, newBins);
 503   delete [] newBins;
 504
 505 }
 506
 507
 508
 509
 510
 511 void AliTPCcalibPID::MakeReport(const char *outputpath) {
 512   //
 513   // Make a standard QA plots
 514   //
 515   for (Int_t ipad=0;ipad<4;ipad++){
 516     DrawRatioTot(ipad,outputpath);
 517     DrawRatioMax(ipad,outputpath);
 518   }
 519   DrawRatiodEdx(0.5,3,outputpath);
 520   DrawResolBGQtot(140,160,1,40,outputpath);
 521   DrawResolBGQmax(140,160,1,40,outputpath);
 522   return;
 523 }
 524
 525 void  AliTPCcalibPID::DrawRatioTot(Int_t ipad, const char* outputpath){
 526   //
 527   // Draw default ratio histogram for given pad type
 528   // ipad - 0 - short pads
 529   //        1 - medium pads
 530   //        2 - long pads
 531   //
 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);
 536   canvas->Divide(3,2);
 537   pid->GetHistRatioQtot()->GetAxis(7)->SetRange(ipad+2,ipad+2);
 538   canvas->cd(1);
 539   TH1 * his0 = pid->GetHistRatioQtot()->Projection(0);
 540   his0->SetXTitle("dEdx_{tracklet}/dEdx_{track}");
 541   his0->SetYTitle("");
 542   his0->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
 543   his0->Draw();
 544   //
 545   canvas->cd(2);
 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);
 553   hprofz->Draw();
 554   //
 555   canvas->cd(3);
 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);
 564   hprofphi->Draw();
 565   //
 566   canvas->cd(4);
 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);
 574   hproftheta->Draw();
 575   //
 576   canvas->cd(5);
 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);
 584   hprofdedx->Draw();
 585
 586   canvas->cd(6);
 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);
 594   hprofdedxL->Draw();
 595
 596
 597
 598   canvas->SaveAs(Form("%s/QtotRatioType%d.eps",outputpath,ipad));
 599   canvas->SaveAs(Form("%s/QtotRatioType%d.gif",outputpath,ipad));
 600 }
 601
 602 void  AliTPCcalibPID::DrawRatioMax(Int_t ipad, const char* outputpath){
 603   //
 604   // Draw default ration histogram for given pad type
 605   // ipad - 0 - short pads
 606   //        1 - medium pads
 607   //        2 - long pads
 608   //
 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);
 613   canvas->Divide(3,2);
 614   pid->GetHistRatioQmax()->GetAxis(7)->SetRange(ipad+2,ipad+2);
 615   canvas->cd(1);
 616   TH1 * his0 = pid->GetHistRatioQmax()->Projection(0);
 617   his0->SetXTitle("dEdx_{tracklet}/dEdx_{track}");
 618   his0->SetYTitle("");
 619   his0->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
 620   his0->Draw();
 621   //
 622   canvas->cd(2);
 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);
 630   hprofz->Draw();
 631   //
 632   canvas->cd(3);
 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);
 641   hprofphi->Draw();
 642   //
 643   canvas->cd(4);
 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);
 651   hproftheta->Draw();
 652
 653   canvas->cd(5);
 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);
 661   hprofdedx->Draw();
 662
 663   canvas->cd(6);
 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);
 671   hprofdedxL->Draw();
 672
 673
 674   canvas->SaveAs(Form("%s/QmaxRatioType%d.eps",outputpath,ipad));
 675   canvas->SaveAs(Form("%s/QmaxRatioType%d.gif",outputpath,ipad));
 676 }
 677
 678
 679
 680 void AliTPCcalibPID::DrawRatiodEdx(Float_t demin, Float_t demax, const char* outputpath){
 681   //
 682   //
 683   //
 684   //
 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);
 689   canvas->Divide(2,4);
 690   canvas->SetLogx(kTRUE);
 691   TH1 * hprofP=0;
 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);
 697
 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);
 706     hprofP->Draw();
 707     // pad  Tot
 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);
 716     hprofP->Draw();
 717   }
 718   //
 719   //
 720   canvas->SaveAs(Form("%s/QratiodEdx.eps",outputpath));
 721   canvas->SaveAs(Form("%s/QratiodEdx.gif",outputpath));
 722 }
 723
 724
 725
 726 void AliTPCcalibPID::DrawResolBGQtot(Int_t minClusters, Int_t maxClusters, Float_t minp, Float_t maxp, const char *outputpath, Bool_t resol){
 727   //
 728   //
 729   //
 730   AliTPCcalibPID * pid = this;
 731   Int_t kmimarkers[10]={21,22,23,24,25,26,27,28,29,30};
 732
 733   TObjArray arr;
 734   TH2 * his   =0;
 735   TH1 * hmean =0;
 736   TH1 * hsigma=0;
 737   //
 738   // set cut
 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);
 743   canvas->Divide(2,3);
 744   //
 745   //
 746   //
 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);
 753     if (resol){
 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);
 759     }else{
 760       hsigma = (TH1*)arr.At(1)->Clone();
 761       hsigma->SetMaximum(2);
 762       hsigma->SetMinimum(0.5);
 763     }
 764     arr.SetOwner(kTRUE);
 765     arr.Delete();
 766     delete his;
 767
 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));
 772     if (ipad==4) {
 773       hsigma->SetTitle(Form("#sigma_{dEdx}/dEdx_{tot} Full"));
 774       hsigma->SetName(Form("#sigma_{dEdx}/dEdx_{tot} Full"));
 775     }
 776     hsigma->SetMarkerStyle(kmimarkers[0]);
 777     hsigma->Draw();
 778   }
 779   if (resol){
 780     canvas->SaveAs(Form("%s/dEdxResolTot.eps",outputpath));
 781     canvas->SaveAs(Form("%s/dEdxResolTot.gif",outputpath));
 782   }else {
 783     canvas->SaveAs(Form("%s/dEdxBGTot.eps",outputpath));
 784     canvas->SaveAs(Form("%s/dEdxBGTot.gif",outputpath));
 785   }
 786 }
 787
 788 void AliTPCcalibPID::DrawResolBGQmax(Int_t minClusters, Int_t maxClusters, Float_t minp, Float_t maxp, const char *outputpath, Bool_t resol){
 789   //
 790   // Int_t minClusters=140;  Int_t maxClusters=200; const char *outputpath="picPID06"
 791   //
 792   Int_t kmimarkers[10]={21,22,23,24,25,26,27,28,29,30};
 793   AliTPCcalibPID * pid = this;
 794   TObjArray arr;
 795   TH2 * his   =0;
 796   TH1 * hmean =0;
 797   TH1 * hsigma=0;
 798   //
 799   // set cut
 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);
 804   canvas->Divide(2,3);
 805   //
 806   //
 807   //
 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);
 814     if (resol){
 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);
 820     }else{
 821       hsigma = (TH1*)arr.At(1)->Clone();
 822       hsigma->SetMaximum(2.);
 823       hsigma->SetMinimum(0.5);
 824     }
 825     arr.SetOwner(kTRUE);
 826     arr.Delete();
 827     delete his;
 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));
 832     if (ipad==4){
 833       hsigma->SetTitle(Form("#sigma_{dEdx}/dEdx_{max} Full"));
 834       hsigma->SetName(Form("#sigma_{dEdx}/dEdx_{max} Full"));
 835     }
 836     hsigma->SetMarkerStyle(kmimarkers[0]);
 837     hsigma->Draw();
 838   }
 839
 840   if (resol){
 841     canvas->SaveAs(Form("%s/dEdxResolMax.eps",outputpath));
 842     canvas->SaveAs(Form("%s/dEdxResolMax.gif",outputpath));
 843   }else {
 844     canvas->SaveAs(Form("%s/dEdxBGMax.eps",outputpath));
 845     canvas->SaveAs(Form("%s/dEdxBGMax.gif",outputpath));
 846   }
 847
 848
 849 }
 850
 851
 852
 853 void AliTPCcalibPID::DumpTree(THnSparse * hndim, const char * outname){
 854   //
 855   // make a tree
 856   // the tree will be written to the file - outname
 857   //
 858   //
 859
 860
 861   TTreeSRedirector *pcstream = new TTreeSRedirector(outname);
 862   //
 863   //
 864   Double_t position[10];
 865   Double_t value;
 866   Int_t bins[1000];
 867   //
 868   //
 869   const Float_t rmsy0=0.5, rmsz0=1;
 870   //
 871   Int_t ndim = hndim->GetNdimensions();
 872   //
 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]);
 878     }
 879     Int_t sector=36;
 880     Int_t row=0;
 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
 883
 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;
 893     //
 894     // transform angular effect to pad units
 895     Double_t pky   = ty*padLength/padWidth;
 896     Double_t pkz   = tz*padLength/zwidth;
 897     //
 898     //
 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);
 902      //
 903     (*pcstream)<<"Dump"<<
 904       "bincont="<<value<<      // bin content
 905       "val="<<position[0]<<    // parameter difference
 906       "ipad="<<ipad<<          // pad type
 907       "dr="<<position[1]<<     //drift
 908       "ty="<<ty<<              //phi
 909       "tz="<<tz<<              //theta
 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
 921       "diy="<<diffT<<
 922       "diz="<<diffL<<
 923       "\n";
 924   }
 925   delete pcstream;
 926 }
 927
 928
 929 void AliTPCcalibPID::DumpTrees(){
 930   //
 931   // dump the content of histogram to the tree
 932   //
 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");
 938 }
 939
 940
 941