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->GetHistLongMediumRatio()->Projection(0)->Fit("gaus")
  29 cal->GetHistShortMediumRatio()->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    fDeDxRatio(0),
 124    fDeDxShortMediumRatio(0),
 125    fDeDxLongMediumRatio(0)
 126 {
 127   //
 128   // Empty default cosntructor
 129   //
 130   AliInfo("Default Constructor");
 131 }
 132
 133
 134 AliTPCcalibPID::AliTPCcalibPID(const Text_t *name, const Text_t *title)
 135   :AliTPCcalibBase(),
 136    fMIP(0),
 137    fLowerTrunc(0),
 138    fUpperTrunc(0),
 139    fUseShapeNorm(0),
 140    fUsePosNorm(0),
 141    fUsePadNorm(0),
 142    fIsCosmic(0),
 143    fHistNTracks(0),
 144    fClusters(0),
 145    fPileUp(0),
 146    fLandau(0),
 147    fDeDxQmax(0),
 148    fDeDxQtot(0),
 149    fDeDxRatio(0),
 150    fDeDxShortMediumRatio(0),
 151    fDeDxLongMediumRatio(0)
 152 {
 153   //
 154   //
 155   //
 156   SetName(name);
 157   SetTitle(title);
 158   //
 159   fMIP = 40.;
 160   fLowerTrunc = 0;
 161   fUpperTrunc = 0.6;
 162   //
 163   fUseShapeNorm = kTRUE;
 164   fUsePosNorm = kFALSE;
 165   fUsePadNorm = 1;
 166   //
 167   fIsCosmic  = kTRUE;
 168   //
 169   //                  dE/dx,  z, phi, theta,    p,  bg, ncls
 170   Int_t binsQA[7]    = {150, 10,  10,    10,   50,   1,  8};
 171   Double_t xminQA[7] = {0.,  0,    0,     0, 0.01, 0.1, 60};
 172   Double_t xmaxQA[7] = {10., 1,  0.6,   1.5,   50, 500, 160};
 173   //
 174   Double_t xminRa[7] = {0.5, 0,    0,     0, 0.01, 0.1, 60};
 175   Double_t xmaxRa[7] = { 2., 1,  0.6,   1.5,   50, 500, 160};
 176
 177   // z,sin(phi),tan(theta),p,betaGamma,ncls
 178   fDeDxQmax  = new THnSparseS("fDeDxQmax","Qmax;(z,sin(phi),tan(theta),p,betaGamma,ncls); TPC signal Qmax (a.u.)",7,binsQA,xminQA,xmaxQA);
 179   fDeDxQtot  = new THnSparseS("fDeDxQtot","Qmax;(z,sin(phi),tan(theta),p,betaGamma,ncls); TPC signal Qmax (a.u.)",7,binsQA,xminQA,xmaxQA);
 180   fDeDxRatio = new THnSparseS("fDeDxRatio","Qmax;(z,sin(phi),tan(theta),p,betaGamma,ncls); TPC signal Qmax (a.u.)",7,binsQA,xminRa,xmaxRa);
 181   fDeDxShortMediumRatio = new THnSparseS("fDeDxQmax","Qmax;(z,sin(phi),tan(theta),p,betaGamma,ncls); TPC signal Qmax (a.u.)",7,binsQA,xminRa,xmaxRa);
 182   fDeDxLongMediumRatio  = new THnSparseS("fDeDxQmax","Qmax;(z,sin(phi),tan(theta),p,betaGamma,ncls); TPC signal Qmax (a.u.)",7,binsQA,xminRa,xmaxRa);
 183   BinLogX(fDeDxQmax,4); BinLogX(fDeDxQmax,5);
 184   BinLogX(fDeDxQtot,4); BinLogX(fDeDxQtot,5);
 185   BinLogX(fDeDxRatio,4); BinLogX(fDeDxRatio,5);
 186   BinLogX(fDeDxShortMediumRatio,4); BinLogX(fDeDxShortMediumRatio,5);
 187   BinLogX(fDeDxLongMediumRatio,4); BinLogX(fDeDxLongMediumRatio,5);
 188   //
 189   fHistNTracks = new TH1F("ntracks","Number of Tracks per Event; number of tracks per event; number of tracks",1001,-0.5,1000.5);
 190   fClusters = new TH1F("signal","Number of Clusters per track; number of clusters per track n_{cl}; counts",40,0,160);
 191   fPileUp = new TH2F("PileUp","timing plots.; #Delta L ; dEdx signal ",400,-1,1,400,0,200);
 192   fLandau = new TH2F("Landau","Landau.; pad type ; cluster charge ",3,-0.5,2.5,400,0,1000);
 193   //
 194   AliInfo("Non Default Constructor");
 195   //
 196 }
 197
 198
 199 AliTPCcalibPID::~AliTPCcalibPID(){
 200   //
 201   //
 202   //
 203 }
 204
 205
 206
 207 void AliTPCcalibPID::Process(AliESDEvent *event) {
 208   //
 209   //
 210   //
 211   if (!event) {
 212     Printf("ERROR: ESD not available");
 213     return;
 214   }
 215
 216   Int_t ntracks=event->GetNumberOfTracks();
 217   fHistNTracks->Fill(ntracks);
 218
 219   AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
 220   if (!ESDfriend) {
 221    Printf("ERROR: ESDfriend not available");
 222    return;
 223   }
 224   //
 225   // track loop
 226   //
 227   for (Int_t i=0;i<ntracks;++i) {
 228
 229     AliESDtrack *track = event->GetTrack(i);
 230     if (!track) continue;
 231
 232
 233     const AliExternalTrackParam * trackIn = track->GetInnerParam();
 234     const AliExternalTrackParam * trackOut = track->GetOuterParam();
 235     if (!trackIn) continue;
 236     if (!trackOut) continue;
 237
 238     // calculate necessary track parameters
 239     //Double_t meanP = 0.5*(trackIn->GetP() + trackOut->GetP());
 240     Double_t meanP = trackIn->GetP();
 241     Double_t d = trackIn->GetLinearD(0,0);
 242     Int_t NclsDeDx = track->GetTPCNcls();
 243
 244     //if (meanP > 0.7 || meanP < 0.2) continue;
 245      if (NclsDeDx < 60) continue;
 246
 247     // exclude tracks which do not look like primaries or are simply too short or on wrong sectors
 248
 249     //if (TMath::Abs(trackIn->GetSnp()) > 3*0.4) continue;
 250     //if (TMath::Abs(trackIn->GetZ()) > 150) continue;
 251     //if (seed->CookShape(1) > 1) continue;
 252     //if (TMath::Abs(trackIn->GetY()) > 20) continue;
 253     //if (TMath::Abs(d)>20) continue;   // distance to the 0,0; select only tracks which cross chambers under proper angle
 254     //if (TMath::Abs(trackIn->GetSnp()) > 0.6) continue;
 255     //if (TMath::Abs(trackOut->GetSnp()) > 0.2) continue;
 256     if (TMath::Abs(trackIn->GetAlpha()+0.872665)<0.01 || TMath::Abs(trackOut->GetAlpha()+0.872665)<0.01) continue;  // Funny sector !
 257
 258     // Get seeds
 259     AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
 260     if (!friendTrack) continue;
 261     TObject *calibObject;
 262     AliTPCseed *seed = 0;
 263     for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
 264       if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
 265     }
 266
 267     if (seed) {
 268       if (meanP > 30 && TMath::Abs(trackIn->GetSnp()) < 0.2) fClusters->Fill(track->GetTPCNcls());
 269       // calculate dEdx
 270       // (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)
 271       Double_t TPCsignalTot    = (1/fMIP)*seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,0,0,159,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
 272       Double_t TPCsignalMax    = (1/fMIP)*seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,0,159,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
 273       //
 274       Double_t TPCsignalShort  = seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,0,63,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
 275       Double_t TPCsignalMedium = seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,63,63+64,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
 276       Double_t TPCsignalLong   = seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,63+64,159,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
 277       //
 278       //Double_t driftMismatch = seed->GetDriftTimeMismatch(0,0.95,0.5);
 279       Double_t driftMismatch = 0;
 280       Double_t drift = 1 - (TMath::Abs(trackIn->GetZ()) + TMath::Abs(trackOut->GetZ()))/500.;
 281
 282       // particle identification
 283       Double_t mass = 0.105658369;// muon
 284
 285       if (meanP > 30) {
 286         fPileUp->Fill(driftMismatch,TPCsignalMax);
 287         //
 288         fLandau->Fill(0.1,0.6);
 289       }
 290       //var. of interest, z,sin(phi),tan(theta),p,betaGamma,ncls
 291       Double_t snp = TMath::Abs(trackIn->GetSnp());
 292       Double_t tgl = TMath::Abs(trackIn->GetTgl());
 293       //
 294       Double_t vecQmax[7] = {TPCsignalMax,drift,snp,tgl,meanP,meanP/mass,NclsDeDx};
 295       Double_t vecQtot[7] = {TPCsignalTot,drift,snp,tgl,meanP,meanP/mass,NclsDeDx};
 296       Double_t vecRatio[7] = {TPCsignalMax/TPCsignalTot,drift,snp,tgl,meanP,meanP/mass,NclsDeDx};
 297       Double_t vecShortMediumRatio[7] = {TPCsignalShort/TPCsignalMedium,drift,snp,tgl,meanP,meanP/mass,NclsDeDx};
 298       Double_t vecLongMediumRatio[7] = {TPCsignalLong/TPCsignalMedium,drift,snp,tgl,meanP,meanP/mass,NclsDeDx};
 299
 300       fDeDxQmax->Fill(vecQmax);
 301       fDeDxQtot->Fill(vecQtot);
 302       fDeDxRatio->Fill(vecRatio);
 303       fDeDxShortMediumRatio->Fill(vecShortMediumRatio);
 304       fDeDxLongMediumRatio->Fill(vecLongMediumRatio);
 305
 306     }
 307
 308   }
 309
 310
 311 }
 312
 313
 314
 315 void AliTPCcalibPID::Analyze() {
 316
 317
 318   return;
 319
 320 }
 321
 322
 323 Long64_t AliTPCcalibPID::Merge(TCollection *li) {
 324
 325   TIterator* iter = li->MakeIterator();
 326   AliTPCcalibPID* cal = 0;
 327
 328   while ((cal = (AliTPCcalibPID*)iter->Next())) {
 329     if (!cal->InheritsFrom(AliTPCcalibPID::Class())) {
 330       Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
 331       return -1;
 332     }
 333
 334     if (cal->GetHistNTracks()) fHistNTracks->Add(cal->GetHistNTracks());
 335     if (cal->GetHistClusters()) fClusters->Add(cal->GetHistClusters());
 336     if (cal->GetHistPileUp()) fPileUp->Add(cal->GetHistPileUp());
 337     if (cal->GetHistLandau()) fLandau->Add(cal->GetHistLandau());
 338     //
 339     if (cal->GetHistQmax()) fDeDxQmax->Add(cal->GetHistQmax());
 340     if (cal->GetHistQtot()) fDeDxQtot->Add(cal->GetHistQtot());
 341     if (cal->GetHistRatio()) fDeDxRatio->Add(cal->GetHistRatio());
 342     if (cal->GetHistShortMediumRatio()) fDeDxShortMediumRatio->Add(cal->GetHistShortMediumRatio());
 343     if (cal->GetHistLongMediumRatio()) fDeDxLongMediumRatio->Add(cal->GetHistLongMediumRatio());
 344   }
 345
 346   return 0;
 347
 348 }
 349
 350
 351
 352
 353 void AliTPCcalibPID::MakeReport() {
 354   //
 355   // 1. standard dEdx picture
 356   TCanvas *cDeDx = new TCanvas("cDeDx","cDeDx");
 357   GetHistQmax()->Projection(0,4)->Draw("colz");
 358   return;
 359 }
 360
 361
 362
 363 void AliTPCcalibPID::BinLogX(THnSparse *h, Int_t axisDim) {
 364
 365   // Method for the correct logarithmic binning of histograms
 366
 367   TAxis *axis = h->GetAxis(axisDim);
 368   int bins = axis->GetNbins();
 369
 370   Double_t from = axis->GetXmin();
 371   Double_t to = axis->GetXmax();
 372   Double_t *new_bins = new Double_t[bins + 1];
 373
 374   new_bins[0] = from;
 375   Double_t factor = pow(to/from, 1./bins);
 376
 377   for (int i = 1; i <= bins; i++) {
 378    new_bins[i] = factor * new_bins[i-1];
 379   }
 380   axis->Set(bins, new_bins);
 381   delete new_bins;
 382
 383 }
 384