remove AliTPCkalmandEdx.cxx AliTPCkalmandEdx.h
authormarian <marian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 14 Apr 2009 16:51:07 +0000 (16:51 +0000)
committermarian <marian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 14 Apr 2009 16:51:07 +0000 (16:51 +0000)
adding AliTPCcalibTimeGain.cxx AliTPCcalibTimeGain.h

AliTPCcalibTracksGain -removing kalmandEdx part+ histograms - created in the AliTPCcalibPID

(Alexander Kalweit)

TPC/AliTPCcalibTimeGain.cxx [new file with mode: 0644]
TPC/AliTPCcalibTimeGain.h [new file with mode: 0644]
TPC/AliTPCcalibTracksGain.cxx
TPC/AliTPCcalibTracksGain.h
TPC/AliTPCkalmandEdx.cxx [deleted file]
TPC/AliTPCkalmandEdx.h [deleted file]
TPC/TPCcalibLinkDef.h
TPC/libTPCcalib.pkg

diff --git a/TPC/AliTPCcalibTimeGain.cxx b/TPC/AliTPCcalibTimeGain.cxx
new file mode 100644 (file)
index 0000000..cc0fe79
--- /dev/null
@@ -0,0 +1,518 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/*
+
+//0.  Libraries to lod
+gSystem->Load("libANALYSIS");
+gSystem->Load("libSTAT");
+gSystem->Load("libTPCcalib");
+
+
+//1. Do calibration ...
+//
+// compare reference
+
+//
+//2. Visulaize results
+//
+TFile fcalib("CalibObjects.root");
+TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
+AliTPCcalibTimeGain * gain = ( AliTPCcalibTimeGain *)array->FindObject("calibTimeGain");
+TGraph * gr = AliTPCcalibBase::FitSlices(gain->GetHistGainTime(),0,1,2000,10);
+gr->Draw("ALPsame")
+
+
+//
+// MakeSlineFit example
+//
+AliSplineFit fit;
+fit.SetGraph(gr)
+fit->SetMinPoints(gr->GetN()+1);
+fit->InitKnots(gr,2,0,0.001)
+fit.SplineFit(0)
+fit.MakeDiffHisto(gr)->Draw();
+TGraph * grfit = fit.MakeGraph(gr->GetX()[0],gr->GetX()[gr->GetN()-1],50000,0);
+
+gr->SetMarkerStyle(25);
+gr->Draw("alp");
+grfit->SetLineColor(2);
+grfit->Draw("lu");
+
+*/
+
+
+#include "Riostream.h"
+#include "TChain.h"
+#include "TTree.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TH3F.h"
+#include "THnSparse.h"
+#include "TList.h"
+#include "TMath.h"
+#include "TCanvas.h"
+#include "TFile.h"
+#include "TF1.h"
+#include "TVectorD.h"
+#include "TProfile.h"
+#include "TGraphErrors.h"
+#include "TCanvas.h"
+
+#include "AliTPCclusterMI.h"
+#include "AliTPCseed.h"
+#include "AliESDVertex.h"
+#include "AliESDEvent.h"
+#include "AliESDfriend.h"
+#include "AliESDInputHandler.h"
+#include "AliAnalysisManager.h"
+
+#include "AliTracker.h"
+#include "AliMagF.h"
+#include "AliTPCCalROC.h"
+
+#include "AliLog.h"
+
+#include "AliTPCcalibTimeGain.h"
+
+#include "TTreeStream.h"
+#include "AliTPCTracklet.h"
+#include "TTimeStamp.h"
+#include "AliTPCcalibDB.h"
+#include "AliTPCcalibLaser.h"
+#include "AliDCSSensorArray.h"
+#include "AliDCSSensor.h"
+
+ClassImp(AliTPCcalibTimeGain)
+
+
+AliTPCcalibTimeGain::AliTPCcalibTimeGain() 
+  :AliTPCcalibBase(), 
+   fHistGainTime(0),
+   fGainVsTime(0),
+   fHistDeDxTotal(0),
+   fIntegrationTimeDeDx(0),
+   fMIP(0),
+   fLowerTrunc(0),
+   fUpperTrunc(0),
+   fUseShapeNorm(0),
+   fUsePosNorm(0),
+   fUsePadNorm(0),
+   fIsCosmic(0)
+{  
+  AliInfo("Default Constructor");  
+}
+
+
+AliTPCcalibTimeGain::AliTPCcalibTimeGain(const Text_t *name, const Text_t *title, UInt_t StartTime, UInt_t EndTime, Int_t deltaIntegrationTimeGain)
+  :AliTPCcalibBase(),
+   fHistGainTime(0),
+   fGainVsTime(0),
+   fHistDeDxTotal(0),
+   fIntegrationTimeDeDx(0),
+   fMIP(0),
+   fLowerTrunc(0),
+   fUpperTrunc(0),
+   fUseShapeNorm(0),
+   fUsePosNorm(0),
+   fUsePadNorm(0),
+   fIsCosmic(0)
+ {
+  
+  SetName(name);
+  SetTitle(title);
+
+  AliInfo("Non Default Constructor");
+
+  fIntegrationTimeDeDx = deltaIntegrationTimeGain;
+  Double_t deltaTime = EndTime - StartTime;
+  
+
+  // 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)
+  Int_t binsGainTime[5] = {100, TMath::Nint(deltaTime/deltaIntegrationTimeGain), 3, 25, 200};
+  Double_t xminGainTime[5] = {0.5, StartTime, -0.5,   0, 0.1};
+  Double_t xmaxGainTime[5] = {  4,   EndTime,  2.5, 250, 50};
+  fHistGainTime = new THnSparseF("HistGainTime","dEdx l;time;dEdx",5,binsGainTime,xminGainTime,xmaxGainTime);
+  BinLogX(fHistGainTime, 4);
+  //
+  fHistDeDxTotal = new TH2F("DeDx","dEdx; momentum p (GeV); TPC signal (a.u.)",250,0.01,100.,1000,0.,1000);
+  BinLogX(fHistDeDxTotal);
+  
+  // default values for dE/dx
+  fMIP = 50.;
+  fLowerTrunc = 0.0;
+  fUpperTrunc = 0.7;
+  fUseShapeNorm = kTRUE;
+  fUsePosNorm = kFALSE;
+  fUsePadNorm = kFALSE;
+  //
+  fIsCosmic = kTRUE;
+  //
+  
+ }
+
+
+
+AliTPCcalibTimeGain::~AliTPCcalibTimeGain(){
+  //
+  //
+  //
+}
+
+
+void AliTPCcalibTimeGain::Process(AliESDEvent *event) {
+  //
+  // main track loop
+  //
+  if (!event) {
+    Printf("ERROR: ESD not available");
+    return;
+  }  
+
+  Int_t ntracks=event->GetNumberOfTracks();
+  
+  AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
+  if (!ESDfriend) {
+   Printf("ERROR: ESDfriend not available");
+   return;
+  }
+  //
+  UInt_t time = event->GetTimeStamp();
+  if (time < 0.1) time = (UInt_t)(fHistGainTime->GetAxis(1)->GetXmin() + 1);
+  //
+  // track loop
+  //
+  for (Int_t i=0;i<ntracks;++i) {
+
+    AliESDtrack *track = event->GetTrack(i);
+    if (!track) continue;
+        
+    const AliExternalTrackParam * trackIn = track->GetInnerParam();
+    const AliExternalTrackParam * trackOut = track->GetOuterParam();
+    if (!trackIn) continue;
+    if (!trackOut) continue;
+
+    // calculate necessary track parameters
+    //Double_t meanP = 0.5*(trackIn->GetP() + trackOut->GetP());
+    Double_t meanP = trackIn->GetP();
+    Double_t meanDrift = 250 - 0.5*TMath::Abs(trackIn->GetZ() + trackOut->GetZ());
+    //Double_t d = trackIn->GetLinearD(0,0);
+    Int_t NclsDeDx = track->GetTPCNcls();
+
+    //if (meanP > 0.7 || meanP < 0.2) continue;
+    if (fIsCosmic && meanP < 20) continue;
+    if (NclsDeDx < 60) continue;     
+
+    // exclude tracks which do not look like primaries or are simply too short or on wrong sectors
+
+    //if (TMath::Abs(trackIn->GetSnp()) > 3*0.4) continue;
+    //if (TMath::Abs(trackIn->GetZ()) > 150) continue;   
+    //if (seed->CookShape(1) > 1) continue;
+    //if (TMath::Abs(trackIn->GetY()) > 20) continue;
+    //if (TMath::Abs(d)>20) continue;   // distance to the 0,0; select only tracks which cross chambers under proper angle
+    if (TMath::Abs(trackIn->GetTgl()) > 1) continue;
+    if (TMath::Abs(trackIn->GetSnp()) > 0.6) continue;
+    
+    // Get seeds
+    AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
+    if (!friendTrack) continue;
+    TObject *calibObject;
+    AliTPCseed *seed = 0;
+    for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
+      if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
+    }    
+
+    if (seed) { 
+      Double_t TPCsignalMax = (1/fMIP)*seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,0,159,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
+      //Double_t TPCsignalMax = (1/fMIP)*track->GetTPCsignal();
+      fHistDeDxTotal->Fill(meanP, TPCsignalMax);
+
+      //
+      //dE/dx, time, type (0-all,1-muon cosmic,2-pion beam data), momenta
+      Double_t vec[5] = {TPCsignalMax,time,0,meanDrift,meanP}; 
+      fHistGainTime->Fill(vec); // avoid this filling if memory consumption is too high
+
+      // only partial filling if memory consumption has to be kept low; for cosmic and beam data
+      if (fIsCosmic) {
+       Double_t vecCos[5] = {TPCsignalMax,time,1,meanDrift,20}; 
+       if (meanP > 20) fHistGainTime->Fill(vecCos);
+      } else {
+       Double_t vecBeam[5] = {TPCsignalMax,time,2,meanDrift,0.5};
+       if (meanP > 0.4 && meanP < 0.5) fHistGainTime->Fill(vecBeam);
+      }
+    
+
+    } else {
+      cout << "ERROR: TPC seed not found" << endl;
+    }
+    
+  }
+  
+  
+}
+
+
+void AliTPCcalibTimeGain::Analyze() {
+  //
+  //
+  //
+  TObjArray arr;
+  if (fIsCosmic) {
+    fHistGainTime->GetAxis(2)->SetRangeUser(0.51,1.49);
+  } else {
+    fHistGainTime->GetAxis(2)->SetRangeUser(1.51,2.49);
+  }
+  fHistGainTime->Projection(0,1)->FitSlicesY(0,0,-1,0,"QNR",&arr);
+  TH1D * fitMean = (TH1D*) arr.At(1);
+  //
+  fGainVsTime = new TGraph(fitMean);
+  //
+  return;
+}
+
+
+Long64_t AliTPCcalibTimeGain::Merge(TCollection *li) {
+
+  TIterator* iter = li->MakeIterator();
+  AliTPCcalibTimeGain* cal = 0;
+
+  while ((cal = (AliTPCcalibTimeGain*)iter->Next())) {
+    if (!cal->InheritsFrom(AliTPCcalibTimeGain::Class())) {
+      Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
+      return -1;
+    }
+
+    // add histograms here...
+    if (cal->GetHistGainTime()) fHistGainTime->Add(cal->GetHistGainTime());
+    if (cal->GetHistDeDxTotal()) fHistDeDxTotal->Add(cal->GetHistDeDxTotal());
+
+  }
+  
+  return 0;
+  
+}
+
+
+TGraph * AliTPCcalibTimeGain::FitSlices(THnSparse *h, Int_t axisDim1, Int_t axisDim2, Int_t minEntries){
+  //
+  // Fitting slices of the projection(axisDim1,axisDim2) of a sparse histogram
+  //
+  TH2D * hist = h->Projection(axisDim1, axisDim2);
+
+  Double_t xvec[10000];
+  Double_t yvec[10000];
+  Int_t counter = 0;
+
+  for(Int_t i=1; i < hist->GetNbinsX(); i++) {
+    Int_t interval = 0;
+    if (hist->Integral(i,i,0,hist->GetNbinsY()) < minEntries) {
+      if (hist->Integral(i,i+1,0,hist->GetNbinsY()) < minEntries) {
+       if (hist->Integral(i,i+2,0,hist->GetNbinsY()) < minEntries) {
+         continue;
+       } else {
+         interval = 2;
+       }
+      } else {
+       interval = 1;
+      }
+    }
+    counter++;
+    i += interval;
+    //
+    Double_t x = hist->GetXaxis()->GetBinCenter(i); 
+    TH1D * projectionHist = hist->ProjectionY("dummy",i,i + interval);
+    TF1 funcGaus("funcGaus","gaus");
+    projectionHist->Fit(&funcGaus,"QN");
+    //  
+    xvec[counter-1] = x;
+    yvec[counter-1] = funcGaus.GetParameter(1);
+    delete projectionHist;
+  }
+  
+  TGraph * graph = new TGraph(counter, xvec, yvec);
+
+  return graph;
+}
+
+
+
+void AliTPCcalibTimeGain::BinLogX(THnSparse *h, Int_t axisDim) {
+
+  // Method for the correct logarithmic binning of histograms
+
+  TAxis *axis = h->GetAxis(axisDim);
+  int bins = axis->GetNbins();
+
+  Double_t from = axis->GetXmin();
+  Double_t to = axis->GetXmax();
+  Double_t *new_bins = new Double_t[bins + 1];
+   
+  new_bins[0] = from;
+  Double_t factor = pow(to/from, 1./bins);
+  
+  for (int i = 1; i <= bins; i++) {
+   new_bins[i] = factor * new_bins[i-1];
+  }
+  axis->Set(bins, new_bins);
+  delete new_bins;
+  
+}
+
+
+void AliTPCcalibTimeGain::BinLogX(TH1 *h) {
+
+  // Method for the correct logarithmic binning of histograms
+
+  TAxis *axis = h->GetXaxis();
+  int bins = axis->GetNbins();
+
+  Double_t from = axis->GetXmin();
+  Double_t to = axis->GetXmax();
+  Double_t *new_bins = new Double_t[bins + 1];
+   
+  new_bins[0] = from;
+  Double_t factor = pow(to/from, 1./bins);
+  
+  for (int i = 1; i <= bins; i++) {
+   new_bins[i] = factor * new_bins[i-1];
+  }
+  axis->Set(bins, new_bins);
+  delete new_bins;
+  
+}
+
+
+
+void AliTPCcalibTimeGain::CalculateBetheAlephParams(TH2F *hist, Double_t * ini) {
+  //{0.0762*MIP,10.632,1.34e-05,1.863,1.948}
+  const Double_t sigma = 0.06;
+
+  TH2F *histBG = new TH2F("histBG","dEdxBg; #beta #gamma; TPC signal (a.u.)",hist->GetNbinsX(),0.1,5000.,hist->GetNbinsY(),0.5,5.);
+  BinLogX(histBG);
+
+  TF1 *foElectron = new TF1("foElectron", "(1.7/1.6)*AliExternalTrackParam::BetheBlochAleph(x/0.000511,[0],[1],[2],[3],[4])",0.01,100);
+  TF1 *foMuon = new TF1("foMuon", "AliExternalTrackParam::BetheBlochAleph(x/0.1056,[0],[1],[2],[3],[4])",0.01,100);
+  TF1 *foPion = new TF1("foPion", "AliExternalTrackParam::BetheBlochAleph(x/0.138,[0],[1],[2],[3],[4])",0.01,100);
+  TF1 *foKaon = new TF1("foKaon", "AliExternalTrackParam::BetheBlochAleph(x/0.498,[0],[1],[2],[3],[4])",0.01,100);
+  TF1 *foProton = new TF1("foProton", "AliExternalTrackParam::BetheBlochAleph(x/0.938,[0],[1],[2],[3],[4])",0.01,100);
+  foElectron->SetParameters(ini);
+  foMuon->SetParameters(ini);
+  foPion->SetParameters(ini);
+  foKaon->SetParameters(ini);
+  foProton->SetParameters(ini);
+  
+  TCanvas *CanvCheck1 = new TCanvas();
+  hist->Draw("colz");
+  foElectron->Draw("same");
+  foMuon->Draw("same");
+  foPion->Draw("same");
+  foKaon->Draw("same");  
+  foProton->Draw("same");
+  // Loop over all points of the input histogram
+  
+  for(Int_t i=1; i < hist->GetNbinsX(); i++) {
+    Double_t x = hist->GetXaxis()->GetBinCenter(i);   
+    for(Int_t j=1; j < hist->GetNbinsY(); j++) {
+      Long64_t n = hist->GetBin(i, j);
+      Double_t y = hist->GetYaxis()->GetBinCenter(j);
+      Double_t entries = hist->GetBinContent(n);
+      Double_t mass = 0;
+
+      // 1. identify protons
+      mass = 0.938;
+      if (TMath::Abs(y - foProton->Eval(x))/foProton->Eval(x) < 4*sigma && x < 0.65 ) {
+       for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
+      }
+
+      // 2. identify electrons
+      mass = 0.000511;
+      if (fIsCosmic) {
+       if (TMath::Abs(y - foElectron->Eval(x))/foElectron->Eval(x) < 4*sigma && (x>0.25&&x<0.7) && fIsCosmic) {
+         for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
+       }
+      } else {
+       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))) {
+         for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
+       }
+      }
+      
+      // 3. identify either muons or pions depending on cosmic or not
+      if (fIsCosmic) {
+       mass = 0.1056;
+       if (TMath::Abs(y - foMuon->Eval(x))/foMuon->Eval(x) < 4*sigma && x > 0.25 && x < 100 ) {
+         for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
+       }
+      } else {
+       mass = 0.1396;
+       if (TMath::Abs(y - foPion->Eval(x))/foPion->Eval(x) < 4*sigma && x > 0.15 && x < 2) {
+         for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
+       }
+      }
+      
+      // 4. for pp also Kaons must be included
+      if (!fIsCosmic) {
+       mass = 0.4936;
+       if (TMath::Abs(y - foKaon->Eval(x))/foKaon->Eval(x) < 4*sigma && x > 0.2 && x < 0.3 ) {
+         for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
+       }
+      }
+    }
+  }
+
+  // Fit Aleph-Parameters to the obtained profile
+  TF1 * funcAlephD = new TF1("AlephParametrizationD", "AliExternalTrackParam::BetheBlochAleph(x,[0],[1],[2],[3],[4])",0.3,10000);
+  funcAlephD->SetParameters(ini);
+
+  TCanvas *CanvCheck2 = new TCanvas();
+  histBG->Draw();
+  
+  //FitSlices
+  TObjArray * arr = new TObjArray();
+  histBG->FitSlicesY(0,0,-1,0,"QN",arr);
+  TH1D * fitMean = (TH1D*) arr->At(1);
+  fitMean->Draw("same");
+
+  funcAlephD->SetParLimits(2,1e-3,1e-7);
+  funcAlephD->SetParLimits(3,0.5,3.5);
+  funcAlephD->SetParLimits(4,0.5,3.5);
+  fitMean->Fit(funcAlephD, "QNR");
+  funcAlephD->Draw("same");
+
+  for(Int_t i=0;i<5;i++) ini[i] = funcAlephD->GetParameter(i);
+
+  foElectron->SetParameters(ini);
+  foMuon->SetParameters(ini);
+  foPion->SetParameters(ini);
+  foKaon->SetParameters(ini);
+  foProton->SetParameters(ini);
+  
+  TCanvas *CanvCheck3 = new TCanvas();
+  hist->Draw("colz");
+  foElectron->Draw("same");
+  foMuon->Draw("same");
+  foPion->Draw("same");
+  foKaon->Draw("same");  
+  foProton->Draw("same");
+  
+  CanvCheck1->Print();
+  CanvCheck2->Print();
+  CanvCheck3->Print();
+
+  return;
+
+
+}
diff --git a/TPC/AliTPCcalibTimeGain.h b/TPC/AliTPCcalibTimeGain.h
new file mode 100644 (file)
index 0000000..5e31fac
--- /dev/null
@@ -0,0 +1,80 @@
+#ifndef ALITPCCALIBTIMEGAIN_H
+#define ALITPCCALIBTIMEGAIN_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+#include "AliTPCcalibBase.h"
+#include "TH2F.h"
+#include "TF1.h"
+#include "TArrayD.h"
+#include "TObjArray.h"
+
+class TH1F;
+class TH3F;
+class TH2F;
+class THnSparse;
+class TList;
+class AliESDEvent;
+class AliESDtrack;
+class AliTPCcalibLaser;
+
+#include "TTreeStream.h"
+
+class AliTPCcalibTimeGain:public AliTPCcalibBase {
+public:
+  AliTPCcalibTimeGain(); 
+  AliTPCcalibTimeGain(const Text_t *name, const Text_t *title, UInt_t StartTime, UInt_t EndTime, Int_t deltaIntegrationTimeGain);
+  virtual ~AliTPCcalibTimeGain();
+  //
+  virtual void           Process(AliESDEvent *event);
+  virtual Long64_t       Merge(TCollection *li);
+  virtual void           Analyze();
+  //
+  void                   CalculateBetheAlephParams(TH2F *hist, Double_t * ini);
+  static void            BinLogX(THnSparse *h, Int_t axisDim);
+  static void            BinLogX(TH1 *h);
+  static TGraph *        FitSlices(THnSparse *h, Int_t axisDim1, Int_t axisDim2, Int_t minEntries);
+  //
+  THnSparse *            GetHistGainTime(){return (THnSparse*) fHistGainTime;};
+  TGraph    *            GetGraphGainVsTime(){return fGainVsTime;};
+  TH2F      *            GetHistDeDxTotal(){return (TH2F*) fHistDeDxTotal;};
+  //
+  void SetMIP(Float_t MIP){fMIP = MIP;};  
+  void SetLowerTrunc(Float_t LowerTrunc){fLowerTrunc = LowerTrunc;};
+  void SetUpperTrunc(Float_t UpperTrunc){fUpperTrunc = UpperTrunc;};
+  void SetUseShapeNorm(Bool_t UseShapeNorm){fUseShapeNorm = UseShapeNorm;};
+  void SetUsePosNorm(Bool_t UsePosNorm){fUsePosNorm = UsePosNorm;};
+  void SetUsePadNorm(Int_t UsePadNorm){fUsePadNorm = UsePadNorm;};
+  void SetIsCosmic(Bool_t IsCosmic){fIsCosmic = IsCosmic;};
+
+private:
+  //
+  THnSparse * fHistGainTime;            // dEdx vs. time, type, Driftlength, momentum P
+  TGraph    * fGainVsTime;              // multiplication factor vs. time
+  TH2F      * fHistDeDxTotal;           // dEdx vs. momentum for quality assurance
+  //
+  Float_t fIntegrationTimeDeDx;         // required statistics for each dEdx time bin
+  //
+  Float_t fMIP;                         // rough MIP position in order to have scalable histograms
+  //
+  Float_t fLowerTrunc;                  // lower truncation of dE/dx ; at most 5%
+  Float_t fUpperTrunc;                  // upper truncation of dE/dx ; ca. 70%
+  Bool_t  fUseShapeNorm;                // use empirical correction of dependencies
+  Bool_t  fUsePosNorm;                  // charge correction (analytical?)
+  Int_t   fUsePadNorm;                  // normalization of pad geometries
+  //
+  Bool_t  fIsCosmic;                    // kTRUE if the analyzed runs are contain cosmic events
+  //
+  AliTPCcalibTimeGain(const AliTPCcalibTimeGain&); 
+  AliTPCcalibTimeGain& operator=(const AliTPCcalibTimeGain&); 
+  void     Process(AliESDtrack *track, Int_t runNo=-1){AliTPCcalibBase::Process(track,runNo);};
+  void     Process(AliTPCseed *track){return AliTPCcalibBase::Process(track);}
+
+  ClassDef(AliTPCcalibTimeGain, 1); 
+};
+
+#endif
+
+
index 232c1f6..33a4298 100644 (file)
 
 
 
+#include "AliTPCcalibTracksGain.h"
 
 
 #include <TPDGCode.h>
 #include "AliTPCClusterParam.h"
 #include "AliTrackPointArray.h"
 #include "TCint.h"
-#include "AliTPCcalibTracksGain.h"
 #include <TH1.h>
 #include <TH3F.h>
 #include <TLinearFitter.h>
 #include "AliTPCCalROC.h"
 #include "AliTPCCalPad.h"
 #include "AliTPCClusterParam.h"
+#include "AliTPCcalibDB.h"
 //
 #include "AliTracker.h"
 #include "AliESD.h"
@@ -244,16 +245,6 @@ AliTPCParamSR* AliTPCcalibTracksGain::fgTPCparam = new AliTPCParamSR();
 AliTPCcalibTracksGain::AliTPCcalibTracksGain() :
   AliTPCcalibBase(),
   fCuts(0),            // cuts that are used for sieving the tracks used for calibration
-  fGainMap(0),                //  gain map to be applied
-  //
-  // Simple Array of histograms
-  // 
-  fArrayQM(0),                // Qmax normalized
-  fArrayQT(0),                // Qtot normalized 
-  fProfileArrayQM(0),         // Qmax normalized  versus local X
-  fProfileArrayQT(0),         // Qtot normalized  versus local X 
-  fProfileArrayQM2D(0),       // Qmax normalized  versus local X and phi
-  fProfileArrayQT2D(0),       // Qtot normalized  versus local X and phi
   //
   // Fitters
   //
@@ -280,9 +271,7 @@ AliTPCcalibTracksGain::AliTPCcalibTracksGain() :
   // Counters
   //
   fTotalTracks(0),         // just for debugging
-  fAcceptedTracks(0),      // just for debugging
-  fDebugCalPadRaw(0),      // just for debugging
-  fDebugCalPadCorr(0)     // just for debugging
+  fAcceptedTracks(0)      // just for debugging
 
 {
    //
@@ -293,16 +282,6 @@ AliTPCcalibTracksGain::AliTPCcalibTracksGain() :
 AliTPCcalibTracksGain::AliTPCcalibTracksGain(const AliTPCcalibTracksGain& obj) :
   AliTPCcalibBase(obj),
   fCuts(obj.fCuts),            // cuts that are used for sieving the tracks used for calibration
-  fGainMap(new AliTPCCalPad(*(obj.fGainMap))),                //  gain map to be applied
-  fArrayQM(0),                // Qmax normalized
-  fArrayQT(0),                // Qtot normalized 
-  //
-  // Simple Histograms
-  //
-  fProfileArrayQM(obj.fProfileArrayQM),         // Qmax normalized  versus local X
-  fProfileArrayQT(obj.fProfileArrayQT),         // Qtot normalized  versus local X 
-  fProfileArrayQM2D(obj.fProfileArrayQM2D),       // Qmax normalized  versus local X and phi
-  fProfileArrayQT2D(obj.fProfileArrayQT2D),       // Qtot normalized  versus local X and phi
   //
   // Fitters
   //
@@ -327,9 +306,7 @@ AliTPCcalibTracksGain::AliTPCcalibTracksGain(const AliTPCcalibTracksGain& obj) :
   // Counters
   //
   fTotalTracks(obj.fTotalTracks),         // just for debugging
-  fAcceptedTracks(obj.fAcceptedTracks),      // just for debugging
-  fDebugCalPadRaw(obj.fDebugCalPadRaw),      // just for debugging
-  fDebugCalPadCorr(obj.fDebugCalPadCorr)     // just for debugging
+  fAcceptedTracks(obj.fAcceptedTracks)      // just for debugging
 
 {
    //
@@ -344,14 +321,11 @@ AliTPCcalibTracksGain& AliTPCcalibTracksGain::operator=(const AliTPCcalibTracksG
 
    if (this != &rhs) {
       TNamed::operator=(rhs);
-      fDebugCalPadRaw = new AliTPCCalPad(*(rhs.fDebugCalPadRaw));
-      fDebugCalPadCorr = new AliTPCCalPad(*(rhs.fDebugCalPadCorr));
       fSimpleFitter = new AliTPCFitPad(*(rhs.fSimpleFitter));
       fSqrtFitter = new AliTPCFitPad(*(rhs.fSqrtFitter));
       fLogFitter = new AliTPCFitPad(*(rhs.fLogFitter));
       fSingleSectorFitter = new AliTPCFitPad(*(rhs.fSingleSectorFitter));
       fCuts = new AliTPCcalibTracksCuts(*(rhs.fCuts));
-      fGainMap = new AliTPCCalPad(*(rhs.fGainMap));
    }
    return *this;
 }
@@ -359,16 +333,6 @@ AliTPCcalibTracksGain& AliTPCcalibTracksGain::operator=(const AliTPCcalibTracksG
 AliTPCcalibTracksGain::AliTPCcalibTracksGain(const char* name, const char* title, AliTPCcalibTracksCuts* cuts) :
   AliTPCcalibBase(),
   fCuts(0),            // cuts that are used for sieving the tracks used for calibration
-  fGainMap(0),                //  gain map to be applied
-  fArrayQM(0),                // Qmax normalized
-  fArrayQT(0),                // Qtot normalized 
-  //
-  // Simple Histograms
-  //
-  fProfileArrayQM(0),         // Qmax normalized  versus local X
-  fProfileArrayQT(0),         // Qtot normalized  versus local X 
-  fProfileArrayQM2D(0),       // Qmax normalized  versus local X and phi
-  fProfileArrayQT2D(0),       // Qtot normalized  versus local X and phi
   //
   // Fitters
   //
@@ -393,9 +357,7 @@ AliTPCcalibTracksGain::AliTPCcalibTracksGain(const char* name, const char* title
   // Counters
   //
   fTotalTracks(0),         // just for debugging
-  fAcceptedTracks(0),      // just for debugging
-  fDebugCalPadRaw(0),      // just for debugging
-  fDebugCalPadCorr(0)     // just for debugging
+  fAcceptedTracks(0)      // just for debugging
   
 {
    //
@@ -441,47 +403,17 @@ AliTPCcalibTracksGain::AliTPCcalibTracksGain(const char* name, const char* title
    fDFitter2T->StoreData(kFALSE);   
    //
    //
-   // Add profile histograms -JUST for visualization - Not used for real calibration
-   // 
-   //
-   fArrayQM=new TObjArray(73);                // Qmax normalized
-   fArrayQT=new TObjArray(73);                // Qtot normalized 
-   fProfileArrayQM   = new TObjArray(37);         // Qmax normalized  versus local X
-   fProfileArrayQT   = new TObjArray(37);         // Qtot normalized  versus local X 
-   fProfileArrayQM2D = new TObjArray(37);         // Qmax normalized  versus local X and phi
-   fProfileArrayQT2D = new TObjArray(37);         // Qtot normalized  versus local X and phi
-   char hname[1000];
-   for (Int_t i=0; i<73; i++){
-     sprintf(hname,"QM_%d",i);
-     fArrayQM->AddAt(new TH1F(hname,hname,200,0,1000),i);
-     sprintf(hname,"QT_%d",i);
-     fArrayQT->AddAt(new TH1F(hname,hname,200,0,1000),i);
-   }
-
-   for (Int_t i=0; i<37;i++){
-     sprintf(hname,"QMvsx_%d",i);
-     fProfileArrayQM->AddAt(new TProfile(hname,hname,50,89,250),i);
-     sprintf(hname,"QTvsx_%d",i);
-     fProfileArrayQT->AddAt(new TProfile(hname,hname,50,89,250),i);
-     sprintf(hname,"QM2D_%d",i);
-     fProfileArrayQM2D->AddAt(new TProfile2D(hname,hname,50,89,250,10,-0.15,0.15),i);
-     sprintf(hname,"QT2D_%d",i);
-     fProfileArrayQT2D->AddAt(new TProfile2D(hname,hname,50,89,250,10,-0.15,0.15),i);
-   }
-   //
    // just for debugging -counters
    //
    fTotalTracks     = 0;
    fAcceptedTracks  = 0;
-   fDebugCalPadRaw  = new AliTPCCalPad("DebugCalPadRaw", "All clusters simply added up before correction");
-   fDebugCalPadCorr = new AliTPCCalPad("DebugCalPadCorr", "All clusters simply added up after correction");   
    // this will be gone for the a new ROOT version > v5-17-05
    for (UInt_t i = 0; i < 36; i++) {
       fNShortClusters[i]  = 0;
       fNMediumClusters[i] = 0;
       fNLongClusters[i]   = 0;
    }
- }
+}
 
 AliTPCcalibTracksGain::~AliTPCcalibTracksGain() {
    //
@@ -494,10 +426,11 @@ AliTPCcalibTracksGain::~AliTPCcalibTracksGain() {
    if (fLogFitter) delete fLogFitter;
    if (fSingleSectorFitter) delete fSingleSectorFitter;
 
-   if (fDebugCalPadRaw) delete fDebugCalPadRaw;
-   if (fDebugCalPadCorr) delete fDebugCalPadCorr;
 }
 
+
+
+
 void AliTPCcalibTracksGain::Terminate(){
    //
    // Evaluate fitters and close the debug stream.
@@ -569,9 +502,6 @@ Long64_t AliTPCcalibTracksGain::Merge(TCollection *list) {
    if (!fSingleSectorFitter) fSingleSectorFitter = new AliTPCFitPad(8, "hyp7", "");
 
 
-   // just for debugging
-   if (!fDebugCalPadRaw) fDebugCalPadRaw  = new AliTPCCalPad("DebugCalPadRaw", "All clusters simply added up before correction");
-    if (!fDebugCalPadCorr) fDebugCalPadCorr = new AliTPCCalPad("DebugCalPadCorr", "All clusters simply added up after correction");
    
    TIterator* iter = list->MakeIterator();
    AliTPCcalibTracksGain* cal = 0;
@@ -587,40 +517,39 @@ Long64_t AliTPCcalibTracksGain::Merge(TCollection *list) {
    return 0;
 }
 
-Float_t  AliTPCcalibTracksGain::GetGain(AliTPCclusterMI* cl){
+Float_t  AliTPCcalibTracksGain::GetGain(AliTPCclusterMI *cluster){
   //
-  // Return local gain at cluster position
+  // get gain
   //
-  Float_t factor = 1;
-  if(!fGainMap) return factor;
-
-  AliTPCCalROC * roc = fGainMap->GetCalROC(cl->GetDetector());
-  Int_t irow = cl->GetRow();
-  if (roc){
-    if (irow < 63) { // IROC
-      factor = roc->GetValue(irow, TMath::Nint(cl->GetPad()));
-    } else {         // OROC
-      factor = roc->GetValue(irow - 63, TMath::Nint(cl->GetPad()));
-    }
+  Float_t gainPad= 1;
+  AliTPCCalPad * gainMap =  AliTPCcalibDB::Instance()->GetDedxGainFactor();
+  if (gainMap) {
+    AliTPCCalROC * roc = gainMap->GetCalROC(cluster->GetDetector());
+    gainPad  = roc->GetValue(cluster->GetRow(), TMath::Nint(cluster->GetPad()));
   }
-  if (factor<0.1) factor=0.1;
-  return factor;
+  return gainPad;
 }
 
-
-Float_t   AliTPCcalibTracksGain::GetMaxNorm(AliTPCclusterMI * cl){
+Float_t   AliTPCcalibTracksGain::GetMaxNorm(AliTPCclusterMI * cluster, Float_t ky, Float_t kz){
   //
   // Get normalized amplituded if the gain map provided
   //
-  return cl->GetMax()/GetGain(cl);
+  AliTPCClusterParam * parcl = AliTPCcalibDB::Instance()->GetClusterParam();
+  Float_t maxNorm =
+    parcl->QmaxCorrection(cluster->GetDetector(), cluster->GetRow(),cluster->GetPad(), 
+                         cluster->GetTimeBin(),ky,kz,0.5,0.5,1.6);
+
+  return GetGain(cluster)*maxNorm;
 }
 
 
-Float_t   AliTPCcalibTracksGain::GetQNorm(AliTPCclusterMI * cl){
+Float_t   AliTPCcalibTracksGain::GetQNorm(AliTPCclusterMI * cluster,  Float_t ky, Float_t kz){
   //
   // Get normalized amplituded if the gain map provided
   //
-  return cl->GetQ()/GetGain(cl);
+  AliTPCClusterParam * parcl = AliTPCcalibDB::Instance()->GetClusterParam();
+  Float_t totNorm = parcl->QtotCorrection(cluster->GetDetector(), cluster->GetRow(),cluster->GetPad(),                               cluster->GetTimeBin(),ky,kz,0.5,0.5,cluster->GetQ(),2.5,1.6);
+  return GetGain(cluster)*totNorm;
 }
 
 
@@ -651,41 +580,6 @@ void AliTPCcalibTracksGain::Add(AliTPCcalibTracksGain* cal) {
   if (cal->fDFitter1T->GetNpoints()>0) fDFitter1T->Add(cal->fDFitter1T);
   if (cal->fDFitter2T->GetNpoints()>0) fDFitter2T->Add(cal->fDFitter2T);
    //
-   //
-   // histograms
-   //
-   for (Int_t i=0; i<73; i++){
-     TH1F *his,*hism; 
-     his = (TH1F*)fArrayQM->At(i);
-     hism =  (TH1F*)cal->fArrayQM->At(i);
-     if (his && hism) his->Add(hism);
-     his =  (TH1F*)fArrayQT->At(i);
-     hism =  (TH1F*)cal->fArrayQT->At(i);
-     if (his && hism) his->Add(hism);
-   }
-   //
-   //
-   for (Int_t i=0; i<37; i++){
-     TProfile *his,*hism; 
-     his = (TProfile*)fProfileArrayQM->At(i);
-     hism = (TProfile*)cal->fProfileArrayQM->At(i);
-     if (his && hism) his->Add(hism);
-     his = (TProfile*)fProfileArrayQT->At(i);
-     hism = (TProfile*)cal->fProfileArrayQT->At(i);
-     if (his && hism) his->Add(hism);
-   }
-   //
-   //
-   for (Int_t i=0; i<37; i++){
-     TProfile2D *his,*hism; 
-     his = (TProfile2D*)fProfileArrayQM2D->At(i);
-     hism = (TProfile2D*)cal->fProfileArrayQM2D->At(i);
-     if (his && hism) his->Add(hism);
-     his = (TProfile2D*)fProfileArrayQT2D->At(i);
-     hism = (TProfile2D*)cal->fProfileArrayQT2D->At(i);
-     if (his && hism) his->Add(hism);
-   }
-   //
    // this will be gone for the a new ROOT version > v5-17-05
    for (UInt_t iSegment = 0; iSegment < 36; iSegment++) {
       fNShortClusters[iSegment] += cal->fNShortClusters[iSegment];
@@ -696,8 +590,6 @@ void AliTPCcalibTracksGain::Add(AliTPCcalibTracksGain* cal) {
    // just for debugging, remove me
    fTotalTracks += cal->fTotalTracks;
    fAcceptedTracks += cal->fAcceptedTracks;
-   fDebugCalPadRaw->Add(cal->fDebugCalPadRaw);
-   fDebugCalPadCorr->Add(cal->fDebugCalPadCorr);
 
 }
 
@@ -708,66 +600,14 @@ void AliTPCcalibTracksGain::AddTrack(AliTPCseed* seed) {
    //
    
    DumpTrack(seed);   
-   //
-   // simple histograming part
-   for (Int_t i=0; i<159; i++){
-     AliTPCclusterMI* cluster = seed->GetClusterPointer(i);
-     if (cluster) AddCluster(cluster);
-   }
 }
    
-void AliTPCcalibTracksGain::AddCluster(AliTPCclusterMI* cluster){
-  //
-  // Adding cluster information to the simple histograms
-  // No correction, fittings are applied   
-
-
-  Float_t kThreshold=5;
-  if (cluster->GetX()<=0) return;
-  if (cluster->GetQ()<=kThreshold) return;
-  //
-
-  Int_t sector = cluster->GetDetector();
-  TH1F  * his;
-  his = GetQT(sector);
-  if (his) his->Fill(GetQNorm(cluster));
-  his = GetQT(-1);
-  if (his) his->Fill(GetQNorm(cluster));
-  his = GetQM(sector);
-  if (his) his->Fill(GetMaxNorm(cluster));
-  his = GetQM(-1);
-  if (his) his->Fill(GetMaxNorm(cluster));
-  //
-  sector = sector%36;
-  TProfile *prof;
-  prof = GetProfileQT(sector);
-  if (prof) prof->Fill(cluster->GetX(),GetQNorm(cluster));
-  prof = GetProfileQT(-1);
-  if (prof) prof->Fill(cluster->GetX(),GetQNorm(cluster));
-  prof = GetProfileQM(sector);
-  if (prof) prof->Fill(cluster->GetX(),GetMaxNorm(cluster));
-  prof = GetProfileQM(-1);
-  if (prof) prof->Fill(cluster->GetX(),GetMaxNorm(cluster));
-  //  
-  Float_t phi = cluster->GetY()/cluster->GetX();
-  TProfile2D *prof2;
-  prof2 = GetProfileQT2D(sector);
-  if (prof2) prof2->Fill(cluster->GetX(),phi,GetQNorm(cluster));
-  prof2 = GetProfileQT2D(-1);
-  if (prof2) prof2->Fill(cluster->GetX(),phi,GetQNorm(cluster));
-  prof2 = GetProfileQM2D(sector);
-  if (prof2) prof2->Fill(cluster->GetX(),phi,GetMaxNorm(cluster));
-  prof2 = GetProfileQM2D(-1);
-  if (prof2) prof2->Fill(cluster->GetX(),phi,GetMaxNorm(cluster));
-
-  //
-}
 
 
 
 void AliTPCcalibTracksGain::AddCluster(AliTPCclusterMI* cluster, Float_t /*momenta*/, Float_t/* mdedx*/, Int_t padType,
                                       Float_t xcenter, TVectorD& dedxQ, TVectorD& /*dedxM*/, Float_t /*fraction*/, Float_t fraction2, Float_t dedge,
-                                      TVectorD& /*parY*/, TVectorD& /*parZ*/, TVectorD& meanPos) {
+                                      TVectorD& parY, TVectorD& parZ, TVectorD& meanPos) {
    //
    // Adds cluster to the appropriate fitter for later analysis.
    // The charge used for the fit is the maximum charge for this specific cluster or the
@@ -825,27 +665,18 @@ void AliTPCcalibTracksGain::AddCluster(AliTPCclusterMI* cluster, Float_t /*momen
    xx[4] = xx[0] * xx[1];
    xx[5] = TMath::Abs(cluster->GetZ()) - TMath::Abs(meanPos[4]);
    xx[6] = xx[5] * xx[5];
-   //
-   // Update profile histograms
-   //
 
    //
    // Update fitters
    //
    Int_t segment = cluster->GetDetector() % 36;
-   Double_t q = fgkUseTotalCharge ? ((Double_t)(GetQNorm(cluster))) : ((Double_t)(GetMaxNorm(cluster)));  // note: no normalization to pad size!
-
-   // just for debugging
-   Int_t row = 0;
-   Int_t pad = 0;
-   GetRowPad(cluster->GetX(), cluster->GetY(), row, pad);
-   fDebugCalPadRaw->GetCalROC(cluster->GetDetector())->SetValue(row, pad, q + fDebugCalPadRaw->GetCalROC(cluster->GetDetector())->GetValue(row, pad));
-   
+   Double_t q = fgkUseTotalCharge ? 
+     ((Double_t)(cluster->GetQ()/GetQNorm(cluster,parY[1], parZ[1]))) : ((Double_t)(cluster->GetMax()/GetMaxNorm(cluster,parY[1], parZ[1])));  
+      
    // correct charge by normalising to mean charge per track
    q /= dedxQ[kinorm];
 
    // just for debugging
-   fDebugCalPadCorr->GetCalROC(cluster->GetDetector())->SetValue(row, pad, q + fDebugCalPadCorr->GetCalROC(cluster->GetDetector())->GetValue(row, pad));
 
    Double_t sqrtQ = TMath::Sqrt(q);
    Double_t logQ = fgkM * TMath::Log(1 + q / fgkM);
@@ -1160,35 +991,6 @@ Int_t AliTPCcalibTracksGain::GetPadType(Double_t lx) {
    return -1;
 }
 
-// ONLY FOR DEBUGGING PURPOSES - REMOVE ME WHEN NOT NEEDED ANYMORE
-Bool_t AliTPCcalibTracksGain::GetRowPad(Double_t lx, Double_t ly, Int_t& row, Int_t& pad) {
-   //
-   // Calculate the row and pad number when the local coordinates are given.
-   // Returns kFALSE if the position is out of range, otherwise return kTRUE.
-   // WARNING: This function is preliminary and probably isn't very accurate!!
-   //
-   
-   Double_t irocLow = fgTPCparam->GetPadRowRadiiLow(0) - fgTPCparam->GetInnerPadPitchLength()/2;
-   //Double_t irocUp = fgTPCparam->GetPadRowRadiiLow(fgTPCparam->GetNRowLow()-1) + fgTPCparam->GetInnerPadPitchLength()/2;
-   Double_t orocLow1 = fgTPCparam->GetPadRowRadiiUp(0) - fgTPCparam->GetOuter1PadPitchLength()/2;
-   //Double_t orocUp1 = fgTPCparam->GetPadRowRadiiUp(fgTPCparam->GetNRowUp1()-1) + fgTPCparam->GetOuter1PadPitchLength()/2;
-   Double_t orocLow2 = fgTPCparam->GetPadRowRadiiUp(fgTPCparam->GetNRowUp1()) - fgTPCparam->GetOuter2PadPitchLength()/2;
-   //Double_t orocUp2 = fgTPCparam->GetPadRowRadiiUp(fgTPCparam->GetNRowUp()-1) + fgTPCparam->GetOuter2PadPitchLength()/2;
-
-   if (GetPadType(lx) == 0) {
-      row = (Int_t)((lx - irocLow) / fgTPCparam->GetInnerPadPitchLength());
-      pad = (Int_t)((ly + fgTPCparam->GetYInner(row)) / fgTPCparam->GetInnerPadPitchWidth());
-   } else if (GetPadType(lx) == 1) {
-      row = (Int_t)((lx - orocLow1) / fgTPCparam->GetOuter1PadPitchLength());
-      pad = (Int_t)((ly + fgTPCparam->GetYOuter(row)) / fgTPCparam->GetOuterPadPitchWidth());
-   } else if (GetPadType(lx) == 2) {
-      row = fgTPCparam->GetNRowUp1() + (Int_t)((lx - orocLow2) / fgTPCparam->GetOuter2PadPitchLength());
-      pad = (Int_t)((ly + fgTPCparam->GetYOuter(row)) / fgTPCparam->GetOuterPadPitchWidth());
-   }
-   else return kFALSE;
-   return kTRUE;
-}
-
 void AliTPCcalibTracksGain::DumpTrack(AliTPCseed* track) {
    //
    //  Dump track information to the debug stream
@@ -1286,34 +1088,34 @@ Bool_t AliTPCcalibTracksGain::GetDedx(AliTPCseed* track, Int_t padType, Int_t* /
   fitZ.StoreData(kFALSE);
   //
   //   
-   Int_t firstRow = 0, lastRow = 0;
-   Int_t minRow = 100;
-   Float_t xcenter = 0;
-   const Float_t ktany = TMath::Tan(TMath::DegToRad() * 10);
-   const Float_t kedgey = 4.;
-   if (padType == 0) {
-      firstRow = 0;
-      lastRow = fgTPCparam->GetNRowLow();
-      xcenter = 108.47;
-   }
-   if (padType == 1) {
-      firstRow = fgTPCparam->GetNRowLow();
-      lastRow = fgTPCparam->GetNRowLow() + fgTPCparam->GetNRowUp1();
-      xcenter = 166.60;
-   }
-   if (padType == 2) {
-      firstRow = fgTPCparam->GetNRowLow() + fgTPCparam->GetNRowUp1();
-      lastRow = fgTPCparam->GetNRowLow() + fgTPCparam->GetNRowUp();
-      xcenter = 222.6;
+  Int_t firstRow = 0, lastRow = 0;
+  Int_t minRow = 100;
+  Float_t xcenter = 0;
+  const Float_t ktany = TMath::Tan(TMath::DegToRad() * 10);
+  const Float_t kedgey = 4.;
+  if (padType == 0) {
+    firstRow = 0;
+    lastRow = fgTPCparam->GetNRowLow();
+    xcenter = 108.47;
+  }
+  if (padType == 1) {
+    firstRow = fgTPCparam->GetNRowLow();
+    lastRow = fgTPCparam->GetNRowLow() + fgTPCparam->GetNRowUp1();
+    xcenter = 166.60;
    }
-   minRow = (lastRow - firstRow) / 2;
-   //
-   //
-   Int_t nclusters = 0;
-   Int_t nclustersNE = 0; // number of not edge clusters
-   Int_t lastSector = -1;
-   Float_t amplitudeQ[100];
-   Float_t amplitudeM[100];
+  if (padType == 2) {
+    firstRow = fgTPCparam->GetNRowLow() + fgTPCparam->GetNRowUp1();
+    lastRow = fgTPCparam->GetNRowLow() + fgTPCparam->GetNRowUp();
+    xcenter = 222.6;
+  }
+  minRow = (lastRow - firstRow) / 2;
+  //
+  //
+  Int_t nclusters = 0;
+  Int_t nclustersNE = 0; // number of not edge clusters
+  Int_t lastSector = -1;
+  Float_t amplitudeQ[100];
+  Float_t amplitudeM[100];
    Int_t rowIn[100];
    Int_t index[100];
    //
@@ -1327,8 +1129,8 @@ Bool_t AliTPCcalibTracksGain::GetDedx(AliTPCseed* track, Int_t padType, Int_t* /
          Int_t detector = cluster->GetDetector() ;
          if (lastSector == -1) lastSector = detector;
          if (lastSector != detector) continue;
-         amplitudeQ[nclusters] = GetQNorm(cluster);
-         amplitudeM[nclusters] = GetMaxNorm(cluster);
+         amplitudeQ[nclusters] = cluster->GetQ()/GetQNorm(cluster,parY[1], parZ[1]);
+         amplitudeM[nclusters] = cluster->GetMax()/GetMaxNorm(cluster,parY[1], parZ[1]);
          rowIn[nclusters] = iCluster;
          nclusters++;
          Double_t dx = cluster->GetX() - xcenter;
@@ -1423,6 +1225,11 @@ Bool_t AliTPCcalibTracksGain::GetDedx(AliTPCseed* track, Int_t padType, Int_t* /
       Float_t fraction2 = Float_t(inonEdge) / Float_t(nclustersNE);
 
       AddCluster(cluster, momenta, mdedx, padType, xcenter, dedxQ, dedxM, fraction, fraction2, dedge, parY, parZ, meanPos);
+      
+      Double_t qNorm = GetQNorm(cluster,parY[1], parZ[1]);
+      Double_t mNorm = GetMaxNorm(cluster,parY[1], parZ[1]);
+
+
       Float_t gain = GetGain(cluster);
       if (cstream) (*cstream) << "dEdx" <<
        "run="<<fRun<<              //  run number
@@ -1431,21 +1238,23 @@ Bool_t AliTPCcalibTracksGain::GetDedx(AliTPCseed* track, Int_t padType, Int_t* /
        "trigger="<<fTrigger<<      //  trigger
        "mag="<<fMagF<<             //  magnetic field        
 
-                    "Cl.=" << cluster <<           // cluster of interest
-                    "gain="<<gain<<                // gain at cluster position
-                    "P=" << momenta <<             // track momenta
-                    "dedx=" << mdedx <<            // mean dedx - corrected for angle
-                    "IPad=" << padType <<          // pad type 0..2
-                    "xc=" << xcenter <<            // x center of chamber
-                    "dedxQ.=" << &dedxQ <<         // dedxQ  - total charge
-                    "dedxM.=" << &dedxM <<         // dedxM  - maximal charge
-                    "fraction=" << fraction <<     // fraction - order in statistic (0,1)
-                    "fraction2=" << fraction2 <<   // fraction - order in statistic (0,1)
-                    "dedge=" << dedge <<           // distance to the edge
-                    "parY.=" << &parY <<           // line fit
-                    "parZ.=" << &parZ <<           // line fit
-                    "meanPos.=" << &meanPos <<     // mean position (dx, dx^2, y,y^2, z, z^2)
-                    "\n";
+       "Cl.=" << cluster <<           // cluster of interest
+       "gain="<<gain<<                // gain at cluster position
+       "mNorm="<<mNorm<<              // Q max normalization
+       "qNorm="<<qNorm<<              // Q tot normalization
+       "P=" << momenta <<             // track momenta
+       "dedx=" << mdedx <<            // mean dedx - corrected for angle
+       "IPad=" << padType <<          // pad type 0..2
+       "xc=" << xcenter <<            // x center of chamber
+       "dedxQ.=" << &dedxQ <<         // dedxQ  - total charge
+       "dedxM.=" << &dedxM <<         // dedxM  - maximal charge
+       "fraction=" << fraction <<     // fraction - order in statistic (0,1)
+       "fraction2=" << fraction2 <<   // fraction - order in statistic (0,1)
+       "dedge=" << dedge <<           // distance to the edge
+       "parY.=" << &parY <<           // line fit
+       "parZ.=" << &parZ <<           // line fit
+       "meanPos.=" << &meanPos <<     // mean position (dx, dx^2, y,y^2, z, z^2)
+       "\n";
    }
    
    if (cstream) (*cstream) << "dEdxT" <<
@@ -1600,194 +1409,3 @@ void   AliTPCcalibTracksGain::Analyze(){
 }
 
 
-
-TVectorD *  AliTPCcalibTracksGain::MakeQPosNorm(TTree * chain0, Int_t ipad, Bool_t isMax, Int_t maxPoints, Int_t verbose){
-  //
-  // Input parameters
-  // chain0 - the tree with information  -Debug stream
-  // ipad   - 0 IROC
-  //        - 1 OROC medium
-  //        - 2 OROC LONG
-  // isMax  - kFALSE - total charge param
-  //          kTRUE  - Max charge param
-  //
-  // maxPoints - number of points for fit
-  //
-  // verbose -
-  //
-  /* e.g 
-    ipad=0
-    isMax=kTRUE;
-    maxPoints=1000000;
-  */
-  // Make Q normalization as function of following parameters
-  // 1 - dp   - relative pad position 
-  // 2 - dt   - relative time position
-  // 3 - di   - drift length (norm to 1);
-  // 4 - dq0  - Tot/Max charge
-  // 5 - dq1  - Max/Tot charge
-  // 6 - sy   - sigma y - shape
-  // 7 - sz   - sigma z - shape
-  //
-  // Coeficient of Taylor expansion fitted
-  // Fit parameters returned as TVectorD
-  // Fit parameters to be used in corresponding correction function
-  // in AliTPCclusterParam 
-  //
-  //
-  TStatToolkit toolkit;
-  Double_t chi2;
-  TVectorD fitParam;
-  TMatrixD covMatrix;
-  Int_t npoints;
-  TCut cutA("dedge>3&&fraction2<0.7");
-  chain0->SetAlias("dp","((Cl.fPad-int(Cl.fPad)-0.5)/0.5)");
-  chain0->SetAlias("dt","((Cl.fTimeBin-int(Cl.fTimeBin)-0.5)/0.5)");
-  chain0->SetAlias("di","(sqrt(1.-abs(Cl.fZ)/250.))");
-  chain0->SetAlias("dq0","(0.2*(Cl.fQ+2)/(Cl.fMax+2))");
-  chain0->SetAlias("dq1","(5*(Cl.fMax+2)/(Cl.fQ+2))");
-  chain0->SetAlias("sy","(0.32/sqrt(0.01^2+Cl.fSigmaY2))");
-  chain0->SetAlias("sz","(0.32/sqrt(0.01^2+Cl.fSigmaZ2))");
-  //
-  TString fstring="";    
-  fstring+="dp++";                               //1
-  fstring+="dt++";                               //2
-  fstring+="dp*dp++";                             //3
-  fstring+="dt*dt++";                             //4
-  fstring+="dt*dt*dt++";                             //5
-  fstring+="dp*dt++";                            //6
-  fstring+="dp*dt*dt++";                          //7
-  fstring+="(dq0)++";                            //8
-  fstring+="(dq1)++";                            //9
-  //
-  //
-  fstring+="dp*dp*(di)++";                        //10
-  fstring+="dt*dt*(di)++";                        //11
-  fstring+="dp*dp*sy++";                          //12
-  fstring+="dt*sz++";                          //13
-  fstring+="dt*dt*sz++";                          //14
-  fstring+="dt*dt*dt*sz++";                          //15
-  //
-  fstring+="dp*dp*1*sy*sz++";                     //16
-  fstring+="dt*sy*sz++";                       //17
-  fstring+="dt*dt*sy*sz++";                       //18
-  fstring+="dt*dt*dt*sy*sz++";                       //19
-  //
-  fstring+="dp*dp*(dq0)++";                       //20
-  fstring+="dt*1*(dq0)++";                       //21
-  fstring+="dt*dt*(dq0)++";                       //22
-  fstring+="dt*dt*dt*(dq0)++";                       //23
-  //
-  fstring+="dp*dp*(dq1)++";                       //24
-  fstring+="dt*(dq1)++";                       //25
-  fstring+="dt*dt*(dq1)++";                       //26
-  fstring+="dt*dt*dt*(dq1)++";                       //27
-  
-  TString var;
-  if (isMax)  var = "Cl.fMax/gain/dedxM.fElements[2]";
-  if (!isMax) var = "Cl.fQ/gain/dedxQ.fElements[2]";
-  TString cutP="IPad==";
-  cutP+=ipad;
-  //
-  TString *strq0 = toolkit.FitPlane(chain0,var.Data(),fstring.Data(), cutP.Data()+cutA, chi2,npoints,fitParam,covMatrix,-1,0,maxPoints);
-  //
-  //
-  if (verbose){
-    printf("Chi2/npoints = %f",TMath::Sqrt(chi2/npoints));
-    printf("\nFit function\n:%s\n",strq0->Data());
-  }
-  TVectorD *vec = new TVectorD(fitParam);
-  return vec;
-}
-
-void  AliTPCcalibTracksGain::MakeQPosNormAll(TTree * chain, AliTPCClusterParam * param, Int_t maxPoints, Int_t verbose){
-  //
-  // Fill the content of the of the AliTPCclusterParam
-  // with fitted values of corrections 
-  //
-  param->fPosQTnorm[0] = MakeQPosNorm(chain,0,kTRUE,maxPoints,verbose);
-  param->fPosQTnorm[1] = MakeQPosNorm(chain,1,kTRUE,maxPoints,verbose);
-  param->fPosQTnorm[2] = MakeQPosNorm(chain,1,kTRUE,maxPoints,verbose);
-  //
-  param->fPosQMnorm[0] = MakeQPosNorm(chain,0,kFALSE,maxPoints,verbose);
-  param->fPosQMnorm[1] = MakeQPosNorm(chain,1,kFALSE,maxPoints,verbose);
-  param->fPosQMnorm[2] = MakeQPosNorm(chain,2,kFALSE,maxPoints,verbose);
-}
-
-
-
-/*
-  
- Position correction fit:
- //
-TStatToolkit toolkit;
-Double_t chi2;
-TVectorD fitParam;
-TMatrixD covMatrix;
-Int_t npoints;
-//
-TCut cutA("dedge>3&&fraction2<0.7");
-chain0->SetAlias("dp","((Cl.fPad-int(Cl.fPad)-0.5)/0.5)");
-chain0->SetAlias("dt","((Cl.fTimeBin-int(Cl.fTimeBin)-0.5)/0.5)");
-chain0->SetAlias("di","(sqrt(1.-abs(Cl.fZ)/250.))");
-chain0->SetAlias("dq0","(0.2*(Cl.fQ+2)/(Cl.fMax+2))");
-chain0->SetAlias("dq1","(5*(Cl.fMax+2)/(Cl.fQ+2))");
-chain0->SetAlias("sy","(0.2/sqrt(0.01^2+Cl.fSigmaY2))");
-chain0->SetAlias("sz","(0.2/sqrt(0.01^2+Cl.fSigmaZ2))");
-
-TString fstring="";  
-
-fstring+="dp++";                               //1
-fstring+="dt++";                               //2
-fstring+="dp*dp++";                             //3
-fstring+="dt*dt++";                             //4
-fstring+="dt*dt*dt++";                             //5
-fstring+="dp*dt++";                            //6
-fstring+="dp*dt*dt++";                          //7
-fstring+="(dq0)++";                            //8
-fstring+="(dq1)++";                            //9
-//
-//
-fstring+="dp*dp*(di)++";                        //10
-fstring+="dt*dt*(di)++";                        //11
-fstring+="dp*dp*sy++";                          //12
-fstring+="dt*sz++";                          //13
-fstring+="dt*dt*sz++";                          //14
-fstring+="dt*dt*dt*sz++";                          //15
-//
-fstring+="dp*dp*1*sy*sz++";                     //16
-fstring+="dt*sy*sz++";                       //17
-fstring+="dt*dt*sy*sz++";                       //18
-fstring+="dt*dt*dt*sy*sz++";                       //19
-//
-fstring+="dp*dp*(dq0)++";                       //20
-fstring+="dt*1*(dq0)++";                       //21
-fstring+="dt*dt*(dq0)++";                       //22
-fstring+="dt*dt*dt*(dq0)++";                       //23
-
-fstring+="dp*dp*(dq1)++";                       //24
-fstring+="dt*(dq1)++";                       //25
-fstring+="dt*dt*(dq1)++";                       //26
-fstring+="dt*dt*dt*(dq1)++";                       //27
-
-
- TString *strq0 = toolkit.FitPlane(chain0,"Cl.fMax/gain/dedxM.fElements[2]",fstring->Data(), "IPad==0"+cutA, chi2,npoints,fitParam,covMatrix,-1,0,200000);
- TString *strqt0 = toolkit.FitPlane(chain0,"Cl.fQ/gain/dedxQ.fElements[2]",fstring->Data(), "IPad==0"+cutA, chi2,npoints,fitParam,covMatrix,-1,0,200000);
-
- chain0->SetAlias("qcorM0",strq0->Data());
- chain0->SetAlias("qcorT0",strqt0->Data());
-//chain0->SetAlias("mmqcorM0","min(max(qcorM0,0.75),1.15)");
- chain0->Draw("(Cl.fMax/gain/dedxM.fElements[2]):min(max(qcorM0,0.75),1.15)","IPad==0"+cutA,"prof",100000)
-
- fraction05 - 
-     sigma                  0.2419
-     sigma fit              0.2302
-     sigma fit with shape   0.2257
- fraction 07    
-     qtot sigma                  0.322
-     qmax sigma                  0.292
-     qmax sigma fit              0.2702
-     qmax sigma fit+ratio        0.2638
-
-*/
-
index 1f09ece..767e047 100644 (file)
@@ -62,13 +62,12 @@ public:
   virtual void            Process(AliTPCseed* seed);
   virtual void            Terminate();
   virtual void            Analyze();
-  void                    SetGainMap(AliTPCCalPad *GainMap){fGainMap = GainMap;};
   //
   // Tracks and cluster manipulation
   //
-  Float_t         GetMaxNorm(AliTPCclusterMI * cl);
-  Float_t         GetQNorm(AliTPCclusterMI * cl);
-  Float_t         GetGain(AliTPCclusterMI* cl);
+  Float_t         GetGain(AliTPCclusterMI *cl);
+  Float_t         GetMaxNorm(AliTPCclusterMI * cl,  Float_t ky, Float_t kz);
+  Float_t         GetQNorm(AliTPCclusterMI * cl,  Float_t ky, Float_t kz);
   void            AddTrack(AliTPCseed* seed);
   void            DumpTrack(AliTPCseed* track);
   Bool_t          GetDedx(AliTPCseed* track, Int_t padType,  Int_t*rows,
@@ -78,26 +77,9 @@ public:
   void            AddCluster(AliTPCclusterMI* cluster, Float_t momenta, Float_t mdedx, Int_t padType, Float_t xcenter, TVectorD &dedxQ, TVectorD & dedxM, Float_t fraction, Float_t fraction2, Float_t dedge, TVectorD& parY, TVectorD& parZ, TVectorD& meanPos);
   void            AddTracklet(UInt_t sector, UInt_t padType,TVectorD &dedxQ, TVectorD &dedxM,TVectorD& parY, TVectorD& parZ, TVectorD& meanPos);
 
-  void            AddCluster(AliTPCclusterMI* cluster);
   void            Add(AliTPCcalibTracksGain* cal);
 
   //
-  // Debug stream analyze part
-  //
-  static TVectorD * MakeQPosNorm(TTree * chain, Int_t ipad, Bool_t isMax, Int_t maxPoints=1000000, Int_t verbose=0);
-  
-  static void MakeQPosNormAll(TTree * chain, AliTPCClusterParam * param, Int_t maxPoints=1000000, Int_t verbose=0);
-
-  //
-  // Histogram part
-  //
-  TH1F  * GetQM(Int_t sector=-1){return (TH1F*)(sector<0 ?  fArrayQM->At(72): fArrayQM->At(sector));}
-  TH1F  * GetQT(Int_t sector=-1){return (TH1F*)(sector<0 ?  fArrayQT->At(72): fArrayQT->At(sector));}
-  TProfile* GetProfileQM(Int_t sector){return (TProfile*)(sector<0 ? fProfileArrayQM->At(36): fProfileArrayQM->At(sector));}
-  TProfile* GetProfileQT(Int_t sector){return (TProfile*)(sector<0 ? fProfileArrayQT->At(36): fProfileArrayQT->At(sector));}
-  TProfile2D* GetProfileQM2D(Int_t sector){return (TProfile2D*)(sector<0 ? fProfileArrayQM2D->At(36): fProfileArrayQM2D->At(sector));}
-  TProfile2D* GetProfileQT2D(Int_t sector){return (TProfile2D*)(sector<0 ? fProfileArrayQT2D->At(36): fProfileArrayQT2D->At(sector));}
-  //
   // Get Derived results - gain maps
   //
   AliTPCCalPad*   CreateFitCalPad(UInt_t fitType, Bool_t undoTransformation = kFALSE, Bool_t normalizeToPadSize = kFALSE);
@@ -126,21 +108,12 @@ public:
   //
   static Double_t GetPadLength(Double_t lx);
   static Int_t    GetPadType(Double_t lx);  
-  static Bool_t   GetRowPad(Double_t lx, Double_t ly, Int_t& row, Int_t& pad); // just for debugging
   //
   //
   AliTPCcalibTracksCuts* fCuts;            // cuts that are used for sieving the tracks used for calibration
-  AliTPCCalPad *fGainMap;                //  gain map to be applied
-  //
   //
-  // Simple Profiles and histograms - per chambers + 1 total
+  // kalman fit - alignment of dEdx
   //
-  TObjArray*        fArrayQM;                // Qmax normalized
-  TObjArray*        fArrayQT;                // Qtot normalized 
-  TObjArray*        fProfileArrayQM;         // Qmax normalized  versus local X
-  TObjArray*        fProfileArrayQT;         // Qtot normalized  versus local X 
-  TObjArray*        fProfileArrayQM2D;       // Qmax normalized  versus local X and phi
-  TObjArray*        fProfileArrayQT2D;       // Qtot normalized  versus local X and phi
   //
   // Fitters
   //
@@ -164,14 +137,12 @@ public:
   TLinearFitter*    fDFitter1T;          // fitting of the atenuation, angular correction
   TLinearFitter*    fDFitter2T;          // fitting of the atenuation, angular correction
   //
-  AliTPCFitPad*     fSingleSectorFitter;   // just for debugging
+  AliTPCFitPad*     fSingleSectorFitter;   // just for debugging  
   //
   // Conters
   //
   UInt_t          fTotalTracks;         // just for debugging
   UInt_t          fAcceptedTracks;      // just for debugging
-  AliTPCCalPad*   fDebugCalPadRaw;      // just for debugging
-  AliTPCCalPad*   fDebugCalPadCorr;     // just for debugging
   UInt_t          fNShortClusters[36];   // number of clusters registered on short pads
   UInt_t          fNMediumClusters[36];  // number of clusters registered on medium pads
   UInt_t          fNLongClusters[36];    // number of clusters registered on long pads
diff --git a/TPC/AliTPCkalmandEdx.cxx b/TPC/AliTPCkalmandEdx.cxx
deleted file mode 100644 (file)
index d9173f7..0000000
+++ /dev/null
@@ -1,394 +0,0 @@
-/**************************************************************************
- * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- *                                                                        *
- * Author: The ALICE Off-line Project.                                    *
- * Contributors are mentioned in the code where appropriate.              *
- *                                                                        *
- * Permission to use, copy, modify and distribute this software and its   *
- * documentation strictly for non-commercial purposes is hereby granted   *
- * without fee, provided that the above copyright notice appears in all   *
- * copies and that both the copyright notice and this permission notice   *
- * appear in the supporting documentation. The authors make no claims     *
- * about the suitability of this software for any purpose. It is          *
- * provided "as is" without express or implied warranty.                  *
- **************************************************************************/
-
-/* 
-   Origin: Marian Ivanov   Marian.Ivanov@cern.ch
-   Not final version yet - fitting using histograms to be implemented
-   Expected to be more robust.                     
-
-
-   Model:
-   TPC Kalman filter implementation for TPC dEdx alignment:
-   The dEdx equalization for 3 types of the pad:
-   Short (0.75 cm) Medium (1 cm) and Long(1.5 cm)
-
-   Model of correction ci:
-   corrected   raw    correction
-   dEdxti    = dEdxri*ci = dEdxri*p01*(1+p1i*kY+p2i*kZ+p3i*dR+p4i/sqrt(1+kY^2+kZ^2)
-   //
-   Matching - update using 2 tracklets  ::UpdatedEdxPair
-   dEdxti-dEdxtj=0 = dEdxri*ci-dEdxrj*cj
-   //
-   Matching - normalization of the signal ::UpdatedEdx 
-   possible only for identified particles
-   dEdxti = dN/dx * sqrt(1+kY^2+kZ^2)
-*/
-
-#include "TMath.h"
-#include "TTreeStream.h"
-#include "TGraph.h"
-#include "TGraphErrors.h"
-#include "TVectorD.h"
-#include "AliTPCROC.h"
-#include "AliTPCkalmandEdx.h"
-
-
-AliTPCkalmandEdx::AliTPCkalmandEdx():
-  TNamed(),
-  fState(0),
-  fCovariance(0), 
-  fMat1(0),                  
-  fNpad(3),                  // number of pad types
-  fNpar(5),                // number of parameters
-  fNelem(3*5),                // number of elements
-  fSampleSize(0),
-  fInit(0)
-{
-  //
-  // Default constructor
-  //
-}
-
-AliTPCkalmandEdx::AliTPCkalmandEdx(const char* name, const char* title, Int_t sampleSize): 
-  TNamed(name,title),
-  fState(0),
-  fCovariance(0),   
-  fMat1(0),                  
-  fNpad(3),                  // number of pad types
-  fNpar(5),                // number of parameters
-  fNelem(3*5),             // number of elements
-  fSampleSize(sampleSize),
-  fInit(0)
-{
-  //
-  // Default constructor
-  //
-  Init();
-}
-
-AliTPCkalmandEdx::AliTPCkalmandEdx(const AliTPCkalmandEdx & kalman):
-  TNamed(kalman),
-  fState(0),
-  fCovariance(0),   
-  fMat1(0),                  
-  fNpad(kalman.fNpad),                  // number of pad types
-  fNpar(kalman.fNpar),                // number of parameters
-  fNelem(kalman.fNelem),                // number of elements
-  fSampleSize(kalman.fSampleSize),
-  fInit(kalman.fInit)
-{
-  //
-  // copy constructor
-  //
-  fState      = new TMatrixD(*(kalman.fState));
-  fCovariance = new TMatrixD(*(kalman.fCovariance));
-  fMat1       = new TMatrixD(*(kalman.fMat1));
-} 
-
-void AliTPCkalmandEdx::Init(){
-  //
-  // Default parameters setting
-  //
-  fState      = new TMatrixD(fNelem,1);
-  fCovariance = new TMatrixD(fNelem,fNelem);
-  fMat1       = new TMatrixD(fNelem,fNelem);
-
-  fInit=0;
-  for (Int_t i=0;i<3; i++) {
-    fSample[i].ResizeTo(fSampleSize);
-    fSampleStat[i].ResizeTo(2);
-    fCounter[i]=0;
-  }
-
-  TMatrixD &vecXk=*fState;
-  TMatrixD &covXk=*fCovariance;
-  TMatrixD &mat1=*fMat1;
-  //
-  //
-  for (Int_t i=0;i<fNelem;i++) 
-    for(Int_t j=0;j<fNelem;j++) {covXk(i,j)=0;  mat1(i,j)=0;}
-  //
-  for (Int_t i=0;i<fNelem;i++)  {covXk(i,i)=1.; mat1(i,i)=1; vecXk(i,0)=0;}
-  for (Int_t ipad=0;ipad<3;ipad++){
-    vecXk(ipad*fNpar,0)=1;
-  }
-//   //
-//   // set reference ipad=0
-//   vecXk(1*fNpar+0,0)=1;
-//   vecXk(1*fNpar+1,0)=0;
-//   vecXk(1*fNpar+2,0)=0;
-//   vecXk(1*fNpar+3,0)=0;
-//   vecXk(1*fNpar+4,0)=0;
-  
-
-}
-
-void AliTPCkalmandEdx::UpdatedEdxPair(Int_t ip0, Int_t ip1,
-                                     Double_t dedx0, Double_t dedx1, 
-                                     Double_t s0, Double_t s1, 
-                                     Double_t kY0, Double_t kY1,
-                                     Double_t kZ0, Double_t kZ1,
-                                     Double_t dR0, Double_t dR1, 
-                                     TTreeSRedirector *debug){
-  //
-  // update relative measurement
-  // 0 = dEdxti-dEdxtj  = dEdxri*ci-dEdxrj*cj
-  //
-  // Model of correction ci:
-  // dEdxti = dEdxri*ci = dEdxri*p01*(1+p1i*kY+p2i*kZ+p3i*dR+p4i/sqrt(1+kY^2+kZ^2)
-  if (fInit<3) return;  // not initialized parameters
-  const Double_t kchi2Cut = 3.*3.;
-
-  Int_t nmeas = 1;
-  TMatrixD vecXk=*fState;           // X vector
-  TMatrixD covXk=*fCovariance;      // X covariance
-  TMatrixD &mat1=*fMat1;            // update covariance matrix
-
-  Double_t length0 = TMath::Sqrt(1+kY0*kY0+kZ0*kZ0);
-  Double_t length1 = TMath::Sqrt(1+kY1*kY1+kZ1*kZ1);
-
-  Double_t corr0 = vecXk(ip0*fNpar,0)*
-    (1+vecXk(ip0*fNpar+1,0)*kY0+
-     vecXk(ip0*fNpar+2,0)*kZ0+vecXk(ip0*fNpar+3,0)*dR0+vecXk(ip0*fNpar+4,0)/length0);
-  Double_t corr1 = vecXk(ip1*fNpar,0)*
-    (1+vecXk(ip1*fNpar+1,0)*kY1+
-     vecXk(ip1*fNpar+2,0)*kZ1+vecXk(ip1*fNpar+3,0)*dR1+vecXk(ip1*fNpar+4,0)/length1);
-  //
-  TMatrixD vecZk(nmeas,nmeas);
-  TMatrixD measR(nmeas,nmeas);
-  TMatrixD matHk(nmeas,fNelem);     // vector to mesurement
-  TMatrixD vecYk(nmeas,nmeas);          // Innovation or measurement residual
-  TMatrixD matHkT(fNelem,nmeas);    // helper matrix Hk transpose
-  TMatrixD matSk(nmeas,nmeas);      // Innovation (or residual) covariance
-  TMatrixD matKk(fNelem,nmeas);     // Optimal Kalman gain
-  TMatrixD covXk2(fNelem,fNelem);   // helper matrix
-  TMatrixD covXk3(fNelem,fNelem);   // helper matrix
-
-  //
-  vecZk(0,0) =dedx1/dedx0;          // dedx ratio
-  measR(0,0) =(s0*s0+s1*s1);
-  //
-  // reset matHk
-  for (Int_t iel=0;iel<fNelem;iel++) 
-    for (Int_t ip=0;ip<nmeas;ip++) matHk(ip,iel)=0; 
-  //
-  matHk(0, ip0*fNpar+0) = corr0/(corr1*vecXk(ip0*fNpar+0,0));
-  matHk(0, ip0*fNpar+1) = (vecXk(ip0*fNpar+0,0)*kY0)/corr1;
-  matHk(0, ip0*fNpar+2) = (vecXk(ip0*fNpar+0,0)*kZ0)/corr1;
-  matHk(0, ip0*fNpar+3) = (vecXk(ip0*fNpar+0,0)*dR0)/corr1;
-  matHk(0, ip0*fNpar+4) = (vecXk(ip0*fNpar+0,0)/length0)/corr1;
-  //
-  //  matHk(0, ip1*fNpar+0) = ;
-  //matHk(0, ip1*fNpar+1) = ;
-  //matHk(0, ip1*fNpar+2) = ;
-  //matHk(0, ip1*fNpar+3) = ;
-  //matHk(0, ip1*fNpar+4) = ;
-  //
-  //
-  vecYk(0,0) = vecZk(0,0)- corr0/corr1;               // Innovation or measurement residual
-  matHkT=matHk.T(); matHk.T();
-  matSk = (matHk*(covXk*matHkT))+measR;    // Innovation (or residual) covariance
-  matSk.Invert();
-  matKk = (covXk*matHkT)*matSk;            //  Optimal Kalman gain
-  vecXk += matKk*vecYk;                    //  updated vector 
-  covXk2= (mat1-(matKk*matHk));
-  covXk3 =  covXk2*covXk;          
-  covXk = covXk3;  
-  Double_t chi2 = vecYk(0,0)* vecYk(0,0)*matSk(0,0);
-  //
-  //
-  //
-  if (debug){ // dump debug info
-    char chname[1000];
-    sprintf(chname,"updatePair_%s",GetName());
-    (*debug)<<chname<<
-      // input parameters                                      
-      "ip0="<<ip0<<
-      "ip1="<<ip1<<
-      "dedx0="<<dedx0<<
-      "dedx1="<<dedx1<<      
-      "s0="<<s0<<
-      "s1="<<s1<<
-      "kY0="<<kY0<<
-      "kY1="<<kY1<<
-      "kZ0="<<kZ0<<
-      "kZ1="<<kZ1<<
-      "dR0="<<dR0<<
-      "dR1="<<dR1<<
-      "chi2="<<chi2<<
-      "corr0="<<corr0<<
-      "corr1="<<corr1<<
-      "length0="<<length0<<
-      "length1="<<length1<<
-      //
-      "vecYk.="<<&vecYk<<
-      "vecZk.="<<&vecZk<<
-      "matHk.="<<&matHk<<
-      "matHkT.="<<&matHkT<<
-      "matSk.="<<&matSk<<
-      "matKk.="<<&matKk<<
-      "covXk2.="<<&covXk2<<
-      "covXk.="<<&covXk<<
-      "cov.="<<fCovariance<<
-      "vecXk.="<<&vecXk<<
-      "vec.="<<fState<<
-      "\n";
-  }
-  if (chi2<kchi2Cut){
-    //    (*fCovariance)=covXk;
-    //(*fState)=vecXk;
-  }
-  else{
-    printf("chi2=%f\n",chi2);
-  }
-}
-
-void AliTPCkalmandEdx::UpdatedEdx(Int_t ip0,
-                                 Double_t dedx0, 
-                                 Double_t dedxRef, 
-                                 Double_t s0,
-                                 Double_t kY0,
-                                 Double_t kZ0,
-                                 Double_t dR0,
-                                 TTreeSRedirector *debug){
-  //
-  // update relative measurement
-  // dEdx  = dEdxti = dEdxri*ci
-  //
-  // Model of correction ci:
-  // dEdxti = dEdxri*ci = dEdxri*p0i*(1+p1i*kY+p2i*kZ+p3i*dR+p4i/sqrt(1+kY^2+kZ^2)
-  const Double_t kchi2Cut = 3.*3.;
-  // removing of "big outliers
-  if (fInit<3) return AdddEdx(ip0,dedx0,dedxRef);
-  if (TMath::Abs(dedxRef/dedx0-fSampleStat[ip0][0])>4*fSampleStat[ip0][1]) return;
-  //
-  //
-  Int_t nmeas = 1;
-  TMatrixD vecXk=*fState;           // X vector
-  TMatrixD covXk=*fCovariance;      // X covariance
-  TMatrixD &mat1=*fMat1;            // update covariance matrix
-
-  Double_t length0 = TMath::Sqrt(1+kY0*kY0+kZ0*kZ0);
-  //
-  Double_t corr0 = vecXk(ip0*fNpar,0)*
-    (1+vecXk(ip0*fNpar+1,0)*kY0+
-     vecXk(ip0*fNpar+2,0)*kZ0+vecXk(ip0*fNpar+3,0)*dR0+vecXk(ip0*fNpar+4,0)/length0);
-
-
-  //
-  TMatrixD vecZk(nmeas,nmeas);
-  TMatrixD measR(nmeas,nmeas);
-  TMatrixD matHk(nmeas,fNelem);     // vector to mesurement
-  TMatrixD vecYk(nmeas,nmeas);          // Innovation or measurement residual
-  TMatrixD matHkT(fNelem,nmeas);    // helper matrix Hk transpose
-  TMatrixD matSk(nmeas,nmeas);      // Innovation (or residual) covariance
-  TMatrixD matKk(fNelem,nmeas);     // Optimal Kalman gain
-  TMatrixD covXk2(fNelem,fNelem);   // helper matrix
-  TMatrixD covXk3(fNelem,fNelem);   // helper matrix
-
-  vecZk(0,0) =dedxRef/dedx0;
-  measR(0,0) =s0*s0*dedxRef*dedxRef/(dedx0*dedx0);
-  //
-  // reset matHk
-  for (Int_t iel=0;iel<fNelem;iel++) 
-    for (Int_t ip=0;ip<nmeas;ip++) matHk(ip,iel)=0; 
-  //
-  //
-  matHk(0, ip0*fNpar+0) = corr0/vecXk(ip0*fNpar,0);
-  matHk(0, ip0*fNpar+1) = kY0*vecXk(ip0*fNpar,0);
-  matHk(0, ip0*fNpar+2) = kZ0*vecXk(ip0*fNpar,0);
-  matHk(0, ip0*fNpar+3) = dR0*vecXk(ip0*fNpar,0);
-  matHk(0, ip0*fNpar+4) = vecXk(ip0*fNpar,0)/length0;
-  //
-  //
-  vecYk(0,0) = vecZk(0,0)-corr0;               // Innovation or measurement residual
-  matHkT=matHk.T(); matHk.T();
-  matSk = (matHk*(covXk*matHkT))+measR;    // Innovation (or residual) covariance
-  matSk.Invert();
-  matKk = (covXk*matHkT)*matSk;            //  Optimal Kalman gain
-  vecXk += matKk*vecYk;                    //  updated vector 
-  covXk2= (mat1-(matKk*matHk));
-  covXk3 =  covXk2*covXk;          
-  covXk = covXk3;  
-  Double_t chi2 = vecYk(0,0)*vecYk(0,0)*matSk(0,0);
-  //
-  //
-  //
-  if (debug){ // dump debug info
-    char chname[1000];
-    sprintf(chname,"update_%s",GetName());
-    (*debug)<<chname<<
-      //
-      // input parameters                                      
-      "ip0="<<ip0<<
-      "dedxRef="<<dedxRef<<
-      "dedx0="<<dedx0<<
-      "s0="<<s0<<
-      "kY0="<<kY0<<
-      "kZ0="<<kZ0<<
-      "dR0="<<dR0<<
-      "chi2="<<chi2<<
-      "corr0="<<corr0<<
-      "length0="<<length0<<
-      //
-      "vecYk.="<<&vecYk<<
-      "vecZk.="<<&vecZk<<
-      "matHk.="<<&matHk<<
-      "matHkT.="<<&matHkT<<
-      "matSk.="<<&matSk<<
-      "matKk.="<<&matKk<<
-      "covXk2.="<<&covXk2<<
-      "covXk.="<<&covXk<<
-      "cov.="<<fCovariance<<
-      "vecXk.="<<&vecXk<<
-      "vec.="<<fState<<
-      "\n";
-  }
-  if (chi2<kchi2Cut){
-    (*fCovariance)=covXk;
-    (*fState)=vecXk;
-  }else{
-    printf("chi2=%f\n",chi2);
-  }
-}
-
-void AliTPCkalmandEdx::AdddEdx(Int_t ip0,Double_t dedx0, Double_t dedxRef){
-  //
-  //
-  //
-  if (fCounter[ip0]>=fSampleSize) return;
-  fSample[ip0][fCounter[ip0]]=dedxRef/dedx0;
-  fCounter[ip0]++;
-  if (fCounter[ip0]==fSampleSize){
-    fSampleStat[ip0].ResizeTo(2);
-    fSampleStat[ip0][0] = TMath::Median(fSampleSize,fSample[ip0].GetMatrixArray());
-    fSampleStat[ip0][1] = TMath::RMS(fSampleSize,fSample[ip0].GetMatrixArray());
-    fInit++;
-    (*fState)(ip0*fNpar,0)    = fSampleStat[ip0][0];
-    (*fCovariance)(ip0*fNpar,ip0*fNpar) = fSampleStat[ip0][1]*fSampleStat[ip0][1]/fSampleSize;
-    Double_t med2  = fSampleStat[ip0][0]*fSampleStat[ip0][0];
-    //
-    //
-    (*fCovariance)(ip0*fNpar+1,ip0*fNpar+1)=0.01*med2;
-    (*fCovariance)(ip0*fNpar+2,ip0*fNpar+2)=0.01*med2;
-    (*fCovariance)(ip0*fNpar+3,ip0*fNpar+3)=0.01*med2;
-    (*fCovariance)(ip0*fNpar+4,ip0*fNpar+4)=0.01*med2;
-  }
-}
-
-
-
diff --git a/TPC/AliTPCkalmandEdx.h b/TPC/AliTPCkalmandEdx.h
deleted file mode 100644 (file)
index d229318..0000000
+++ /dev/null
@@ -1,60 +0,0 @@
-#ifndef ALITPCKALMANDEDX_H
-#define ALITPCKALMANDEDX_H
-
-/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- * See cxx source for full Copyright notice                               */
-
-#include "TNamed.h"
-#include "TMatrixD.h"
-#include "TVectorD.h"
-class TTreeSRedirector;
-class TObjArray;
-
-class AliTPCkalmandEdx: public TNamed{
-public:
-  AliTPCkalmandEdx();
-  AliTPCkalmandEdx(const char* name, const char* title, Int_t sampleSize=50);
-  AliTPCkalmandEdx(const AliTPCkalmandEdx & kalman);
-  void UpdatedEdxPair(Int_t ip0, Int_t ip1,
-                     Double_t dedx0, Double_t dedx1, 
-                     Double_t s0, Double_t s1, 
-                     Double_t kY0, Double_t kY1,
-                     Double_t kZ0, Double_t kZ1,
-                     Double_t dR0, Double_t dR1, 
-                     TTreeSRedirector *debug=0);
-  void UpdatedEdx(Int_t ip0,
-                 Double_t dedx0, 
-                 Double_t dedxRef, 
-                 Double_t s0,
-                 Double_t kY0,
-                 Double_t kZ0,
-                 Double_t dR0,
-                 TTreeSRedirector *debug=0);
-  static Int_t GetIndex(Int_t i, Int_t j) { return 15*i+j;}
-  //
-public:
-  void AdddEdx(Int_t ip0,Double_t dedx0, Double_t dedxRef);
-  void Init();
-  TMatrixD * fState;           // state vector
-  TMatrixD * fCovariance;      // covariance
-  TMatrixD * fMat1;            // helper unit matrix
-  Int_t      fNpad;            // number of pad types
-  Int_t      fNpar;            // number of parameters
-  Int_t      fNelem;           // number of elements
-  //
-  // initial parameters estimate
-  //
-  Int_t      fSampleSize;       // size of  starting  sample
-  Int_t      fInit;             // number of initialized estimators
-  //
-  TVectorD   fSample[3];        // !training sample for robust estimate of initial parameters
-  TVectorD   fSampleStat[3];    // sample statistic
-  Int_t      fCounter[3];       // counter of samples
-  //
-private:
-  AliTPCkalmandEdx&  operator=(const AliTPCkalmandEdx&);// not implemented
-  ClassDef(AliTPCkalmandEdx,1);
-};
-
-#endif
-
index b1c1ae2..74031fe 100644 (file)
@@ -28,8 +28,8 @@
 #pragma link C++ class AliTPCcalibCosmic+;
 #pragma link C++ class AliTPCLaserTrack+;
 #pragma link C++ class AliTPCcalibTime+;
+#pragma link C++ class AliTPCcalibTimeGain+;
 #pragma link C++ class AliTPCcalibUnlinearity+;
-#pragma link C++ class AliTPCkalmandEdx+;
 #pragma link C++ class AliTPCcalibPID+;
 
 #endif
index 940c09f..9b9d01e 100644 (file)
@@ -5,7 +5,7 @@ SRCS = AliTPCcalibTracksCuts.cxx   AliTPCcalibTracks.cxx   AliTPCcalibTracksGain
        AliTPCcalibV0.cxx AliTPCCalibKr.cxx AliTPCcalibBase.cxx AliTPCAnalysisTaskcalib.cxx  \
        AliTPCCalibKrTask.cxx AliTPCcalibLaser.cxx AliTPCcalibCosmic.cxx AliTPCLaserTrack.cxx \
        AliTPCcalibCalib.cxx   AliTPCcalibTime.cxx AliTPCcalibUnlinearity.cxx \
-       AliTPCcalibPID.cxx AliTPCkalmandEdx.cxx
+       AliTPCcalibPID.cxx AliTPCcalibTimeGain.cxx
 
 HDRS:= $(SRCS:.cxx=.h)