#include <AliCDBStorage.h>
#include <AliCDBId.h>
#include <AliCDBMetaData.h>
+#include "TVectorD.h"
+
+#include "TRandom.h"
+#include "AliTPCTransform.h"
+#include "AliTPCcalibDB.h"
+#include "AliTPCExB.h"
+#include "AliTPCCorrection.h"
+#include "AliTPCRecoParam.h"
-#include "TRandom.h"
#include "AliExternalTrackParam.h"
#include "AliTrackPointArray.h"
#include "TDatabasePDG.h"
#include "AliTrackerBase.h"
#include "AliTPCROC.h"
#include "THnSparse.h"
-
+#include "AliTPCLaserTrack.h"
#include "AliTPCCorrection.h"
for (Double_t radius=kRTPC0; radius<kRTPC1; radius++){
AliTrackerBase::PropagateTrackToBxByBz(&track,radius,kMass,3,kTRUE,kMaxSnp);
track.GetXYZ(xyz);
- xyz[0]+=gRandom->Gaus(0,0.001);
- xyz[1]+=gRandom->Gaus(0,0.001);
+ xyz[0]+=gRandom->Gaus(0,0.005);
+ xyz[1]+=gRandom->Gaus(0,0.005);
+ xyz[2]+=gRandom->Gaus(0,0.005);
AliTrackPoint pIn0; // space point
AliTrackPoint pIn1;
Int_t sector= (xyz[2]>0)? 0:18;
npoints++;
if (npoints>=npoints0) break;
}
+ if (npoints<npoints0/2) return 0;
//
// refit track
//
track0 = AliTrackerBase::MakeSeed(point1, point2, point3);
track1 = AliTrackerBase::MakeSeed(point1, point2, point3);
-
for (Int_t jpoint=0; jpoint<npoints; jpoint++){
Int_t ipoint= (dir>0) ? jpoint: npoints-1-jpoint;
//
Double_t deltaYX=deltaX*TMath::Tan(TMath::ASin(track1->GetSnp())); // deltaY due delta X
Double_t deltaZX=deltaX*track1->GetTgl(); // deltaZ due delta X
- pointPos[0]=prot1.GetY()-deltaYX;//local y
- pointPos[1]=prot1.GetZ()-deltaZX;//local z
+ pointPos[0]=prot1.GetY()+deltaYX;//local y is sign correct?
+ pointPos[1]=prot1.GetZ()+deltaZX;//local z is sign correct?
pointCov[0]=prot1.GetCov()[3];//simay^2
pointCov[1]=prot1.GetCov()[4];//sigmayz
pointCov[2]=prot1.GetCov()[5];//sigmaz^2
const Double_t kMaxSnp = 0.85;
const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
// const Double_t kB2C=-0.299792458e-3;
- const Int_t kMinEntries=50;
+ const Int_t kMinEntries=50;
Double_t phi,theta, snp, mean,rms, entries;
tinput->SetBranchAddress("theta",&theta);
tinput->SetBranchAddress("phi", &phi);
//
Int_t nentries=tinput->GetEntries();
Int_t ncorr=corrArray->GetEntries();
- Double_t corrections[100]; //
+ 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};
Double_t refX=0;
Int_t dir=0;
- if (dtype==0) {refX=85; dir=-1;}
- if (dtype==1) {refX=275; dir=1;}
- if (dtype==2) {refX=85; dir=-1;}
+ if (dtype==0) {refX=85.; dir=-1;}
+ if (dtype==1) {refX=275.; dir=1;}
+ if (dtype==2) {refX=85.; dir=-1;}
+ if (dtype==3) {refX=360.; dir=-1;}
//
for (Int_t ientry=0; ientry<nentries; ientry+=step){
tinput->GetEntry(ientry);
+ if (TMath::Abs(snp)>kMaxSnp) continue;
tPar[0]=0;
tPar[1]=theta*refX;
tPar[2]=snp;
Double_t bz=AliTrackerBase::GetBz();
if (refX>10. && TMath::Abs(bz)>0.1 ) tPar[4]=snp/(refX*bz*kB2C*2);
tPar[4]+=(gRandom->Rndm()-0.5)*0.02;
- if (TMath::Abs(snp)>0.251) continue;
+ AliExternalTrackParam track(refX,phi,tPar,cov);
+ Double_t xyz[3];
+ track.GetXYZ(xyz);
+ Int_t id=0;
+ Double_t dRrec=0; // dummy value - needed for points - e.g for laser
+
(*pcstream)<<"fit"<<
"bz="<<bz<< // magnetic filed used
"dtype="<<dtype<< // detector match type
"snp="<<snp<< // snp
"mean="<<mean<< // mean dist value
"rms="<<rms<< // rms
+ "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
+ "id="<<id<< // track id
"entries="<<entries;// number of entries in bin
//
for (Int_t icorr=0; icorr<ncorr; icorr++) {
AliExternalTrackParam *trackOut = 0;
if (debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,pcstream);
if (!debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,0);
- if (dtype==0) {refX=85; dir=-1;}
- if (dtype==1) {refX=275; dir=1;}
+ if (dtype==0) {refX=85.; dir=-1;}
+ if (dtype==1) {refX=275.; dir=1;}
if (dtype==2) {refX=0; dir=-1;}
+ if (dtype==3) {refX=360.; dir=-1;}
//
- AliTrackerBase::PropagateTrackToBxByBz(&trackIn,refX,kMass,3,kTRUE,kMaxSnp);
- AliTrackerBase::PropagateTrackToBxByBz(trackOut,refX,kMass,3,kTRUE,kMaxSnp);
- //
- corrections[icorr]= trackOut->GetParameter()[ptype]-trackIn.GetParameter()[ptype];
- delete trackOut;
+ if (trackOut){
+ AliTrackerBase::PropagateTrackToBxByBz(&trackIn,refX,kMass,3,kTRUE,kMaxSnp);
+ trackOut->Rotate(trackIn.GetAlpha());
+ trackOut->PropagateTo(trackIn.GetX(),AliTrackerBase::GetBz());
+ //
+ corrections[icorr]= trackOut->GetParameter()[ptype]-trackIn.GetParameter()[ptype];
+ delete trackOut;
+ }else{
+ corrections[icorr]=0;
+ }
}
+ Double_t dRdummy=0;
(*pcstream)<<"fit"<<
- Form("%s=",corr->GetName())<<corrections[icorr]; // dump correction value
+ Form("%s=",corr->GetName())<<corrections[icorr]<< // dump correction value
+ Form("dR%s=",corr->GetName())<<dRdummy; // dump dummy correction value not needed for tracks
+ // for points it is neccessary
}
(*pcstream)<<"fit"<<"\n";
}
+void AliTPCCorrection::MakeLaserDistortionTree(TTree* tree, TObjArray *corrArray, Int_t itype){
+ //
+ // Make a laser fit tree for global minimization
+ //
+ const Double_t cutErrY=0.1;
+ const Double_t cutErrZ=0.1;
+ const Double_t kEpsilon=0.00000001;
+ 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("distortion4_0.root");
+ Double_t bz=AliTrackerBase::GetBz();
+ //
+
+ for (Int_t ientry=0; ientry<entries; ientry++){
+ tree->GetEntry(ientry);
+ if (!ltr->GetVecGX()){
+ ltr->UpdatePoints();
+ }
+ TVectorD * delta= (itype==0)? vecdY:vecdZ;
+ TVectorD * err= (itype==0)? veceY:veceZ;
+
+ for (Int_t irow=0; irow<159; irow++){
+ Int_t nentries = 1000;
+ if (veceY->GetMatrixArray()[irow]>cutErrY||veceZ->GetMatrixArray()[irow]>cutErrZ) nentries=0;
+ if (veceY->GetMatrixArray()[irow]<kEpsilon||veceZ->GetMatrixArray()[irow]<kEpsilon) nentries=0;
+ Int_t dtype=4;
+ Double_t phi =(*ltr->GetVecPhi())[irow];
+ Double_t theta =ltr->GetTgl();
+ Double_t mean=delta->GetMatrixArray()[irow];
+ Double_t gx=0,gy=0,gz=0;
+ Double_t snp = (*ltr->GetVecP2())[irow];
+ Double_t rms = 0.1+err->GetMatrixArray()[irow];
+ gx = (*ltr->GetVecGX())[irow];
+ gy = (*ltr->GetVecGY())[irow];
+ gz = (*ltr->GetVecGZ())[irow];
+ Int_t bundle= ltr->GetBundle();
+ Double_t dRrec=0;
+ //
+ // get delta R used in reconstruction
+ AliTPCcalibDB* calib=AliTPCcalibDB::Instance();
+ AliTPCCorrection * correction = calib->GetTPCComposedCorrection();
+ const AliTPCRecoParam * recoParam = calib->GetTransform()->GetCurrentRecoParam();
+ Double_t xyz0[3]={gx,gy,gz};
+ Double_t oldR=TMath::Sqrt(gx*gx+gy*gy);
+ //
+ // old ExB correction
+ //
+ if(recoParam&&recoParam->GetUseExBCorrection()) {
+ Double_t xyz1[3]={gx,gy,gz};
+ calib->GetExB()->Correct(xyz0,xyz1);
+ Double_t newR=TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
+ dRrec=oldR-newR;
+ }
+ if(recoParam&&recoParam->GetUseComposedCorrection()&&correction) {
+ Float_t xyz1[3]={gx,gy,gz};
+ Int_t sector=(gz>0)?0:18;
+ correction->CorrectPoint(xyz1, sector);
+ Double_t newR=TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
+ dRrec=oldR-newR;
+ }
+
+
+ (*pcstream)<<"fit"<<
+ "bz="<<bz<< // magnetic filed used
+ "dtype="<<dtype<< // detector match type
+ "ptype="<<itype<< // parameter type
+ "theta="<<theta<< // theta
+ "phi="<<phi<< // phi
+ "snp="<<snp<< // snp
+ "mean="<<mean<< // mean dist value
+ "rms="<<rms<< // rms
+ "gx="<<gx<< // global position
+ "gy="<<gy<< // global position
+ "gz="<<gz<< // global position
+ "dRrec="<<dRrec<< // delta Radius in reconstruction
+ "id="<<bundle<< //bundle
+ "entries="<<nentries;// number of entries in bin
+ //
+ //
+ Double_t ky = TMath::Tan(TMath::ASin(snp));
+ Int_t ncorr = corrArray->GetEntries();
+ Double_t r0 = TMath::Sqrt(gx*gx+gy*gy);
+ Double_t phi0 = TMath::ATan2(gy,gx);
+ Double_t distortions[1000]={0};
+ Double_t distortionsR[1000]={0};
+ for (Int_t icorr=0; icorr<ncorr; icorr++) {
+ AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
+ Float_t distPoint[3]={gx,gy,gz};
+ Int_t sector= (gz>0)? 0:18;
+ if (r0>80){
+ corr->DistortPoint(distPoint, sector);
+ }
+ Double_t value=distPoint[2]-gz;
+ if (itype==0){
+ Double_t r1 = TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
+ Double_t phi1 = TMath::ATan2(distPoint[1],distPoint[0]);
+ Double_t drphi= r0*(phi1-phi0);
+ Double_t dr = r1-r0;
+ distortions[icorr] = drphi-ky*dr;
+ distortionsR[icorr] = dr;
+ }
+ (*pcstream)<<"fit"<<
+ Form("%s=",corr->GetName())<<distortions[icorr]<< // dump correction value
+ Form("dR%s=",corr->GetName())<<distortionsR[icorr]; // dump correction R value
+ }
+ (*pcstream)<<"fit"<<"\n";
+ }
+ }
+ delete pcstream;
+}
+
+
void AliTPCCorrection::MakeDistortionMap(THnSparse * his0, TTreeSRedirector * const pcstream, const char* hname, Int_t run){
//
--- /dev/null
+ /*
+ gROOT->Macro("~/rootlogon.C")
+ gSystem->AddIncludePath("-I$ALICE_ROOT/STAT")
+ gSystem->AddIncludePath("-I$ALICE_ROOT/TPC")
+ gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros")
+ gSystem->Load("libANALYSIS");
+ gSystem->Load("libTPCcalib");
+
+ .L $ALICE_ROOT/TPC/CalibMacros/MakeGlobalFit.C+
+ MakeGlobalFit();
+
+ rm matching.ps
+ aliroot -b -q $ALICE_ROOT/TPC/CalibMacros/MakeGlobalFit.C | tee globalFit.log
+
+*/
+
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include "THnSparse.h"
+#include "TLatex.h"
+#include "TCanvas.h"
+#include "TLegend.h"
+#include "TSystem.h"
+#include "TFile.h"
+#include "TChain.h"
+#include "TCut.h"
+#include "TH3.h"
+#include "TH2F.h"
+#include "TProfile3D.h"
+#include "TMath.h"
+#include "TVectorD.h"
+#include "TMatrixD.h"
+#include "TStatToolkit.h"
+#include "TTreeStream.h"
+#include "AliExternalTrackParam.h"
+#include "AliESDfriend.h"
+#include "AliTPCcalibTime.h"
+#include "TROOT.h"
+#include "AliXRDPROOFtoolkit.h"
+#include "AliTPCCorrection.h"
+#include "AliTPCExBTwist.h"
+#include "AliTPCGGVoltError.h"
+#include "AliTPCComposedCorrection.h"
+#include "AliTPCExBConical.h"
+#include "TPostScript.h"
+#include "TStyle.h"
+#include "AliTrackerBase.h"
+#include "AliTPCCalibGlobalMisalignment.h"
+#include "AliTPCExBEffective.h"
+#include "TEntryList.h"
+#endif
+
+const char* chDetName[5]={"ITS","TRD", "Vertex", "TOF","Laser"};
+const char* chParName[5]={"rphi","z", "snp", "tan","1/pt"};
+const char* chSideName[2]={"A side","C side"};
+Bool_t enableDet[5]={1,1,1,1,1}; // detector for fit ITS=0, TRD=1, Vertex=2, Laser=4
+Bool_t enableParam[5]={1,0,1,0,1}; //
+Bool_t enableSign=kFALSE;
+Bool_t useEff0=kFALSE;
+Bool_t useEffD=kFALSE;
+Bool_t useEffR=kFALSE;
+//
+TChain *chain = 0;
+TChain *chainRef = 0;
+Bool_t printMatrix=kFALSE;
+TString *fitString=0;
+//
+// cuts
+//
+TCut cut="1";
+TCut cutD="1";
+
+void MakeChain();
+void MakeAliases();
+void MakeCuts();
+void MakeFit(TCut cutCustom);
+TCanvas* DrawFitITS(const char *name);
+TCanvas* DrawFitVertex(const char *name);
+TCanvas* DrawFitLaser(const char *cname);
+TCanvas* DrawCorrdY();
+TCanvas* DrawCorrdSnp();
+TCanvas * DrawFitdY(const char *name);
+TCanvas * DrawFitdSnp(const char *name);
+void PrintMatch();
+TCanvas * MakeComposedCorrection(const char * name);
+
+void MakeAliases(){
+ chain->SetAlias("isITS","(dtype==0||dtype==2)");
+ chain->SetAlias("isTRD","(dtype==1)");
+ chain->SetAlias("isLaser","(dtype==4)");
+ chain->SetAlias("r","sqrt(gx*gx+gy*gy)");
+ chain->SetAlias("r250","(sqrt(gx*gx+gy*gy)/250.)");
+ chain->SetAlias("weight","((dtype==4)*rms*10+rms)"); // downscale laser
+ chain->SetAlias("side","(1+(theta>0)*2)");
+ chain->SetAlias("mdelta","(mean-R.mean-isLaser*((dRrec-R.dRrec)*tan(asin(snp))))");
+ chain->SetAlias("drift","(1-abs(gz/250))");
+ chain->SetAlias("r0","(r-160)/80");
+ chain->SetAlias("delta","(0*gx)");
+}
+
+
+void MakeGlobalFit(){
+ //
+ //
+ //
+ gROOT->Macro("~/rootlogon.C");
+ gROOT->Macro("NimStyle.C");
+ gSystem->AddIncludePath("-I$ALICE_ROOT/STAT");
+ gSystem->AddIncludePath("-I$ALICE_ROOT/TPC");
+ gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
+ gSystem->Load("libANALYSIS");
+ gSystem->Load("libTPCcalib");
+
+ MakeChain();
+ PrintMatch();
+ TPostScript *ps = new TPostScript("matching.ps", 112);
+ gStyle->SetOptTitle(1);
+ gStyle->SetOptStat(0);
+
+ useEff0=kFALSE; useEffD=kFALSE; useEffR=kFALSE;
+ MakeFit("1");
+ ps->NewPage();
+ DrawFitdY("dY-Physical")->Draw();
+ ps->NewPage();
+ DrawFitdSnp("dSnp-Physical")->Draw();
+ ps->NewPage();
+ DrawFitITS("ITS Physical")->Draw();
+ ps->NewPage();
+ DrawFitVertex("Vertex Physical")->Draw();
+ ps->NewPage();
+ DrawFitLaser("Laser Physical")->Draw();
+ ps->NewPage();
+ MakeComposedCorrection("Correction physical")->Draw();
+ //
+ //
+ //
+ //
+ useEff0=kTRUE; useEffD=kTRUE; useEffR=kTRUE; enableSign=kFALSE;
+ ps->NewPage();
+ MakeFit("1");
+ ps->NewPage();
+ DrawFitdY("dY-Physical+Effective ")->Draw();
+ ps->NewPage();
+ DrawFitdSnp("dSnp-Physical+Effective ")->Draw();
+ ps->NewPage();
+ DrawFitITS("ITS Physical+Effective")->Draw();
+ ps->NewPage();
+ DrawFitVertex("Vertex Physical+Effective")->Draw();
+ ps->NewPage();
+ DrawFitLaser("Laser Physical +Effective ")->Draw();
+ ps->NewPage();
+ MakeComposedCorrection("Correction physical+Effective")->Draw();
+ //
+ useEff0=kTRUE; useEffD=kTRUE; useEffR=kTRUE; enableSign=kTRUE;
+ ps->NewPage();
+ MakeFit("1");
+ ps->NewPage();
+ DrawFitdY("dY-Physical+Effective Sign")->Draw();
+ ps->NewPage();
+ DrawFitdSnp("dSnp-Physical+Effective Sign")->Draw();
+ ps->NewPage();
+ DrawFitITS("ITS Physical+Effective Sign")->Draw();
+ ps->NewPage();
+ DrawFitVertex("Vertex Physical+Effective Sign")->Draw();
+ ps->NewPage();
+ DrawFitLaser("Laser Physical +Effective Sign")->Draw();
+ ps->NewPage();
+ MakeComposedCorrection("Correction physical+Effective Sign")->Draw();
+ //
+ ps->Close();
+ delete ps;
+
+}
+
+void MakeChain(){
+ //
+ //
+ TH1::AddDirectory(0);
+ TFile * f =0;
+ TTree * tree=0;
+ TTree * treeF=0;
+ //
+ chain = new TChain("fit","fit");
+ chainRef = new TChain("fit","fit");
+ for (Int_t idet=0; idet<5; idet++){
+ for (Int_t ipar=0; ipar<5; ipar++){
+ char fname[1000];
+ f= TFile::Open(Form("../mergePlus/distortion%d_%d.root",idet,ipar));
+ if (f){
+ tree = (TTree*)f->Get("fit");
+ f= TFile::Open(Form("../mergeField0/distortion%d_%d.root",idet,ipar));
+ if (f && ipar!=4 &&ipar!=3 &&ipar!=1){
+ treeF = (TTree*)f->Get("fit");
+ if (tree && treeF)
+ if (tree->GetEntries()>0 &&(tree->GetEntries()==treeF->GetEntries())){
+ printf("%d\t%d\n",Int_t(tree->GetEntries()),Int_t(treeF->GetEntries()));
+ chain->Add(Form("../mergePlus/distortion%d_%d.root",idet,ipar));
+ chainRef->Add(Form("../mergeField0/distortion%d_%d.root",idet,ipar));
+ }
+ }else{
+ f= TFile::Open(Form("../mergeMinus/distortion%d_%d.root",idet,ipar));
+ if (f){
+ treeF = (TTree*)f->Get("fit");
+ if (tree && treeF)
+ if (tree->GetEntries()>0 &&(tree->GetEntries()==treeF->GetEntries())){
+ printf("%d\t%d\n",Int_t(tree->GetEntries()),Int_t(treeF->GetEntries()));
+ chain->Add(Form("../mergePlus/distortion%d_%d.root",idet,ipar));
+ chainRef->Add(Form("../mergeMinus/distortion%d_%d.root",idet,ipar));
+ }
+ }
+ }
+ }
+ f= TFile::Open(Form("../mergeMinus/distortion%d_%d.root",idet,ipar));
+ if (f){
+ tree = (TTree*)f->Get("fit");
+ f= TFile::Open(Form("../mergeField0/distortion%d_%d.root",idet,ipar));
+ if (f){
+ treeF = (TTree*)f->Get("fit");
+ if (tree && treeF)
+ if (tree->GetEntries()>0 &&(tree->GetEntries()==treeF->GetEntries())){
+ printf("%d\t%d\n",Int_t(tree->GetEntries()),Int_t(treeF->GetEntries()));
+ chain->Add(Form("../mergeMinus/distortion%d_%d.root",idet,ipar));
+ chainRef->Add(Form("../mergeField0/distortion%d_%d.root",idet,ipar));
+ }
+ }
+ else{
+ f= TFile::Open(Form("../mergeMinus/distortion%d_%d.root",idet,ipar));
+ if (f && ipar!=4&&ipar!=3 &&ipar!=1){
+ treeF = (TTree*)f->Get("fit");
+ if (tree && treeF)
+ if (tree->GetEntries()>0 &&(tree->GetEntries()==treeF->GetEntries())){
+ printf("%d\t%d\n",Int_t(tree->GetEntries()),Int_t(treeF->GetEntries()));
+ chain->Add(Form("../mergePlus/distortion%d_%d.root",idet,ipar));
+ chainRef->Add(Form("../mergeMinus/distortion%d_%d.root",idet,ipar));
+ }
+ }
+ }
+ }
+ }
+ }
+ chain->AddFriend(chainRef,"R");
+ MakeAliases();
+ MakeCuts();
+}
+
+
+void MakeCuts(){
+ //
+ //
+ //
+ TCut cutS="((rms>0&&R.rms>0&&entries>0&&R.entries>0))"; // statistic cuts
+ TCut cutType="((dtype==R.dtype)&&(ptype==R.ptype))"; // corresponding types
+ TCut cutOut="(ptype==0)*abs(mdelta)<(0.3+rms)||(ptype==0&&abs(mdelta*85)<(0.3+rms*85))"; // corresponding types
+ //
+
+ chain->Draw(">>outList",cutS+cutType+cutOut+"abs(snp)<0.5","entryList");
+ TEntryList *elistOut = (TEntryList*)gDirectory->Get("outList");
+ chain->SetEntryList(elistOut);
+
+
+}
+
+
+TMatrixD * MakeCorrelation(TMatrixD &matrix){
+ //
+ //
+ //
+ Int_t nrows = matrix.GetNrows();
+ TMatrixD * mat = new TMatrixD(nrows,nrows);
+ for (Int_t irow=0; irow<nrows; irow++)
+ for (Int_t icol=0; icol<nrows; icol++){
+ (*mat)(irow,icol)= matrix(irow,icol)/TMath::Sqrt(matrix(irow,irow)*matrix(icol,icol));
+ }
+ return mat;
+}
+
+void MakeDetCut(){
+ cutD=Form("((dtype==%d)||(dtype==%d)||(dtype==%d)||(dtype==%d)||(dtype==%d))",enableDet[0]?0:0,enableDet[1]?1:0,enableDet[2]?2:0,enableDet[3]?3:0,enableDet[4]?4:0);
+ cutD+=Form("((ptype==%d)||(ptype==%d)||(ptype==%d)||(ptype==%d)||(ptype==%d))",enableParam[0]?0:0,enableParam[1]?1:0,enableParam[2]?2:0,enableParam[3]?3:0,enableParam[4]?4:0);
+}
+
+void MakeFit(TCut cutCustom){
+ MakeDetCut();
+ Int_t npointsMax=30000000;
+ TStatToolkit toolkit;
+ Double_t chi2=0;
+ Int_t npoints=0;
+ TVectorD param;
+ TMatrixD covar;
+ //
+ TString fstring=""; // magnetic part
+ fstring+="(tX1-R.tX1)++"; // twist X
+ fstring+="(tY1-R.tY1)++"; // twist Y
+ fstring+="(sign(bz)*(tX1-R.tX1))++"; // twist X
+ fstring+="(sign(bz)*(tY1-R.tY1))++"; // twist Y
+ fstring+="(exbT1-exb11-(R.exbT1-R.exb11))++"; // T1 adjustment
+ fstring+="(exbT2-exb11-(R.exbT2-R.exb11))++"; // T2 adjustment
+
+ {if (enableDet[0]){
+ fstring+="(isITS*shiftX)++"; // shift X - ITS
+ fstring+="(isITS*shiftY)++"; // shift Y - ITS
+ fstring+="(isITS*sign(bz)*shiftX)++"; // shift X - ITS
+ fstring+="(isITS*sign(bz)*shiftY)++"; // shift Y - ITS
+ }}
+ {if (enableDet[1]){
+ fstring+="(isTRD*shiftX)++"; // shift X - TRD
+ fstring+="(isTRD*shiftY)++"; // shift Y - TRD
+ fstring+="(isTRD*sign(bz)*shiftX)++"; // shift X - TRD
+ fstring+="(isTRD*sign(bz)*shiftY)++"; // shift Y - TRD
+ }}
+ //
+ if (enableDet[4]){
+ }
+ TString fstringCustom="";
+ if (useEff0){
+ fstringCustom+="(eff0_0_0-R.eff0_0_0)++"; // effective correction constant part
+ fstringCustom+="(eff1_0_0-R.eff1_0_0)++"; //
+ }
+ if (useEffD){
+ fstringCustom+="(eff0_0_1-R.eff0_0_1)++"; // effective correction drift part
+ fstringCustom+="(eff1_0_1-R.eff1_0_1)++"; //
+ fstringCustom+="(eff0_0_2-R.eff0_0_2)++"; // effective correction drift part
+ fstringCustom+="(eff1_0_2-R.eff1_0_2)++"; //
+ }
+ if (useEffR){
+ fstringCustom+="(eff0_1_0-R.eff0_1_0)++"; // effective correction radial part
+ fstringCustom+="(eff1_1_0-R.eff1_1_0)++"; //
+ fstringCustom+="(eff0_1_1-R.eff0_1_1)++"; // effective correction radial part
+ fstringCustom+="(eff1_1_1-R.eff1_1_1)++"; //
+ fstringCustom+="(eff0_1_2-R.eff0_1_2)++"; // effective correction radial part
+ fstringCustom+="(eff1_1_2-R.eff1_1_2)++"; //
+ fstringCustom+="(eff0_2_0-R.eff0_2_0)++"; // effective correction radial part
+ fstringCustom+="(eff1_2_0-R.eff1_2_0)++"; //
+ fstringCustom+="(eff0_2_1-R.eff0_2_1)++"; // effective correction radial part
+ fstringCustom+="(eff1_2_1-R.eff1_2_1)++"; //
+ fstringCustom+="(eff0_2_2-R.eff0_2_2)++"; // effective correction radial part
+ fstringCustom+="(eff1_2_2-R.eff1_2_2)++"; //
+ }
+ if (enableSign){
+ if (useEff0){
+ fstringCustom+="sign(bz)*(eff0_0_0-R.eff0_0_0)++"; // effective correction constant part
+ fstringCustom+="sign(bz)*(eff1_0_0-R.eff1_0_0)++"; //
+ }
+ if (useEffD){
+ fstringCustom+="sign(bz)*(eff0_0_1-R.eff0_0_1)++"; // effective correction drift part
+ fstringCustom+="sign(bz)*(eff1_0_1-R.eff1_0_1)++"; //
+ fstringCustom+="sign(bz)*(eff0_0_2-R.eff0_0_2)++"; // effective correction drift part
+ fstringCustom+="sign(bz)*(eff1_0_2-R.eff1_0_2)++"; //
+ }
+ if (useEffR){
+ fstringCustom+="sign(bz)*(eff0_1_0-R.eff0_1_0)++"; // effective correction radial part
+ fstringCustom+="sign(bz)*(eff1_1_0-R.eff1_1_0)++"; //
+ fstringCustom+="sign(bz)*(eff0_1_1-R.eff0_1_1)++"; // effective correction radial part
+ fstringCustom+="sign(bz)*(eff1_1_1-R.eff1_1_1)++"; //
+ fstringCustom+="sign(bz)*(eff0_1_2-R.eff0_1_2)++"; // effective correction radial part
+ fstringCustom+="sign(bz)*(eff1_1_2-R.eff1_1_2)++"; //
+ fstringCustom+="sign(bz)*(eff0_2_0-R.eff0_2_0)++"; // effective correction radial part
+ fstringCustom+="sign(bz)*(eff1_2_0-R.eff1_2_0)++"; //
+ fstringCustom+="sign(bz)*(eff0_2_1-R.eff0_2_1)++"; // effective correction radial part
+ fstringCustom+="sign(bz)*(eff1_2_1-R.eff1_2_1)++"; //
+ fstringCustom+="sign(bz)*(eff0_2_2-R.eff0_2_2)++"; // effective correction radial part
+ fstringCustom+="sign(bz)*(eff1_2_2-R.eff1_2_2)++"; //
+ }
+ }
+ //
+ fstring+=fstringCustom;
+ //
+ TString *strDelta = TStatToolkit::FitPlane(chain,"mdelta:weight", fstring.Data(),cut+cutD+cutCustom, chi2,npoints,param,covar,-1,0, npointsMax, kTRUE);
+ if (printMatrix) MakeCorrelation(covar)->Print();
+ chain->SetAlias("delta",strDelta->Data());
+ strDelta->Tokenize("++")->Print();
+ fitString = strDelta;
+ PrintMatch();
+}
+
+
+
+
+
+void PrintMatch(){
+ //
+ // Print detector matching info
+ //
+ for (Int_t ipar=0; ipar<5; ipar++){
+ for (Int_t idet=0; idet<5; idet++){
+ Double_t mean0,rms0,mean1,rms1;
+ Int_t entries = chain->Draw("delta-mdelta",cut+Form("dtype==%d&&ptype==%d",idet,ipar),"");
+ if (entries==0) continue;
+ TH1 * his = (TH1*)(chain->GetHistogram()->Clone());
+ mean1=his->GetMean();
+ rms1 =his->GetRMS();
+ delete his;
+ //
+ entries = chain->Draw("mdelta",cut+Form("dtype==%d&&ptype==%d",idet,ipar),"");
+ if (entries==0) continue;
+ his = (TH1*)(chain->GetHistogram()->Clone());
+ //
+ mean0=his->GetMean();
+ rms0 =his->GetRMS();
+
+ printf("ptype==%s\tdtype==%s\tMean=%f -> %f\tRMS=%f -> %f\n",chParName[ipar],chDetName[idet], mean0,mean1,rms0,rms1);
+ delete his;
+ }
+ }
+}
+
+
+
+
+TCanvas* DrawFitITS(const char *name){
+ //
+ //
+ //
+ TLegend *legend=0;
+ TCanvas *canvas = new TCanvas(name,name,800,800);
+ canvas->Divide(1,2);
+ Int_t entries=0;
+ const char * grname[10]={"A side (B=0.5T)","A side (B=-0.5T)", "C side (B=0.5T)", "C side (B=-0.5T)",0};
+
+ TGraph *graphsdyITS[10];
+ TGraph *graphsdyITSc[10];
+ entries=chain->Draw("mdelta:phi",cut+"ptype==0&&dtype==0&&bz>0&&theta>0","goff");
+ graphsdyITS[0]=new TGraph(entries,chain->GetV2(), chain->GetV1());
+ entries=chain->Draw("mdelta:phi",cut+cut+"ptype==0&&dtype==0&&bz<0&&theta>0","goff");
+ graphsdyITS[1]=new TGraph(entries,chain->GetV2(), chain->GetV1());
+ entries=chain->Draw("mdelta:phi",cut+cut+"ptype==0&&dtype==0&&bz>0&&theta<0","goff");
+ graphsdyITS[2]=new TGraph(entries,chain->GetV2(), chain->GetV1());
+ entries=chain->Draw("mdelta:phi",cut+cut+"ptype==0&&dtype==0&&bz<0&&theta<0","goff");
+ graphsdyITS[3]=new TGraph(entries,chain->GetV2(), chain->GetV1());
+ //
+ entries=chain->Draw("mdelta-delta:phi",cut+cut+"ptype==0&&dtype==0&&bz>0&&theta>0","goff");
+ graphsdyITSc[0]=new TGraph(entries,chain->GetV2(), chain->GetV1());
+ entries=chain->Draw("mdelta-delta:phi",cut+cut+"ptype==0&&dtype==0&&bz<0&&theta>0","goff");
+ graphsdyITSc[1]=new TGraph(entries,chain->GetV2(), chain->GetV1());
+ entries=chain->Draw("mdelta-delta:phi",cut+cut+"ptype==0&&dtype==0&&bz>0&&theta<0","goff");
+ graphsdyITSc[2]=new TGraph(entries,chain->GetV2(), chain->GetV1());
+ entries=chain->Draw("mdelta-delta:phi",cut+cut+"ptype==0&&dtype==0&&bz<0&&theta<0","goff");
+ graphsdyITSc[3]=new TGraph(entries,chain->GetV2(), chain->GetV1());
+
+ canvas->cd(1);
+ legend = new TLegend(0.6,0.6,1.0,1.0,"#Delta R#phi TPC-ITS");
+ for (Int_t i=0; i<4; i++){
+ graphsdyITS[i]->SetMaximum(0.5);
+ graphsdyITS[i]->SetMinimum(-0.5);
+ graphsdyITS[i]->SetMarkerColor(i+1);
+ graphsdyITS[i]->SetMarkerStyle(25+i);
+ // graphsdyITS[i]->SetName(grname[i]);
+ graphsdyITS[i]->GetXaxis()->SetTitle("#phi");
+ graphsdyITS[i]->GetYaxis()->SetTitle("#Delta R#Phi (cm)");
+ if (i==0) graphsdyITS[i]->Draw("ap");
+ graphsdyITS[i]->Draw("p");
+ legend->AddEntry(graphsdyITS[i],grname[i]);
+ }
+ legend->Draw();
+ canvas->cd(2);
+ legend = new TLegend(0.6,0.6,1.0,1.0,"#Delta R#phi TPC-ITS corrected");
+ for (Int_t i=0; i<4; i++){
+ graphsdyITSc[i]->SetMaximum(0.5);
+ graphsdyITSc[i]->SetMinimum(-0.5);
+ graphsdyITSc[i]->SetMarkerColor(i+1);
+ graphsdyITSc[i]->SetMarkerStyle(25+i);
+ // graphsdyITS[i]->SetName(grname[i]);
+ graphsdyITSc[i]->GetXaxis()->SetTitle("#phi");
+ graphsdyITSc[i]->GetYaxis()->SetTitle("#Delta R#Phi (cm)");
+ if (i==0) graphsdyITS[i]->Draw("ap");
+ graphsdyITSc[i]->Draw("p");
+ legend->AddEntry(graphsdyITSc[i],grname[i]);
+ }
+ legend->Draw();
+ return canvas;
+}
+
+
+TCanvas* DrawFitLaser(const char *cname){
+ //
+ //
+ //
+ TH1::AddDirectory(0);
+ TCut cutLaser=cut+"isLaser&&bz<0";
+ TCanvas *canvas= new TCanvas(cname, cname,800,800);
+ canvas->Divide(2,2);
+ canvas->Draw();
+ TLegend*legend=0;
+ //
+ TH1 *his[16]={0};
+ {for (Int_t icorr=0; icorr<2; icorr++){
+ for (Int_t iside=0; iside<2; iside++){
+ canvas->cd(iside+1);
+ for (Int_t ib=0; ib<4; ib++){
+ TCut cutB="1";
+ if (iside==0) cutB=Form("(id==%d&&gz>0)",ib);
+ if (iside==1) cutB=Form("(id==%d&&gz<0)",ib);
+ //cutB.Print();
+ if (icorr==0) chain->Draw("10*mdelta:r",cutLaser+cutB,"prof");
+ if (icorr==1) chain->Draw("10*(mdelta-delta):r",cutLaser+cutB,"prof");
+ his[icorr*8+iside*4+ib] = (TH1*)(chain->GetHistogram()->Clone());
+ his[icorr*8+iside*4+ib]->SetName(Form("B%d%d%d",icorr,iside,ib));
+ his[icorr*8+iside*4+ib]->SetTitle(Form("Bundle %d",ib));
+ his[icorr*8+iside*4+ib]->SetMarkerColor(ib+1);
+ his[icorr*8+iside*4+ib]->SetMarkerStyle(ib+25);
+ his[icorr*8+iside*4+ib]->SetMarkerSize(0.4);
+ his[icorr*8+iside*4+ib]->SetMaximum(3);
+ his[icorr*8+iside*4+ib]->SetMinimum(-3);
+ his[icorr*8+iside*4+ib]->GetXaxis()->SetTitle("r (cm)");
+ his[icorr*8+iside*4+ib]->GetYaxis()->SetTitle("#Delta r#phi (mm)");
+ }
+ }
+ }}
+ //
+ for (Int_t icorr=0; icorr<2; icorr++){
+ for (Int_t iside=0; iside<2; iside++){
+ canvas->cd(icorr*2+iside+1);
+ legend = new TLegend(0.6,0.6,1.0,1.0,Form("#Delta R#phi Laser-%s",chSideName[iside]));
+ for (Int_t ib=0; ib<4; ib++){
+ if (ib==0) his[icorr*2+iside*4+ib]->Draw();
+ his[icorr*8+iside*4+ib]->Draw("same");
+ legend->AddEntry(his[icorr*8+iside*4+ib]);
+ }
+ legend->Draw();
+ }
+ }
+ return canvas;
+}
+
+
+
+
+
+
+
+
+
+
+
+TCanvas* DrawFitVertex(const char *name){
+ //
+ //
+ //
+ TLegend *legend=0;
+ TCanvas *canvas = new TCanvas(name,name,800,800);
+ canvas->Divide(1,2);
+ Int_t entries=0;
+ const char * grname[10]={"A side (B=0.5T)","A side (B=-0.5T)", "C side (B=0.5T)", "C side (B=-0.5T)",0};
+
+ TGraph *graphsdyVertex[10];
+ TGraph *graphsdyVertexc[10];
+ entries=chain->Draw("mdelta:phi",cut+cut+"ptype==0&&dtype==2&&bz>0&&theta>0","goff");
+ graphsdyVertex[0]=new TGraph(entries,chain->GetV2(), chain->GetV1());
+ entries=chain->Draw("mdelta:phi",cut+cut+"ptype==0&&dtype==2&&bz<0&&theta>0","goff");
+ graphsdyVertex[1]=new TGraph(entries,chain->GetV2(), chain->GetV1());
+ entries=chain->Draw("mdelta:phi",cut+cut+"ptype==0&&dtype==2&&bz>0&&theta<0","goff");
+ graphsdyVertex[2]=new TGraph(entries,chain->GetV2(), chain->GetV1());
+ entries=chain->Draw("mdelta:phi",cut+cut+"ptype==0&&dtype==2&&bz<0&&theta<0","goff");
+ graphsdyVertex[3]=new TGraph(entries,chain->GetV2(), chain->GetV1());
+ //
+ entries=chain->Draw("mdelta-delta:phi",cut+cut+"ptype==0&&dtype==2&&bz>0&&theta>0","goff");
+ graphsdyVertexc[0]=new TGraph(entries,chain->GetV2(), chain->GetV1());
+ entries=chain->Draw("mdelta-delta:phi",cut+cut+"ptype==0&&dtype==2&&bz<0&&theta>0","goff");
+ graphsdyVertexc[1]=new TGraph(entries,chain->GetV2(), chain->GetV1());
+ entries=chain->Draw("mdelta-delta:phi",cut+cut+"ptype==0&&dtype==2&&bz>0&&theta<0","goff");
+ graphsdyVertexc[2]=new TGraph(entries,chain->GetV2(), chain->GetV1());
+ entries=chain->Draw("mdelta-delta:phi",cut+cut+"ptype==0&&dtype==2&&bz<0&&theta<0","goff");
+ graphsdyVertexc[3]=new TGraph(entries,chain->GetV2(), chain->GetV1());
+
+ canvas->cd(1);
+ legend = new TLegend(0.6,0.6,1.0,1.0,"#Delta R#phi TPC-Vertex");
+ for (Int_t i=0; i<4; i++){
+ graphsdyVertex[i]->SetMaximum(1);
+ graphsdyVertex[i]->SetMinimum(-1);
+ graphsdyVertex[i]->SetMarkerColor(i+1);
+ graphsdyVertex[i]->SetMarkerStyle(25+i);
+ // graphsdyVertex[i]->SetName(grname[i]);
+ graphsdyVertex[i]->GetXaxis()->SetTitle("#phi");
+ graphsdyVertex[i]->GetYaxis()->SetTitle("#Delta R#Phi (cm)");
+ if (i==0) graphsdyVertex[i]->Draw("ap");
+ graphsdyVertex[i]->Draw("p");
+ legend->AddEntry(graphsdyVertex[i],grname[i]);
+ }
+ legend->Draw();
+ canvas->cd(2);
+ legend = new TLegend(0.6,0.6,1.0,1.0,"#Delta R#phi TPC-Vertex corrected");
+ for (Int_t i=0; i<4; i++){
+ graphsdyVertexc[i]->SetMaximum(1);
+ graphsdyVertexc[i]->SetMinimum(-1);
+ graphsdyVertexc[i]->SetMarkerColor(i+1);
+ graphsdyVertexc[i]->SetMarkerStyle(25+i);
+ // graphsdyVertex[i]->SetName(grname[i]);
+ graphsdyVertexc[i]->GetXaxis()->SetTitle("#phi");
+ graphsdyVertexc[i]->GetYaxis()->SetTitle("#Delta R#Phi (cm)");
+ if (i==0) graphsdyVertex[i]->Draw("ap");
+ graphsdyVertexc[i]->Draw("p");
+ legend->AddEntry(graphsdyVertexc[i],grname[i]);
+ }
+ legend->Draw();
+ return canvas;
+}
+
+
+
+
+
+
+TCanvas* DrawCorrdY(){
+
+ TLegend *legend=0;
+ TCanvas *canvas = new TCanvas("Corrections","Corrections",800,800);
+ canvas->Divide(1,3);
+
+ TGraph *graphsITS[10];
+ TGraph *graphsTRD[10];
+ TGraph *graphsVertex[10];
+ Int_t entries=0;
+ const char * grname[10]={"X twist (1 mrad)","Y twist (1 mrad)", "X shift (1 mm)", "Y shift (1 mm)", "drPhi (1 mm)"};
+ //
+ entries=chain->Draw("tX1:phi",cut+cut+"ptype==0&&dtype==0&&bz>0","goff");
+ graphsITS[0]=new TGraph(entries,chain->GetV2(), chain->GetV1());
+ entries=chain->Draw("tY1:phi",cut+cut+"ptype==0&&dtype==0&&bz>0","goff");
+ graphsITS[1]=new TGraph(entries,chain->GetV2(), chain->GetV1());
+ entries=chain->Draw("shiftX:phi",cut+cut+"ptype==0&&dtype==0&&bz>0","goff");
+ graphsITS[2]=new TGraph(entries,chain->GetV2(), chain->GetV1());
+ entries=chain->Draw("shiftY:phi",cut+cut+"ptype==0&&dtype==0&&bz>0","goff");
+ graphsITS[3]=new TGraph(entries,chain->GetV2(), chain->GetV1());
+ entries=chain->Draw("drPhiA+drPhiC:phi",cut+cut+"ptype==0&&dtype==0&&bz>0","goff");
+ graphsITS[4]=new TGraph(entries,chain->GetV2(), chain->GetV1());
+ //
+ entries=chain->Draw("tX1:phi",cut+cut+"ptype==0&&dtype==1&&bz>0","goff");
+ graphsTRD[0]=new TGraph(entries,chain->GetV2(), chain->GetV1());
+ entries=chain->Draw("tY1:phi",cut+cut+"ptype==0&&dtype==1&&bz>0","goff");
+ graphsTRD[1]=new TGraph(entries,chain->GetV2(), chain->GetV1());
+ entries=chain->Draw("shiftX:phi",cut+cut+"ptype==0&&dtype==1&&bz>0","goff");
+ graphsTRD[2]=new TGraph(entries,chain->GetV2(), chain->GetV1());
+ entries=chain->Draw("shiftY:phi",cut+cut+"ptype==0&&dtype==1&&bz>0","goff");
+ graphsTRD[3]=new TGraph(entries,chain->GetV2(), chain->GetV1());
+ entries=chain->Draw("drPhiA+drPhiC:phi",cut+cut+"ptype==0&&dtype==1&&bz>0","goff");
+ graphsTRD[4]=new TGraph(entries,chain->GetV2(), chain->GetV1());
+ //
+ entries=chain->Draw("tX1:phi",cut+cut+"ptype==0&&dtype==2&&bz>0","goff");
+ graphsVertex[0]=new TGraph(entries,chain->GetV2(), chain->GetV1());
+ entries=chain->Draw("tY1:phi",cut+cut+"ptype==0&&dtype==2&&bz>0","goff");
+ graphsVertex[1]=new TGraph(entries,chain->GetV2(), chain->GetV1());
+ entries=chain->Draw("shiftX:phi",cut+cut+"ptype==0&&dtype==2&&bz>0","goff");
+ graphsVertex[2]=new TGraph(entries,chain->GetV2(), chain->GetV1());
+ entries=chain->Draw("shiftY:phi",cut+cut+"ptype==0&&dtype==2&&bz>0","goff");
+ graphsVertex[3]=new TGraph(entries,chain->GetV2(), chain->GetV1());
+ entries=chain->Draw("drPhiA+drPhiC:phi",cut+cut+"ptype==0&&dtype==2&&bz>0","goff");
+ graphsVertex[4]=new TGraph(entries,chain->GetV2(), chain->GetV1());
+
+ canvas->cd(1);
+ //
+ legend = new TLegend(0.6,0.6,1.0,1.0,"#Delta R#Phi TPC-Vertex");
+ for (Int_t i=0; i<5; i++){
+ graphsVertex[i]->SetMaximum(0.2);
+ graphsVertex[i]->SetMinimum(-0.2);
+ graphsVertex[i]->SetMarkerColor(i+1);
+ graphsVertex[i]->SetMarkerStyle(25+i);
+ graphsVertex[i]->SetName(grname[i]);
+ graphsVertex[i]->GetXaxis()->SetTitle("#phi");
+ graphsVertex[i]->GetYaxis()->SetTitle("#DeltaR#Phi (cm)");
+ if (i==0) graphsVertex[i]->Draw("ap");
+ graphsVertex[i]->Draw("p");
+ legend->AddEntry(graphsITS[i],grname[i]);
+ }
+ legend->Draw();
+
+ canvas->cd(2);
+ legend = new TLegend(0.6,0.6,1.0,1.0,"#Delta R#Phi ITS-TPC");
+ for (Int_t i=0; i<5; i++){
+ graphsITS[i]->SetMaximum(0.2);
+ graphsITS[i]->SetMinimum(-0.2);
+ graphsITS[i]->SetMarkerColor(i+1);
+ graphsITS[i]->SetMarkerStyle(25+i);
+ graphsITS[i]->SetName(grname[i]);
+ graphsITS[i]->GetXaxis()->SetTitle("#phi");
+ graphsITS[i]->GetYaxis()->SetTitle("#Delta R#Phi (cm)");
+ if (i==0) graphsITS[i]->Draw("ap");
+ graphsITS[i]->Draw("p");
+ legend->AddEntry(graphsITS[i],grname[i]);
+ }
+ legend->Draw();
+ //
+ canvas->cd(3);
+ legend = new TLegend(0.6,0.6,1.0,1.0,"#Delta R#Phi TPC-TRD");
+ for (Int_t i=0; i<5; i++){
+ graphsTRD[i]->SetMaximum(0.2);
+ graphsTRD[i]->SetMinimum(-0.2);
+ graphsTRD[i]->SetMarkerColor(i+1);
+ graphsTRD[i]->SetMarkerStyle(25+i);
+ graphsTRD[i]->SetName(grname[i]);
+ graphsTRD[i]->GetXaxis()->SetTitle("#phi");
+ graphsTRD[i]->GetYaxis()->SetTitle("#DeltaR#Phi (cm)");
+ if (i==0) graphsTRD[i]->Draw("ap");
+ graphsTRD[i]->Draw("p");
+ legend->AddEntry(graphsITS[i],grname[i]);
+ }
+ legend->Draw();
+ return canvas;
+}
+
+
+TCanvas * DrawCorrdSnp(){
+
+ TLegend *legend=0;
+ TCanvas *canvas = new TCanvas("Corrections dSnp","Corrections dSnp",800,800);
+ canvas->Divide(1,3);
+
+ TGraph *graphsITS[10];
+ TGraph *graphsTRD[10];
+ TGraph *graphsVertex[10];
+ Int_t entries=0;
+ const char * grname[10]={"X twist (1 mrad)","Y twist (1 mrad)", "X shift (1 mm)", "Y shift (1 mm)",0};
+ //
+ entries=chain->Draw("1000*tX1:phi",cut+cut+"ptype==2&&dtype==0&&bz>0","goff");
+ graphsITS[0]=new TGraph(entries,chain->GetV2(), chain->GetV1());
+ entries=chain->Draw("1000*tY1:phi",cut+cut+"ptype==2&&dtype==0&&bz>0","goff");
+ graphsITS[1]=new TGraph(entries,chain->GetV2(), chain->GetV1());
+ entries=chain->Draw("1000*shiftX:phi",cut+cut+"ptype==2&&dtype==0&&bz>0","goff");
+ graphsITS[2]=new TGraph(entries,chain->GetV2(), chain->GetV1());
+ entries=chain->Draw("1000*shiftY:phi",cut+cut+"ptype==2&&dtype==0&&bz>0","goff");
+ graphsITS[3]=new TGraph(entries,chain->GetV2(), chain->GetV1());
+ //
+ entries=chain->Draw("1000*tX1:phi",cut+cut+"ptype==2&&dtype==1&&bz>0","goff");
+ graphsTRD[0]=new TGraph(entries,chain->GetV2(), chain->GetV1());
+ entries=chain->Draw("1000*tY1:phi",cut+cut+"ptype==2&&dtype==1&&bz>0","goff");
+ graphsTRD[1]=new TGraph(entries,chain->GetV2(), chain->GetV1());
+ entries=chain->Draw("1000*shiftX:phi",cut+cut+"ptype==2&&dtype==1&&bz>0","goff");
+ graphsTRD[2]=new TGraph(entries,chain->GetV2(), chain->GetV1());
+ entries=chain->Draw("1000*shiftY:phi",cut+cut+"ptype==2&&dtype==1&&bz>0","goff");
+ graphsTRD[3]=new TGraph(entries,chain->GetV2(), chain->GetV1());
+ //
+ entries=chain->Draw("1000*tX1:phi",cut+cut+"ptype==2&&dtype==2&&bz>0","goff");
+ graphsVertex[0]=new TGraph(entries,chain->GetV2(), chain->GetV1());
+ entries=chain->Draw("1000*tY1:phi",cut+cut+"ptype==2&&dtype==2&&bz>0","goff");
+ graphsVertex[1]=new TGraph(entries,chain->GetV2(), chain->GetV1());
+ entries=chain->Draw("1000*shiftX:phi",cut+cut+"ptype==2&&dtype==2&&bz>0","goff");
+ graphsVertex[2]=new TGraph(entries,chain->GetV2(), chain->GetV1());
+ entries=chain->Draw("1000*shiftY:phi",cut+cut+"ptype==2&&dtype==2&&bz>0","goff");
+ graphsVertex[3]=new TGraph(entries,chain->GetV2(), chain->GetV1());
+
+ canvas->cd(1);
+ //
+ legend = new TLegend(0.6,0.6,1.0,1.0,"#Delta sin(#phi) TPC-Vertex");
+ for (Int_t i=0; i<4; i++){
+ graphsVertex[i]->SetMaximum(2);
+ graphsVertex[i]->SetMinimum(2);
+ graphsVertex[i]->SetMarkerColor(i+1);
+ graphsVertex[i]->SetMarkerStyle(25+i);
+ graphsVertex[i]->SetName(grname[i]);
+ graphsVertex[i]->GetXaxis()->SetTitle("#phi");
+ graphsVertex[i]->GetYaxis()->SetTitle("#Deltasin(#phi) (mrad)");
+ if (i==0) graphsVertex[i]->Draw("ap");
+ graphsVertex[i]->Draw("p");
+ legend->AddEntry(graphsITS[i],grname[i]);
+ }
+ legend->Draw();
+
+ canvas->cd(2);
+ legend = new TLegend(0.6,0.6,1.0,1.0,"#Delta sin(#phi) ITS-TPC");
+ for (Int_t i=0; i<4; i++){
+ graphsITS[i]->SetMaximum(1);
+ graphsITS[i]->SetMinimum(-1);
+ graphsITS[i]->SetMarkerColor(i+1);
+ graphsITS[i]->SetMarkerStyle(25+i);
+ graphsITS[i]->SetName(grname[i]);
+ graphsITS[i]->GetXaxis()->SetTitle("#phi");
+ graphsITS[i]->GetYaxis()->SetTitle("#Delta sin(#phi) (mrad)");
+ if (i==0) graphsITS[i]->Draw("ap");
+ graphsITS[i]->Draw("p");
+ legend->AddEntry(graphsITS[i],grname[i]);
+ }
+ legend->Draw();
+ //
+ canvas->cd(3);
+ legend = new TLegend(0.6,0.6,1.0,1.0,"#Delta sin(#phi) TPC-TRD");
+ for (Int_t i=0; i<4; i++){
+ graphsTRD[i]->SetMaximum(1);
+ graphsTRD[i]->SetMinimum(-1);
+ graphsTRD[i]->SetMarkerColor(i+1);
+ graphsTRD[i]->SetMarkerStyle(25+i);
+ graphsTRD[i]->SetName(grname[i]);
+ graphsTRD[i]->GetXaxis()->SetTitle("#phi");
+ graphsTRD[i]->GetYaxis()->SetTitle("#Deltasin(#phi) (mrad)");
+ if (i==0) graphsTRD[i]->Draw("ap");
+ graphsTRD[i]->Draw("p");
+ legend->AddEntry(graphsITS[i],grname[i]);
+ }
+ legend->Draw();
+ return canvas;
+ }
+
+
+TCanvas * DrawFitdY(const char *name){
+ //
+ //
+ //
+ TH1::AddDirectory(0);
+ TCanvas *canvas = new TCanvas(name,name,800,800);
+ canvas->Divide(3,5);
+ for (Int_t idet=0; idet<5; idet++){
+ chain->SetMarkerColor(4);
+ chain->SetMarkerStyle(25);
+ chain->SetMarkerSize(0.3);
+ chain->SetLineColor(2);
+ //
+ canvas->cd(idet*3+1);
+ chain->Draw("mdelta:delta",cut+Form("ptype==0&&dtype==%d",idet),"");
+ //
+ canvas->cd(idet*3+2);
+ chain->SetLineColor(2);
+ chain->Draw("mdelta-delta",cut+Form("ptype==0&&dtype==%d",idet),"");
+ chain->SetLineColor(1);
+ chain->Draw("mdelta",cut+Form("ptype==0&&dtype==%d",idet),"same");
+ //
+ canvas->cd(idet*3+3);
+ chain->SetMarkerColor(1);
+ chain->Draw("mdelta:phi",cut+Form("ptype==0&&dtype==%d&&theta>0&&bz>0",idet),"");
+ chain->SetMarkerColor(2);
+ chain->Draw("mdelta-delta:phi",cut+Form("ptype==0&&dtype==%d&&theta>0&&bz>0",idet),"same");
+
+ }
+ return canvas;
+}
+
+TCanvas * DrawFitdSnp(const char *name){
+ //
+ //
+ //
+ TH1::AddDirectory(0);
+ TCanvas *canvas = new TCanvas(name,name,800,800);
+ canvas->Divide(3,5);
+ {for (Int_t idet=0; idet<5; idet++){
+ chain->SetMarkerColor(4);
+ chain->SetMarkerStyle(25);
+ chain->SetMarkerSize(0.3);
+ chain->SetLineColor(2);
+ //
+ canvas->cd(idet*3+1);
+ chain->Draw("1000*mdelta:1000*delta",cut+Form("ptype==2&&dtype==%d",idet),"");
+ //
+ canvas->cd(idet*3+2);
+ chain->SetLineColor(2);
+ chain->Draw("1000*(mdelta-delta)",cut+Form("ptype==2&&dtype==%d",idet),"");
+ chain->SetLineColor(1);
+ chain->Draw("1000*mdelta",cut+Form("ptype==2&&dtype==%d",idet),"same");
+ //
+ canvas->cd(idet*3+3);
+ chain->SetMarkerColor(1);
+ chain->Draw("1000*(mdelta):phi",cut+Form("ptype==2&&dtype==%d&&theta>0&&bz>0",idet),"");
+ chain->SetMarkerColor(2);
+ chain->Draw("1000*(mdelta-delta):phi",cut+Form("ptype==2&&dtype==%d&&theta>0&&bz>0",idet),"same");
+ }}
+ return canvas;
+}
+
+TCanvas * MakeComposedCorrection(const char *name){
+
+ TString fit = chain->GetAlias("delta");
+ TObjArray * array = fit.Tokenize("++");
+ Int_t nfun=array->GetEntries();
+ Double_t wt = 0.3 ;
+ Double_t T1 = 0.9;
+ Double_t T2 = 1.5;
+ //
+ // sign independent correction
+ //
+ AliTPCExBEffective *eff = new AliTPCExBEffective;
+ AliTPCCalibGlobalMisalignment *shiftITS = new AliTPCCalibGlobalMisalignment;
+ AliTPCExBTwist *twist = new AliTPCExBTwist;
+ //
+ // sign dependent correction
+ //
+ AliTPCExBEffective *effS = new AliTPCExBEffective;
+ AliTPCCalibGlobalMisalignment *shiftITSS = new AliTPCCalibGlobalMisalignment;
+ AliTPCExBTwist *twistS = new AliTPCExBTwist;
+
+ TMatrixD polA(100,4);
+ TMatrixD polC(100,4);
+ TMatrixD valA(100,1);
+ TMatrixD valC(100,1);
+ Int_t counterA=0;
+ Int_t counterC=0;
+
+ {
+ for (Int_t i=1; i<nfun;i++){
+ TObjString fitName=array->At(i)->GetName();
+ TObjString fitVal= &(fitName.String()[1+fitName.String().Last('(')]);
+ Double_t value= fitVal.String().Atof();
+ if (fitName.String().Contains("sign")) continue;
+ if (fitName.String().Contains("tX1")){
+ twist->SetXTwist(value*0.001);
+ }
+ if (fitName.String().Contains("tY1")){
+ twist->SetYTwist(value*0.001);
+ }
+
+ if (fitName.String().Contains("shiftX")&&fitName.String().Contains("isITS")){
+ shiftITS->SetXShift(value*0.1);
+ }
+ if (fitName.String().Contains("shiftY")&&fitName.String().Contains("isITS")){
+ shiftITS->SetYShift(value*0.1);
+ }
+
+ if (fitName.String().Contains("eff")){
+ Int_t index=fitName.String().First("_")-1;
+ Int_t side=atoi(&(fitName.String()[index]));
+ Int_t px =atoi(&(fitName.String()[index+2]));
+ Int_t pd =atoi(&(fitName.String()[index+4]));
+ Int_t pp =0;
+ //printf("%s\t%d\t%d\t%d\t%d\t%f\n",fitName.GetName(),side,px,pd,pp, value);
+ if (side==0){
+ polA(counterA,0)=0;
+ polA(counterA,1)=px;
+ polA(counterA,2)=pd;
+ polA(counterA,3)=pp;
+ valA(counterA,0)=value*0.1;
+ counterA++;
+ }
+ if (side==1){
+ polC(counterC,0)=0;
+ polC(counterC,1)=px;
+ polC(counterC,2)=pd;
+ polC(counterC,3)=pp;
+ valC(counterC,0)=value*0.1;
+ counterC++;
+ }
+ }
+ }
+ }
+ polA.ResizeTo(counterA,4);
+ polC.ResizeTo(counterC,4);
+ valA.ResizeTo(counterA,1);
+ valC.ResizeTo(counterC,1);
+ eff->SetPolynoms(&polA,&polC);
+ eff->SetCoeficients(&valA,&valC);
+ eff->SetOmegaTauT1T2(wt,T1,T2);
+ shiftITS->SetOmegaTauT1T2(wt,T1,T2);
+ twist->SetOmegaTauT1T2(wt,T1,T2);
+ //
+ //
+ //
+ {
+ counterA=0;
+ counterC=0;
+ for (Int_t i=1; i<nfun;i++){
+ TObjString fitName=array->At(i)->GetName();
+ TObjString fitVal= &(fitName.String()[1+fitName.String().Last('(')]);
+ if (!fitName.String().Contains("sign")) continue;
+ Double_t value= fitVal.String().Atof();
+ if (fitName.String().Contains("tX1")){
+ twistS->SetXTwist(value*0.001);
+ }
+ if (fitName.String().Contains("tY1")){
+ twistS->SetYTwist(value*0.001);
+ }
+
+ if (fitName.String().Contains("shiftX")&&fitName.String().Contains("isITS")){
+ shiftITSS->SetXShift(value*0.1);
+ }
+ if (fitName.String().Contains("shiftY")&&fitName.String().Contains("isITS")){
+ shiftITSS->SetYShift(value*0.1);
+ }
+
+ if (fitName.String().Contains("eff")){
+ Int_t index=fitName.String().First("_")-1;
+ Int_t side=atoi(&(fitName.String()[index]));
+ Int_t px =atoi(&(fitName.String()[index+2]));
+ Int_t pd =atoi(&(fitName.String()[index+4]));
+ Int_t pp =0;
+ //printf("%s\t%d\t%d\t%d\t%d\t%f\n",fitName.GetName(),side,px,pd,pp, value);
+ if (side==0){
+ polA(counterA,0)=0;
+ polA(counterA,1)=px;
+ polA(counterA,2)=pd;
+ polA(counterA,3)=pp;
+ valA(counterA,0)=value*0.1;
+ counterA++;
+ }
+ if (side==1){
+ polC(counterC,0)=0;
+ polC(counterC,1)=px;
+ polC(counterC,2)=pd;
+ polC(counterC,3)=pp;
+ valC(counterC,0)=value*0.1;
+ counterC++;
+ }
+ }
+ }
+ }
+ polA.ResizeTo(counterA,4);
+ polC.ResizeTo(counterC,4);
+ valA.ResizeTo(counterA,1);
+ valC.ResizeTo(counterC,1);
+ effS->SetPolynoms(&polA,&polC);
+ effS->SetCoeficients(&valA,&valC);
+ effS->SetOmegaTauT1T2(wt,T1,T2);
+ shiftITSS->SetOmegaTauT1T2(wt,T1,T2);
+ twistS->SetOmegaTauT1T2(wt,T1,T2);
+
+
+
+ //
+ TCanvas *c = new TCanvas(name,name,800,800);
+ c->Divide(3,2);
+ c->cd(1);
+ shiftITS->CreateHistoDRPhiinXY()->Draw("surf2");
+ c->cd(2);
+ twist->CreateHistoDRPhiinXY()->Draw("surf2");
+ c->cd(3);
+ eff->CreateHistoDRPhiinZR()->Draw("surf2");
+ c->cd(4);
+ shiftITSS->CreateHistoDRPhiinXY()->Draw("surf2");
+ c->cd(5);
+ twistS->CreateHistoDRPhiinXY()->Draw("surf2");
+ c->cd(6);
+ effS->CreateHistoDRPhiinZR()->Draw("surf2");
+ TObjArray * corr = new TObjArray;
+
+ return c;
+}
/*
Example usage:
-
+ gROOT->Macro("~/rootlogon.C")
gSystem->AddIncludePath("-I$ALICE_ROOT/STAT")
gSystem->AddIncludePath("-I$ALICE_ROOT/TPC")
gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros")
gSystem->Load("libANALYSIS");
gSystem->Load("libTPCcalib");
.L $ALICE_ROOT/TPC/CalibMacros/MakeLookup.C+
- Int_t run=11111
+ Int_t run=115888
+ // .x $ALICE_ROOT/ANALYSIS/CalibMacros/Pass0/ConfigCalibTrain.C(run,"local:///lustre/alice/alien/alice/data/2010/OCDB")
+
MakeLookup(run,0);
//Local check of procedure
#include "AliTPCGGVoltError.h"
#include "AliTPCComposedCorrection.h"
#include "AliTPCExBConical.h"
+#include "TPostScript.h"
+#include "TStyle.h"
#include "AliTrackerBase.h"
+#include "AliTPCCalibGlobalMisalignment.h"
+#include "AliTPCExBEffective.h"
+#include "AliTPCExBBShape.h"
#endif
TChain * chain=0;
-void MakeLookup(Int_t run, Int_t mode=0);
+void MakeLookup(Int_t run, Int_t mode);
//void MakeLookupHisto(THnSparse * his, TTreeSRedirector *pcstream, const char* hname, Int_t run);
void FitLookup(TChain *chainIn, const char *prefix, TVectorD &vecA, TVectorD &vecC, TVectorD& vecStatA, TVectorD &vecStatD, TCut cut, TObjArray *picArray);
void MakeFits(Int_t run);
void MakeFitTree();
void DrawDistortionDy(TCut cutUser="", Double_t ymin=-0.6, Double_t ymax=0.6);
void DrawDistortionMaps(const char *fname="mean.root");
-
+TCanvas * DrawDistortionMaps(TTree * treedy, TTree * treedsnp, TTree* treed1Pt, const char * name, const char *title);
+AliTPCComposedCorrection * MakeComposedCorrection();
+void MakeLaserTree();
+void AddEffectiveCorrection(AliTPCComposedCorrection* comp);
void MakeLookup(Int_t run, Int_t mode){
//
// make a lookup tree with mean values
- //
+ // 5. make laser tree
+ // 4.
gSystem->AddIncludePath("-I$ALICE_ROOT/STAT");
gSystem->Load("libANALYSIS");
gSystem->Load("libTPCcalib");
+ if (mode==5) {MakeLaserTree();return;};
+ if (mode==4) {DrawDistortionMaps();return;}
if (mode==1){
DrawDistortionDy("abs(snp)<0.25");
DrawDistortionMaps();
return;
}
if (mode==2) return MakeFitTree();
-
- TFile f("CalibObjects.root");
+ TFile f("TPCTimeObjects.root");
AliTPCcalibTime *calibTime= (AliTPCcalibTime*)f.Get("calibTime");
+ if (!calibTime){
+ TFile* f2 = new TFile("CalibObjects.root");
+ calibTime= (AliTPCcalibTime*)f2->Get("calibTime");
+ if (!calibTime){
+ TFile* f3 = new TFile("../CalibObjects.root");
+ calibTime= (AliTPCcalibTime*)f3->Get("calibTime");
+ }
+ }
//
TTreeSRedirector * pcstream = new TTreeSRedirector("mean.root");
THnSparse * his = 0;
his->GetAxis(2)->SetRange(0,1000000);
his->GetAxis(3)->SetRangeUser(-0.5,0.5);
AliTPCCorrection::MakeDistortionMap(his,pcstream, Form("TRD%s",hname[ihis]),run);
+ his = calibTime->GetResHistoTPCTOF(ihis);
+ if (his){
+ his->GetAxis(1)->SetRangeUser(-1.1,1.1);
+ his->GetAxis(2)->SetRange(0,1000000);
+ his->GetAxis(3)->SetRangeUser(-0.5,0.5);
+ AliTPCCorrection::MakeDistortionMap(his,pcstream, Form("TOF%s",hname[ihis]),run);
+ }
}
}
delete pcstream;
}
-void DrawDistortionMaps(const char *fname){
- //
- //
- //
- TFile f(fname);
- TTree *ITSdy=(TTree*)f.Get("ITSdy");
- TTree *ITSdsnp=(TTree*)f.Get("ITSdsnp");
- TTree *ITSd1pt=(TTree*)f.Get("ITSd1pt");
- TCanvas *cdist = new TCanvas("distITS","distITS",1200,800);
- cdist->Divide(3,2);
- cdist->Draw("");
- TH1 * his=0;
- cdist->cd(1);
- ITSdy->Draw("10*mean","rms>0&&abs(snp)<0.25&&entries>100","");
- his = (TH1*)ITSdy->GetHistogram()->Clone();
- his->SetName("dY"); his->SetXTitle("#Delta_{r#phi} (cm)");
- his->Draw();
- //
- cdist->cd(2);
- ITSdsnp->Draw("1000*mean","rms>0&&abs(snp)<0.25&&entries>200","");
- his = (TH1*)ITSdsnp->GetHistogram()->Clone();
- his->SetName("dsnp"); his->SetXTitle("#Delta_{sin(#phi)}");
- his->Draw();
- //
- cdist->cd(3);
- ITSd1pt->Draw("mean","rms>0&&abs(snp)<0.15&&entries>400","");
- his = (TH1*)ITSd1pt->GetHistogram()->Clone();
- his->SetName("d1pt"); his->SetXTitle("#Delta_{1/pt} (1/GeV)");
- his->Draw();
- //
- ITSdy->SetMarkerSize(0.3);
- ITSdsnp->SetMarkerSize(0.3);
- ITSd1pt->SetMarkerSize(0.3);
- {for (Int_t type=0; type<3; type++){
- cdist->cd(4+type);
- TTree * tree =ITSdy;
- if (type==1) tree=ITSdsnp;
- if (type==2) tree=ITSd1pt;
- Int_t counter=0;
- for (Double_t theta=-1; theta<=1; theta+=0.2){
- TCut cut=Form("rms>0&&abs(snp)<0.25&&entries>20&&abs(theta-%f)<0.05",theta);
- tree->SetMarkerStyle(20+counter);
- Option_t *option=(counter==0)? "prof": "profsame";
- if (theta>0) tree->SetMarkerColor(2);
- if (theta<0) tree->SetMarkerColor(4);
- if (type==0) tree->Draw("10*mean:phi>>his(45,-pi,pi)",cut,option);
- if (type==1) tree->Draw("1000*mean:phi>>his(45,-pi,pi)",cut,option);
- if (type==2) tree->Draw("mean:phi>>his(45,-pi,pi)",cut,option);
- his = (TH1*)tree->GetHistogram()->Clone();
- his->SetName(Form("%d_%d",counter,type));
- his->SetXTitle("#phi");
- his->SetMaximum(4);
- his->SetMinimum(-4);
- if (type==2){
- his->SetMaximum(0.06);
- his->SetMinimum(-0.06);
- }
- his->Draw(option);
- counter++;
- }
- }
- }
- cdist->SaveAs("itstpcdistrotion.pdf");
-}
+
void MakeFitTree(){
// .x $ALICE_ROOT/ANALYSIS/CalibMacros/Pass0/ConfigCalibTrain.C(114972,"local:///lustre/alice/alien/alice/data/2010/OCDB")
//
//
- Double_t bzField=AliTrackerBase::GetBz();
- Double_t vdrift = 2.6; // [cm/us] // to be updated: per second (ideally)
- Double_t ezField = 400; // [V/cm] // to be updated: never (hopefully)
- Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ;
- Double_t T1 = 0.9;
- Double_t T2 = 1.5;
- //
- AliTPCExBTwist *twistX1= new AliTPCExBTwist;
- twistX1->SetXTwist(0.001); // 1 mrad twist in x
- twistX1->SetName("tX1");
- //
- AliTPCExBTwist *twistY1= new AliTPCExBTwist;
- twistY1->SetYTwist(0.001); // 1 mrad twist in Y
- twistY1->SetName("tY1");
- //
- AliTPCGGVoltError* ggErrorA = new AliTPCGGVoltError;
- ggErrorA->SetDeltaVGGA(40.); // delta GG - 1 mm
- ggErrorA->SetName("ggA");
- //
- AliTPCGGVoltError* ggErrorC = new AliTPCGGVoltError;
- ggErrorC->SetDeltaVGGC(40.); // delta GG - 1 mm
- ggErrorC->SetName("ggC");
- // conical free param
- //
- AliTPCExBConical *exbConA = new AliTPCExBConical;
- Float_t conicalAA[3]={1., 0.0,0.00};
- Float_t conicalCA[3]={0., 0.0,0.00};
- exbConA->SetConicalA(conicalAA);
- exbConA->SetConicalC(conicalCA);
- exbConA->SetConicalFactor(0.25);
- exbConA->SetName("ExBConA"); // conical factor - A side
-
- //conical free param
- AliTPCExBConical *exbConC = new AliTPCExBConical;
- Float_t conicalAC[3]={0., 0.0,0.00};
- Float_t conicalCC[3]={1., 0.0,0.00};
- exbConC->SetConicalA(conicalAC);
- exbConC->SetConicalC(conicalCC);
- exbConC->SetConicalFactor(0.25);
- exbConC->SetName("ExBConC"); // conical factor - c side
- // conical as measured
- AliTPCExBConical *exbCon = new AliTPCExBConical;
- Float_t conicalC[3]={1., 0.6,-0.08};
- Float_t conicalA[3]={0.9, 0.3,0.04};
- exbCon->SetConicalA(conicalA);
- exbCon->SetConicalC(conicalC);
- exbCon->SetConicalFactor(0.25);
- exbCon->SetName("ExBCon"); // conical factor
-
- TObjArray * corr = new TObjArray;
- corr->AddLast(twistX1);
- corr->AddLast(twistY1);
- corr->AddLast(ggErrorA);
- corr->AddLast(ggErrorC);
- corr->AddLast(exbConA);
- corr->AddLast(exbConC);
- corr->AddLast(exbCon);
- //
- AliTPCComposedCorrection *cc= new AliTPCComposedCorrection ;
- cc->SetCorrections(corr);
- cc->SetOmegaTauT1T2(wt,T1,T2);
- cc->Init();
- // cc->SetMode(1);
- cc->Print("DA"); // Print used correction classes
- cc->SetName("Comp");
- //corr->AddLast(cc);
+ AliTPCComposedCorrection *cc= MakeComposedCorrection();
+ TObjArray * corr = (TObjArray*)(cc->GetCorrections());
+ // corr->AddLast(cc);
+ const Int_t kskip=3;
//
TFile f("mean.root");
TTree * tree= 0;
tree = (TTree*)f.Get("ITSdy");
- AliTPCCorrection::MakeTrackDistortionTree(tree,0,0,corr,7);
+ printf("ITSdy\n");
+ AliTPCCorrection::MakeTrackDistortionTree(tree,0,0,corr,kskip);
tree = (TTree*)f.Get("TRDdy");
- AliTPCCorrection::MakeTrackDistortionTree(tree,1,0,corr,7);
+ printf("TRDdy\n");
+ AliTPCCorrection::MakeTrackDistortionTree(tree,1,0,corr,kskip);
tree = (TTree*)f.Get("Vertexdy");
- AliTPCCorrection::MakeTrackDistortionTree(tree,2,0,corr,7);
+ printf("Vertexdy\n");
+ AliTPCCorrection::MakeTrackDistortionTree(tree,2,0,corr,kskip);
+ tree = (TTree*)f.Get("TOFdy");
+ printf("TOFdy\n");
+ if (tree) AliTPCCorrection::MakeTrackDistortionTree(tree,3,0,corr,kskip);
+
//
tree = (TTree*)f.Get("ITSdsnp");
- AliTPCCorrection::MakeTrackDistortionTree(tree,0,2,corr,7);
+ printf("ITSdsnp\n");
+ AliTPCCorrection::MakeTrackDistortionTree(tree,0,2,corr,kskip);
tree = (TTree*)f.Get("TRDdsnp");
- AliTPCCorrection::MakeTrackDistortionTree(tree,1,2,corr,7);
+ printf("TRDdsnp\n");
+ AliTPCCorrection::MakeTrackDistortionTree(tree,1,2,corr,kskip);
tree = (TTree*)f.Get("Vertexdsnp");
- AliTPCCorrection::MakeTrackDistortionTree(tree,2,2,corr,7);
+ printf("Vertexdsnp\n");
+ AliTPCCorrection::MakeTrackDistortionTree(tree,2,2,corr,kskip);
//
+
+ tree = (TTree*)f.Get("ITSd1pt");
+ printf("ITSd1pt\n");
+ if (tree) AliTPCCorrection::MakeTrackDistortionTree(tree,0,4,corr,kskip);
+
+ tree = (TTree*)f.Get("TOFdz");
+ printf("TOFdz\n");
+ if (tree) AliTPCCorrection::MakeTrackDistortionTree(tree,3,0,corr,kskip);
+
+
tree = (TTree*)f.Get("ITSdz");
- AliTPCCorrection::MakeTrackDistortionTree(tree,0,1,corr,7);
- //tree = (TTree*)f.Get("TRDdz");
- //AliTPCCorrection::MakeTrackDistortionTree(tree,1,1,corr,3);
+ printf("ITSdz\n");
+ AliTPCCorrection::MakeTrackDistortionTree(tree,0,1,corr,kskip);
tree = (TTree*)f.Get("Vertexdz");
- AliTPCCorrection::MakeTrackDistortionTree(tree,2,1,corr,7);
+ printf("Vertexdz\n");
+ AliTPCCorrection::MakeTrackDistortionTree(tree,2,1,corr,kskip);
//
tree = (TTree*)f.Get("ITSdtheta");
- AliTPCCorrection::MakeTrackDistortionTree(tree,0,3,corr,7);
- //tree = (TTree*)f.Get("TRDdtheta");
- // AliTPCCorrection::MakeTrackDistortionTree(tree,1,3,corr,3);
+ printf("ITSdtheta\n");
+ AliTPCCorrection::MakeTrackDistortionTree(tree,0,3,corr,kskip);
//
tree = (TTree*)f.Get("Vertexdtheta");
- AliTPCCorrection::MakeTrackDistortionTree(tree,2,3,corr,7);
-
- tree = (TTree*)f.Get("ITSd1pt");
- AliTPCCorrection::MakeTrackDistortionTree(tree,0,4,corr,7);
- tree = (TTree*)f.Get("TRDd1pt");
- AliTPCCorrection::MakeTrackDistortionTree(tree,1,4,corr,7);
- tree = (TTree*)f.Get("Vertexd1pt");
- AliTPCCorrection::MakeTrackDistortionTree(tree,2,4,corr,7);
+ printf("Vertexdtheta\n");
+ AliTPCCorrection::MakeTrackDistortionTree(tree,2,3,corr,kskip);
}
//
AliXRDPROOFtoolkit tool;
- TChain *chain = tool.MakeChain("distortion.txt","fit",0,100000);
+ // TChain *chain = tool.MakeChain("distortion.txt","fit",0,100000);
TChain *chainPlus = tool.MakeChain("distortionPlus.txt","fit",0,100000);
TChain *chainMinus = tool.MakeChain("distortionMinus.txt","fit",0,100000);
+ TChain *chain0 = tool.MakeChain("distortionField0.txt","fit",0,100000);
+
chainPlus->AddFriend(chainMinus,"M");
+ chainPlus->AddFriend(chain0,"F0");
TCut cut="rms>0&&M.rms>0&&entries>100&&dtype==0&&(ptype==0||ptype==2||ptype==4)";
chainPlus->SetAlias("deltaM","mean-M.mean");
canvasdSnp->SaveAs(Form("fitdsnp%d.pdf",highFrequency));
delete fpic;
}
+
+void DrawDistortionMaps(const char *fname){
+ TFile f(fname);
+ TTree *ITSdy=(TTree*)f.Get("ITSdy");
+ TTree *ITSdsnp=(TTree*)f.Get("ITSdsnp");
+ TTree *ITSd1pt=(TTree*)f.Get("ITSd1pt");
+ TTree *TRDdy=(TTree*)f.Get("TRDdy");
+ TTree *TRDdsnp=(TTree*)f.Get("TRDdsnp");
+ TTree *TRDd1pt=(TTree*)f.Get("TRDd1pt");
+ TTree *Vertexdy=(TTree*)f.Get("Vertexdy");
+ TTree *Vertexdsnp=(TTree*)f.Get("Vertexdsnp");
+ TTree *Vertexd1pt=(TTree*)f.Get("Vertexd1pt");
+ TCanvas *cits = DrawDistortionMaps(ITSdy,ITSdsnp,ITSd1pt,"ITS","ITS");
+ cits->SaveAs("matchingITS.pdf");
+ TCanvas *ctrd = DrawDistortionMaps(TRDdy,TRDdsnp,TRDd1pt,"TRD","TRD");
+ ctrd->SaveAs("matchingTRD.pdf");
+ TCanvas *cvertex = DrawDistortionMaps(Vertexdy,Vertexdsnp,Vertexd1pt,"Vertex","Vertex");
+ cvertex->SaveAs("matchingVertex.pdf");
+ //
+ //
+ TPostScript *ps = new TPostScript("matching.ps", 112);
+ gStyle->SetOptTitle(1);
+ gStyle->SetOptStat(0);
+ ps->NewPage();
+ cits->Draw();
+ ps->NewPage();
+ ctrd->Draw();
+ ps->NewPage();
+ cvertex->Draw();
+ ps->Close();
+ delete ps;
+}
+
+
+TCanvas * DrawDistortionMaps(TTree * treedy, TTree * treedsnp, TTree* treed1pt, const char * name, const char *title){
+ //
+ //
+ //
+ TH1::AddDirectory(0);
+ TCanvas *cdist = new TCanvas(name,name,1200,800);
+ cdist->Divide(3,2);
+ cdist->Draw("");
+ TH1 * his=0;
+
+ cdist->cd(1);
+ treedy->Draw("10*mean","rms>0&&abs(snp)<0.25&&entries>100","");
+ his = (TH1*)treedy->GetHistogram()->Clone();
+ his->SetName("dY"); his->SetXTitle("#Delta_{r#phi} (cm)");
+ his->Draw();
+ //
+ cdist->cd(2);
+ treedsnp->Draw("1000*mean","rms>0&&abs(snp)<0.25&&entries>200","");
+ his = (TH1*)treedsnp->GetHistogram()->Clone();
+ his->SetName("dsnp"); his->SetXTitle("#Delta_{sin(#phi)}");
+ his->Draw();
+ //
+ cdist->cd(3);
+ treed1pt->Draw("mean","rms>0&&abs(snp)<0.15&&entries>400","");
+ his = (TH1*)treed1pt->GetHistogram()->Clone();
+ his->SetName("d1pt"); his->SetXTitle("#Delta_{1/pt} (1/GeV)");
+ his->Draw();
+ //
+ treedy->SetMarkerSize(1);
+ treedsnp->SetMarkerSize(1);
+ treed1pt->SetMarkerSize(1);
+ {for (Int_t type=0; type<3; type++){
+ cdist->cd(4+type);
+ TTree * tree =treedy;
+ if (type==1) tree=treedsnp;
+ if (type==2) tree=treed1pt;
+ Int_t counter=0;
+ for (Double_t theta=-1; theta<=1; theta+=0.2){
+ if (TMath::Abs(theta)<0.11) continue;
+ TCut cut=Form("rms>0&&abs(snp)<0.25&&entries>20&&abs(theta-%f)<0.1",theta);
+ tree->SetMarkerStyle(20+counter);
+ Option_t *option=(counter==0)? "": "same";
+ if (theta>0) tree->SetMarkerColor(2);
+ if (theta<0) tree->SetMarkerColor(4);
+ if (type==0) tree->Draw("10*mean:phi>>his(45,-pi,pi)",cut,"profgoff");
+ if (type==1) tree->Draw("1000*mean:phi>>his(45,-pi,pi)",cut,"profgoff");
+ if (type==2) tree->Draw("mean:phi>>his(45,-pi,pi)",cut,"profgoff");
+ his = (TH1*)tree->GetHistogram()->Clone();
+ his->SetName(Form("%d_%d",counter,type));
+ his->SetXTitle("#phi");
+ his->SetMaximum(4);
+ his->SetMinimum(-4);
+ if (type==2){
+ his->SetMaximum(0.06);
+ his->SetMinimum(-0.06);
+ }
+ his->Draw(option);
+ counter++;
+ }
+ }
+ }
+ // cdist->SaveAs(Form("%s.pdf"),name);
+ return cdist;
+}
+
+
+
+
+
+AliTPCComposedCorrection * MakeComposedCorrection(){
+
+ //
+ Double_t bzField=AliTrackerBase::GetBz();
+ AliMagF* magF= (AliMagF*)(TGeoGlobalMagField::Instance()->GetField());
+ Double_t vdrift = 2.6; // [cm/us] // to be updated: per second (ideally)
+ Double_t ezField = 400; // [V/cm] // to be updated: never (hopefully)
+ Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ;
+ Double_t T1 = 0.9;
+ Double_t T2 = 1.5;
+ //
+ //
+
+
+ AliTPCExBTwist *twistX1= new AliTPCExBTwist;
+ twistX1->SetXTwist(0.001); // 1 mrad twist in x
+ twistX1->SetName("tX1");
+ //
+ AliTPCExBTwist *twistY1= new AliTPCExBTwist;
+ twistY1->SetYTwist(0.001); // 1 mrad twist in Y
+ twistY1->SetName("tY1");
+ //
+ AliTPCGGVoltError* ggErrorA = new AliTPCGGVoltError;
+ ggErrorA->SetDeltaVGGA(40.); // delta GG - 1 mm
+ ggErrorA->SetName("ggA");
+ //
+ AliTPCGGVoltError* ggErrorC = new AliTPCGGVoltError;
+ ggErrorC->SetDeltaVGGC(40.); // delta GG - 1 mm
+ ggErrorC->SetName("ggC");
+ // conical free param
+ //
+ AliTPCCalibGlobalMisalignment *shiftX=new AliTPCCalibGlobalMisalignment;
+ shiftX->SetXShift(0.1);
+ shiftX->SetName("shiftX");
+ AliTPCCalibGlobalMisalignment *shiftY=new AliTPCCalibGlobalMisalignment;
+ shiftY->SetYShift(0.1);
+ shiftY->SetName("shiftY");
+
+ AliTPCExBBShape *exb11 = new AliTPCExBBShape;
+ exb11->SetBField(magF);
+ exb11->SetName("exb11");
+ AliTPCExBBShape *exbT1 = new AliTPCExBBShape;
+ exbT1->SetBField(magF);
+ exbT1->SetName("exbT1");
+ AliTPCExBBShape *exbT2 = new AliTPCExBBShape;
+ exbT2->SetBField(magF);
+ exbT2->SetName("exbT2");
+
+ TObjArray * corr = new TObjArray;
+ corr->AddLast(twistX1);
+ corr->AddLast(twistY1);
+ corr->AddLast(ggErrorA);
+ corr->AddLast(ggErrorC);
+ corr->AddLast(shiftX);
+ corr->AddLast(shiftY);
+ corr->AddLast(exb11);
+ corr->AddLast(exbT1);
+ corr->AddLast(exbT2);
+ //
+ AliTPCComposedCorrection *cc= new AliTPCComposedCorrection ;
+ cc->SetCorrections((TObjArray*)(corr));
+ AddEffectiveCorrection(cc);
+ cc->SetOmegaTauT1T2(wt,T1,T2);
+ exbT1->SetOmegaTauT1T2(wt,T1+0.1,T2);
+ exbT2->SetOmegaTauT1T2(wt,T1,T2+0.1);
+ cc->Init();
+ // cc->SetMode(1);
+ cc->Print("DA"); // Print used correction classes
+ cc->SetName("Comp");
+ return cc;
+}
+
+
+
+void AddEffectiveCorrection(AliTPCComposedCorrection* comp){
+ //
+ //
+ //
+ TMatrixD polA(18,4);
+ TMatrixD polC(18,4);
+ TMatrixD valA(18,1);
+ TMatrixD valC(18,1);
+ Int_t counter=0;
+ { for (Int_t iside=0; iside<=1; iside++){
+ {for (Int_t ir=0; ir<3; ir++)
+ for (Int_t id=0; id<3; id++) {
+ polA(counter,0)=0;
+ polA(counter,1)=ir;
+ polA(counter,2)=id;
+ polC(counter,0)=0;
+ polC(counter,1)=ir;
+ polC(counter,2)=id;
+ valA(counter,0)=0;
+ valC(counter,0)=0;
+ counter++;
+ }
+ }
+ }
+ }
+ counter=0;
+ TObjArray * array = (TObjArray*) comp->GetCorrections();
+ { for (Int_t iside=0; iside<=1; iside++)
+ for (Int_t ir=0; ir<3; ir++)
+ for (Int_t id=0; id<3; id++) {
+ AliTPCExBEffective *eff=new AliTPCExBEffective;
+ valA*=0;
+ valC*=0;
+ eff->SetName(Form("eff%d_%d_%d",iside,ir,id));
+ if (iside==0) valA(counter,0)=0.1;
+ if (iside==1) valC(counter,0)=0.1;
+ eff->SetPolynoms(&polA,&polC);
+ eff->SetCoeficients(&valA,&valC);
+ valA.Print();
+ valC.Print();
+ array->AddLast(eff);
+ counter++;
+ }
+ }
+}
+
+
+
+void MakeLaserTree(){
+ //
+ //
+ //
+ AliTPCComposedCorrection *cc= MakeComposedCorrection();
+ TObjArray * corr = (TObjArray*)(cc->GetCorrections());
+ // corr->AddLast(cc);
+ TFile f("laserMean.root");
+ TTree* tree =(TTree*)f.Get("Mean");
+ AliTPCCorrection::MakeLaserDistortionTree(tree,corr,0);
+
+}