--- /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. *
+ **************************************************************************/
+
+
+
+/*
+ Responsible: marian.ivanov@cern.ch
+ Tools for fitting of the space point distortion parameters.
+ Functionality
+
+
+1. Creation of the distortion maps from the residual histograms
+2. Making fit trees
+3. Parameters to calculate- dO/dp
+
+ Some functions, for the moment function present in the AliTPCPreprocesorOffline, some will be
+ extracted from the old macros
+
+
+*/
+
+#include "Riostream.h"
+#include <fstream>
+#include "TMap.h"
+#include "TGraphErrors.h"
+#include "AliExternalTrackParam.h"
+#include "TFile.h"
+#include "TGraph.h"
+#include "TMultiGraph.h"
+#include "TCanvas.h"
+#include "THnSparse.h"
+#include "THnBase.h"
+#include "TProfile.h"
+#include "TROOT.h"
+#include "TLegend.h"
+#include "TPad.h"
+#include "TH2D.h"
+#include "TH3D.h"
+#include "AliTPCROC.h"
+#include "AliTPCCalROC.h"
+#include "AliESDfriend.h"
+#include "AliTPCcalibTime.h"
+#include "AliSplineFit.h"
+#include "AliCDBMetaData.h"
+#include "AliCDBId.h"
+#include "AliCDBManager.h"
+#include "AliCDBStorage.h"
+#include "AliTPCcalibBase.h"
+#include "AliTPCcalibDB.h"
+#include "AliTPCcalibDButil.h"
+#include "AliRelAlignerKalman.h"
+#include "AliTPCParamSR.h"
+#include "AliTPCcalibTimeGain.h"
+#include "AliTPCcalibGainMult.h"
+#include "AliSplineFit.h"
+#include "AliTPCComposedCorrection.h"
+#include "AliTPCExBTwist.h"
+#include "AliTPCCalibGlobalMisalignment.h"
+#include "TStatToolkit.h"
+#include "TChain.h"
+#include "TCut.h"
+#include "AliTrackerBase.h"
+#include "AliTPCCorrectionFit.h"
+#include "AliTPCLaserTrack.h"
+#include "TDatabasePDG.h"
+#include "AliTPCcalibAlign.h"
+
+ClassImp(AliTPCCorrectionFit)
+
+AliTPCCorrectionFit::AliTPCCorrectionFit():
+ TNamed("TPCCorrectionFit","TPCCorrectionFit")
+{
+ //
+ // default constructor
+ //
+}
+
+AliTPCCorrectionFit::~AliTPCCorrectionFit() {
+ //
+ // Destructor
+ //
+}
+
+
+Double_t AliTPCCorrectionFit::EvalAt(Double_t phi, Double_t refX, Double_t theta, Int_t corr, Int_t ptype){
+ //
+ //
+ //
+ Double_t sector = 9*phi/TMath::Pi();
+ if (sector<0) sector+=18;
+ Double_t y85=AliTPCCorrection::GetCorrSector(sector,85,theta,1,corr);
+ Double_t y245=AliTPCCorrection::GetCorrSector(sector,245,theta,1,corr);
+ if (ptype==0) return y85+(y245-y85)*(refX-85.)/(245.-85.);
+ if (ptype==2) return (y245-y85)/(245.-85.);
+ return 0;
+}
+
+
+
+
+Double_t AliTPCCorrectionFit::EvalAtPar(Double_t phi0, Double_t snp, Double_t refX, Double_t theta, Int_t corr, Int_t ptype, Int_t nsteps){
+ //
+ // Fit the distortion along the line with the parabolic model
+ // Parameters:
+ // phi0 - phi at the entrance of the TPC
+ // snp - local inclination angle at the entrance of the TPC
+ // refX - ref X where the distortion is evanluated
+ // theta
+ //
+ static TLinearFitter fitter(3,"pol2");
+ fitter.ClearPoints();
+ if (nsteps<3) nsteps=3;
+ Double_t deltaX=(245-85)/(nsteps);
+ for (Int_t istep=0; istep<(nsteps+1); istep++){
+ //
+ Double_t localX =85.+deltaX*istep;
+ Double_t localPhi=phi0+deltaX*snp*istep;
+ Double_t sector = 9*localPhi/TMath::Pi();
+ if (sector<0) sector+=18;
+ Double_t y=AliTPCCorrection::GetCorrSector(sector,localX,theta,1,corr);
+ Double_t dlocalX=AliTPCCorrection::GetCorrSector(sector,localX,theta,0,corr);
+ Double_t x[1]={localX-dlocalX};
+ fitter.AddPoint(x,y);
+ }
+ fitter.Eval();
+ Double_t par[3];
+ par[0]=fitter.GetParameter(0);
+ par[1]=fitter.GetParameter(1);
+ par[2]=fitter.GetParameter(2);
+
+ if (ptype==0) return par[0]+par[1]*refX+par[2]*refX*refX;
+ if (ptype==2) return par[1]+2*par[2]*refX;
+ if (ptype==4) return par[2];
+ return 0;
+}
+
+
+
+
+void AliTPCCorrectionFit::CreateAlignMaps(Double_t bz, Int_t run){
+ //
+ // Create cluster distortion map
+ //
+ TFile *falign = TFile::Open("CalibObjects.root");
+ TObjArray * arrayAlign = (TObjArray *)falign->Get("TPCAlign");
+ AliTPCcalibAlign * align = (AliTPCcalibAlign *)arrayAlign->FindObject("alignTPC");
+ TTreeSRedirector * pcstream = new TTreeSRedirector("TPCAlign.root");
+
+ THnBase * hdY = (THnBase*)align->GetClusterDelta(0);
+ //THnBase * hdZ = (THnBase*)align->GetClusterDelta(1);
+ AliTPCCorrectionFit::MakeClusterDistortionMap(hdY,pcstream,0, bz);
+ // AliTPCCorrectionFit::MakeClusterDistortionMap(hdZ,pcstream,1, bz);
+
+ const char * hname[5]={"dy","dz","dsnp","dtheta","d1pt"};
+ for (Int_t ihis=0; ihis<4; ihis++){
+ THnSparse * hisAlign =align->GetTrackletDelta(ihis);
+ AliTPCCorrectionFit::MakeDistortionMapSector(hisAlign, pcstream, hname[ihis], run, ihis,bz);
+ }
+ delete pcstream;
+ delete falign;
+}
+
+
+
+
+
+void AliTPCCorrectionFit::MakeClusterDistortionMap(THnBase * hisInput,TTreeSRedirector *pcstream , Int_t ptype, Int_t dtype){
+ //
+ // Make cluster residual map from the n-dimensional histogram
+ // hisInput supposed to have given format:
+ // 4 Dim:
+ // delta,
+ // sector
+ // localX
+ // kZ
+ // Vertex position assumed to be at (0,0,0)
+ //
+ //TTreeSRedirector *pcstream=new TTreeSRedirector(sname);
+ //
+ Int_t nbins1=hisInput->GetAxis(1)->GetNbins();
+ Int_t nbins2=hisInput->GetAxis(2)->GetNbins();
+ Int_t nbins3=hisInput->GetAxis(3)->GetNbins();
+ TF1 *fgaus=0;
+ TH3F * hisResMap3D =
+ new TH3F("his3D","his3D",
+ nbins1,hisInput->GetAxis(1)->GetXmin(), hisInput->GetAxis(1)->GetXmax(),
+ nbins2,hisInput->GetAxis(2)->GetXmin(), hisInput->GetAxis(2)->GetXmax(),
+ nbins3,hisInput->GetAxis(3)->GetXmin(), hisInput->GetAxis(3)->GetXmax());
+ hisResMap3D->GetXaxis()->SetTitle("sector");
+ hisResMap3D->GetYaxis()->SetTitle("localX");
+ hisResMap3D->GetZaxis()->SetTitle("kZ");
+
+ TH2F * hisResMap2D[4] ={0,0,0,0};
+ for (Int_t i=0; i<4; i++){
+ hisResMap2D[i]=
+ new TH2F(Form("his2D_0%d",i),Form("his2D_0%d",i),
+ nbins1,hisInput->GetAxis(1)->GetXmin(), hisInput->GetAxis(1)->GetXmax(),
+ nbins2,hisInput->GetAxis(2)->GetXmin(), hisInput->GetAxis(2)->GetXmax());
+ hisResMap2D[i]->GetXaxis()->SetTitle("sector");
+ hisResMap2D[i]->GetYaxis()->SetTitle("localX");
+ }
+ //
+ //
+ //
+ TF1 * f1= 0;
+ Int_t axis0[4]={0,1,2,3};
+ Int_t axis1[4]={0,1,2,3};
+ Int_t counter=0;
+ for (Int_t ibin1=1; ibin1<nbins1; ibin1+=1){
+ // phi- sector range
+ hisInput->GetAxis(1)->SetRange(ibin1-1,ibin1+1);
+ THnBase *his1=(THnBase *)hisInput->ProjectionND(4,axis0);
+ Double_t sector=hisInput->GetAxis(1)->GetBinCenter(ibin1);
+ //
+ for (Int_t ibin2=1; ibin2<nbins2; ibin2+=1){
+ // local x range
+ // kz fits
+ his1->GetAxis(2)->SetRange(ibin2-1,ibin2+1);
+ THnBase *his2=(THnBase *)his1->ProjectionND(4,axis1);
+ Double_t localX=hisInput->GetAxis(2)->GetBinCenter(ibin2);
+ //
+ //A side
+ his2->GetAxis(3)->SetRangeUser(0.01,0.3);
+ TH1 * hisA = his2->Projection(0);
+ Double_t meanA= hisA->GetMean();
+ Double_t rmsA= hisA->GetRMS();
+ Double_t entriesA= hisA->GetEntries();
+ delete hisA;
+ //C side
+ his2->GetAxis(3)->SetRangeUser(0.01,0.3);
+ TH1 * hisC = his2->Projection(0);
+ Double_t meanC= hisC->GetMean();
+ Double_t rmsC= hisC->GetRMS();
+ Double_t entriesC= hisC->GetEntries();
+ delete hisC;
+ his2->GetAxis(3)->SetRangeUser(-1.2,1.2);
+ TH2 * hisAC = his2->Projection(0,3);
+ TProfile *profAC = hisAC->ProfileX();
+ delete hisAC;
+ profAC->Fit("pol1","QNR","QNR",0.05,1);
+ if (!f1) f1=(TF1*)gROOT->FindObject("pol1");
+ Double_t offsetA=f1->GetParameter(0);
+ Double_t slopeA=f1->GetParameter(1);
+ Double_t offsetAE=f1->GetParError(0);
+ Double_t slopeAE=f1->GetParError(1);
+ Double_t chi2A=f1->GetChisquare()/f1->GetNumberFreeParameters();
+ profAC->Fit("pol1","QNR","QNR",-1.1,-0.1);
+ f1=(TF1*)gROOT->FindObject("pol1");
+ Double_t offsetC=f1->GetParameter(0);
+ Double_t slopeC=f1->GetParameter(1);
+ Double_t offsetCE=f1->GetParError(0);
+ Double_t slopeCE=f1->GetParError(1);
+ Double_t chi2C=f1->GetChisquare()/f1->GetNumberFreeParameters();
+ if (counter%50==0) printf("%f\t%f\t%f\t%f\t%f\t%f\t%f\n", sector,localX, entriesA+entriesC, slopeA,slopeC, chi2A, chi2C);
+ counter++;
+ (*pcstream)<<"deltaFit"<<
+ "sector="<<sector<<
+ "localX="<<localX<<
+ "meanA="<<meanA<<
+ "rmsA="<<rmsA<<
+ "entriesA="<<entriesA<<
+ "meanC="<<meanC<<
+ "rmsC="<<rmsC<<
+ "entriesC="<<entriesC<<
+ "offsetA="<<offsetA<<
+ "slopeA="<<slopeA<<
+ "offsetAE="<<offsetAE<<
+ "slopeAE="<<slopeAE<<
+ "chi2A="<<chi2A<<
+ "offsetC="<<offsetC<<
+ "slopeC="<<slopeC<<
+ "offsetCE="<<offsetCE<<
+ "slopeCE="<<slopeCE<<
+ "chi2C="<<chi2C<<
+ "\n";
+ //
+ hisResMap2D[0]->SetBinContent(ibin1,ibin2, offsetA);
+ hisResMap2D[1]->SetBinContent(ibin1,ibin2, slopeA);
+ hisResMap2D[2]->SetBinContent(ibin1,ibin2, offsetC);
+ hisResMap2D[3]->SetBinContent(ibin1,ibin2, slopeC);
+
+ for (Int_t ibin3=1; ibin3<nbins3; ibin3++){
+ Double_t kZ=hisInput->GetAxis(3)->GetBinCenter(ibin3);
+ if (TMath::Abs(kZ)<0.05) continue; // crossing
+ his2->GetAxis(3)->SetRange(ibin3,ibin3);
+ if (TMath::Abs(kZ)>0.15){
+ his2->GetAxis(3)->SetRange(ibin3,ibin3);
+ }
+ TH1 * his = his2->Projection(0);
+ Double_t mean= his->GetMean();
+ Double_t rms= his->GetRMS();
+ Double_t entries= his->GetEntries();
+ //printf("%f\t%f\t%f\t%f\t%f\t%f\n", sector,localX,kZ, entries, mean,rms);
+ hisResMap3D->SetBinContent(ibin1,ibin2,ibin3, mean);
+ Double_t phi=TMath::Pi()*sector/9;
+ if (phi>TMath::Pi()) phi+=TMath::Pi();
+ Double_t meanG=0;
+ Double_t rmsG=0;
+ if (entries>50){
+ if (!fgaus) {
+ his->Fit("gaus","Q","goff");
+ fgaus= (TF1*)((his->GetListOfFunctions()->FindObject("gaus"))->Clone());
+ }
+ if (fgaus) {
+ his->Fit(fgaus,"Q","goff");
+ meanG=fgaus->GetParameter(1);
+ rmsG=fgaus->GetParameter(2);
+ }
+ }
+ Double_t dsec=sector-Int_t(sector)-0.5;
+ Double_t snp=dsec*TMath::Pi()/9.;
+ (*pcstream)<<"delta"<<
+ "ptype="<<ptype<<
+ "dtype="<<dtype<<
+ "sector="<<sector<<
+ "dsec="<<dsec<<
+ "snp="<<snp<<
+ "phi="<<phi<<
+ "localX="<<localX<<
+ "kZ="<<kZ<<
+ "theta="<<kZ<<
+ "mean="<<mean<<
+ "rms="<<rms<<
+ "meanG="<<meanG<<
+ "rmsG="<<rmsG<<
+ "entries="<<entries<<
+ "meanA="<<meanA<<
+ "rmsA="<<rmsA<<
+ "entriesA="<<entriesA<<
+ "meanC="<<meanC<<
+ "rmsC="<<rmsC<<
+ "entriesC="<<entriesC<<
+ "offsetA="<<offsetA<<
+ "slopeA="<<slopeA<<
+ "chi2A="<<chi2A<<
+ "offsetC="<<offsetC<<
+ "slopeC="<<slopeC<<
+ "chi2C="<<chi2C<<
+ "\n";
+ delete his;
+ }
+ delete his2;
+ }
+ delete his1;
+ }
+ hisResMap3D->Write();
+ hisResMap2D[0]->Write();
+ hisResMap2D[1]->Write();
+ hisResMap2D[2]->Write();
+ hisResMap2D[3]->Write();
+ // delete pcstream;
+}
+
+
+
+void AliTPCCorrectionFit::MakeDistortionMap(THnSparse * his0, TTreeSRedirector * const pcstream, const char* hname, Int_t run, Float_t refX, Int_t type, Int_t integ, Double_t bz){
+ //
+ // make a distortion map out ou fthe residual histogram
+ // Results are written to the debug streamer - pcstream
+ // Parameters:
+ // his0 - input (4D) residual histogram
+ // pcstream - file to write the tree
+ // run - run number
+ // refX - track matching reference X
+ // type - 0- y 1-z,2 -snp, 3-theta, 4=1/pt
+ // THnSparse axes:
+ // OBJ: TAxis #Delta #Delta
+ // OBJ: TAxis tanTheta tan(#Theta)
+ // OBJ: TAxis phi #phi
+ // OBJ: TAxis snp snp
+
+ // marian.ivanov@cern.ch
+ const Int_t kMinEntries=10;
+ Int_t idim[4]={0,1,2,3};
+ //
+ //
+ //
+ Int_t nbins3=his0->GetAxis(3)->GetNbins();
+ Int_t first3=his0->GetAxis(3)->GetFirst();
+ Int_t last3 =his0->GetAxis(3)->GetLast();
+ //
+ for (Int_t ibin3=first3; ibin3<last3; ibin3+=1){ // axis 3 - local angle
+ his0->GetAxis(3)->SetRange(TMath::Max(ibin3-integ,1),TMath::Min(ibin3+integ,nbins3));
+ Double_t x3= his0->GetAxis(3)->GetBinCenter(ibin3);
+ THnSparse * his3= his0->Projection(3,idim); //projected histogram according selection 3
+ //
+ Int_t nbins2 = his3->GetAxis(2)->GetNbins();
+ Int_t first2 = his3->GetAxis(2)->GetFirst();
+ Int_t last2 = his3->GetAxis(2)->GetLast();
+ //
+ for (Int_t ibin2=first2; ibin2<last2; ibin2+=1){ // axis 2 - phi
+ his3->GetAxis(2)->SetRange(TMath::Max(ibin2-integ,1),TMath::Min(ibin2+integ,nbins2));
+ Double_t x2= his3->GetAxis(2)->GetBinCenter(ibin2);
+ THnSparse * his2= his3->Projection(2,idim); //projected histogram according selection 2
+ Int_t nbins1 = his2->GetAxis(1)->GetNbins();
+ Int_t first1 = his2->GetAxis(1)->GetFirst();
+ Int_t last1 = his2->GetAxis(1)->GetLast();
+ for (Int_t ibin1=first1; ibin1<last1; ibin1++){ //axis 1 - theta
+ //
+ Double_t x1= his2->GetAxis(1)->GetBinCenter(ibin1);
+ his2->GetAxis(1)->SetRange(TMath::Max(ibin1-1,1),TMath::Min(ibin1+1,nbins1));
+ if (TMath::Abs(x1)<0.1){
+ if (x1<0) his2->GetAxis(1)->SetRange(TMath::Max(ibin1-1,1),TMath::Min(ibin1,nbins1));
+ if (x1>0) his2->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1+1,nbins1));
+ }
+ if (TMath::Abs(x1)<0.06){
+ his2->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1,nbins1));
+ }
+ TH1 * hisDelta = his2->Projection(0);
+ //
+ Double_t entries = hisDelta->GetEntries();
+ Double_t mean=0, rms=0;
+ if (entries>kMinEntries){
+ mean = hisDelta->GetMean();
+ rms = hisDelta->GetRMS();
+ }
+ Double_t sector = 9.*x2/TMath::Pi();
+ if (sector<0) sector+=18;
+ Double_t dsec = sector-Int_t(sector)-0.5;
+ Double_t z=refX*x1;
+ (*pcstream)<<hname<<
+ "run="<<run<<
+ "bz="<<bz<<
+ "theta="<<x1<<
+ "phi="<<x2<<
+ "z="<<z<< // dummy z
+ "snp="<<x3<<
+ "entries="<<entries<<
+ "mean="<<mean<<
+ "rms="<<rms<<
+ "refX="<<refX<< // track matching refernce plane
+ "type="<<type<< //
+ "sector="<<sector<<
+ "dsec="<<dsec<<
+ "\n";
+ delete hisDelta;
+ //printf("%f\t%f\t%f\t%f\t%f\n",x3,x2,x1, entries,mean);
+ }
+ delete his2;
+ }
+ delete his3;
+ }
+}
+
+
+
+void AliTPCCorrectionFit::MakeDistortionMapCosmic(THnSparse * hisInput, TTreeSRedirector * const pcstream, const char* hname, Int_t run, Float_t refX, Int_t type){
+ //
+ // make a distortion map out ou fthe residual histogram
+ // Results are written to the debug streamer - pcstream
+ // Parameters:
+ // his0 - input (4D) residual histogram
+ // pcstream - file to write the tree
+ // run - run number
+ // refX - track matching reference X
+ // type - 0- y 1-z,2 -snp, 3-theta, 4=1/pt
+ // marian.ivanov@cern.ch
+ //
+ // Histo axeses
+ // Collection name='TObjArray', class='TObjArray', size=16
+ // 0. OBJ: TAxis #Delta #Delta
+ // 1. OBJ: TAxis N_{cl} N_{cl}
+ // 2. OBJ: TAxis dca_{r} (cm) dca_{r} (cm)
+ // 3. OBJ: TAxis z (cm) z (cm)
+ // 4. OBJ: TAxis sin(#phi) sin(#phi)
+ // 5. OBJ: TAxis tan(#theta) tan(#theta)
+ // 6. OBJ: TAxis 1/pt (1/GeV) 1/pt (1/GeV)
+ // 7. OBJ: TAxis pt (GeV) pt (GeV)
+ // 8. OBJ: TAxis alpha alpha
+ const Int_t kMinEntries=10;
+ //
+ // 1. make default selections
+ //
+ TH1 * hisDelta=0;
+ Int_t idim0[4]={0 , 5, 8, 3}; // delta, theta, alpha, z
+ hisInput->GetAxis(1)->SetRangeUser(110,190); //long tracks
+ hisInput->GetAxis(2)->SetRangeUser(-10,35); //tracks close to beam pipe
+ hisInput->GetAxis(4)->SetRangeUser(-0.3,0.3); //small snp at TPC entrance
+ hisInput->GetAxis(7)->SetRangeUser(3,100); //"high pt tracks"
+ hisDelta= hisInput->Projection(0);
+ hisInput->GetAxis(0)->SetRangeUser(-6.*hisDelta->GetRMS(), +6.*hisDelta->GetRMS());
+ delete hisDelta;
+ THnSparse *his0= hisInput->Projection(4,idim0);
+ //
+ // 2. Get mean in diferent bins
+ //
+ Int_t nbins1=his0->GetAxis(1)->GetNbins();
+ Int_t first1=his0->GetAxis(1)->GetFirst();
+ Int_t last1 =his0->GetAxis(1)->GetLast();
+ //
+ Double_t bz=AliTrackerBase::GetBz();
+ Int_t idim[4]={0,1, 2, 3}; // delta, theta,alpha,z
+ //
+ for (Int_t ibin1=first1; ibin1<=last1; ibin1++){ //axis 1 - theta
+ //
+ Double_t x1= his0->GetAxis(1)->GetBinCenter(ibin1);
+ his0->GetAxis(1)->SetRange(TMath::Max(ibin1-1,1),TMath::Min(ibin1+1,nbins1));
+ //
+ THnSparse * his1 = his0->Projection(4,idim); // projected histogram according range1
+ Int_t nbins3 = his1->GetAxis(3)->GetNbins();
+ Int_t first3 = his1->GetAxis(3)->GetFirst();
+ Int_t last3 = his1->GetAxis(3)->GetLast();
+ //
+ for (Int_t ibin3=first3-1; ibin3<=last3; ibin3+=1){ // axis 3 - z at "vertex"
+ his1->GetAxis(3)->SetRange(TMath::Max(ibin3-1,1),TMath::Min(ibin3+1,nbins3));
+ Double_t x3= his1->GetAxis(3)->GetBinCenter(ibin3);
+ if (ibin3<first3) {
+ his1->GetAxis(3)->SetRangeUser(-1,1);
+ x3=0;
+ }
+ THnSparse * his3= his1->Projection(4,idim); //projected histogram according selection 3
+ Int_t nbins2 = his3->GetAxis(2)->GetNbins();
+ Int_t first2 = his3->GetAxis(2)->GetFirst();
+ Int_t last2 = his3->GetAxis(2)->GetLast();
+ //
+ for (Int_t ibin2=first2; ibin2<=last2; ibin2+=1){
+ his3->GetAxis(2)->SetRange(TMath::Max(ibin2-1,1),TMath::Min(ibin2+1,nbins2));
+ Double_t x2= his3->GetAxis(2)->GetBinCenter(ibin2);
+ hisDelta = his3->Projection(0);
+ //
+ Double_t entries = hisDelta->GetEntries();
+ Double_t mean=0, rms=0;
+ if (entries>kMinEntries){
+ mean = hisDelta->GetMean();
+ rms = hisDelta->GetRMS();
+ }
+ Double_t sector = 9.*x2/TMath::Pi();
+ if (sector<0) sector+=18;
+ Double_t dsec = sector-Int_t(sector)-0.5;
+ Double_t snp=0; // dummy snp - equal 0
+ (*pcstream)<<hname<<
+ "run="<<run<<
+ "bz="<<bz<< // magnetic field
+ "theta="<<x1<< // theta
+ "phi="<<x2<< // phi (alpha)
+ "z="<<x3<< // z at "vertex"
+ "snp="<<snp<< // dummy snp
+ "entries="<<entries<< // entries in bin
+ "mean="<<mean<< // mean
+ "rms="<<rms<<
+ "refX="<<refX<< // track matching refernce plane
+ "type="<<type<< // parameter type
+ "sector="<<sector<< // sector
+ "dsec="<<dsec<< // dummy delta sector
+ "\n";
+ delete hisDelta;
+ printf("%f\t%f\t%f\t%f\t%f\n",x1,x3,x2, entries,mean);
+ }
+ delete his3;
+ }
+ delete his1;
+ }
+ delete his0;
+}
+
+
+
+void AliTPCCorrectionFit::MakeDistortionMapSector(THnSparse * hisInput, TTreeSRedirector * const pcstream, const char* hname, Int_t run, Int_t type, Double_t bz){
+ //
+ // make a distortion map out of the residual histogram
+ // Results are written to the debug streamer - pcstream
+ // Parameters:
+ // his0 - input (4D) residual histogram
+ // pcstream - file to write the tree
+ // run - run number
+ // type - 0- y 1-z,2 -snp, 3-theta
+ // bz - magnetic field
+ // marian.ivanov@cern.ch
+
+ //Collection name='TObjArray', class='TObjArray', size=16
+ //0 OBJ: TAxis delta delta
+ //1 OBJ: TAxis phi phi
+ //2 OBJ: TAxis localX localX
+ //3 OBJ: TAxis kY kY
+ //4 OBJ: TAxis kZ kZ
+ //5 OBJ: TAxis is1 is1
+ //6 OBJ: TAxis is0 is0
+ //7. OBJ: TAxis z z
+ //8. OBJ: TAxis IsPrimary IsPrimary
+
+ const Int_t kMinEntries=10;
+ THnSparse * hisSector0=0;
+ TH1 * htemp=0; // histogram to calculate mean value of parameter
+ // Double_t bz=AliTrackerBase::GetBz();
+
+ //
+ // Loop over pair of sector:
+ // isPrim - 8 ==> 8
+ // isec0 - 6 ==> 7
+ // isec1 - 5 ==> 6
+ // refX - 2 ==> 5
+ //
+ // phi - 1 ==> 4
+ // z - 7 ==> 3
+ // snp - 3 ==> 2
+ // theta- 4 ==> 1
+ // 0 ==> 0;
+ for (Int_t isec0=0; isec0<72; isec0++){
+ Int_t index0[9]={0, 4, 3, 7, 1, 2, 5, 6,8}; //regroup indeces
+ //
+ //hisInput->GetAxis(8)->SetRangeUser(-0.1,0.4); // select secondaries only ? - get out later ?
+ hisInput->GetAxis(6)->SetRangeUser(isec0-0.1,isec0+0.1);
+ hisSector0=hisInput->Projection(7,index0);
+ //
+ //
+ for (Int_t isec1=isec0+1; isec1<72; isec1++){
+ //if (isec1!=isec0+36) continue;
+ if ( TMath::Abs((isec0%18)-(isec1%18))>1.5 && TMath::Abs((isec0%18)-(isec1%18))<16.5) continue;
+ printf("Sectors %d\t%d\n",isec1,isec0);
+ hisSector0->GetAxis(6)->SetRangeUser(isec1-0.1,isec1+0.1);
+ TH1 * hisX=hisSector0->Projection(5);
+ Double_t refX= hisX->GetMean();
+ delete hisX;
+ TH1 *hisDelta=hisSector0->Projection(0);
+ Double_t dmean = hisDelta->GetMean();
+ Double_t drms = hisDelta->GetRMS();
+ hisSector0->GetAxis(0)->SetRangeUser(dmean-5.*drms, dmean+5.*drms);
+ delete hisDelta;
+ //
+ // 1. make default selections
+ //
+ Int_t idim0[5]={0 , 1, 2, 3, 4}; // {delta, theta, snp, z, phi }
+ THnBase *hisSector1= hisSector0->ProjectionND(5,idim0);
+ //
+ // 2. Get mean in diferent bins
+ //
+ Int_t idim[5]={0, 1, 2, 3, 4}; // {delta, theta-1,snp-2 ,z-3, phi-4}
+ //
+ // Int_t nbinsPhi=hisSector1->GetAxis(4)->GetNbins();
+ Int_t firstPhi=hisSector1->GetAxis(4)->GetFirst();
+ Int_t lastPhi =hisSector1->GetAxis(4)->GetLast();
+ //
+ for (Int_t ibinPhi=firstPhi; ibinPhi<=lastPhi; ibinPhi+=2){ //axis 4 - phi
+ //
+ // Phi loop
+ //
+ Double_t xPhi= hisSector1->GetAxis(4)->GetBinCenter(ibinPhi);
+ Double_t psec = (9*xPhi/TMath::Pi());
+ if (psec<0) psec+=18;
+ Bool_t isOK0=kFALSE;
+ Bool_t isOK1=kFALSE;
+ if (TMath::Abs(psec-isec0%18-0.5)<1. || TMath::Abs(psec-isec0%18-17.5)<1.) isOK0=kTRUE;
+ if (TMath::Abs(psec-isec1%18-0.5)<1. || TMath::Abs(psec-isec1%18-17.5)<1.) isOK1=kTRUE;
+ if (!isOK0) continue;
+ if (!isOK1) continue;
+ //
+ hisSector1->GetAxis(4)->SetRange(TMath::Max(ibinPhi-2,firstPhi),TMath::Min(ibinPhi+2,lastPhi));
+ if (isec1!=isec0+36) {
+ hisSector1->GetAxis(4)->SetRange(TMath::Max(ibinPhi-3,firstPhi),TMath::Min(ibinPhi+3,lastPhi));
+ }
+ //
+ htemp = hisSector1->Projection(4);
+ xPhi=htemp->GetMean();
+ delete htemp;
+ THnBase * hisPhi = hisSector1->ProjectionND(4,idim);
+ //Int_t nbinsZ = hisPhi->GetAxis(3)->GetNbins();
+ Int_t firstZ = hisPhi->GetAxis(3)->GetFirst();
+ Int_t lastZ = hisPhi->GetAxis(3)->GetLast();
+ //
+ for (Int_t ibinZ=firstZ; ibinZ<=lastZ; ibinZ+=2){ // axis 3 - z
+ //
+ // Z loop
+ //
+ hisPhi->GetAxis(3)->SetRange(TMath::Max(ibinZ,firstZ),TMath::Min(ibinZ,lastZ));
+ if (isec1!=isec0+36) {
+ hisPhi->GetAxis(3)->SetRange(TMath::Max(ibinZ-1,firstZ),TMath::Min(ibinZ-1,lastZ));
+ }
+ htemp = hisPhi->Projection(3);
+ Double_t xZ= htemp->GetMean();
+ delete htemp;
+ THnBase * hisZ= hisPhi->ProjectionND(3,idim);
+ //projected histogram according selection 3 -z
+ //
+ //
+ //Int_t nbinsSnp = hisZ->GetAxis(2)->GetNbins();
+ Int_t firstSnp = hisZ->GetAxis(2)->GetFirst();
+ Int_t lastSnp = hisZ->GetAxis(2)->GetLast();
+ for (Int_t ibinSnp=firstSnp; ibinSnp<=lastSnp; ibinSnp+=2){ // axis 2 - snp
+ //
+ // Snp loop
+ //
+ hisZ->GetAxis(2)->SetRange(TMath::Max(ibinSnp-1,firstSnp),TMath::Min(ibinSnp+1,lastSnp));
+ if (isec1!=isec0+36) {
+ hisZ->GetAxis(2)->SetRange(TMath::Max(ibinSnp-2,firstSnp),TMath::Min(ibinSnp+2,lastSnp));
+ }
+ htemp = hisZ->Projection(2);
+ Double_t xSnp= htemp->GetMean();
+ delete htemp;
+ THnBase * hisSnp= hisZ->ProjectionND(2,idim);
+ //projected histogram according selection 2 - snp
+
+ //Int_t nbinsTheta = hisSnp->GetAxis(1)->GetNbins();
+ Int_t firstTheta = hisSnp->GetAxis(1)->GetFirst();
+ Int_t lastTheta = hisSnp->GetAxis(1)->GetLast();
+ //
+ for (Int_t ibinTheta=firstTheta; ibinTheta<=lastTheta; ibinTheta+=2){ // axis1 theta
+
+
+ hisSnp->GetAxis(1)->SetRange(TMath::Max(ibinTheta-2,firstTheta),TMath::Min(ibinTheta+2,lastTheta));
+ if (isec1!=isec0+36) {
+ hisSnp->GetAxis(1)->SetRange(TMath::Max(ibinTheta-3,firstTheta),TMath::Min(ibinTheta+3,lastTheta));
+ }
+ htemp = hisSnp->Projection(1);
+ Double_t xTheta=htemp->GetMean();
+ delete htemp;
+ hisDelta = hisSnp->Projection(0);
+ //
+ Double_t entries = hisDelta->GetEntries();
+ Double_t mean=0, rms=0;
+ if (entries>kMinEntries){
+ mean = hisDelta->GetMean();
+ rms = hisDelta->GetRMS();
+ }
+ Double_t sector = 9.*xPhi/TMath::Pi();
+ if (sector<0) sector+=18;
+ Double_t dsec = sector-Int_t(sector)-0.5;
+ Int_t dtype=1; // TPC alignment type
+ (*pcstream)<<hname<<
+ "run="<<run<<
+ "bz="<<bz<< // magnetic field
+ "ptype="<<type<< // parameter type
+ "dtype="<<dtype<< // parameter type
+ "isec0="<<isec0<< // sector 0
+ "isec1="<<isec1<< // sector 1
+ "sector="<<sector<< // sector as float
+ "dsec="<<dsec<< // delta sector
+ //
+ "theta="<<xTheta<< // theta
+ "phi="<<xPhi<< // phi (alpha)
+ "z="<<xZ<< // z
+ "snp="<<xSnp<< // snp
+
+ //
+ "entries="<<entries<< // entries in bin
+ "mean="<<mean<< // mean
+ "rms="<<rms<< // rms
+ "refX="<<refX<< // track matching reference plane
+ "\n";
+ delete hisDelta;
+ //printf("%d\t%d\t%f\t%f\t%f\t%f\t%f\t%f\n",isec0, isec1, xPhi,xZ,xSnp, xTheta, entries,mean);
+ //
+ }//ibinTheta
+ delete hisSnp;
+ } //ibinSnp
+ delete hisZ;
+ }//ibinZ
+ delete hisPhi;
+ }//ibinPhi
+ delete hisSector1;
+ }//isec1
+ delete hisSector0;
+ }//isec0
+}
+
+
+
+
+
+// void AliTPCCorrectionFit::MakeLaserDistortionTree(TTree* tree, TObjArray */*corrArray*/, Int_t /*itype*/){
+// //
+// // Make a laser fit tree for global minimization
+// //
+// AliTPCcalibDB* calib=AliTPCcalibDB::Instance();
+// AliTPCCorrection * correction = calib->GetTPCComposedCorrection();
+// if (!correction) correction = calib->GetTPCComposedCorrection(AliTrackerBase::GetBz());
+// correction->AddVisualCorrection(correction,0); //register correction
+
+// // AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
+// //AliTPCParam *param = AliTPCcalibDB::Instance()->GetParameters();
+// //
+// const Double_t cutErrY=0.05;
+// const Double_t kSigmaCut=4;
+// // const Double_t cutErrZ=0.03;
+// const Double_t kEpsilon=0.00000001;
+// // const Double_t kMaxDist=1.; // max distance - space correction
+// TVectorD *vecdY=0;
+// TVectorD *vecdZ=0;
+// TVectorD *veceY=0;
+// TVectorD *veceZ=0;
+// AliTPCLaserTrack *ltr=0;
+// AliTPCLaserTrack::LoadTracks();
+// tree->SetBranchAddress("dY.",&vecdY);
+// tree->SetBranchAddress("dZ.",&vecdZ);
+// tree->SetBranchAddress("eY.",&veceY);
+// tree->SetBranchAddress("eZ.",&veceZ);
+// tree->SetBranchAddress("LTr.",<r);
+// Int_t entries= tree->GetEntries();
+// TTreeSRedirector *pcstream= new TTreeSRedirector("distortionLaser_0.root");
+// Double_t bz=AliTrackerBase::GetBz();
+// //
+// // Double_t globalXYZ[3];
+// //Double_t globalXYZCorr[3];
+// for (Int_t ientry=0; ientry<entries; ientry++){
+// tree->GetEntry(ientry);
+// if (!ltr->GetVecGX()){
+// ltr->UpdatePoints();
+// }
+// //
+// TVectorD fit10(5);
+// TVectorD fit5(5);
+// printf("Entry\t%d\n",ientry);
+// for (Int_t irow0=0; irow0<158; irow0+=1){
+// //
+// TLinearFitter fitter10(4,"hyp3");
+// TLinearFitter fitter5(2,"hyp1");
+// Int_t sector= (Int_t)(*ltr->GetVecSec())[irow0];
+// if (sector<0) continue;
+// //if (TMath::Abs(vecdY->GetMatrixArray()[irow0])<kEpsilon) continue;
+
+// Double_t refX= (*ltr->GetVecLX())[irow0];
+// Int_t firstRow1 = TMath::Max(irow0-10,0);
+// Int_t lastRow1 = TMath::Min(irow0+10,158);
+// Double_t padWidth=(irow0<64)?0.4:0.6;
+// // make long range fit
+// for (Int_t irow1=firstRow1; irow1<=lastRow1; irow1++){
+// if (TMath::Abs((*ltr->GetVecSec())[irow1]-sector)>kEpsilon) continue;
+// if (veceY->GetMatrixArray()[irow1]>cutErrY) continue;
+// if (TMath::Abs(vecdY->GetMatrixArray()[irow1])<kEpsilon) continue;
+// Double_t idealX= (*ltr->GetVecLX())[irow1];
+// Double_t idealY= (*ltr->GetVecLY())[irow1];
+// // Double_t idealZ= (*ltr->GetVecLZ())[irow1];
+// Double_t gx= (*ltr->GetVecGX())[irow1];
+// Double_t gy= (*ltr->GetVecGY())[irow1];
+// Double_t gz= (*ltr->GetVecGZ())[irow1];
+// Double_t measY=(*vecdY)[irow1]+idealY;
+// Double_t deltaR = GetCorrXYZ(gx, gy, gz, 0,0);
+// // deltaR = R distorted -R ideal
+// Double_t xxx[4]={idealX+deltaR-refX,TMath::Cos(idealY/padWidth), TMath::Sin(idealY/padWidth)};
+// fitter10.AddPoint(xxx,measY,1);
+// }
+// Bool_t isOK=kTRUE;
+// Double_t rms10=0;//TMath::Sqrt(fitter10.GetChisquare()/(fitter10.GetNpoints()-4));
+// Double_t mean10 =0;// fitter10.GetParameter(0);
+// Double_t slope10 =0;// fitter10.GetParameter(0);
+// Double_t cosPart10 = 0;// fitter10.GetParameter(2);
+// Double_t sinPart10 =0;// fitter10.GetParameter(3);
+
+// if (fitter10.GetNpoints()>10){
+// fitter10.Eval();
+// rms10=TMath::Sqrt(fitter10.GetChisquare()/(fitter10.GetNpoints()-4));
+// mean10 = fitter10.GetParameter(0);
+// slope10 = fitter10.GetParameter(1);
+// cosPart10 = fitter10.GetParameter(2);
+// sinPart10 = fitter10.GetParameter(3);
+// //
+// // make short range fit
+// //
+// for (Int_t irow1=firstRow1+5; irow1<=lastRow1-5; irow1++){
+// if (TMath::Abs((*ltr->GetVecSec())[irow1]-sector)>kEpsilon) continue;
+// if (veceY->GetMatrixArray()[irow1]>cutErrY) continue;
+// if (TMath::Abs(vecdY->GetMatrixArray()[irow1])<kEpsilon) continue;
+// Double_t idealX= (*ltr->GetVecLX())[irow1];
+// Double_t idealY= (*ltr->GetVecLY())[irow1];
+// // Double_t idealZ= (*ltr->GetVecLZ())[irow1];
+// Double_t gx= (*ltr->GetVecGX())[irow1];
+// Double_t gy= (*ltr->GetVecGY())[irow1];
+// Double_t gz= (*ltr->GetVecGZ())[irow1];
+// Double_t measY=(*vecdY)[irow1]+idealY;
+// Double_t deltaR = GetCorrXYZ(gx, gy, gz, 0,0);
+// // deltaR = R distorted -R ideal
+// Double_t expY= mean10+slope10*(idealX+deltaR-refX);
+// if (TMath::Abs(measY-expY)>kSigmaCut*rms10) continue;
+// //
+// Double_t corr=cosPart10*TMath::Cos(idealY/padWidth)+sinPart10*TMath::Sin(idealY/padWidth);
+// Double_t xxx[4]={idealX+deltaR-refX,TMath::Cos(idealY/padWidth), TMath::Sin(idealY/padWidth)};
+// fitter5.AddPoint(xxx,measY-corr,1);
+// }
+// }else{
+// isOK=kFALSE;
+// }
+// if (fitter5.GetNpoints()<8) isOK=kFALSE;
+
+// Double_t rms5=0;//TMath::Sqrt(fitter5.GetChisquare()/(fitter5.GetNpoints()-4));
+// Double_t offset5 =0;// fitter5.GetParameter(0);
+// Double_t slope5 =0;// fitter5.GetParameter(0);
+// if (isOK){
+// fitter5.Eval();
+// rms5=TMath::Sqrt(fitter5.GetChisquare()/(fitter5.GetNpoints()-4));
+// offset5 = fitter5.GetParameter(0);
+// slope5 = fitter5.GetParameter(0);
+// }
+// //
+// Double_t dtype=5;
+// Double_t ptype=0;
+// Double_t phi =(*ltr->GetVecPhi())[irow0];
+// Double_t theta =ltr->GetTgl();
+// Double_t mean=(vecdY)->GetMatrixArray()[irow0];
+// Double_t gx=0,gy=0,gz=0;
+// Double_t snp = (*ltr->GetVecP2())[irow0];
+// Int_t bundle= ltr->GetBundle();
+// Int_t id= ltr->GetId();
+// // Double_t rms = err->GetMatrixArray()[irow];
+// //
+// gx = (*ltr->GetVecGX())[irow0];
+// gy = (*ltr->GetVecGY())[irow0];
+// gz = (*ltr->GetVecGZ())[irow0];
+// Double_t dRrec = GetCorrXYZ(gx, gy, gz, 0,0);
+// fitter10.GetParameters(fit10);
+// fitter5.GetParameters(fit5);
+// Double_t idealY= (*ltr->GetVecLY())[irow0];
+// Double_t measY=(*vecdY)[irow0]+idealY;
+// Double_t corr=cosPart10*TMath::Cos(idealY/padWidth)+sinPart10*TMath::Sin(idealY/padWidth);
+// if (TMath::Max(rms5,rms10)>0.06) isOK=kFALSE;
+// //
+// (*pcstream)<<"fitFull"<< // dumpe also intermediate results
+// "bz="<<bz<< // magnetic filed used
+// "dtype="<<dtype<< // detector match type
+// "ptype="<<ptype<< // parameter type
+// "theta="<<theta<< // theta
+// "phi="<<phi<< // phi
+// "snp="<<snp<< // snp
+// "sector="<<sector<<
+// "bundle="<<bundle<<
+// // // "dsec="<<dsec<<
+// "refX="<<refX<< // reference radius
+// "gx="<<gx<< // global position
+// "gy="<<gy<< // global position
+// "gz="<<gz<< // global position
+// "dRrec="<<dRrec<< // delta Radius in reconstruction
+// "id="<<id<< //bundle
+// "rms10="<<rms10<<
+// "rms5="<<rms5<<
+// "fit10.="<<&fit10<<
+// "fit5.="<<&fit5<<
+// "measY="<<measY<<
+// "mean="<<mean<<
+// "idealY="<<idealY<<
+// "corr="<<corr<<
+// "isOK="<<isOK<<
+// "\n";
+// }
+// }
+// delete pcstream;
+// }
+
+
+void AliTPCCorrectionFit::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Int_t offset, Bool_t debug ){
+ //
+ // Make a fit tree:
+ // For each partial correction (specified in array) and given track topology (phi, theta, snp, refX)
+ // calculates partial distortions
+ // Partial distortion is stored in the resulting tree
+ // Output is storred in the file distortion_<dettype>_<partype>.root
+ // Partial distortion is stored with the name given by correction name
+ //
+ //
+ // Parameters of function:
+ // input - input tree
+ // dtype - distortion type 0 - ITSTPC, 1 -TPCTRD, 2 - TPCvertex , 3 - TPC-TOF, 4 - TPCTPC track crossing
+ // ppype - parameter type
+ // corrArray - array with partial corrections
+ // step - skipe entries - if 1 all entries processed - it is slow
+ // debug 0 if debug on also space points dumped - it is slow
+
+ const Double_t kMaxSnp = 0.85;
+ const Double_t kcutSnp=0.25;
+ const Double_t kcutTheta=1.;
+ const Double_t kRadiusTPC=85;
+ // AliTPCROC *tpcRoc =AliTPCROC::Instance();
+ //
+ const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
+ // const Double_t kB2C=-0.299792458e-3;
+ const Int_t kMinEntries=20;
+ Double_t phi,theta, snp, mean,rms, entries,sector,dsec;
+ Float_t refX;
+ Int_t run;
+ tinput->SetBranchAddress("run",&run);
+ tinput->SetBranchAddress("theta",&theta);
+ tinput->SetBranchAddress("phi", &phi);
+ tinput->SetBranchAddress("snp",&snp);
+ tinput->SetBranchAddress("mean",&mean);
+ tinput->SetBranchAddress("rms",&rms);
+ tinput->SetBranchAddress("entries",&entries);
+ tinput->SetBranchAddress("sector",§or);
+ tinput->SetBranchAddress("dsec",&dsec);
+ tinput->SetBranchAddress("refX",&refX);
+ TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortion%d_%d_%d.root",dtype,ptype,offset));
+ //
+ Int_t nentries=tinput->GetEntries();
+ Int_t ncorr=corrArray->GetEntries();
+ Double_t corrections[100]={0}; //
+ Double_t tPar[5];
+ Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
+ Int_t dir=0;
+ if (dtype==5 || dtype==6) dtype=4;
+ if (dtype==0) { dir=-1;}
+ if (dtype==1) { dir=1;}
+ if (dtype==2) { dir=-1;}
+ if (dtype==3) { dir=1;}
+ if (dtype==4) { dir=-1;}
+ //
+ for (Int_t ientry=offset; ientry<nentries; ientry+=step){
+ tinput->GetEntry(ientry);
+ if (TMath::Abs(snp)>kMaxSnp) continue;
+ tPar[0]=0;
+ tPar[1]=theta*refX;
+ if (dtype==2) tPar[1]=theta*kRadiusTPC;
+ tPar[2]=snp;
+ tPar[3]=theta;
+ tPar[4]=(gRandom->Rndm()-0.5)*0.02; // should be calculated - non equal to 0
+ if (dtype==4){
+ // tracks crossing CE
+ tPar[1]=0; // track at the CE
+ //if (TMath::Abs(theta) <0.05) continue; // deep cross
+ }
+
+ if (TMath::Abs(snp) >kcutSnp) continue;
+ if (TMath::Abs(theta) >kcutTheta) continue;
+ printf("%f\t%f\t%f\t%f\t%f\t%f\n",entries, sector,theta,snp, mean,rms);
+ Double_t bz=AliTrackerBase::GetBz();
+ if (dtype !=4) { //exclude TPC - for TPC mainly non primary tracks
+ if (dtype!=2 && TMath::Abs(bz)>0.1 ) tPar[4]=snp/(refX*bz*kB2C*2);
+
+ if (dtype==2 && TMath::Abs(bz)>0.1 ) {
+ tPar[4]=snp/(kRadiusTPC*bz*kB2C*2);//
+ // snp at the TPC inner radius in case the vertex match used
+ }
+ }
+ //
+ tPar[4]+=(gRandom->Rndm()-0.5)*0.02;
+ AliExternalTrackParam track(refX,phi,tPar,cov);
+ Double_t xyz[3];
+ track.GetXYZ(xyz);
+ Int_t id=0;
+ Double_t pt=1./tPar[4];
+ Double_t dRrec=0; // dummy value - needed for points - e.g for laser
+ //if (ptype==4 &&bz<0) mean*=-1; // interpret as curvature -- COMMENTED out - in lookup signed 1/pt used
+ Double_t refXD=refX;
+ (*pcstream)<<"fit"<<
+ "run="<<run<< // run number
+ "bz="<<bz<< // magnetic filed used
+ "dtype="<<dtype<< // detector match type
+ "ptype="<<ptype<< // parameter type
+ "theta="<<theta<< // theta
+ "phi="<<phi<< // phi
+ "snp="<<snp<< // snp
+ "mean="<<mean<< // mean dist value
+ "rms="<<rms<< // rms
+ "sector="<<sector<<
+ "dsec="<<dsec<<
+ "refX="<<refXD<< // referece X as double
+ "gx="<<xyz[0]<< // global position at reference
+ "gy="<<xyz[1]<< // global position at reference
+ "gz="<<xyz[2]<< // global position at reference
+ "dRrec="<<dRrec<< // delta Radius in reconstruction
+ "pt="<<pt<< // pt
+ "id="<<id<< // track id
+ "entries="<<entries;// number of entries in bin
+ //
+ Bool_t isOK=kTRUE;
+ if (entries<kMinEntries) isOK=kFALSE;
+ //
+ if (dtype!=4) for (Int_t icorr=0; icorr<ncorr; icorr++) {
+ AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
+ corrections[icorr]=0;
+ if (entries>kMinEntries){
+ AliExternalTrackParam trackIn(refX,phi,tPar,cov);
+ AliExternalTrackParam *trackOut = 0;
+ if (debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,pcstream);
+ if (!debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,0);
+ if (dtype==0) {dir= -1;}
+ if (dtype==1) {dir= 1;}
+ if (dtype==2) {dir= -1;}
+ if (dtype==3) {dir= 1;}
+ //
+ if (trackOut){
+ if (!AliTrackerBase::PropagateTrackTo(&trackIn,refX,kMass,5,kTRUE,kMaxSnp)) isOK=kFALSE;
+ if (!trackOut->Rotate(trackIn.GetAlpha())) isOK=kFALSE;
+ if (!AliTrackerBase::PropagateTrackTo(trackOut,trackIn.GetX(),kMass,5,kFALSE,kMaxSnp)) isOK=kFALSE;
+ // trackOut->PropagateTo(trackIn.GetX(),AliTrackerBase::GetBz());
+ //
+ corrections[icorr]= trackOut->GetParameter()[ptype]-trackIn.GetParameter()[ptype];
+ delete trackOut;
+ }else{
+ corrections[icorr]=0;
+ isOK=kFALSE;
+ }
+ //if (ptype==4 &&bz<0) corrections[icorr]*=-1; // interpret as curvature - commented out
+ }
+ (*pcstream)<<"fit"<<
+ Form("%s=",corr->GetName())<<corrections[icorr]; // dump correction value
+ }
+
+ if (dtype==4) for (Int_t icorr=0; icorr<ncorr; icorr++) {
+ //
+ // special case of the TPC tracks crossing the CE
+ //
+ AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
+ corrections[icorr]=0;
+ if (entries>kMinEntries){
+ AliExternalTrackParam trackIn0(refX,phi,tPar,cov); //Outer - direction to vertex
+ AliExternalTrackParam trackIn1(refX,phi,tPar,cov); //Inner - direction magnet
+ AliExternalTrackParam *trackOut0 = 0;
+ AliExternalTrackParam *trackOut1 = 0;
+ //
+ if (debug) trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,pcstream);
+ if (!debug) trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,0);
+ if (debug) trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,pcstream);
+ if (!debug) trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,0);
+ //
+ if (trackOut0 && trackOut1){
+ if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,5,kTRUE,kMaxSnp)) isOK=kFALSE;
+ if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
+ if (!trackOut0->Rotate(trackIn0.GetAlpha())) isOK=kFALSE;
+ if (!AliTrackerBase::PropagateTrackTo(trackOut0,trackIn0.GetX(),kMass,5,kFALSE,kMaxSnp)) isOK=kFALSE;
+ //
+ if (!AliTrackerBase::PropagateTrackTo(&trackIn1,refX,kMass,5,kTRUE,kMaxSnp)) isOK=kFALSE;
+ if (!trackIn1.Rotate(trackIn0.GetAlpha())) isOK=kFALSE;
+ if (!AliTrackerBase::PropagateTrackTo(&trackIn1,trackIn0.GetX(),kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
+ if (!trackOut1->Rotate(trackIn1.GetAlpha())) isOK=kFALSE;
+ if (!AliTrackerBase::PropagateTrackTo(trackOut1,trackIn1.GetX(),kMass,5,kFALSE,kMaxSnp)) isOK=kFALSE;
+ //
+ corrections[icorr] = (trackOut0->GetParameter()[ptype]-trackIn0.GetParameter()[ptype]);
+ corrections[icorr]-= (trackOut1->GetParameter()[ptype]-trackIn1.GetParameter()[ptype]);
+ if (isOK)
+ if ((TMath::Abs(trackOut0->GetX()-trackOut1->GetX())>0.1)||
+ (TMath::Abs(trackOut0->GetX()-trackIn1.GetX())>0.1)||
+ (TMath::Abs(trackOut0->GetAlpha()-trackOut1->GetAlpha())>0.00001)||
+ (TMath::Abs(trackOut0->GetAlpha()-trackIn1.GetAlpha())>0.00001)||
+ (TMath::Abs(trackIn0.GetTgl()-trackIn1.GetTgl())>0.0001)||
+ (TMath::Abs(trackIn0.GetSnp()-trackIn1.GetSnp())>0.0001)
+ ){
+ isOK=kFALSE;
+ }
+ delete trackOut0;
+ delete trackOut1;
+ }else{
+ corrections[icorr]=0;
+ isOK=kFALSE;
+ }
+ //
+ //if (ptype==4 &&bz<0) corrections[icorr]*=-1; // interpret as curvature - commented out no in lookup
+ }
+ (*pcstream)<<"fit"<<
+ Form("%s=",corr->GetName())<<corrections[icorr]; // dump correction value
+ }
+ //
+ (*pcstream)<<"fit"<<"isOK="<<isOK<<"\n";
+ }
+
+
+ delete pcstream;
+}
+
+
+
+void AliTPCCorrectionFit::MakeSectorDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Int_t offset, Bool_t debug ){
+ //
+ // Make a fit tree:
+ // For each partial correction (specified in array) and given track topology (phi, theta, snp, refX)
+ // calculates partial distortions
+ // Partial distortion is stored in the resulting tree
+ // Output is storred in the file distortion_<dettype>_<partype>.root
+ // Partial distortion is stored with the name given by correction name
+ //
+ //
+ // Parameters of function:
+ // input - input tree
+ // dtype - distortion type 10 - IROC-OROC
+ // ppype - parameter type
+ // corrArray - array with partial corrections
+ // step - skipe entries - if 1 all entries processed - it is slow
+ // debug 0 if debug on also space points dumped - it is slow
+
+ const Double_t kMaxSnp = 0.8;
+ const Int_t kMinEntries=200;
+ // AliTPCROC *tpcRoc =AliTPCROC::Instance();
+ //
+ const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
+ // const Double_t kB2C=-0.299792458e-3;
+ Double_t phi,theta, snp, mean,rms, entries,sector,dsec,globalZ;
+ Int_t isec1, isec0;
+ Double_t refXD;
+ Float_t refX;
+ Int_t run;
+ tinput->SetBranchAddress("run",&run);
+ tinput->SetBranchAddress("theta",&theta);
+ tinput->SetBranchAddress("phi", &phi);
+ tinput->SetBranchAddress("snp",&snp);
+ tinput->SetBranchAddress("mean",&mean);
+ tinput->SetBranchAddress("rms",&rms);
+ tinput->SetBranchAddress("entries",&entries);
+ tinput->SetBranchAddress("sector",§or);
+ tinput->SetBranchAddress("dsec",&dsec);
+ tinput->SetBranchAddress("refX",&refXD);
+ tinput->SetBranchAddress("z",&globalZ);
+ tinput->SetBranchAddress("isec0",&isec0);
+ tinput->SetBranchAddress("isec1",&isec1);
+ TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortionSector%d_%d_%d.root",dtype,ptype,offset));
+ //
+ Int_t nentries=tinput->GetEntries();
+ Int_t ncorr=corrArray->GetEntries();
+ Double_t corrections[100]={0}; //
+ Double_t tPar[5];
+ Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
+ Int_t dir=0;
+ //
+ for (Int_t ientry=offset; ientry<nentries; ientry+=step){
+ tinput->GetEntry(ientry);
+ refX=refXD;
+ Int_t id=-1;
+ if (TMath::Abs(TMath::Abs(isec0%18)-TMath::Abs(isec1%18))==0) id=1; // IROC-OROC - opposite side
+ if (TMath::Abs(TMath::Abs(isec0%36)-TMath::Abs(isec1%36))==0) id=2; // IROC-OROC - same side
+ if (dtype==10 && id==-1) continue;
+ //
+ dir=-1;
+ tPar[0]=0;
+ tPar[1]=globalZ;
+ tPar[2]=snp;
+ tPar[3]=theta;
+ tPar[4]=(gRandom->Rndm()-0.1)*0.2; //
+ Double_t pt=1./tPar[4];
+ //
+ printf("%f\t%f\t%f\t%f\t%f\t%f\n",entries, sector,theta,snp, mean,rms);
+ Double_t bz=AliTrackerBase::GetBz();
+ AliExternalTrackParam track(refX,phi,tPar,cov);
+ Double_t xyz[3],xyzIn[3],xyzOut[3];
+ track.GetXYZ(xyz);
+ track.GetXYZAt(85,bz,xyzIn);
+ track.GetXYZAt(245,bz,xyzOut);
+ Double_t phiIn = TMath::ATan2(xyzIn[1],xyzIn[0]);
+ Double_t phiOut = TMath::ATan2(xyzOut[1],xyzOut[0]);
+ Double_t phiRef = TMath::ATan2(xyz[1],xyz[0]);
+ Int_t sectorRef = TMath::Nint(9.*phiRef/TMath::Pi()-0.5);
+ Int_t sectorIn = TMath::Nint(9.*phiIn/TMath::Pi()-0.5);
+ Int_t sectorOut = TMath::Nint(9.*phiOut/TMath::Pi()-0.5);
+ //
+ Bool_t isOK=kTRUE;
+ if (sectorIn!=sectorOut) isOK=kFALSE; // requironment - cluster in the same sector
+ if (sectorIn!=sectorRef) isOK=kFALSE; // requironment - cluster in the same sector
+ if (entries<kMinEntries/(1+TMath::Abs(globalZ/100.))) isOK=kFALSE; // requironment - minimal amount of tracks in bin
+ // Do downscale
+ if (TMath::Abs(theta)>1) isOK=kFALSE;
+ //
+ Double_t dRrec=0; // dummy value - needed for points - e.g for laser
+ //
+ (*pcstream)<<"fit"<<
+ "run="<<run<< //run
+ "bz="<<bz<< // magnetic filed used
+ "dtype="<<dtype<< // detector match type
+ "ptype="<<ptype<< // parameter type
+ "theta="<<theta<< // theta
+ "phi="<<phi<< // phi
+ "snp="<<snp<< // snp
+ "mean="<<mean<< // mean dist value
+ "rms="<<rms<< // rms
+ "sector="<<sector<<
+ "dsec="<<dsec<<
+ "refX="<<refXD<< // referece X
+ "gx="<<xyz[0]<< // global position at reference
+ "gy="<<xyz[1]<< // global position at reference
+ "gz="<<xyz[2]<< // global position at reference
+ "dRrec="<<dRrec<< // delta Radius in reconstruction
+ "pt="<<pt<< //pt
+ "id="<<id<< // track id
+ "entries="<<entries;// number of entries in bin
+ //
+ AliExternalTrackParam *trackOut0 = 0;
+ AliExternalTrackParam *trackOut1 = 0;
+ AliExternalTrackParam *ptrackIn0 = 0;
+ AliExternalTrackParam *ptrackIn1 = 0;
+
+ for (Int_t icorr=0; icorr<ncorr; icorr++) {
+ //
+ // special case of the TPC tracks crossing the CE
+ //
+ AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
+ corrections[icorr]=0;
+ if (entries>kMinEntries &&isOK){
+ AliExternalTrackParam trackIn0(refX,phi,tPar,cov);
+ AliExternalTrackParam trackIn1(refX,phi,tPar,cov);
+ ptrackIn1=&trackIn0;
+ ptrackIn0=&trackIn1;
+ //
+ if (debug) trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,pcstream);
+ if (!debug) trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,0);
+ if (debug) trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,pcstream);
+ if (!debug) trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,0);
+ //
+ if (trackOut0 && trackOut1){
+ //
+ if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,1,kTRUE,kMaxSnp)) isOK=kFALSE;
+ if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
+ // rotate all tracks to the same frame
+ if (!trackOut0->Rotate(trackIn0.GetAlpha())) isOK=kFALSE;
+ if (!trackIn1.Rotate(trackIn0.GetAlpha())) isOK=kFALSE;
+ if (!trackOut1->Rotate(trackIn0.GetAlpha())) isOK=kFALSE;
+ //
+ if (!AliTrackerBase::PropagateTrackTo(trackOut0,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
+ if (!AliTrackerBase::PropagateTrackTo(&trackIn1,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
+ if (!AliTrackerBase::PropagateTrackTo(trackOut1,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
+ //
+ corrections[icorr] = (trackOut0->GetParameter()[ptype]-trackIn0.GetParameter()[ptype]);
+ corrections[icorr]-= (trackOut1->GetParameter()[ptype]-trackIn1.GetParameter()[ptype]);
+ (*pcstream)<<"fitDebug"<< // just to debug the correction
+ "mean="<<mean<<
+ "pIn0.="<<ptrackIn0<<
+ "pIn1.="<<ptrackIn1<<
+ "pOut0.="<<trackOut0<<
+ "pOut1.="<<trackOut1<<
+ "refX="<<refXD<<
+ "\n";
+ delete trackOut0;
+ delete trackOut1;
+ }else{
+ corrections[icorr]=0;
+ isOK=kFALSE;
+ }
+ }
+ (*pcstream)<<"fit"<<
+ Form("%s=",corr->GetName())<<corrections[icorr]; // dump correction value
+ }
+ //
+ (*pcstream)<<"fit"<<"isOK="<<isOK<<"\n";
+ }
+ delete pcstream;
+}
+
+
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff