dy = +dr(r,z)*tan(phi)
-
- //
+ .x ~/NimStyle.C
+ gSystem->Load("libANALYSIS");
+ gSystem->Load("libTPCcalib");
+ TFile fcalib("CalibObjects.root");
+ TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
+ AliTPCcalibUnlinearity * calibUnlin = ( AliTPCcalibUnlinearity *)array->FindObject("calibUnlinearity");
+ //
+
*/
#include "AliESDInputHandler.h"
#include "AliAnalysisManager.h"
#include "AliTracker.h"
-#include "AliMagFMaps.h"
+#include "AliMagF.h"
#include "AliTPCCalROC.h"
#include "AliLog.h"
#include "AliTPCTracklet.h"
#include "TTimeStamp.h"
#include "AliTPCcalibDB.h"
-#include "AliTPCcalibLaser.h"
#include "AliDCSSensorArray.h"
#include "AliDCSSensor.h"
#include "TLinearFitter.h"
#include "AliTPCcalibUnlinearity.h"
+#include "AliTPCPointCorrection.h"
-AliTPCcalibUnlinearity* AliTPCcalibUnlinearity::fgInstance = 0;
ClassImp(AliTPCcalibUnlinearity)
AliTPCcalibUnlinearity::AliTPCcalibUnlinearity():
AliTPCcalibBase(),
- fDiffHistoLine(0),
- fDiffHistoPar(0),
+ fInit(kFALSE),
+ fDiffHistoLineY(0),
+ fDiffHistoParY(0),
+ fDiffHistoLineZ(0),
+ fDiffHistoParZ(0),
fFittersOutR(38),
fFittersOutZ(38),
fParamsOutR(38),
fParamsOutZ(38),
fErrorsOutR(38),
- fErrorsOutZ(38)
+ fErrorsOutZ(38),
+ fDistRPHIPlus(74),
+ fDistRPHIMinus(74),
+ fFitterQuadrantY(16*38),
+ fFitterQuadrantPhi(16*38)
{
//
// Defualt constructor
AliTPCcalibUnlinearity::AliTPCcalibUnlinearity(const Text_t *name, const Text_t *title):
AliTPCcalibBase(name,title),
- fDiffHistoLine(0),
- fDiffHistoPar(0),
+ fInit(kFALSE),
+ fDiffHistoLineY(0),
+ fDiffHistoParY(0),
+ fDiffHistoLineZ(0),
+ fDiffHistoParZ(0),
fFittersOutR(38),
fFittersOutZ(38),
fParamsOutR(38),
fParamsOutZ(38),
fErrorsOutR(38),
- fErrorsOutZ(38)
+ fErrorsOutZ(38),
+ fDistRPHIPlus(74),
+ fDistRPHIMinus(74),
+ fFitterQuadrantY(16*38),
+ fFitterQuadrantPhi(16*38)
{
//
// Non default constructor
//
- MakeFitters();
+ MakeHisto();
}
AliTPCcalibUnlinearity::~AliTPCcalibUnlinearity(){
//
//
//
- if (fDiffHistoLine) delete fDiffHistoLine;
- if (fDiffHistoPar) delete fDiffHistoPar;
+ if (fDiffHistoLineZ) delete fDiffHistoLineY;
+ if (fDiffHistoParZ) delete fDiffHistoParY;
+ if (fDiffHistoLineZ) delete fDiffHistoLineZ;
+ if (fDiffHistoParZ) delete fDiffHistoParZ;
}
-AliTPCcalibUnlinearity* AliTPCcalibUnlinearity::Instance()
-{
- //
- // Singleton implementation
- // Returns an instance of this class, it is created if neccessary
- //
- if (fgInstance == 0){
- fgInstance = new AliTPCcalibUnlinearity();
- }
- return fgInstance;
-}
-
//
//
//
+ if (!fInit) {
+ Init(); // work around for PROOF - initialize firs time used
+ }
const Int_t kdrow=10;
const Int_t kMinCluster=35;
- if (!fDiffHistoLine) MakeHisto();
if (!seed) return;
if (TMath::Abs(fMagF)>0.01) return; // use only data without field
Int_t counter[72];
//
for (Int_t is=0; is<72;is++){
if (counter[is]>kMinCluster) ProcessDiff(seed, is);
+ if (counter[is]>kMinCluster) {
+ AlignOROC(seed, is);
+ }
}
}
//
const Double_t kXmean=134;
const Int_t kdrow=10;
+ const Float_t kSagitaCut = 1;
static TLinearFitter fy1(2,"pol1");
static TLinearFitter fy2(3,"pol2");
static TLinearFitter fz1(2,"pol1");
AliTPCclusterMI *c=track->GetClusterPointer(irow);
if (!c) continue;
if (c->GetDetector()!=isec) continue;
- if (c->GetType()<0) continue;
Double_t dx = c->GetX()-kXmean;
Double_t y1 = py1[0]+py1[1]*dx;
Double_t y2 = py2[0]+py2[1]*dx+py2[2]*dx*dx;
Double_t z1 = pz1[0]+pz1[1]*dx;
Double_t z2 = pz2[0]+pz2[1]*dx+pz2[2]*dx*dx;
//
+ Double_t edgeMinus= -y1+c->GetX()*TMath::Tan(TMath::Pi()/18.);
+ Double_t edgePlus = y1+c->GetX()*TMath::Tan(TMath::Pi()/18.);
+ Int_t npads = AliTPCROC::Instance()->GetNPads(isec,c->GetRow());
+ Float_t dpad = TMath::Min(c->GetPad(), npads-1- c->GetPad());
+ Float_t dy =c->GetY()-y1;
+ //
+ // Corrections
+ //
+ AliTPCPointCorrection * corr = AliTPCPointCorrection::Instance();
+ Double_t corrclY=0, corrtrY=0, corrR=0;
+ corrclY = corr->RPhiCOGCorrection(isec,c->GetRow(), c->GetPad(),
+ c->GetY(),c->GetY(), c->GetZ(), py1[1],c->GetMax(),2.5);
+ corrtrY = corr->RPhiCOGCorrection(isec,c->GetRow(), c->GetPad(),
+ c->GetY(),y1, c->GetZ(), py1[1], c->GetMax(),2.5);
+ corrR = corr->CorrectionOutR0(kFALSE,kFALSE,c->GetX(),c->GetY(),c->GetZ(),isec);
+ //
//
- Double_t x[10];
- x[0]=isec;
- x[1]=irow;
- x[2]=c->GetY();
- x[3]=c->GetZ();
- x[4]=py1[1];
- x[5]=pz1[1];
- x[6]=py2[2]*150*150/4; // sagita 150 cm
//
- x[7]=c->GetY()-y1;
- x[8]=c->GetZ()-z1;
- fDiffHistoLine->Fill(x);
- x[7]=c->GetY()-y2;
- x[8]=c->GetZ()-z2;
- fDiffHistoPar->Fill(x);
-
- if (TMath::Abs(py2[2]*150*150/4)<1 && TMath::Abs(pz2[2]*150*150/4)<1){
- // Apply sagita cut
- AddPoint(isec,irow,c->GetZ()-z1, c->GetY()-y1,
- py1[1],pz1[1],1-TMath::Abs(c->GetZ()/250.),1);
- }
-
if (fStreamLevel>1){
TTreeSRedirector *cstream = GetDebugStreamer();
if (cstream){
"time="<<fTime<< // time stamp of event
"trigger="<<fTrigger<< // trigger
"mag="<<fMagF<< // magnetic field
- "isec="<<isec<<
+ "isec="<<isec<< // current sector
+ "npads="<<npads<< // number of pads at given sector
+ "dpad="<<dpad<< // distance to edge pad
+ //
+ "crR="<<corrR<< // radial correction
+ "cclY="<<corrclY<< // edge effect correction cl
+ "ctrY="<<corrtrY<< // edge effect correction using track
+ //
"Cl.="<<c<<
"y1="<<y1<<
"y2="<<y2<<
"py2.="<<&py2<<
"pz1.="<<&pz1<<
"pz2.="<<&pz2<<
+ "eP="<<edgePlus<<
+ "eM="<<edgeMinus<<
"\n";
}
}
+ if (TMath::Min(edgeMinus,edgePlus)<6){
+ AddPointRPHI(isec,c->GetX(),c->GetY(),c->GetZ(),y1,z1,py1[1],pz1[1],1);
+ TTreeSRedirector *cstream = GetDebugStreamer();
+ if (TMath::Abs(dy)<2 && cstream){
+ (*cstream)<<"rphi"<<
+ "isec="<<isec<< // current sector
+ "npads="<<npads<< // number of pads at given sector
+ "dpad="<<dpad<< // distance to edge pad
+ "eP="<<edgePlus<< // distance to edge
+ "eM="<<edgeMinus<< // distance to edge minus
+ //
+ "crR="<<corrR<< // radial correction
+ "cclY="<<corrclY<< // edge effect correction cl
+ "ctrY="<<corrtrY<< // edge effect correction using track
+ //
+ "dy="<<dy<< // dy
+ "Cl.="<<c<<
+ "y1="<<y1<<
+ "y2="<<y2<<
+ "z1="<<z1<<
+ "z2="<<z2<<
+ "py1.="<<&py1<<
+ "pz1.="<<&pz1<<
+ "\n";
+ }
+ }
+ if (c->GetType()<0) continue; // don't use edge rphi cluster
+ //
+ //
+ Double_t x[10];
+ x[0]=c->GetDetector();
+ x[1]=c->GetRow();
+ x[2]=c->GetY()/c->GetX();
+ x[3]=c->GetZ();
+ //
+ // apply sagita cut
+ //
+ if (TMath::Abs(py2[2]*150*150/4)<kSagitaCut &&
+ TMath::Abs(pz2[2]*150*150/4)<kSagitaCut){
+ //
+ x[4]=py1[1];
+ x[5]=c->GetY()-y1;
+ fDiffHistoLineY->Fill(x);
+ x[5]=c->GetY()-y2;
+ //fDiffHistoParY->Fill(x);
+ //
+ x[4]=pz1[1];
+ x[5]=c->GetY()-z1;
+ fDiffHistoLineZ->Fill(x);
+ x[5]=c->GetY()-z2;
+ //fDiffHistoParZ->Fill(x);
+ // Apply sagita cut
+ AddPoint(isec,c->GetX(),c->GetY(),c->GetZ(),y1,z1,py1[1],pz1[1],1);
+ }
+
}
}
+void AliTPCcalibUnlinearity::AlignOROC(AliTPCseed *track, Int_t isec){
+ //
+ // fit the tracklet in OROC sepatately in Quadrant
+ //
+ //
+ if (isec<36) return;
+ const Int_t kMinClusterF=35;
+ const Int_t kMinClusterQ=10;
+ //
+ const Int_t kdrow1 =8; // rows to skip at the end
+ const Int_t kdrow0 =2; // rows to skip at beginning
+ const Float_t kedgey=3;
+ const Float_t kMaxDist=0.5;
+ const Float_t kMaxCorrY=0.1;
+ const Float_t kPRFWidth = 0.6; //cut 2 sigma of PRF
+ //
+ //
+ AliTPCROC * roc = AliTPCROC::Instance();
+ AliTPCPointCorrection * corr = AliTPCPointCorrection::Instance();
+ const Double_t kXmean=roc->GetPadRowRadii(isec,53);
+ //
+ // full track fit parameters
+ //
+ TLinearFitter fyf(2,"pol1");
+ TLinearFitter fzf(2,"pol1");
+ TVectorD pyf(2), peyf(2),pzf(2), pezf(2);
+ Int_t nf=0;
+ //
+ // make full fit as reference
+ //
+ for (Int_t iter=0; iter<2; iter++){
+ fyf.ClearPoints();
+ for (Int_t irow=kdrow0;irow<159-kdrow1;irow++) {
+ AliTPCclusterMI *c=track->GetClusterPointer(irow);
+ if (!c) continue;
+ if (c->GetDetector()!=isec) continue;
+ if (c->GetRow()<kdrow0) continue;
+ Double_t dx = c->GetX()-kXmean;
+ Double_t x[2]={dx, dx*dx};
+ if (iter==0 &&c->GetType()<0) continue;
+ if (iter==1){
+ Double_t yfit = pyf[0]+pyf[1]*dx;
+ Double_t dedge = c->GetX()*TMath::Tan(TMath::Pi()/18.)-TMath::Abs(yfit);
+ if (TMath::Abs(c->GetY()-yfit)>kMaxDist) continue;
+ if (dedge<kedgey) continue;
+ Double_t corrtrY =
+ corr->RPhiCOGCorrection(isec,c->GetRow(), c->GetPad(),
+ c->GetY(),yfit, c->GetZ(), pyf[1], c->GetMax(),2.5);
+ if (TMath::Abs(corrtrY)>kMaxCorrY) continue;
+ }
+ fyf.AddPoint(x,c->GetY(),0.1);
+ fzf.AddPoint(x,c->GetZ(),0.1);
+ }
+ nf = fyf.GetNpoints();
+ if (nf<kMinClusterF) return; // not enough points - skip
+ fyf.Eval();
+ fyf.GetParameters(pyf);
+ fyf.GetErrors(peyf);
+ fzf.Eval();
+ fzf.GetParameters(pzf);
+ fzf.GetErrors(pezf);
+ }
+ //
+ // Make Fitters and params for 3 fitters
+ //
+ TLinearFitter *fitters[4];
+ Int_t npoints[4];
+ TVectorD params[4];
+ TVectorD errors[4];
+ Double_t chi2C[4];
+ for (Int_t i=0;i<4;i++) {
+ fitters[i] = new TLinearFitter(2,"pol1");
+ npoints[i]=0;
+ params[i].ResizeTo(2);
+ errors[i].ResizeTo(2);
+ }
+ //
+ // Update fitters
+ //
+ for (Int_t irow=kdrow0;irow<159-kdrow1;irow++) {
+ AliTPCclusterMI *c=track->GetClusterPointer(irow);
+ if (!c) continue;
+ if (c->GetDetector()!=isec) continue;
+ if (c->GetRow()<kdrow0) continue;
+ Double_t dx = c->GetX()-kXmean;
+ Double_t x[2]={dx, dx*dx};
+ Double_t yfit = pyf[0]+pyf[1]*dx;
+ Double_t dedge = c->GetX()*TMath::Tan(TMath::Pi()/18.)-TMath::Abs(yfit);
+ if (TMath::Abs(c->GetY()-yfit)>kMaxDist) continue;
+ if (dedge<kedgey) continue;
+ Double_t corrtrY =
+ corr->RPhiCOGCorrection(isec,c->GetRow(), c->GetPad(),
+ c->GetY(),yfit, c->GetZ(), pyf[1], c->GetMax(),2.5);
+ if (TMath::Abs(corrtrY)>kMaxCorrY) continue;
+ if (dx<0){
+ if (yfit<-kPRFWidth) fitters[0]->AddPoint(x,c->GetY(),0.1);
+ if (yfit>kPRFWidth) fitters[1]->AddPoint(x,c->GetY(),0.1);
+ }
+ if (dx>0){
+ if (yfit<-kPRFWidth) fitters[2]->AddPoint(x,c->GetY(),0.1);
+ if (yfit>kPRFWidth) fitters[3]->AddPoint(x,c->GetY(),0.1);
+ }
+ }
+
+ for (Int_t i=0;i<4;i++) {
+ npoints[i] = fitters[i]->GetNpoints();
+ if (npoints[i]>=kMinClusterQ){
+ fitters[i]->Eval();
+ Double_t chi2Fac = TMath::Sqrt(fitters[i]->GetChisquare()/(fitters[i]->GetNpoints()-2));
+ chi2C[i]=chi2Fac;
+ fitters[i]->GetParameters(params[i]);
+ fitters[i]->GetErrors(errors[i]);
+ // renormalize errors
+ errors[i][0]*=chi2Fac;
+ errors[i][1]*=chi2Fac;
+ }
+ }
+ //
+ // Fill fitters
+ //
+ for (Int_t i0=0;i0<4;i0++){
+ for (Int_t i1=0;i1<4;i1++){
+ if (i0==i1) continue;
+ if(npoints[i0]<kMinClusterQ) continue;
+ if(npoints[i1]<kMinClusterQ) continue;
+ Int_t index0=i0*4+i1;
+ Double_t xy[1];
+ Double_t xfi[1];
+ xy[0] = pyf[1];
+ xfi[0] = pyf[1];
+ //
+ Double_t dy = params[i1][0]-params[i0][0];
+ Double_t sy = TMath::Sqrt(errors[i1][0]*errors[i1][0]+errors[i0][0]*errors[i0][0]);
+ Double_t dphi = params[i1][1]-params[i0][1];
+ Double_t sphi = TMath::Sqrt(errors[i1][1]*errors[i1][1]+errors[i0][1]*errors[i0][1]);
+ //
+ Int_t sector = isec-36;
+ //
+ if (TMath::Abs(dy/(sy+0.1))>5.) continue; // 5 sigma cut
+ if (TMath::Abs(dphi/(sphi+0.004))>5.) continue; // 5 sigma cut
+ TLinearFitter * fitterY,*fitterPhi;
+ fitterY = (TLinearFitter*) fFitterQuadrantY.At(16*sector+index0);
+ if (TMath::Abs(dy/(sy+0.1))<5.) fitterY->AddPoint(xy,dy,sy);
+ fitterY = (TLinearFitter*) fFitterQuadrantY.At(16*36+index0);
+ if (TMath::Abs(dy/(sy+0.1))<5.) fitterY->AddPoint(xy,dy,sy);
+ //
+ fitterPhi = (TLinearFitter*) fFitterQuadrantPhi.At(16*sector+index0);
+ if (TMath::Abs(dphi/(sphi+0.001))<5.) fitterPhi->AddPoint(xfi,dphi,sphi);
+ fitterPhi = (TLinearFitter*) fFitterQuadrantPhi.At(16*36+index0);
+ if (TMath::Abs(dphi/(sphi+0.001))<5.) fitterPhi->AddPoint(xfi,dphi,sphi);
+ }
+ }
+ //
+ // Dump debug information
+ //
+ if (fStreamLevel>0){
+ TTreeSRedirector *cstream = GetDebugStreamer();
+ Int_t sector = isec-36;
+ if (cstream){
+ for (Int_t i0=0;i0<4;i0++){
+ for (Int_t i1=i0;i1<4;i1++){
+ if (i0==i1) continue;
+ if(npoints[i0]<kMinClusterQ) continue;
+ if(npoints[i1]<kMinClusterQ) continue;
+ (*cstream)<<"fitQ"<<
+ "run="<<fRun<< // run number
+ "event="<<fEvent<< // event number
+ "time="<<fTime<< // time stamp of event
+ "trigger="<<fTrigger<< // trigger
+ "mag="<<fMagF<< // magnetic field
+ "sec="<<sector<< // current sector from 0 to 36
+ "isec="<<isec<< // current sector
+ // full fit
+ "nf="<<nf<< // total number of points
+ "pyf.="<<&pyf<< // full OROC fit y
+ "pzf.="<<&pzf<< // full OROC fit z
+ // quadrant fit
+ "i0="<<i0<< // quadrant number
+ "i1="<<i1<<
+ "n0="<<npoints[i0]<< // number of points
+ "n1="<<npoints[i1]<<
+ "p0.="<<¶ms[i0]<< // parameters
+ "p1.="<<¶ms[i1]<<
+ "e0.="<<&errors[i0]<< // errors
+ "e1.="<<&errors[i1]<<
+ "chi0="<<chi2C[i0]<< // chi2s
+ "chi1="<<chi2C[i1]<<
+
+ "\n";
+ }
+ }
+ }
+ }
+}
+
+
+
+
+
void AliTPCcalibUnlinearity::MakeHisto(){
//
//
//
axisName[0]="sector";
xmin[0]=0; xmax[0]=72; nbins[0]=72;
- //
+ //
axisName[1]="row";
- xmin[1]=0; xmax[1]=159; nbins[1]=159;
- //
- axisName[2]="ly";
- xmin[2]=-50; xmax[2]=50; nbins[2]=10;
- //
- axisName[3]="lz";
- xmin[3]=-250; xmax[3]=250; nbins[3]=50;
- //
- axisName[4]="p2";
- xmin[4]=-0.6; xmax[4]=0.6; nbins[4]=12;
- //
- axisName[5]="p3";
- xmin[5]=-0.6; xmax[5]=0.6; nbins[5]=12;
+ xmin[1]=0; xmax[1]=96; nbins[1]=96;
//
- axisName[6]="p4";
- xmin[6]=-2; xmax[6]=2; nbins[6]=20;
- //
- axisName[7]="dy";
- xmin[7]=-0.5; xmax[7]=0.5; nbins[7]=100;
+ axisName[2]="tphi"; // tan phi - ly/lx
+ xmin[2]=-0.17; xmax[2]=0.17; nbins[2]=5;
+ //
+ axisName[3]="lz"; // global z
+ xmin[3]=-250; xmax[3]=250; nbins[3]=10;
+ //
+ axisName[4]="k"; // tangent - in angle
+ xmin[4]=-0.5; xmax[4]=0.5; nbins[4]=10;
+ //
+ //
+ axisName[5]="delta"; // delta
+ xmin[5]=-0.5; xmax[5]=0.5; nbins[5]=100;
+ //
+ //
+ fDiffHistoLineY = new THnSparseF("hnDistHistoLineY","hnDistHistoLineY",6,nbins,xmin,xmax);
+ fDiffHistoParY = new THnSparseF("hnDistHistoParY","hnDistHistoParY",6,nbins,xmin,xmax);
+ fDiffHistoLineZ = new THnSparseF("hnDistHistoLineZ","hnDistHistoLineZ",6,nbins,xmin,xmax);
+ fDiffHistoParZ = new THnSparseF("hnDistHistoParZ","hnDistHistoParZ",6,nbins,xmin,xmax);
+ // for (Int_t i=0; i<8;i++){
+// fDiffHistoLineY->GetAxis(i)->SetName(axisName[i].Data());
+// fDiffHistoParY->GetAxis(i)->SetName(axisName[i].Data());
+// fDiffHistoLineY->GetAxis(i)->SetTitle(axisName[i].Data());
+// fDiffHistoParY->GetAxis(i)->SetTitle(axisName[i].Data());
+// fDiffHistoLineZ->GetAxis(i)->SetName(axisName[i].Data());
+// fDiffHistoParZ->GetAxis(i)->SetName(axisName[i].Data());
+// fDiffHistoLineZ->GetAxis(i)->SetTitle(axisName[i].Data());
+// fDiffHistoParZ->GetAxis(i)->SetTitle(axisName[i].Data());
+// }
+
//
- axisName[8]="dz";
- xmin[8]=-0.5; xmax[8]=0.5; nbins[8]=100;
+ //
//
- fDiffHistoLine = new THnSparseF("DistHistoLine","DistHistoLine",9,nbins,xmin,xmax);
- fDiffHistoPar = new THnSparseF("DistHistoPar","DistHistoPar",9,nbins,xmin,xmax);
- for (Int_t i=0; i<9;i++){
- fDiffHistoLine->GetAxis(i)->SetName(axisName[i].Data());
- fDiffHistoPar->GetAxis(i)->SetName(axisName[i].Data());
- fDiffHistoLine->GetAxis(i)->SetTitle(axisName[i].Data());
- fDiffHistoPar->GetAxis(i)->SetTitle(axisName[i].Data());
+ char hname[1000];
+ TH2F * his=0;
+ for (Int_t isec=0;isec<74;isec++){
+ sprintf(hname,"DeltaRPhiPlus_Sector%d",isec);
+ his = new TH2F(hname,hname,100,1,4,100,-0.5,0.5);
+ his->SetDirectory(0);
+ fDistRPHIPlus.AddAt(his,isec);
+ sprintf(hname,"DeltaRPhiMinus_Sector%d",isec);
+ his = new TH2F(hname,hname,100,1,4,100,-0.5,0.5);
+ his->SetDirectory(0);
+ fDistRPHIMinus.AddAt(his,isec);
}
}
AliTPCcalibBase::Terminate();
}
+Long64_t AliTPCcalibUnlinearity::Merge(TCollection *list) {
+ //
+ // merge results
+ //
+ TIterator* iter = list->MakeIterator();
+ AliTPCcalibUnlinearity* cal = 0;
+ //
+ while ((cal = (AliTPCcalibUnlinearity*)iter->Next())) {
+ if (!cal->InheritsFrom(AliTPCcalibUnlinearity::Class())) {
+ return -1;
+ }
+ Add(cal);
+ }
+ return 0;
+}
+
+void AliTPCcalibUnlinearity::Add(AliTPCcalibUnlinearity * calib){
+ //
+ //
+ //
+ if (!fInit) Init();
+ if (fDiffHistoLineY && calib->fDiffHistoLineY){
+ fDiffHistoLineY->Add(calib->fDiffHistoLineY);
+ }
+ if (fDiffHistoParY && calib->fDiffHistoParY){
+ fDiffHistoParY->Add(calib->fDiffHistoParY);
+ }
+ if (fDiffHistoLineZ && calib->fDiffHistoLineZ){
+ fDiffHistoLineZ->Add(calib->fDiffHistoLineZ);
+ }
+ if (fDiffHistoParZ && calib->fDiffHistoParZ){
+ fDiffHistoParZ->Add(calib->fDiffHistoParZ);
+ }
+ for (Int_t isec=0;isec<38;isec++){
+ TLinearFitter *f0r = (TLinearFitter*)fFittersOutR.At(isec);
+ TLinearFitter *f1r = (TLinearFitter*)(calib->fFittersOutR.At(isec));
+ if (f0r&&f1r) f0r->Add(f1r);
+ TLinearFitter *f0z = (TLinearFitter*)fFittersOutZ.At(isec);
+ TLinearFitter *f1z = (TLinearFitter*)(calib->fFittersOutZ.At(isec));
+ if (f0z&&f1z) f0z->Add(f1z);
+ }
+
+ for (Int_t isec=0;isec<16*38;isec++){
+ TLinearFitter *f0y = (TLinearFitter*)fFitterQuadrantY.At(isec);
+ TLinearFitter *f1y = (TLinearFitter*)(calib->fFitterQuadrantY.At(isec));
+ if (f0y&&f1y) f0y->Add(f1y);
+ TLinearFitter *f0phi = (TLinearFitter*)fFitterQuadrantPhi.At(isec);
+ TLinearFitter *f1phi = (TLinearFitter*)(calib->fFitterQuadrantPhi.At(isec));
+ if (f0phi&&f1phi) f0phi->Add(f1phi);
+ }
+
+
+ for (Int_t isec=0;isec<74;isec++){
+ TH2F * h0p= (TH2F*)(fDistRPHIPlus.At(isec));
+ TH2F * h1p= (TH2F*)(calib->fDistRPHIPlus.At(isec));
+ if (h0p&&h1p) h0p->Add(h1p);
+ TH2F * h0m= (TH2F*)(fDistRPHIMinus.At(isec));
+ TH2F * h1m= (TH2F*)(calib->fDistRPHIMinus.At(isec));
+ if (h0m&&h1m) h0m->Add(h1m);
+ }
+
+
+}
+
+
-void AliTPCcalibUnlinearity::DumpTree(){
+void AliTPCcalibUnlinearity::DumpTree(const char *fname){
//
//
//
- if (fStreamLevel==0) return;
- TTreeSRedirector *cstream = GetDebugStreamer();
+ TTreeSRedirector *cstream =new TTreeSRedirector(fname);
if (!cstream) return;
//
THnSparse *his=0;
Int_t *bins = new Int_t[10];
//
for (Int_t ihis=0; ihis<2; ihis++){
- his = (ihis==0)? fDiffHistoLine:fDiffHistoPar;
+ his = (ihis==0)? fDiffHistoLineY:fDiffHistoParY;
if (!his) continue;
//
Int_t ndim = his->GetNdimensions();
"\n";
}
}
+ delete cstream;
}
-void AliTPCcalibUnlinearity::AddPoint(Int_t sector, Int_t row, Double_t dz, Double_t dy, Double_t p2, Double_t p3, Double_t dr, Int_t npoints){
+void AliTPCcalibUnlinearity::AddPoint(Int_t sector, Double_t cx, Double_t cy, Double_t cz, Double_t ty, Double_t tz, Double_t ky, Double_t kz, Int_t npoints){
//
//
// Simple distortion fit in outer sectors
// - Decay length 10 and 5 cm
// - Decay amplitude - Z dependent
//
- /*
- chainUnlin->SetAlias("side","(-1+((sector%36)<18)*2)");
- chainUnlin->SetAlias("sa","sin(pi*sector*20/180)");
- chainUnlin->SetAlias("ca","cos(pi*sector*20/180)");
- chainUnlin->SetAlias("dzt","dz*side");
- chainUnlin->SetAlias("dr","(1-abs(lz/250))");
- chainUnlin->SetAlias("dout","(159-row)*1.5");
- chainUnlin->SetAlias("din","row*0.75");
- chainUnlin->SetAlias("eout10","exp(-(dout)/10.)");
- chainUnlin->SetAlias("eout5","exp(-(dout)/5.)");
- */
-
+
Double_t side = (-1+((sector%36)<18)*2); // A side +1 C side -1
- Double_t dzt = dz*side;
- Double_t dout = (159-row)*1.5; // distance to the last pad row - (valid only for long pads)
- if (dout>30) return; // use only edge pad rows
- dr-=0.5; // dr shifted to the middle - reduce numerical instabilities
-
+ Double_t dzt = (cz-tz)*side;
+ Double_t radius = TMath::Sqrt(cx*cx+cy*cy);
+ AliTPCROC * roc = AliTPCROC::Instance();
+ Double_t xout = roc->GetPadRowRadiiUp(roc->GetNRows(37)-1);
+ Double_t dout = xout-radius;
+ if (dout>30) return;
+ //drift
+ Double_t dr = 0.5 - TMath::Abs(cz/250.);
+ Double_t dy = cy-ty;
Double_t eout10 = TMath::Exp(-dout/10.);
Double_t eout5 = TMath::Exp(-dout/5.);
//
xxxZ[4]=dr*dr*eoutp; //fstring+="dr*dr*eoutp++";
xxxZ[5]=dr*dr*eoutm; //fstring+="dr*dr*eoutm++";
//
- xxxR[0]=p2*eoutp; //fstring+="eoutp++";
- xxxR[1]=p2*eoutm; //fstring+="eoutm++";
- xxxR[2]=p2*dr*eoutp; //fstring+="dr*eoutp++";
- xxxR[3]=p2*dr*eoutm; //fstring+="dr*eoutm++";
- xxxR[4]=p2*dr*dr*eoutp; //fstring+="dr*dr*eoutp++";
- xxxR[5]=p2*dr*dr*eoutm; //fstring+="dr*dr*eoutm++";
+ xxxR[0]=ky*eoutp; //fstring+="eoutp++";
+ xxxR[1]=ky*eoutm; //fstring+="eoutm++";
+ xxxR[2]=ky*dr*eoutp; //fstring+="dr*eoutp++";
+ xxxR[3]=ky*dr*eoutm; //fstring+="dr*eoutm++";
+ xxxR[4]=ky*dr*dr*eoutp; //fstring+="dr*dr*eoutp++";
+ xxxR[5]=ky*dr*dr*eoutm; //fstring+="dr*dr*eoutm++";
//
- xxxRZ[0]=p3*eoutp; //fstring+="eoutp++";
- xxxRZ[1]=p3*eoutm; //fstring+="eoutm++";
- xxxRZ[2]=p3*dr*eoutp; //fstring+="dr*eoutp++";
- xxxRZ[3]=p3*dr*eoutm; //fstring+="dr*eoutm++";
- xxxRZ[4]=p3*dr*dr*eoutp; //fstring+="dr*dr*eoutp++";
- xxxRZ[5]=p3*dr*dr*eoutm; //fstring+="dr*dr*eoutm++";
+ xxxRZ[0]=kz*eoutp; //fstring+="eoutp++";
+ xxxRZ[1]=kz*eoutm; //fstring+="eoutm++";
+ xxxRZ[2]=kz*dr*eoutp; //fstring+="dr*eoutp++";
+ xxxRZ[3]=kz*dr*eoutm; //fstring+="dr*eoutm++";
+ xxxRZ[4]=kz*dr*dr*eoutp; //fstring+="dr*dr*eoutp++";
+ xxxRZ[5]=kz*dr*dr*eoutm; //fstring+="dr*dr*eoutm++";
//
TLinearFitter * fitter=0;
Double_t err=0.1;
}
}
}
+void AliTPCcalibUnlinearity::AddPointRPHI(Int_t sector, Double_t cx, Double_t cy, Double_t cz, Double_t ty, Double_t tz, Double_t /*ky*/, Double_t /*kz*/, Int_t npoints){
+ //
+ //
+ //
+ Float_t kMaxDz = 0.5; // cut on z in cm
+ Double_t dy = cy-ty;
+ Double_t dz = cz-tz;
+ if (TMath::Abs(dz)>kMaxDz) return;
+ Double_t edgeMinus= -ty+cx*TMath::Tan(TMath::Pi()/18.);
+ Double_t edgePlus = ty+cx*TMath::Tan(TMath::Pi()/18.);
+ //
+ TH2F * his =0;
+ his = (TH2F*)fDistRPHIPlus.At(sector);
+ his->Fill(edgePlus,dy,npoints);
+ his = (TH2F*)((sector<36)? fDistRPHIPlus.At(72):fDistRPHIPlus.At(73));
+ his->Fill(edgePlus,dy,npoints);
+ his = (TH2F*)fDistRPHIMinus.At(sector);
+ his->Fill(edgeMinus,dy,npoints);
+ his = (TH2F*)((sector<36)? fDistRPHIMinus.At(72):fDistRPHIMinus.At(73));
+ his->Fill(edgeMinus,dy,npoints);
+}
+
-void AliTPCcalibUnlinearity::MakeFitters(){
+void AliTPCcalibUnlinearity::Init(){
//
//
//
+ //
+ MakeHisto();
// Make outer fitters
TLinearFitter * fitter = 0;
for (Int_t ifit=0; ifit<38; ifit++){
fitter->StoreData(kFALSE);
fFittersOutZ.AddAt(fitter,ifit);
}
+ for (Int_t ifit=0; ifit<16*38;ifit++){
+ fitter = new TLinearFitter(2,"hyp1");
+ fitter->StoreData(kFALSE);
+ fFitterQuadrantY.AddAt(fitter,ifit);
+ fitter = new TLinearFitter(2,"hyp1");
+ fitter->StoreData(kFALSE);
+ fFitterQuadrantPhi.AddAt(fitter,ifit);
+ }
+ fInit= kTRUE;
}
void AliTPCcalibUnlinearity::EvalFitters(){
}
}
-Double_t AliTPCcalibUnlinearity::GetDr(Int_t sector, Double_t dout, Double_t dr){
- //
- //
- //
- TVectorD * vec = GetParamOutR(sector);
- if (!vec) return 0;
- dr-=0.5; // dr=0 at the middle drift length
- Double_t eout10 = TMath::Exp(-dout/10.);
- Double_t eout5 = TMath::Exp(-dout/5.);
- Double_t eoutp = (eout10+eout5)*0.5;
- Double_t eoutm = (eout10-eout5)*0.5;
-
- Double_t result=0;
- result+=(*vec)[1]*eoutp;
- result+=(*vec)[2]*eoutm;
- result+=(*vec)[3]*eoutp*dr;
- result+=(*vec)[4]*eoutm*dr;
- result+=(*vec)[5]*eoutp*dr*dr;
- result+=(*vec)[6]*eoutm*dr*dr;
- return result;
-}
-
-
-Double_t AliTPCcalibUnlinearity::GetGDr(Int_t stype,Float_t gx, Float_t gy,Float_t gz){
- //
- //
- //
- Double_t alpha = TMath::ATan2(gy,gx);
- if (alpha<0) alpha+=TMath::Pi()*2;
- Int_t lsector = Int_t(18*alpha/(TMath::Pi()*2));
- alpha = (lsector+0.5)*TMath::Pi()/9.;
- //
- Double_t r[3];
- r[0]=gx; r[1]=gy;
- Double_t cs=TMath::Cos(-alpha), sn=TMath::Sin(-alpha), x=r[0];
- r[0]=x*cs - r[1]*sn; r[1]=x*sn + r[1]*cs;
- //
- AliTPCROC * roc = AliTPCROC::Instance();
- Double_t xout = roc->GetPadRowRadiiUp(roc->GetNRows(37)-1);
- Double_t dout = xout-r[0];
- if (dout<0) return 0;
- Double_t dr = 1-TMath::Abs(gz/250.);
- if (gz<0) lsector+=18;
- if (stype==0) lsector = (gz>0) ? 36:37;
- if (stype<0) return lsector; //
- return GetDr(lsector,dout,dr);
-}
-
-
-Double_t AliTPCcalibUnlinearity::GetDz(Int_t sector, Double_t dout, Double_t dr){
- //
- //
- //
- TVectorD * vec = GetParamOutZ(sector);
- if (!vec) return 0;
- dr-=0.5; // 0 at the middle
- Double_t eout10 = TMath::Exp(-dout/10.);
- Double_t eout5 = TMath::Exp(-dout/5.);
- Double_t eoutp = (eout10+eout5)*0.5;
- Double_t eoutm = (eout10-eout5)*0.5;
-
-
- Double_t result=0;
- result+=(*vec)[1]*eoutp;
- result+=(*vec)[2]*eoutm;
- result+=(*vec)[3]*eoutp*dr;
- result+=(*vec)[4]*eoutm*dr;
- result+=(*vec)[5]*eoutp*dr*dr;
- result+=(*vec)[6]*eoutm*dr*dr;
- return result;
-}
-
-
-Double_t AliTPCcalibUnlinearity::SGetDr(Int_t sector, Double_t dout, Double_t dr){
- //
- //
- //
- return fgInstance->GetDr(sector,dout,dr);
-}
-Double_t AliTPCcalibUnlinearity::SGetGDr(Int_t stype,Float_t gx, Float_t gy,Float_t gz){
- //
- //
- //
- return fgInstance->GetGDr(stype,gx,gy,gz);
-}
-Double_t AliTPCcalibUnlinearity::SGetDz(Int_t sector, Double_t dout, Double_t dr){
- //
- //
- //
- return fgInstance->GetDz(sector,dout,dr);
-}
if (i%10000==0) printf("%d\n",(Int_t)i);
Int_t row= cl->GetRow();
if (cl->GetDetector()>36) row+=64;
+ calib->AddPointRPHI(cl->GetDetector(),cl->GetX(),cl->GetY(),cl->GetZ(),
+ ty,tz,(*vy)[1],(*vz)[1],1);
+
if (cl->GetType()<0) continue;
Double_t dy = cl->GetY()-ty;
Double_t dz = cl->GetZ()-tz;
- Double_t dr = 1-TMath::Abs(cl->GetZ()/250);
//
//
if (TMath::Abs(dy)>0.25) continue;
if (TMath::Abs((*vy2)[2]*150*150/4)>1) continue;
if (TMath::Abs((*vz2)[2]*150*150/4)>1) continue;
// Apply sagita cut
- calib->AddPoint(cl->GetDetector(), row, dz, dy,
- (*vy)[1],(*vz)[1], dr, 1);
+ if (cl->GetType()>=0) {
+ calib->AddPoint(cl->GetDetector(),cl->GetX(),cl->GetY(),cl->GetZ(),
+ ty,tz,(*vy)[1],(*vz)[1],1);
+ }
}
}
}
//1. Load Parameters to be modified:
//e.g:
- AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT");
+ AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
AliCDBManager::Instance()->SetRun(0)
AliTPCClusterParam * param = AliTPCcalibDB::Instance()->GetClusterParam();
//
TString *strZ0 = toolkit.FitPlane(chainUnlinD,"(Cl.fZ-z2):1",fstringZ.Data(), "Cl.fDetector<36"+cutAll, chi2,npoints,fitParamZ0,covMatrix,-1,0,maxPoints);
printf("Z0 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
strZ0->Tokenize("++")->Print();
- param->fPosZcor[0] = (TVectorD*) fitParamZ0.Clone();
+ param->PosZcor(0) = (TVectorD*) fitParamZ0.Clone();
//
TString *strZ1 = toolkit.FitPlane(chainUnlinD,"(Cl.fZ-z2):1",fstringZ.Data(), "Cl.fDetector>36"+cutAll, chi2,npoints,fitParamZ1,covMatrix,-1,0,maxPoints);
//
printf("Z1 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
strZ1->Tokenize("++")->Print();
- param->fPosZcor[1] = (TVectorD*) fitParamZ1.Clone();
- param->fPosZcor[2] = (TVectorD*) fitParamZ1.Clone();
+ param->PosZcor(1) = (TVectorD*) fitParamZ1.Clone();
+ param->PosZcor(2) = (TVectorD*) fitParamZ1.Clone();
//
// Y corrections
//
TString *strY0 = toolkit.FitPlane(chainUnlinD,"(Cl.fY-y2):1",fstringY.Data(), "Cl.fDetector<36"+cutAll, chi2,npoints,fitParamY0,covMatrix,-1,0,maxPoints);
printf("Y0 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
strY0->Tokenize("++")->Print();
- param->fPosYcor[0] = (TVectorD*) fitParamY0.Clone();
+ param->PosYcor(0) = (TVectorD*) fitParamY0.Clone();
TString *strY1 = toolkit.FitPlane(chainUnlinD,"(Cl.fY-y2):1",fstringY.Data(), "Cl.fDetector>36"+cutAll, chi2,npoints,fitParamY1,covMatrix,-1,0,maxPoints);
//
printf("Y1 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
strY1->Tokenize("++")->Print();
- param->fPosYcor[1] = (TVectorD*) fitParamY1.Clone();
- param->fPosYcor[2] = (TVectorD*) fitParamY1.Clone();
+ param->PosYcor(1) = (TVectorD*) fitParamY1.Clone();
+ param->PosYcor(2) = (TVectorD*) fitParamY1.Clone();
//
//
//
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff