TPC/AliTPCcalibTimeGain.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 //0.  Libraries to lod
  19 gSystem->Load("libANALYSIS");
  20 gSystem->Load("libSTAT");
  21 gSystem->Load("libTPCcalib");
  22
  23
  24 //1. Do calibration ...
  25 //
  26 // compare reference
  27
  28 //
  29 //2. Visualize results
  30 //
  31 TFile fcalib("CalibObjects.root");
  32 TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
  33 AliTPCcalibTimeGain * gain = ( AliTPCcalibTimeGain *)array->FindObject("calibTimeGain");
  34 TGraphErrors * gr = AliTPCcalibBase::FitSlices(gain->GetHistGainTime(),0,1,2000,10,0.2,0.8);
  35 gain->GetHistGainTime()->Projection(0,1)->Draw("colz")
  36 gr->SetMarkerStyle(25);
  37 gr->Draw("lp")
  38
  39
  40 //
  41 // MakeSlineFit example
  42 //
  43 AliSplineFit fit;
  44 fit.SetGraph(gr)
  45 fit->SetMinPoints(gr->GetN()+1);
  46 fit->InitKnots(gr,2,0,0.001)
  47 fit.SplineFit(0)
  48 fit.MakeDiffHisto(gr)->Draw();
  49 TGraph * grfit = fit.MakeGraph(gr->GetX()[0],gr->GetX()[gr->GetN()-1],50000,0);
  50
  51 gr->SetMarkerStyle(25);
  52 gr->Draw("alp");
  53 grfit->SetLineColor(2);
  54 grfit->Draw("lu");
  55
  56 */
  57
  58
  59 #include "Riostream.h"
  60 #include "TChain.h"
  61 #include "TTree.h"
  62 #include "TH1F.h"
  63 #include "TH2F.h"
  64 #include "TH3F.h"
  65 #include "THnSparse.h"
  66 #include "TGraphErrors.h"
  67 #include "TList.h"
  68 #include "TMath.h"
  69 #include "TCanvas.h"
  70 #include "TFile.h"
  71 #include "TF1.h"
  72 #include "TVectorD.h"
  73 #include "TProfile.h"
  74 #include "TGraphErrors.h"
  75 #include "TCanvas.h"
  76
  77 #include "AliTPCclusterMI.h"
  78 #include "AliTPCseed.h"
  79 #include "AliESDVertex.h"
  80 #include "AliESDEvent.h"
  81 #include "AliESDfriend.h"
  82 #include "AliESDInputHandler.h"
  83 #include "AliAnalysisManager.h"
  84
  85 #include "AliTracker.h"
  86 #include "AliMagF.h"
  87 #include "AliTPCCalROC.h"
  88
  89 #include "AliLog.h"
  90
  91 #include "AliTPCcalibTimeGain.h"
  92
  93 #include "TTreeStream.h"
  94 #include "AliTPCTracklet.h"
  95 #include "TTimeStamp.h"
  96 #include "AliTPCcalibDB.h"
  97 #include "AliTPCcalibLaser.h"
  98 #include "AliDCSSensorArray.h"
  99 #include "AliDCSSensor.h"
 100
 101 ClassImp(AliTPCcalibTimeGain)
 102
 103
 104 AliTPCcalibTimeGain::AliTPCcalibTimeGain()
 105   :AliTPCcalibBase(),
 106    fHistGainTime(0),
 107    fGainVsTime(0),
 108    fHistDeDxTotal(0),
 109    fIntegrationTimeDeDx(0),
 110    fMIP(0),
 111    fLowerTrunc(0),
 112    fUpperTrunc(0),
 113    fUseShapeNorm(0),
 114    fUsePosNorm(0),
 115    fUsePadNorm(0),
 116    fIsCosmic(0),
 117    fLowMemoryConsumption(0)
 118 {
 119   AliInfo("Default Constructor");
 120 }
 121
 122
 123 AliTPCcalibTimeGain::AliTPCcalibTimeGain(const Text_t *name, const Text_t *title, UInt_t StartTime, UInt_t EndTime, Int_t deltaIntegrationTimeGain)
 124   :AliTPCcalibBase(),
 125    fHistGainTime(0),
 126    fGainVsTime(0),
 127    fHistDeDxTotal(0),
 128    fIntegrationTimeDeDx(0),
 129    fMIP(0),
 130    fLowerTrunc(0),
 131    fUpperTrunc(0),
 132    fUseShapeNorm(0),
 133    fUsePosNorm(0),
 134    fUsePadNorm(0),
 135    fIsCosmic(0),
 136    fLowMemoryConsumption(0)
 137  {
 138
 139   SetName(name);
 140   SetTitle(title);
 141
 142   AliInfo("Non Default Constructor");
 143
 144   fIntegrationTimeDeDx = deltaIntegrationTimeGain;
 145
 146   Double_t deltaTime = EndTime - StartTime;
 147
 148
 149   // main histogram for time dependence: dE/dx, time, type (1-muon cosmic,2-pion beam data), meanDriftlength, momenta (only filled if enough space is available)
 150   Int_t timeBins = TMath::Nint(deltaTime/deltaIntegrationTimeGain);
 151   Int_t binsGainTime[5]    = {100,  timeBins,    2,  25, 200};
 152   Double_t xminGainTime[5] = {0.5, StartTime,  0.5,   0, 0.1};
 153   Double_t xmaxGainTime[5] = {  4,   EndTime,  2.5, 250, 50};
 154   fHistGainTime = new THnSparseF("HistGainTime","dEdx time dep.;dEdx;dEdx,time,type,driftlength,momenta",5,binsGainTime,xminGainTime,xmaxGainTime);
 155   BinLogX(fHistGainTime, 4);
 156   //
 157   fHistDeDxTotal = new TH2F("DeDx","dEdx; momentum p (GeV); TPC signal (a.u.)",250,0.01,100.,1000,0.,1000);
 158   BinLogX(fHistDeDxTotal);
 159
 160   // default values for dE/dx
 161   fMIP = 50.;
 162   fLowerTrunc = 0.0;
 163   fUpperTrunc = 0.7;
 164   fUseShapeNorm = kTRUE;
 165   fUsePosNorm = kFALSE;
 166   fUsePadNorm = kFALSE;
 167   //
 168   fIsCosmic = kTRUE;
 169   fLowMemoryConsumption = kTRUE;
 170   //
 171
 172  }
 173
 174
 175
 176 AliTPCcalibTimeGain::~AliTPCcalibTimeGain(){
 177   //
 178   //
 179   //
 180 }
 181
 182
 183 void AliTPCcalibTimeGain::Process(AliESDEvent *event) {
 184   //
 185   // main track loop
 186   //
 187   if (!event) {
 188     Printf("ERROR: ESD not available");
 189     return;
 190   }
 191
 192   if (fIsCosmic) { // this should be removed at some point based on trigger mask !?
 193     ProcessCosmicEvent(event);
 194   } else {
 195     ProcessBeamEvent(event);
 196   }
 197
 198
 199
 200
 201 }
 202
 203
 204 void AliTPCcalibTimeGain::ProcessCosmicEvent(AliESDEvent *event) {
 205
 206   AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
 207   if (!ESDfriend) {
 208    Printf("ERROR: ESDfriend not available");
 209    return;
 210   }
 211   //
 212   UInt_t time = event->GetTimeStamp();
 213   Int_t ntracks = event->GetNumberOfTracks();
 214   //
 215   // track loop
 216   //
 217   for (Int_t i=0;i<ntracks;++i) {
 218
 219     AliESDtrack *track = event->GetTrack(i);
 220     if (!track) continue;
 221
 222     const AliExternalTrackParam * trackIn = track->GetInnerParam();
 223     const AliExternalTrackParam * trackOut = track->GetOuterParam();
 224     if (!trackIn) continue;
 225     if (!trackOut) continue;
 226
 227     // calculate necessary track parameters
 228     Double_t meanP = trackIn->GetP();
 229     Double_t meanDrift = 250 - 0.5*TMath::Abs(trackIn->GetZ() + trackOut->GetZ());
 230     Int_t NclsDeDx = track->GetTPCNcls();
 231
 232     // exclude tracks which do not look like primaries or are simply too short or on wrong sectors
 233     if (NclsDeDx < 60) continue;
 234     if (TMath::Abs(trackIn->GetTgl()) > 1) continue;
 235     if (TMath::Abs(trackIn->GetSnp()) > 0.6) continue;
 236
 237     // Get seeds
 238     AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
 239     if (!friendTrack) continue;
 240     TObject *calibObject;
 241     AliTPCseed *seed = 0;
 242     for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
 243       if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
 244     }
 245
 246     if (seed) {
 247       Double_t TPCsignalMax = (1/fMIP)*seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,0,159,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
 248       fHistDeDxTotal->Fill(meanP, TPCsignalMax);
 249       //
 250       if (fLowMemoryConsumption) {
 251         if (meanP < 20) continue;
 252         meanP = 30; // set all momenta to one in order to save memory
 253       }
 254       //dE/dx, time, type (1-muon cosmic,2-pion beam data), momenta
 255       Double_t vec[5] = {TPCsignalMax,time,1,meanDrift,meanP};
 256       fHistGainTime->Fill(vec);
 257     }
 258
 259   }
 260
 261 }
 262
 263
 264
 265 void AliTPCcalibTimeGain::ProcessBeamEvent(AliESDEvent *event) {
 266
 267   AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
 268   if (!ESDfriend) {
 269    Printf("ERROR: ESDfriend not available");
 270    return;
 271   }
 272   //
 273   UInt_t time = event->GetTimeStamp();
 274   Int_t ntracks = event->GetNumberOfTracks();
 275   //
 276   // track loop
 277   //
 278   for (Int_t i=0;i<ntracks;++i) {
 279
 280     AliESDtrack *track = event->GetTrack(i);
 281     if (!track) continue;
 282
 283     const AliExternalTrackParam * trackIn = track->GetInnerParam();
 284     const AliExternalTrackParam * trackOut = track->GetOuterParam();
 285     if (!trackIn) continue;
 286     if (!trackOut) continue;
 287
 288     // calculate necessary track parameters
 289     Double_t meanP = trackIn->GetP();
 290     Double_t meanDrift = 250 - 0.5*TMath::Abs(trackIn->GetZ() + trackOut->GetZ());
 291     Int_t NclsDeDx = track->GetTPCNcls();
 292
 293     // exclude tracks which do not look like primaries or are simply too short or on wrong sectors
 294     if (NclsDeDx < 60) continue;
 295     if (TMath::Abs(trackIn->GetTgl()) > 1) continue;
 296     if (TMath::Abs(trackIn->GetSnp()) > 0.6) continue;
 297
 298     // Get seeds
 299     AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
 300     if (!friendTrack) continue;
 301     TObject *calibObject;
 302     AliTPCseed *seed = 0;
 303     for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
 304       if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
 305     }
 306
 307     if (seed) {
 308       Double_t TPCsignalMax = (1/fMIP)*seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,0,159,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
 309       fHistDeDxTotal->Fill(meanP, TPCsignalMax);
 310       //
 311       if (fLowMemoryConsumption) {
 312         if (meanP > 0.5 || meanP < 0.4) continue;
 313         meanP = 0.45; // set all momenta to one in order to save memory
 314       }
 315       //dE/dx, time, type (1-muon cosmic,2-pion beam data), momenta
 316       Double_t vec[5] = {TPCsignalMax,time,1,meanDrift,meanP};
 317       fHistGainTime->Fill(vec);
 318     }
 319
 320   }
 321
 322 }
 323
 324
 325 void AliTPCcalibTimeGain::Analyze() {
 326   //
 327   //
 328   //
 329   if (fIsCosmic) {
 330     fHistGainTime->GetAxis(2)->SetRangeUser(0.51,1.49); // only cosmics
 331     fHistGainTime->GetAxis(4)->SetRangeUser(20,100);    // only Fermi-Plateau muons
 332   } else {
 333     fHistGainTime->GetAxis(2)->SetRangeUser(1.51,2.49); // only beam data
 334     fHistGainTime->GetAxis(4)->SetRangeUser(0.39,0.51); // only MIP pions
 335   }
 336   //
 337   fGainVsTime = AliTPCcalibBase::FitSlices(fHistGainTime,0,1,2000,10);
 338   //
 339   return;
 340 }
 341
 342
 343 Long64_t AliTPCcalibTimeGain::Merge(TCollection *li) {
 344
 345   TIterator* iter = li->MakeIterator();
 346   AliTPCcalibTimeGain* cal = 0;
 347
 348   while ((cal = (AliTPCcalibTimeGain*)iter->Next())) {
 349     if (!cal->InheritsFrom(AliTPCcalibTimeGain::Class())) {
 350       Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
 351       return -1;
 352     }
 353
 354     // add histograms here...
 355     if (cal->GetHistGainTime()) fHistGainTime->Add(cal->GetHistGainTime());
 356     if (cal->GetHistDeDxTotal()) fHistDeDxTotal->Add(cal->GetHistDeDxTotal());
 357
 358   }
 359
 360   return 0;
 361
 362 }
 363
 364
 365
 366 void AliTPCcalibTimeGain::BinLogX(THnSparse *h, Int_t axisDim) {
 367
 368   // Method for the correct logarithmic binning of histograms
 369
 370   TAxis *axis = h->GetAxis(axisDim);
 371   int bins = axis->GetNbins();
 372
 373   Double_t from = axis->GetXmin();
 374   Double_t to = axis->GetXmax();
 375   Double_t *new_bins = new Double_t[bins + 1];
 376
 377   new_bins[0] = from;
 378   Double_t factor = pow(to/from, 1./bins);
 379
 380   for (int i = 1; i <= bins; i++) {
 381    new_bins[i] = factor * new_bins[i-1];
 382   }
 383   axis->Set(bins, new_bins);
 384   delete new_bins;
 385
 386 }
 387
 388
 389 void AliTPCcalibTimeGain::BinLogX(TH1 *h) {
 390
 391   // Method for the correct logarithmic binning of histograms
 392
 393   TAxis *axis = h->GetXaxis();
 394   int bins = axis->GetNbins();
 395
 396   Double_t from = axis->GetXmin();
 397   Double_t to = axis->GetXmax();
 398   Double_t *new_bins = new Double_t[bins + 1];
 399
 400   new_bins[0] = from;
 401   Double_t factor = pow(to/from, 1./bins);
 402
 403   for (int i = 1; i <= bins; i++) {
 404    new_bins[i] = factor * new_bins[i-1];
 405   }
 406   axis->Set(bins, new_bins);
 407   delete new_bins;
 408
 409 }
 410
 411
 412
 413 void AliTPCcalibTimeGain::CalculateBetheAlephParams(TH2F *hist, Double_t * ini) {
 414   //{0.0762*MIP,10.632,1.34e-05,1.863,1.948}
 415   const Double_t sigma = 0.06;
 416
 417   TH2F *histBG = new TH2F("histBG","dEdxBg; #beta #gamma; TPC signal (a.u.)",TMath::Nint(0.5*hist->GetNbinsX()),0.1,5000.,hist->GetNbinsY(),hist->GetYaxis()->GetXmin(),hist->GetYaxis()->GetXmax());
 418   BinLogX(histBG);
 419
 420   TF1 *foElectron = new TF1("foElectron", "AliExternalTrackParam::BetheBlochAleph(x/0.000511,[0],[1],[2],[3],[4])",0.01,100);
 421   TF1 *foMuon = new TF1("foMuon", "AliExternalTrackParam::BetheBlochAleph(x/0.1056,[0],[1],[2],[3],[4])",0.01,100);
 422   TF1 *foPion = new TF1("foPion", "AliExternalTrackParam::BetheBlochAleph(x/0.138,[0],[1],[2],[3],[4])",0.01,100);
 423   TF1 *foKaon = new TF1("foKaon", "AliExternalTrackParam::BetheBlochAleph(x/0.498,[0],[1],[2],[3],[4])",0.01,100);
 424   TF1 *foProton = new TF1("foProton", "AliExternalTrackParam::BetheBlochAleph(x/0.938,[0],[1],[2],[3],[4])",0.01,100);
 425   foElectron->SetParameters(ini);
 426   foMuon->SetParameters(ini);
 427   foPion->SetParameters(ini);
 428   foKaon->SetParameters(ini);
 429   foProton->SetParameters(ini);
 430
 431   TCanvas *CanvCheck1 = new TCanvas();
 432   hist->Draw("colz");
 433   foElectron->Draw("same");
 434   foMuon->Draw("same");
 435   foPion->Draw("same");
 436   foKaon->Draw("same");
 437   foProton->Draw("same");
 438
 439   // Loop over all points of the input histogram
 440
 441   for(Int_t i=1; i < hist->GetNbinsX(); i++) {
 442     Double_t x = hist->GetXaxis()->GetBinCenter(i);
 443     for(Int_t j=1; j < hist->GetNbinsY(); j++) {
 444       Long64_t n = hist->GetBin(i, j);
 445       Double_t y = hist->GetYaxis()->GetBinCenter(j);
 446       Double_t entries = hist->GetBinContent(n);
 447       Double_t mass = 0;
 448
 449       // 1. identify protons
 450       mass = 0.938;
 451       if (TMath::Abs(y - foProton->Eval(x))/foProton->Eval(x) < 4*sigma && x < 0.65 ) {
 452         for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
 453       }
 454
 455       // 2. identify electrons
 456       mass = 0.000511;
 457       if (fIsCosmic) {
 458         if (TMath::Abs(y - foElectron->Eval(x))/foElectron->Eval(x) < 4*sigma && (x>0.25&&x<0.7) && fIsCosmic) {
 459           for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
 460         }
 461       } else {
 462         if (TMath::Abs(y - foElectron->Eval(x))/foElectron->Eval(x) < 2*sigma && ((x>0.25&&x<0.35) || (x>1.5&&x<1.8) || (x>0.65&&x<0.7))) {
 463           for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
 464         }
 465       }
 466
 467       // 3. identify either muons or pions depending on cosmic or not
 468       if (fIsCosmic) {
 469         mass = 0.1056;
 470         if (TMath::Abs(y - foMuon->Eval(x))/foMuon->Eval(x) < 4*sigma && x > 0.25 && x < 100 ) {
 471           for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
 472         }
 473       } else {
 474         mass = 0.1396;
 475         if (TMath::Abs(y - foPion->Eval(x))/foPion->Eval(x) < 4*sigma && x > 0.15 && x < 2) {
 476           for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
 477         }
 478       }
 479
 480       // 4. for pp also Kaons must be included
 481       if (!fIsCosmic) {
 482         mass = 0.4936;
 483         if (TMath::Abs(y - foKaon->Eval(x))/foKaon->Eval(x) < 4*sigma && x > 0.2 && x < 0.3 ) {
 484           for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
 485         }
 486       }
 487     }
 488   }
 489
 490   // Fit Aleph-Parameters to the obtained profile
 491   TF1 * funcAlephD = new TF1("AlephParametrizationD", "AliExternalTrackParam::BetheBlochAleph(x,[0],[1],[2],[3],[4])",0.3,10000);
 492   funcAlephD->SetParameters(ini);
 493
 494   TCanvas *CanvCheck2 = new TCanvas();
 495   histBG->Draw();
 496
 497   //FitSlices
 498   TObjArray * arr = new TObjArray();
 499   histBG->FitSlicesY(0,0,-1,0,"QN",arr);
 500   TH1D * fitMean = (TH1D*) arr->At(1);
 501   fitMean->Draw("same");
 502
 503   funcAlephD->SetParLimits(2,1e-3,1e-7);
 504   funcAlephD->SetParLimits(3,0.5,3.5);
 505   funcAlephD->SetParLimits(4,0.5,3.5);
 506   fitMean->Fit(funcAlephD, "QNR");
 507   funcAlephD->Draw("same");
 508
 509   for(Int_t i=0;i<5;i++) ini[i] = funcAlephD->GetParameter(i);
 510
 511   foElectron->SetParameters(ini);
 512   foMuon->SetParameters(ini);
 513   foPion->SetParameters(ini);
 514   foKaon->SetParameters(ini);
 515   foProton->SetParameters(ini);
 516
 517   TCanvas *CanvCheck3 = new TCanvas();
 518   hist->Draw("colz");
 519   foElectron->Draw("same");
 520   foMuon->Draw("same");
 521   foPion->Draw("same");
 522   foKaon->Draw("same");
 523   foProton->Draw("same");
 524
 525   CanvCheck1->Print();
 526   CanvCheck2->Print();
 527   CanvCheck3->Print();
 528
 529   return;
 530
 531
 532 }