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