--- /dev/null
+/**************************************************************************
+ * 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;
+
+
+}
--- /dev/null
+#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
+
+
+#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"
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
//
// 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
{
//
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
//
// 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
{
//
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;
}
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
//
// 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
{
//
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() {
//
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.
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;
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;
}
if (cal->fDFitter0T->GetNpoints()>0) fDFitter0T->Add(cal->fDFitter0T);
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++) {
// just for debugging, remove me
fTotalTracks += cal->fTotalTracks;
fAcceptedTracks += cal->fAcceptedTracks;
- fDebugCalPadRaw->Add(cal->fDebugCalPadRaw);
- fDebugCalPadCorr->Add(cal->fDebugCalPadCorr);
}
//
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
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);
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
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];
//
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;
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
"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" <<
}
-
-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
-
-*/
-
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,
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
//
//
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
//
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
+++ /dev/null
-/**************************************************************************
- * 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;
- }
-}
-
-
-
+++ /dev/null
-#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
-
#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
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)
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff