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]    = {150, 10,  10,    10,   50, 50,  8,    5};
 175   Double_t xminQA[8] = {0.,  0,    0,     0, 0.01, 0.1, 60,   0};
 176   Double_t xmaxQA[8] = {10., 1,  0.6,   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, 10,  10,    10,   25,  25,      4,    5};
 182   Double_t xminRa[9] = {0.5, 0,    0,     0,  0.2, 0.2,     60,    0};
 183   Double_t xmaxRa[9] = {4.0, 1,  0.6,   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   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);
 196
 197
 198   BinLogX(fDeDxQmax,4); BinLogX(fDeDxQmax,5);
 199   BinLogX(fDeDxQtot,4); BinLogX(fDeDxQtot,5);
 200   //
 201   BinLogX(fDeDxRatioMaxTot,4); BinLogX(fDeDxRatioMaxTot,5);
 202   BinLogX(fDeDxRatioQmax,4); BinLogX(fDeDxRatioQmax,5);
 203   BinLogX(fDeDxRatioQtot,4); BinLogX(fDeDxRatioQtot,5);
 204   BinLogX(fDeDxRatioTruncQmax,4); BinLogX(fDeDxRatioTruncQmax,5);
 205   BinLogX(fDeDxRatioTruncQtot,4); BinLogX(fDeDxRatioTruncQtot,5);
 206   //
 207   fHistNTracks = new TH1F("ntracks","Number of Tracks per Event; number of tracks per event; number of tracks",1001,-0.5,1000.5);
 208   fClusters = new TH1F("signal","Number of Clusters per track; number of clusters per track n_{cl}; counts",40,0,160);
 209   fPileUp = new TH2F("PileUp","timing plots.; #Delta L ; dEdx signal ",400,-1,1,400,0,200);
 210   fLandau = new TH2F("Landau","Landau.; pad type ; cluster charge ",3,-0.5,2.5,400,0,1000);
 211   //
 212   AliInfo("Non Default Constructor");
 213   //
 214 }
 215
 216
 217 AliTPCcalibPID::~AliTPCcalibPID(){
 218   //
 219   //
 220   //
 221
 222 }
 223
 224
 225
 226 void AliTPCcalibPID::Process(AliESDEvent *event) {
 227   //
 228   //
 229   //
 230   if (!event) {
 231     Printf("ERROR: ESD not available");
 232     return;
 233   }
 234   const Int_t kMinClustersTracklet = 25;
 235   Int_t ntracks=event->GetNumberOfTracks();
 236   fHistNTracks->Fill(ntracks);
 237
 238   AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
 239   if (!ESDfriend) {
 240    Printf("ERROR: ESDfriend not available");
 241    return;
 242   }
 243   //
 244   // track loop
 245   //
 246   for (Int_t i=0;i<ntracks;++i) {
 247
 248     AliESDtrack *track = event->GetTrack(i);
 249     if (!track) continue;
 250
 251
 252     const AliExternalTrackParam * trackIn = track->GetInnerParam();
 253     const AliExternalTrackParam * trackOut = track->GetOuterParam();
 254     if (!trackIn) continue;
 255     if (!trackOut) continue;
 256
 257     // calculate necessary track parameters
 258     //Double_t meanP = 0.5*(trackIn->GetP() + trackOut->GetP());
 259     Double_t meanP = trackIn->GetP();
 260     //Double_t d = trackIn->GetLinearD(0,0);
 261     Int_t NclsDeDx = track->GetTPCNcls();
 262
 263     //if (meanP > 0.7 || meanP < 0.2) continue;
 264      if (NclsDeDx < 60) continue;
 265
 266     // exclude tracks which do not look like primaries or are simply too short or on wrong sectors
 267
 268     //if (TMath::Abs(trackIn->GetSnp()) > 3*0.4) continue;
 269     //if (TMath::Abs(trackIn->GetZ()) > 150) continue;
 270     //if (seed->CookShape(1) > 1) continue;
 271     //if (TMath::Abs(trackIn->GetY()) > 20) continue;
 272     //if (TMath::Abs(d)>20) continue;   // distance to the 0,0; select only tracks which cross chambers under proper angle
 273     //if (TMath::Abs(trackIn->GetSnp()) > 0.6) continue;
 274     //if (TMath::Abs(trackOut->GetSnp()) > 0.2) continue;
 275     if (TMath::Abs(trackIn->GetAlpha()+0.872665)<0.01 || TMath::Abs(trackOut->GetAlpha()+0.872665)<0.01) continue;  // Funny sector !
 276
 277     // Get seeds
 278     AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
 279     if (!friendTrack) continue;
 280     TObject *calibObject;
 281     AliTPCseed *seed = 0;
 282     for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
 283       if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
 284     }
 285
 286     if (seed) {
 287       if (meanP > 30 && TMath::Abs(trackIn->GetSnp()) < 0.2) fClusters->Fill(track->GetTPCNcls());
 288       // calculate dEdx
 289       // (Double_t low, Double_t up, Int_t type, Int_t i1, Int_t i2, Bool_t shapeNorm,Bool_t posNorm, Int_t padNorm, Int_t returnVal)
 290       //Double_t TPCsignalTot    = (1/fMIP)*seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,0,0,159,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
 291       Double_t TPCsignalMax    = (1/fMIP)*seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,0,159,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
 292       //
 293       //Double_t TPCsignalShort  = seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,0,63,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
 294       //Double_t TPCsignalMedium = seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,63,63+64,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
 295       //Double_t TPCsignalLong   = seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,63+64,159,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
 296       //
 297       //Double_t driftMismatch = seed->GetDriftTimeMismatch(0,0.95,0.5);
 298       Double_t driftMismatch = 0;
 299       //      Double_t drift = 1 - (TMath::Abs(trackIn->GetZ()) + TMath::Abs(trackOut->GetZ()))/500.;
 300
 301       // particle identification
 302       Double_t mass = 0.105658369;// muon
 303
 304       if (meanP > 30) {
 305         fPileUp->Fill(driftMismatch,TPCsignalMax);
 306         //
 307         fLandau->Fill(0.1,0.6);
 308       }
 309       //var. of interest, z,sin(phi),tan(theta),p,betaGamma,ncls
 310       Double_t snpIn = TMath::Abs(trackIn->GetSnp());
 311       Double_t snpOut = TMath::Abs(trackOut->GetSnp());
 312       Double_t tglIn = TMath::Abs(trackIn->GetTgl());
 313       Double_t tglOut = TMath::Abs(trackOut->GetTgl());
 314       Double_t driftIn = 1 - (TMath::Abs(trackIn->GetZ()))/250.;
 315       Double_t driftOut = 1 - (TMath::Abs(trackIn->GetZ()))/250.;
 316       //
 317       // dEdx in region
 318       //
 319       Float_t dEdxTot[5],dEdxTotFull[5];
 320       Float_t dEdxMax[5],dEdxMaxFull[5];
 321       Int_t   ncl[5];
 322       for (Int_t itype=0; itype<5;itype++){
 323         Int_t row0=0, row1 =159;
 324         if (itype==1) {row0=0;      row1 = 63;};
 325         if (itype==2) {row0= 64;    row1=63+64;}
 326         if (itype==3) {row0= 63+64; row1=159;}
 327         dEdxTot[itype]= (1/fMIP)*seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,0,row0,row1,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
 328         dEdxMax[itype]= (1/fMIP)*seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,row0,row1,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
 329         // non trucated dEdx
 330         dEdxTotFull[itype]= (1/fMIP)*seed->CookdEdxNorm(0.0,0.99,0,row0,row1,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
 331         dEdxMaxFull[itype]= (1/fMIP)*seed->CookdEdxNorm(0.0,0.99,1,row0,row1,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
 332
 333
 334         ncl[itype]=TMath::Nint(seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,row0,row1,fUseShapeNorm,fUsePosNorm,fUsePadNorm,2));
 335       }
 336       //
 337       //
 338       //
 339       Float_t wmeanTot=0, wmeanMax=0, sumW=0;
 340       Double_t length[3] = {0.75,1,1.5};
 341
 342       for (Int_t ipad=0; ipad<3; ipad++){
 343         if (ncl[1+ipad]<3) continue;
 344         Double_t weight = Double_t(ncl[1+ipad])*TMath::Sqrt(length[ipad]);
 345         wmeanTot+=dEdxTot[1+ipad]*weight;
 346         wmeanMax+=dEdxMax[1+ipad]*weight;
 347         sumW+=weight;
 348       }
 349       if (sumW>0){
 350         dEdxTot[4]= wmeanTot/sumW;
 351         dEdxMax[4]= wmeanMax/sumW;
 352       }
 353       for (Int_t itype=0;itype<5;itype++){
 354         //
 355         Float_t snp=(TMath::Abs(snpIn)+TMath::Abs(snpOut))*0.5;
 356         Float_t tgl=(TMath::Abs(tglIn)+TMath::Abs(tglOut))*0.5;
 357         Float_t drift = (driftIn+driftOut)*0.5;
 358         if (itype==1) {snp = TMath::Abs(snpIn); tgl = TMath::Abs(tglIn); drift= driftIn;};
 359         if (itype==3) {snp = TMath::Abs(snpOut); tgl = TMath::Abs(tglOut);drift=driftOut;};
 360         if (ncl[itype]<kMinClustersTracklet) continue;
 361         Float_t deltaL = TMath::Sqrt(1+snp*snp+tgl*tgl);
 362         //
 363         Double_t vecQmax[8] = {dEdxMax[itype],drift,snp,tgl,meanP,meanP/mass,NclsDeDx, itype};
 364         Double_t vecQtot[8] = {dEdxTot[itype],drift,snp,tgl,meanP,meanP/mass,NclsDeDx, itype};
 365         //
 366         //
 367         //
 368         Double_t ratioMaxTot           = (dEdxTot[itype]>0)  ? dEdxMax[itype]/dEdxTot[itype]:0;
 369         Double_t ratioTrackletTot      = (dEdxTot[0]>0)      ? dEdxTot[itype]/dEdxTot[0]:0;
 370         Double_t ratioTrackletMax      = (dEdxMax[0]>0)      ? dEdxMax[itype]/dEdxMax[0]:0;
 371         Double_t ratioTrackletTruncTot = (dEdxTot[0]>0)      ? dEdxTotFull[itype]/dEdxTot[itype]:0;
 372         Double_t ratioTrackletTruncMax = (dEdxMax[0]>0)      ? dEdxMaxFull[itype]/dEdxMax[itype]:0;
 373
 374         Double_t vecRatioMaxTot[8]      = {ratioMaxTot,      drift,snp,tgl,dEdxTot[0],  dEdxTot[0]*deltaL,NclsDeDx,itype};
 375         Double_t vecRatioTrackletTot[8] = {ratioTrackletTot, drift,snp,tgl,dEdxTot[0],  dEdxTot[0]*deltaL,NclsDeDx,itype};
 376         Double_t vecRatioTrackletMax[8] = {ratioTrackletMax, drift,snp,tgl,dEdxMax[0],  dEdxMax[0]*deltaL,NclsDeDx,itype};
 377         Double_t vecRatioTrackletTruncTot[8] = {ratioTrackletTruncTot, drift,snp,tgl,dEdxTot[0],  dEdxTot[0]*deltaL,NclsDeDx,itype};
 378         Double_t vecRatioTrackletTruncMax[8] = {ratioTrackletTruncMax, drift,snp,tgl,dEdxMax[0],  dEdxMax[0]*deltaL,NclsDeDx,itype};
 379         fDeDxQmax->Fill(vecQmax);
 380         fDeDxQtot->Fill(vecQtot);
 381         fDeDxRatioMaxTot->Fill(vecRatioMaxTot);
 382         fDeDxRatioQmax->Fill(vecRatioTrackletTot);
 383         fDeDxRatioQtot->Fill(vecRatioTrackletMax);
 384         fDeDxRatioTruncQmax->Fill(vecRatioTrackletTruncTot);
 385         fDeDxRatioTruncQtot->Fill(vecRatioTrackletTruncMax);
 386         //
 387         TTreeSRedirector * cstream =  GetDebugStreamer();
 388         if (cstream){
 389           TVectorD vQT(9,vecQtot);
 390           TVectorD vQM(9,vecQmax);
 391           TVectorD vQMT(9, vecRatioMaxTot);
 392           TVectorD vQRT(9,vecRatioTrackletTot);
 393           TVectorD vQRM(9,vecRatioTrackletMax);
 394           TVectorD vQRTT(9,vecRatioTrackletTruncTot);
 395           TVectorD vQRTM(9,vecRatioTrackletTruncMax);
 396           (*cstream) << "dEdx" <<
 397             "run="<<fRun<<              //  run number
 398             "event="<<fEvent<<          //  event number
 399             "time="<<fTime<<            //  time stamp of event
 400             "trigger="<<fTrigger<<      //  trigger
 401             "mag="<<fMagF<<             //  magnetic field
 402             "track.="<<seed<<           //  original seed
 403             //
 404             "vQT.="<<&vQT<<             // trunc mean - total charge
 405             "vQM.="<<&vQM<<             // trunc mean - max charge
 406             //
 407             "vQMT.="<<&vQMT<<           // ratio max/tot
 408             "vQRT.="<<&vQRT<<           // ratio tracklet/track - total charge
 409             "vQRM.="<<&vQRM<<           // ratio tracklet/track - max charge
 410             "vQRTT.="<<&vQRTT<<         // ratio trunc/full     - total charge
 411             "vQRTM.="<<&vQRTM<<         // ratio trunc/full     - max charge
 412             "\n";
 413         }
 414       }
 415     }
 416   }
 417 }
 418
 419
 420
 421 void AliTPCcalibPID::Analyze() {
 422
 423
 424   return;
 425
 426 }
 427
 428
 429 Long64_t AliTPCcalibPID::Merge(TCollection *li) {
 430
 431   TIterator* iter = li->MakeIterator();
 432   AliTPCcalibPID* cal = 0;
 433
 434   while ((cal = (AliTPCcalibPID*)iter->Next())) {
 435     if (!cal->InheritsFrom(AliTPCcalibPID::Class())) {
 436       Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
 437       return -1;
 438     }
 439
 440     if (cal->GetHistNTracks()) fHistNTracks->Add(cal->GetHistNTracks());
 441     if (cal->GetHistClusters()) fClusters->Add(cal->GetHistClusters());
 442     if (cal->GetHistPileUp()) fPileUp->Add(cal->GetHistPileUp());
 443     if (cal->GetHistLandau()) fLandau->Add(cal->GetHistLandau());
 444     //
 445     if (cal->GetHistQmax()) fDeDxQmax->Add(cal->GetHistQmax());
 446     if (cal->GetHistQtot()) fDeDxQtot->Add(cal->GetHistQtot());
 447     if (cal->GetHistRatioMaxTot()) fDeDxRatioMaxTot->Add(cal->GetHistRatioMaxTot());
 448     if (cal->GetHistRatioQmax()) fDeDxRatioQmax->Add(cal->GetHistRatioQmax());
 449     if (cal->GetHistRatioQtot()) fDeDxRatioQtot->Add(cal->GetHistRatioQtot());
 450     if (cal->GetHistRatioTruncQmax()) fDeDxRatioTruncQmax->Add(cal->GetHistRatioTruncQmax());
 451     if (cal->GetHistRatioTruncQtot()) fDeDxRatioTruncQtot->Add(cal->GetHistRatioTruncQtot());
 452   }
 453
 454   return 0;
 455
 456 }
 457
 458
 459
 460 void AliTPCcalibPID::BinLogX(THnSparse *h, Int_t axisDim) {
 461
 462   // Method for the correct logarithmic binning of histograms
 463
 464   TAxis *axis = h->GetAxis(axisDim);
 465   int bins = axis->GetNbins();
 466
 467   Double_t from = axis->GetXmin();
 468   Double_t to = axis->GetXmax();
 469   Double_t *new_bins = new Double_t[bins + 1];
 470
 471   new_bins[0] = from;
 472   Double_t factor = pow(to/from, 1./bins);
 473
 474   for (int i = 1; i <= bins; i++) {
 475    new_bins[i] = factor * new_bins[i-1];
 476   }
 477   axis->Set(bins, new_bins);
 478   delete new_bins;
 479
 480 }
 481
 482
 483
 484
 485
 486 void AliTPCcalibPID::MakeReport(const char *outputpath) {
 487   //
 488   // Make a standard QA plots
 489   //
 490   for (Int_t ipad=0;ipad<4;ipad++){
 491     DrawRatioTot(ipad,outputpath);
 492     DrawRatioMax(ipad,outputpath);
 493   }
 494   DrawRatiodEdx(0.5,3,outputpath);
 495   DrawResolBGQtot(140,160,outputpath);
 496   DrawResolBGQmax(140,160,outputpath);
 497   return;
 498 }
 499
 500 void  AliTPCcalibPID::DrawRatioTot(Int_t ipad, const char* outputpath){
 501   //
 502   // Draw default ratio histogram for given pad type
 503   // ipad - 0 - short pads
 504   //        1 - medium pads
 505   //        2 - long pads
 506   //
 507   //  Int_t kmicolors[10]={1,2,3,4,6,7,8,9,10,11};
 508   Int_t kmimarkers[10]={21,22,23,24,25,26,27,28,29,30};
 509   AliTPCcalibPID * pid = this;
 510   TCanvas *canvas= new TCanvas(Form("QtotRatio%d)",ipad),Form("QtotRatio%d)",ipad));
 511   canvas->Divide(3,2);
 512   pid->GetHistRatioQtot()->GetAxis(7)->SetRange(ipad+2,ipad+2);
 513   canvas->cd(1);
 514   TH1 * his0 = pid->GetHistRatioQtot()->Projection(0);
 515   his0->SetXTitle("dEdx_{tracklet}/dEdx_{track}");
 516   his0->SetYTitle("");
 517   his0->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
 518   his0->Draw();
 519   //
 520   canvas->cd(2);
 521   TH1 * hprofz = (TH1*) (pid->GetHistRatioQtot()->Projection(0,1)->ProfileX());
 522   hprofz->SetXTitle("drift length");
 523   hprofz->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
 524   hprofz->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
 525   hprofz->SetMarkerStyle(kmimarkers[0]);
 526   hprofz->SetMaximum(1.1);
 527   hprofz->SetMinimum(0.9);
 528   hprofz->Draw();
 529   //
 530   canvas->cd(3);
 531   TH1 * hprofphi = (TH1*) (pid->GetHistRatioQtot()->Projection(0,2)->ProfileX());
 532   hprofphi->SetXTitle("sin(#phi)");
 533   hprofphi->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
 534   hprofphi->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
 535   hprofphi->SetMarkerStyle(kmimarkers[0]);
 536   hprofz->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
 537   hprofphi->SetMaximum(1.1);
 538   hprofphi->SetMinimum(0.9);
 539   hprofphi->Draw();
 540   //
 541   canvas->cd(4);
 542   TH1 * hproftheta = (TH1*) (pid->GetHistRatioQtot()->Projection(0,3)->ProfileX());
 543   hproftheta->SetXTitle("tan(#theta)");
 544   hproftheta->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
 545   hproftheta->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
 546   hproftheta->SetMarkerStyle(kmimarkers[0]);
 547   hproftheta->SetMaximum(1.1);
 548   hproftheta->SetMinimum(0.9);
 549   hproftheta->Draw();
 550   //
 551   canvas->cd(5);
 552   TH1 * hprofdedx = (TH1*) (pid->GetHistRatioQtot()->Projection(0,4)->ProfileX());
 553   hprofdedx->SetXTitle("dEdx_{track}");
 554   hprofdedx->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
 555   hprofdedx->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
 556   hprofdedx->SetMarkerStyle(kmimarkers[0]);
 557   hprofdedx->SetMaximum(1.1);
 558   hprofdedx->SetMinimum(0.9);
 559   hprofdedx->Draw();
 560
 561   canvas->cd(6);
 562   TH1 * hprofdedxL = (TH1*) (pid->GetHistRatioQtot()->Projection(0,5)->ProfileX());
 563   hprofdedxL->SetXTitle("dEdx_{track}#Delta_{x}");
 564   hprofdedxL->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
 565   hprofdedxL->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
 566   hprofdedxL->SetMarkerStyle(kmimarkers[0]);
 567   hprofdedxL->SetMaximum(1.1);
 568   hprofdedxL->SetMinimum(0.9);
 569   hprofdedxL->Draw();
 570
 571
 572
 573   canvas->SaveAs(Form("%s/QtotRatioType%d.eps",outputpath,ipad));
 574   canvas->SaveAs(Form("%s/QtotRatioType%d.gif",outputpath,ipad));
 575 }
 576
 577 void  AliTPCcalibPID::DrawRatioMax(Int_t ipad, const char* outputpath){
 578   //
 579   // Draw default ration histogram for given pad type
 580   // ipad - 0 - short pads
 581   //        1 - medium pads
 582   //        2 - long pads
 583   //
 584   //  Int_t kmicolors[10]={1,2,3,4,6,7,8,9,10,11};
 585   Int_t kmimarkers[10]={21,22,23,24,25,26,27,28,29,30};
 586   AliTPCcalibPID * pid = this;
 587   TCanvas *canvas= new TCanvas(Form("QmaxRatio%d)",ipad),Form("QmaxRatio%d)",ipad));
 588   canvas->Divide(3,2);
 589   pid->GetHistRatioQmax()->GetAxis(7)->SetRange(ipad+2,ipad+2);
 590   canvas->cd(1);
 591   TH1 * his0 = pid->GetHistRatioQmax()->Projection(0);
 592   his0->SetXTitle("dEdx_{tracklet}/dEdx_{track}");
 593   his0->SetYTitle("");
 594   his0->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
 595   his0->Draw();
 596   //
 597   canvas->cd(2);
 598   TH1 * hprofz = (TH1*) (pid->GetHistRatioQmax()->Projection(0,1)->ProfileX());
 599   hprofz->SetXTitle("drift length");
 600   hprofz->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
 601   hprofz->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
 602   hprofz->SetMarkerStyle(kmimarkers[0]);
 603   hprofz->SetMaximum(1.1);
 604   hprofz->SetMinimum(0.9);
 605   hprofz->Draw();
 606   //
 607   canvas->cd(3);
 608   TH1 * hprofphi = (TH1*) (pid->GetHistRatioQmax()->Projection(0,2)->ProfileX());
 609   hprofphi->SetXTitle("sin(#phi)");
 610   hprofphi->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
 611   hprofphi->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
 612   hprofphi->SetMarkerStyle(kmimarkers[0]);
 613   hprofphi->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
 614   hprofphi->SetMaximum(1.1);
 615   hprofphi->SetMinimum(0.9);
 616   hprofphi->Draw();
 617   //
 618   canvas->cd(4);
 619   TH1 * hproftheta = (TH1*) (pid->GetHistRatioQmax()->Projection(0,3)->ProfileX());
 620   hproftheta->SetXTitle("tan(#theta)");
 621   hproftheta->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
 622   hproftheta->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
 623   hproftheta->SetMarkerStyle(kmimarkers[0]);
 624   hproftheta->SetMaximum(1.1);
 625   hproftheta->SetMinimum(0.9);
 626   hproftheta->Draw();
 627
 628   canvas->cd(5);
 629   TH1 * hprofdedx = (TH1*) (pid->GetHistRatioQmax()->Projection(0,4)->ProfileX());
 630   hprofdedx->SetXTitle("dEdx_{track}");
 631   hprofdedx->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
 632   hprofdedx->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
 633   hprofdedx->SetMarkerStyle(kmimarkers[0]);
 634   hprofdedx->SetMaximum(1.1);
 635   hprofdedx->SetMinimum(0.9);
 636   hprofdedx->Draw();
 637
 638   canvas->cd(6);
 639   TH1 * hprofdedxL = (TH1*) (pid->GetHistRatioQmax()->Projection(0,5)->ProfileX());
 640   hprofdedxL->SetXTitle("dEdx_{track}#Delta_{x}");
 641   hprofdedxL->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
 642   hprofdedxL->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
 643   hprofdedxL->SetMarkerStyle(kmimarkers[0]);
 644   hprofdedxL->SetMaximum(1.1);
 645   hprofdedxL->SetMinimum(0.9);
 646   hprofdedxL->Draw();
 647
 648
 649   canvas->SaveAs(Form("%s/QmaxRatioType%d.eps",outputpath,ipad));
 650   canvas->SaveAs(Form("%s/QmaxRatioType%d.gif",outputpath,ipad));
 651 }
 652
 653
 654
 655 void AliTPCcalibPID::DrawRatiodEdx(Float_t demin, Float_t demax, const char* outputpath){
 656   //
 657   //
 658   //
 659   //
 660   //  Int_t kmicolors[10]={1,2,3,4,6,7,8,9,10,11};
 661   Int_t kmimarkers[10]={21,22,23,24,25,26,27,28,29,30};
 662   AliTPCcalibPID * pid = this;
 663   TCanvas *canvas= new TCanvas("QRatiodEdx","QRatiodEdx");
 664   canvas->Divide(2,4);
 665   canvas->SetLogx(kTRUE);
 666   TH1 * hprofP=0;
 667   for (Int_t ipad=0;ipad<4;ipad++){
 668     pid->GetHistRatioQmax()->GetAxis(7)->SetRange(ipad+2,ipad+2);
 669     pid->GetHistRatioQtot()->GetAxis(7)->SetRange(ipad+2,ipad+2);
 670     pid->GetHistRatioQmax()->GetAxis(5)->SetRangeUser(demin,demax);
 671     pid->GetHistRatioQtot()->GetAxis(5)->SetRangeUser(demin,demax);
 672
 673     canvas->cd(ipad*2+1)->SetLogx(kFALSE);
 674     hprofP  = (TH1*) (pid->GetHistRatioQmax()->Projection(0,5)->ProfileX());
 675     hprofP->SetXTitle("dEdx_{track}");
 676     hprofP->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
 677     hprofP->SetTitle(Form("Q_{max} dEdx_{tracklet}/dEdx_{track} type %d",0));
 678     hprofP->SetMarkerStyle(kmimarkers[0]);
 679     hprofP->SetMaximum(1.1);
 680     hprofP->SetMinimum(0.9);
 681     hprofP->Draw();
 682     // pad  Tot
 683     canvas->cd(ipad*2+2)->SetLogx(kFALSE);
 684     hprofP  = (TH1*) (pid->GetHistRatioQtot()->Projection(0,5)->ProfileX());
 685     hprofP->SetXTitle("dEdx_{track}");
 686     hprofP->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
 687     hprofP->SetTitle(Form("Q_{tot} dEdx_{tracklet}/dEdx_{track} type %d",0));
 688     hprofP->SetMarkerStyle(kmimarkers[0]);
 689     hprofP->SetMaximum(1.1);
 690     hprofP->SetMinimum(0.9);
 691     hprofP->Draw();
 692   }
 693   //
 694   //
 695   canvas->SaveAs(Form("%s/QratiodEdx%d.eps",outputpath));
 696   canvas->SaveAs(Form("%s/QratiodEdx%d.gif",outputpath));
 697 }
 698
 699
 700
 701 void AliTPCcalibPID::DrawResolBGQtot(Int_t minClusters, Int_t maxClusters, const char *outputpath){
 702   //
 703   //
 704   //
 705   AliTPCcalibPID * pid = this;
 706   Int_t kmimarkers[10]={21,22,23,24,25,26,27,28,29,30};
 707
 708   TObjArray arr;
 709   TH2 * his   =0;
 710   TH1 * hmean =0;
 711   TH1 * hsigma=0;
 712   //
 713   // set cut
 714   pid->GetHistQtot()->GetAxis(6)->SetRangeUser(minClusters,maxClusters);
 715   pid->GetHistQtot()->GetAxis(5)->SetRangeUser(1,10000);
 716   TCanvas *canvas= new TCanvas("dEdxResolQ_{Tot}","dEdxResolQ_{Tot}");
 717   canvas->Divide(2,3);
 718   //
 719   //
 720   //
 721   for (Int_t ipad=0;ipad<5;ipad++){
 722     canvas->cd(1+ipad)->SetLogx(kTRUE);
 723     if (ipad<4) pid->GetHistQtot()->GetAxis(7)->SetRange(ipad+2,ipad+2);
 724     if (ipad==4) pid->GetHistQtot()->GetAxis(7)->SetRange(1,1);
 725     his = (TH2*)(pid->GetHistQtot()->Projection(0,5));
 726     his->FitSlicesY(0,0,-1,0,"QNR",&arr);
 727     hmean = (TH1*)arr.At(1);
 728     hsigma = (TH1*)arr.At(2)->Clone();
 729     hsigma->Divide(hmean);
 730     arr.SetOwner(kTRUE);
 731     arr.Delete();
 732     delete his;
 733
 734     hsigma->SetXTitle("#beta#gamma");
 735     hsigma->SetYTitle("#sigma_{dEdx}/dEdx");
 736     hsigma->SetTitle(Form("#sigma_{dEdx}/dEdx_{tot} Pad %d",ipad));
 737     hsigma->SetName(Form("#sigma_{dEdx}/dEdx_{tot} Pad %d",ipad));
 738     if (ipad==4) {
 739       hsigma->SetTitle(Form("#sigma_{dEdx}/dEdx_{tot} Full"));
 740       hsigma->SetName(Form("#sigma_{dEdx}/dEdx_{tot} Full"));
 741     }
 742     hsigma->SetMarkerStyle(kmimarkers[0]);
 743     hsigma->SetMaximum(0.15);
 744     hsigma->SetMinimum(0.0);
 745     hsigma->Draw();
 746   }
 747   canvas->SaveAs(Form("%s/dEdxResolMax.eps",outputpath));
 748   canvas->SaveAs(Form("%s/dEdxResolMax.gif",outputpath));
 749 }
 750
 751 void AliTPCcalibPID::DrawResolBGQmax(Int_t minClusters, Int_t maxClusters, const char *outputpath){
 752   //
 753   // Int_t minClusters=140;  Int_t maxClusters=200; const char *outputpath="picPID06"
 754   //
 755   Int_t kmimarkers[10]={21,22,23,24,25,26,27,28,29,30};
 756   AliTPCcalibPID * pid = this;
 757   TObjArray arr;
 758   TH2 * his   =0;
 759   TH1 * hmean =0;
 760   TH1 * hsigma=0;
 761   //
 762   // set cut
 763   pid->GetHistQmax()->GetAxis(6)->SetRangeUser(minClusters,maxClusters);
 764   pid->GetHistQmax()->GetAxis(5)->SetRangeUser(1,10000);
 765   TCanvas *canvas= new TCanvas("dEdxResolQ_{Max}","dEdxResolQ_{Max}");
 766   canvas->Divide(2,3);
 767   //
 768   //
 769   //
 770   for (Int_t ipad=0;ipad<5;ipad++){
 771     canvas->cd(1+ipad)->SetLogx(kTRUE);
 772     if (ipad<4) pid->GetHistQmax()->GetAxis(7)->SetRange(ipad+2,ipad+2);
 773     if (ipad==4) pid->GetHistQmax()->GetAxis(7)->SetRange(1,1);
 774     his = (TH2*)(pid->GetHistQmax()->Projection(0,5));
 775     his->FitSlicesY(0,0,-1,0,"QNR",&arr);
 776     hmean = (TH1*)arr.At(1);
 777     hsigma = (TH1*)arr.At(2)->Clone();
 778     hsigma->Divide(hmean);
 779     arr.SetOwner(kTRUE);
 780     arr.Delete();
 781     delete his;
 782     hsigma->SetXTitle("#beta#gamma");
 783     hsigma->SetYTitle("#sigma_{dEdx}/dEdx");
 784     hsigma->SetTitle(Form("#sigma_{dEdx}/dEdx_{max} Pad %d",ipad));
 785     hsigma->SetName(Form("#sigma_{dEdx}/dEdx_{max} Pad %d",ipad));
 786     if (ipad==4){
 787       hsigma->SetTitle(Form("#sigma_{dEdx}/dEdx_{max} Full"));
 788       hsigma->SetName(Form("#sigma_{dEdx}/dEdx_{max} Full"));
 789     }
 790     hsigma->SetMarkerStyle(kmimarkers[0]);
 791     hsigma->SetMaximum(0.15);
 792     hsigma->SetMinimum(0.0);
 793     hsigma->Draw();
 794   }
 795   canvas->SaveAs(Form("%s/dEdxResolMax.eps",outputpath));
 796   canvas->SaveAs(Form("%s/dEdxResolMax.gif",outputpath));
 797 }
 798
 799
 800
 801 void AliTPCcalibPID::DumpTree(THnSparse * hndim, const char * outname){
 802   //
 803   // make a tree
 804   // the tree will be written to the file - outname
 805   //
 806   TTreeSRedirector *pcstream = new TTreeSRedirector(outname);
 807   //
 808   //
 809   Double_t position[10];
 810   Double_t value;
 811   Int_t *bins = new Int_t[10];
 812   //
 813   //
 814   Int_t ndim = hndim->GetNdimensions();
 815   //
 816   for (Long64_t i = 0; i < hndim->GetNbins(); ++i) {
 817     value = hndim->GetBinContent(i, bins);
 818     for (Int_t idim = 0; idim < ndim; idim++) {
 819       position[idim] = hndim->GetAxis(idim)->GetBinCenter(bins[idim]);
 820     }
 821     Double_t ty = TMath::Tan(TMath::ASin(position[2]));
 822     Double_t tz = position[3]*TMath::Sqrt(1+ty*ty);
 823     //
 824     (*pcstream)<<"Dump"<<
 825       "bincont="<<value<<      // bin content
 826       "val="<<position[0]<<    // parameter difference
 827       "dr="<<position[1]<<     //drift
 828       "ty="<<ty<<              //phi
 829       "tz="<<tz<<              //theta
 830       "p2="<<position[2]<<      //p2
 831       "p3="<<position[3]<<     //p3
 832       "p="<<position[4]<<      //p
 833       "bg="<<position[5]<<     //bg
 834       "ncl="<<position[6]<<    //ncl
 835       "type="<<position[7]<<   //tracklet
 836       "tot="<<position[8]<<    //is tot
 837       "\n";
 838   }
 839   delete pcstream;
 840 }
 841
 842
 843 void AliTPCcalibPID::DumpTrees(){
 844   //
 845   // dump the content of histogram to the tree
 846   //
 847   AliTPCcalibPID *pid = this;
 848   DumpTree(pid->GetHistQtot(),"dumpQtot.root");
 849   DumpTree(pid->GetHistQmax(),"dumpQmax.root");
 850   DumpTree(pid->GetHistRatioQtot(),"dumpRatioQtot.root");
 851   DumpTree(pid->GetHistRatioQmax(),"dumpRatioQmax.root");
 852 }
 853
 854
 855