From: marian Date: Fri, 11 May 2012 21:48:11 +0000 (+0000) Subject: Adding first version of class for the TPC correction fit. X-Git-Url: http://git.uio.no/git/?a=commitdiff_plain;h=93096a078cd1c9a5220e0f362b34769ab3ae004a;p=u%2Fmrichter%2FAliRoot.git Adding first version of class for the TPC correction fit. 1. Parameters to calculate- dO/dp 2. Creation of the distortion maps from the residual histograms 3. Making fit trees Some functions, for the moment function present in the AliTPCPreprocesorOffline. --- diff --git a/TPC/AliTPCCorrectionFit.cxx b/TPC/AliTPCCorrectionFit.cxx new file mode 100644 index 00000000000..5fee6d0d9c2 --- /dev/null +++ b/TPC/AliTPCCorrectionFit.cxx @@ -0,0 +1,1331 @@ +/************************************************************************** + * 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 +#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; ibin1GetAxis(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; ibin2GetAxis(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="<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; ibin3GetAxis(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="<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; ibin3GetAxis(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; ibin2GetAxis(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; ibin1GetAxis(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)<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 (ibin3GetAxis(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)< 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)<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; ientryGetEntry(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])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])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])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="<_.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; ientryGetEntry(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="<_.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; ientryGetEntry(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 (entries1) isOK=kFALSE; + // + Double_t dRrec=0; // dummy value - needed for points - e.g for laser + // + (*pcstream)<<"fit"<< + "run="<