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 *new_bins = new Double_t[bins + 1];
 494
 495   new_bins[0] = from;
 496   Double_t factor = pow(to/from, 1./bins);
 497
 498   for (int i = 1; i <= bins; i++) {
 499    new_bins[i] = factor * new_bins[i-1];
 500   }
 501   axis->Set(bins, new_bins);
 502   delete new_bins;
 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