M AliTPCcalibLaser.cxx - Adding Integral of B field
authormarian <marian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 8 Mar 2010 12:55:50 +0000 (12:55 +0000)
committermarian <marian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 8 Mar 2010 12:55:50 +0000 (12:55 +0000)
M      AliTPCcalibTime.cxx   - Protection in case of missing friends
                               In case of filtered ESDs

CalibExB.C                     - Obsolete macro removed
CalibLaserExBscan.C            - Replaced by another

Marian

TPC/AliTPCcalibLaser.cxx
TPC/AliTPCcalibTime.cxx
TPC/CalibMacros/CalibExB.C [deleted file]
TPC/CalibMacros/CalibLaserExBscan.C [new file with mode: 0644]

index bc999cc..ba72a71 100644 (file)
@@ -611,7 +611,7 @@ void AliTPCcalibLaser::Process(AliESDEvent * event) {
   if (counter<kMinTracks) return;
 
   //FitDriftV();
-  FitDriftV(0.3);
+  FitDriftV(0.2);
   if (!fFullCalib) return;
   static Bool_t init=kFALSE;
   if (!init){
@@ -961,8 +961,8 @@ Bool_t  AliTPCcalibLaser::FitDriftV(Float_t minFraction){
   const Float_t kSaturCut    = 0.05;    // remove saturated lasers - cut on fraction of saturated 
   const Float_t kDistCut     = 3.;      // distance sigma cut - 3 sigma
   const Float_t kDistCutAbs  = 1.;      // absolute cut 1 cm
-  const Float_t kMinClusters = 60.;      // minimal amount of the clusters
-  const Float_t kMinSignal   = 10.;      // minimal mean height of the signal
+  const Float_t kMinClusters = 40.;      // minimal amount of the clusters
+  const Float_t kMinSignal   = 2.5;      // minimal mean height of the signal
   const Float_t kChi2Cut     = 1.0;     // chi2 cut to accept drift fit
   //
   static TLinearFitter fdriftA(3,"hyp2");
@@ -2249,8 +2249,9 @@ void AliTPCcalibLaser::DumpMeanInfo(Int_t run){
   const Float_t kmultiCut=2;
   const Float_t kcutP0=0.002;
   AliMagF* magF=  dynamic_cast<AliMagF*> (TGeoGlobalMagField::Instance()->GetField());
-  Double_t xyz[3]={90,0,10};
-  Double_t bxyz[3]={90,0,10};
+  Double_t xyz[3]={90,0,10};         // tmp. global position 
+  Double_t bxyz[3]={90,0,10};        // tmp. mag field  integral - cylindrical
+  Double_t bgxyz[3]={90,0,10};       // tmp. mag field  integral - local
   //
   AliTPCcalibLaser *laser = this;
   TTreeSRedirector *pcstream = new TTreeSRedirector("laserMean.root");
@@ -2556,21 +2557,42 @@ void AliTPCcalibLaser::DumpMeanInfo(Int_t run){
     // magnetic field integrals
     TVectorD vecIBR(159);        // radial
     TVectorD vecIBRPhi(159);     // r-phi
+    TVectorD vecIBLX(159);       // local x
+    TVectorD vecIBLY(159);       // local y
+    TVectorD vecIBGX(159);       // local x
+    TVectorD vecIBGY(159);       // local y
     TVectorD vecIBZ(159);        // z
     //
     for (Int_t irow=0;irow<159;irow++){
       vecIBR[irow]=0;
       vecIBRPhi[irow]=0;
+      vecIBLX[irow]=0;
+      vecIBLY[irow]=0;
+      vecIBGX[irow]=0;
+      vecIBGY[irow]=0;
       vecIBZ[irow]=0;
-      Double_t gx=(*(ltrp->fVecGX))[irow];
-      Double_t gy=(*(ltrp->fVecGY))[irow];
+      Double_t gx    =(*(ltrp->fVecGX))[irow];
+      Double_t gy    =(*(ltrp->fVecGY))[irow];
+      Int_t    lsec  =TMath::Nint((*(ltrp->fVecSec))[irow]);
+      Double_t   ca  =TMath::Cos(TMath::Pi()*(lsec+0.5)/9.);
+      Double_t   sa  =TMath::Sin(TMath::Pi()*(lsec+0.5)/9.);
       xyz[2]=(*(ltrp->fVecGZ))[irow];
       xyz[0]=TMath::Sqrt(gx*gx+gy*gy);
       xyz[1]=TMath::ATan2(gy,gx);
+      Double_t gxyz[3]={gx,gy,(*(ltrp->fVecGZ))[irow]};
       if (magF){
        magF->GetTPCIntCyl(xyz,bxyz);
+       magF->GetTPCInt(gxyz,bgxyz);
        vecIBR[irow]=bxyz[0];
        vecIBRPhi[irow]=bxyz[1];
+       //
+       vecIBGX[irow]=bgxyz[0];
+       vecIBGY[irow]=bgxyz[1];
+       //
+       vecIBLX[irow]=  bgxyz[0]*ca+bgxyz[1]*sa;
+       vecIBLY[irow]= -bgxyz[0]*sa+bgxyz[1]*ca;
+       //
+
        vecIBZ[irow]=bxyz[2];
       }
     }
@@ -2914,7 +2936,14 @@ void AliTPCcalibLaser::DumpMeanInfo(Int_t run){
       "lasTanPhiLocOut="<<lasTanPhiLocOut<<// laser tan phi in local frame (outer)
       "ibr.="<<&vecIBR<<   // radial filed integral
       "ibrphi.="<<&vecIBRPhi<<   // r=phifiled integral
-      "ibz.="<<&vecIBZ<<   // radial filed integral
+      "ibr.="<<&vecIBR<<   // radial filed integral
+      "ibz.="<<&vecIBZ<<   // z filed integral
+      //
+      "iblx.="<<&vecIBLX<<   // local bx  integral
+      "ibly.="<<&vecIBLY<<   // local by integral
+      "ibgx.="<<&vecIBGX<<   // global bx  integral
+      "ibgy.="<<&vecIBGY<<   // global by integral
+      //
       "X.="<<&vecX<<       // local x 
       "Y.="<<&vecY<<       // local y 
       "R.="<<&vecR<<       // radius 
index 6bcfb93..79d585c 100644 (file)
@@ -566,6 +566,7 @@ void AliTPCcalibTime::ProcessCosmic(AliESDEvent *event){
     if (!trackOut) continue;
     
     AliESDfriendTrack *friendTrack = esdFriend->GetTrack(i);
+    if (!friendTrack) continue;
     if (friendTrack) ProcessSame(track,friendTrack,event);
     if (friendTrack) ProcessAlignITS(track,friendTrack,event,esdFriend);
     if (friendTrack) ProcessAlignTRD(track,friendTrack);
diff --git a/TPC/CalibMacros/CalibExB.C b/TPC/CalibMacros/CalibExB.C
deleted file mode 100644 (file)
index 28ff8c0..0000000
+++ /dev/null
@@ -1,277 +0,0 @@
-void Init(){
-  //
-  // Initialize
-  //
-  AliMagF* field = new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG);
-  AliTPCExB::RegisterField(0,field);
-  AliMagF* fieldC0 = new AliMagF("Maps","Maps", 1, 1, AliMagF::k5kG);
-  AliTPCExB::RegisterField(1,fieldC0);
-  AliMagF* fieldC1 = new AliMagF("Maps","Maps", 1, 1, AliMagF::k5kG);
-  AliTPCExB::RegisterField(2,fieldC1);
-
-  gSystem->Load("libSTAT.so");
-  AliTPCExBFirst *exbfirst1  = new  AliTPCExBFirst(fieldC1,0.88*2.6400e+04,50,50,50);
-  AliTPCExB::SetInstance(exbfirst1);
-
-}
-
-void FitField(){
-  //
-  //
-  //
-  AliTPCExB * exb =  AliTPCExB::Instance();
-  exb->TestExB("field.root");
-  TFile f("field.root");
-  TTree * tree = (TTree*)f.Get("positions");
-  //
-  TStatToolkit toolkit;
-  Double_t chi2;
-  TVectorD fitParam;
-  TMatrixD covMatrix;
-  Int_t npoints;
-  //
-  // SetAliases
-  //
-  tree->SetAlias("sa","sin(phi+0.0)");
-  tree->SetAlias("ca","cos(phi+0.0)");
-  tree->SetAlias("sa2","sin(phi*2+0.0)");
-  tree->SetAlias("ca2","cos(phi*2+0.0)");
-  tree->SetAlias("zn","(x2/250.)");
-  tree->SetAlias("rn","(r/250.)");
-    
-  TString fstringSym="";
-  //  
-  fstringSym+="zn++";
-  fstringSym+="rn++";
-  fstringSym+="zn*rn++";
-  fstringSym+="zn*zn++";
-  fstringSym+="zn*zn*rn++";
-  fstringSym+="zn*rn*rn++";
-  //
-  fstringSym+="sa++";
-  fstringSym+="ca++";  
-  fstringSym+="ca2++";
-  fstringSym+="sa2++";
-  fstringSym+="ca*zn++";
-  fstringSym+="sa*zn++";
-  fstringSym+="ca2*zn++";
-  fstringSym+="sa2*zn++";
-  fstringSym+="ca*zn*zn++";
-  fstringSym+="sa*zn*zn++";
-  fstringSym+="ca*zn*rn++";
-  fstringSym+="sa*zn*rn++";
-
-  // BrBz
-  TString *strBrBz = toolkit.FitPlane(tree,"br/bz",fstringSym->Data(), "abs(r)<250&&abs(r)>90", chi2,npoints,fitParam,covMatrix);
-  strBrBz->Tokenize("+")->Print();
-  printf("Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
-  tree->SetAlias("fitBrBz",strBrBz->Data());
-  exb->fMatBrBz = new TVectorD(fitParam);
-  // BrfiBz
-  TString *strBrfiBz = toolkit.FitPlane(tree,"brfi/bz",fstringSym->Data(), "abs(r)<250&&abs(r)>90", chi2,npoints,fitParam,covMatrix);
-  strBrfiBz->Tokenize("+")->Print();
-  printf("Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
-  tree->SetAlias("fitBrfiBz",strBrfiBz->Data());
-  exb->fMatBrfiBz = new TVectorD(fitParam);
-
-  
-  //
-  // BR integral parameterization
-  //
-  TString *strBRI0 = toolkit.FitPlane(tree,"bri",fstringSym->Data(), "abs(r)<250&&abs(r)>90&&x2>0&&x2<250", chi2,npoints,fitParam,covMatrix);
-  strBRI0->Tokenize("+")->Print();
-  printf("Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
-  tree->SetAlias("fitBRI0",strBRI0->Data());
-  exb->fMatBrBzI0 = new TVectorD(fitParam);
-
-  TString *strBRI1 = toolkit.FitPlane(tree,"bri",fstringSym->Data(), "abs(r)<250&&abs(r)>90&&x2<-10&&x2>-250", chi2,npoints,fitParam,covMatrix);
-  strBRI1->Tokenize("+")->Print();
-  printf("Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
-  tree->SetAlias("fitBRI1",strBRI1->Data());
-  exb->fMatBrBzI1 = new TVectorD(fitParam);
-  //
-
-  TString *strBRFII0 = toolkit.FitPlane(tree,"brfii",fstringSym->Data(), "abs(r)<250&&abs(r)>90&&x2>0&&x2<250", chi2,npoints,fitParam,covMatrix);
-  strBRFII0->Tokenize("+")->Print();
-  printf("Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
-  tree->SetAlias("fitBRFII0",strBRFII0->Data());
-  exb->fMatBrfiBzI0 = new TVectorD(fitParam);
-  TString *strBRFII1 = toolkit.FitPlane(tree,"brfii",fstringSym->Data(), "abs(r)<250&&abs(r)>90&&x2<-10&&x2>-240", chi2,npoints,fitParam,covMatrix);
-  strBRFII1->Tokenize("+")->Print();
-  printf("Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
-  tree->SetAlias("fitBRFII1",strBRFII1->Data());
-  exb->fMatBrfiBzI1 = new TVectorD(fitParam);
-
-  TFile * fout = new TFile("bfit.root","recreate");
-  exb->Write("bfit");
-  fout->Close();
-
-}
-
-
-
-
-
-
-void CalibExB(){
-  // 
-  // Macro to test ExB correction versus laser B filed scan data
-  //
-  // Before running it be sure the file laserScan.root exist
-  //
-  // Create this file from scan - See AliTPCCalibLaser.C
-  //
-  TFile fscan("laserScan.root");
-  TTree * treeT = (TTree*)fscan.Get("Mean");
-  TStatToolkit toolkit;
-  Double_t chi2;
-  TVectorD fitParam;
-  TMatrixD covMatrix;
-  Int_t npoints;
-  
-  TCut cutF0("abs(gphi1-(pphi0+pphi1*bz))<0.05");
-  TCut cutF1("abs(gphiP1-(pphiP0+pphiP1*bz))<0.0005");
-  TCut cutN("entries>2");
-  
-  TCut cutA = cutF0+cutF1+cutN;
-  
-  treeT->SetAlias("side","(-1+(LTr.fP[1]>0)*2)");       // side
-  treeT->SetAlias("dr","(abs(LTr.fP[1]/250.))");
-  treeT->SetAlias("sa","sin(atan2(lx1+0.0,lx0+0.0))");
-  treeT->SetAlias("ca","cos(atan2(lx1+0.0,lx0+0.0))");
-  treeT->SetAlias("ta","tan(asin(LTr.fP[2]+0.0))");
-
-  TString fstring="";
-  //
-  fstring+="((dr)^3-1)*bz++";             //1
-  fstring+="((dr)^3-1)*bz*sa++";        //3
-  fstring+="((dr)^3-1)*bz*ca++";        //5
-  //
-  fstring+="((dr)^1-1)*bz++";           //1
-  fstring+="((dr)^1-1)*bz*sa++";        //3
-  fstring+="((dr)^1-1)*bz*ca++";        //5
-  //
-  fstring+="side*((dr)^3-1)*bz++";             //1
-  fstring+="side*((dr)^3-1)*bz*sa++";        //3
-  fstring+="side*((dr)^3-1)*bz*ca++";        //5
-  //
-  fstring+="side*((dr)^1-1)*bz++";           //7
-  fstring+="side*((dr)^1-1)*bz*sa++";        //9
-  fstring+="side*((dr)^1-1)*bz*ca++";        //11
-  //
-  fstring+="((dr)^3-1)*bz^2*ta++";        //2
-  fstring+="((dr)^3-1)*bz^2*sa*ta++";     //4
-  fstring+="((dr)^3-1)*bz^2*ca*ta++";     //6
-
-  fstring+="((dr)^1-1)*bz^2*ta++";        //2
-  fstring+="((dr)^1-1)*bz^2*sa*ta++";     //4
-  fstring+="((dr)^1-1)*bz^2*ca*ta++";     //6
-  //  
-  fstring+="side*((dr)^1-1)*bz^2*ta++";        //8 
-  fstring+="side*((dr)^1-1)*bz^2*sa*ta++";     //10
-  fstring+="side*((dr)^1-1)*bz^2*ca*ta++";     //12
-
-
-
-
-  TString *strq0 = toolkit.FitPlane(treeT,"gphi1-pphi0",fstring->Data(), cutA, chi2,npoints,fitParam,covMatrix);
-  strq0->Tokenize("+")->Print();
-  treeT->SetAlias("fit",strq0->Data());
-  printf("Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
-
-
-  
-  treeT->SetAlias("bir1",  "-AliTPCExB::GetBrI(254,atan2(lx1,lx0+0),LTr.fP[1],2)");
-  treeT->SetAlias("birfi1","-AliTPCExB::GetBrfiI(254,atan2(lx1,lx0+0),LTr.fP[1],2)");
-
-  treeT->SetAlias("bir0",  "AliTPCExB::GetBrI(254,atan2(lx1,lx0+0),LTr.fP[1],2)");
-  treeT->SetAlias("birfi0","AliTPCExB::GetBrfiI(254,atan2(lx1,lx0+0),LTr.fP[1],2)");
-
-  treeT->SetAlias("fbz00", "(bir0+birfi0*ta)");
-  treeT->SetAlias("fbz02", "(birfi0+bir0*ta)");
-  treeT->SetAlias("fbz10", "(bir1+birfi1*ta)");
-  treeT->SetAlias("fbz12", "(birfi1+bir1*ta)");
-  //
-  treeT->SetAlias("fbz0", "((fSide==0)*fbz00+(fSide==1)*fbz10)");
-  treeT->SetAlias("fbz2", "((fSide==0)*fbz02+(fSide==1)*fbz12)");
-
-  //
-  TString fstringeb="";
-  //
-  fstringeb+="bz*fbz0++";                //1
-  fstringeb+="sign(bz)*bz^2*fbz2++";     //2
-  fstringeb+="((dr)-1)*bz*sa++";        //9
-  fstringeb+="side*((dr)-1)*bz*sa++";        //9
-
-
-  fstringeb+="((dr)-1)*bz*ca++";        //9
-  fstringeb+="side*((dr)^3-1)*bz*sa++";        //9
-  fstringeb+="side*((dr)^3-1)*bz++";           //7
-
-  fstringeb+="((dr)^3-1)*bz++";           //7
-  fstringeb+="side*((dr)^3-1)*bz*ca++";        //11
-  
-  //  
-  TString *strExB = toolkit.FitPlane(treeT,"gphi1-pphi0",fstringeb->Data(), "abs(gphi1-pphi0-fit)<0.06&&abs(bz)>0.1"+cutA, chi2,npoints,fitParam,covMatrix);
-  strExB->Tokenize("+")->Print();
-  printf("Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
-  treeT->SetAlias("fitEB",strExB->Data());
-
-
-
-  TString fstringeb="";
-  //
-  fstringeb+="AliTPCExB::GetDrphi(254,atan2(lx1,lx0),LTr.fP[1],bz*10)++";        //1
-  fstringeb+="AliTPCExB::GetDr(254,atan2(lx1,lx0),LTr.fP[1],-bz*10)*ta++";        //2
-  //
-  fstringeb+="side*((dr)^1-1)*bz++";        //9
-  fstringeb+="side*((dr)^1-1)*bz*sa++";        //9
-  //
-  fstringeb+="side*((dr)^1-1)*bz*ta++";        //9
-  fstringeb+="side*((dr)^1-1)*bz*sa*ta++";        //9
-  fstringeb+="side*((dr)^1-1)*bz*ca*ta++";        //9
-
-  TString *strExB = toolkit.FitPlane(treeT,"fit",fstringeb->Data(), "abs(gphi1-pphi0-fit)<0.06&&abs(bz)>0.1"+cutA, chi2,npoints,fitParam,covMatrix);
-  strExB->Tokenize("+")->Print();
-  printf("Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
-  treeT->SetAlias("fitEB",strExB->Data());
-
-
-
-
-
-
-  {
-    Float_t otcor0    =0.7;
-    Float_t omegatau0 =0.35;
-    Float_t rms0      =10;
-    for (Float_t otcor=0.6; otcor<1.1;otcor+=0.1)
-      for (Float_t ometau=0.4; ometau<0.45;ometau+=0.01){
-       char hname[100];
-       sprintf(hname,"hist_ometau_%f(50,-0.10,0.10)",ometau);
-       char expr[100];
-       sprintf(expr,"gphi1-pphi0-bz*fbz0*%f-sign(bz)*(%f*%f*bz)^2*fbz2>>%s",ometau,ometau,otcor,hname);
-       treeT->Draw(expr,"abs(gphi1-pphi0-fit)<0.06&&abs(bz)>0.1"+cutA);    
-       printf("Ometau=%f\tCor=%f\tRMS=%f\n",ometau,otcor,treeT->GetHistogram()->GetRMS());
-       if (rms0>treeT->GetHistogram()->GetRMS()){
-         otcor0=otcor;
-         omegatau0=ometau;
-         rms0=treeT->GetHistogram()->GetRMS();
-       }
-      }
-  }
-}
-
-
-
-
-void MakePic(){
-  //
-  //
-  //
-
-}
-
diff --git a/TPC/CalibMacros/CalibLaserExBscan.C b/TPC/CalibMacros/CalibLaserExBscan.C
new file mode 100644 (file)
index 0000000..abe66fa
--- /dev/null
@@ -0,0 +1,329 @@
+/* 
+//
+   .x ~/rootlogon.C
+   .L $ALICE_ROOT/TPC/CalibMacros/CalibLaserExBscan.C
+
+   
+   // 0. Make a calibration
+   // 1. Make a laser scan list
+   //    e.g in TPC workscape
+   
+   // 2. Define a reference data  e.g:
+   // for a in `cat laserScan.txt`; do echo `pwd`/../mergerunMag0.list/laserMean.root; done >laserScanRef.txt
+  
+   Init();         // init chain
+   MakeAliases();  // make alaises for variables
+
+
+*/
+
+
+
+
+TCut  cutE="eY.fElements<0.02&&Rr.eY.fElements<0.02";         // error cut 200 microns 
+TCut  cutN="nCl.fElements>10&&Rr.nCl.fElements>10";           // number of clusters cut
+TCut  cutEd="(abs(Rr.Y.fElements)-Rr.X.fElements*0.155<-1)";   // edge cut
+TCut  cutB="LTr.fVecLX.fElements>0";                          // local x cut
+TCut  cutR="(LTr.fVecLX.fElements%3)==1&&((abs(LTr.fVecLZ.fElements)+abs(LTr.fVecLY.fElements))%3)==1";   // cut - skip points
+TCut  cutA=cutE+cutN+cutEd+cutB;
+//
+TCut  cutAM5  =cutA+"LTr.fP[1]>0&&abs(bz+5)<0.1";
+TCut  cutAP2  =cutA+"LTr.fP[1]>0&&abs(bz-2)<0.1";
+TCut  cutAP5  =cutA+"LTr.fP[1]>0&&abs(bz-5)<0.1";
+//
+TCut  cutCM5  =cutA+"LTr.fP[1]<0&&abs(bz+5)<0.1";
+TCut  cutCP2  =cutA+"LTr.fP[1]<0&&abs(bz-2)<0.1";
+TCut  cutCP5  =cutA+"LTr.fP[1]<0&&abs(bz-5)<0.1";
+
+
+TChain * chain =0;
+
+TVectorD fitParamB[6];      // magnetic fit param
+TVectorD fitParamBT[6];     // + electric field tilt
+TVectorD fitParamBT0[6];    // + electric field close to ROC              
+TVectorD fitParamA[6];      // + electric field rotation
+//
+TMatrixD covarB[6];         // covariance matrices
+TMatrixD covarBT[6];        //
+TMatrixD covarBT0[6];       //
+TMatrixD covarA[6];         //
+//
+TMatrixD *errB[6];         // covariance matrices
+TMatrixD *errBT[6];        //
+TMatrixD *errBT0[6];       //
+TMatrixD *errA[6];         //
+//
+Double_t chi2B[6];          // chi2
+Double_t chi2BT[6];         // chi2
+Double_t chi2BT0[6];        //
+Double_t chi2A[6];          //
+
+
+void Init(){
+  //
+  // 
+  //
+  gSystem->Load("libANALYSIS");
+  gSystem->Load("libTPCcalib"); 
+  gSystem->Load("libSTAT.so");
+  gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
+  gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+");
+  AliXRDPROOFtoolkit tool;
+  chain = tool.MakeChainRandom("laserScan.txt","Mean",0,10200);
+  chain->Lookup();
+  chainRef = tool.MakeChain("laserScanRef.txt","Mean",0,10200);
+  chainRef->Lookup();
+  chain->AddFriend(chainRef,"Rr");
+}
+
+
+
+void MakeAliases(){
+  //
+  // shortcuts for variables to fit
+  //
+  // bserved distrotions
+  chain->SetAlias("dy","(dY.fElements-Rr.dY.fElements)");              // y ~ (r-phi) distortion
+  chain->SetAlias("ctany","(dY.fElements-Rr.dY.fElements)/(250*dr)");  // mean distortion (angle) over drift length
+  //
+  chain->SetAlias("dr","(250.-abs(LTr.fP[1]))");                       // drift length 0-250
+  chain->SetAlias("dr1","(1.-abs(LTr.fP[1]/250.))");                       // drift length 0-250
+  chain->SetAlias("phiL","(LTr.fVecLY.fElements/LTr.fVecLX.fElements)");       // tann of local phi
+  chain->SetAlias("Esign","(1.-LTr.fSide*2.)");          // E field sign:    +1 for  A side : -1 for c side 
+  //
+  // Br and Brfi  
+  //
+  chain->SetAlias("ablx","(iblx.fElements/(ibz.fElements))");
+  chain->SetAlias("ably","(ibly.fElements/(ibz.fElements))");
+  chain->SetAlias("abr","(ibr.fElements/(ibz.fElements))");
+  chain->SetAlias("abrf","(ibrphi.fElements/(ibz.fElements))");
+
+
+  chain->SetAlias("r","(R.fElements-165.5)/80.3");          // local X - 0 at middle   -1 first +1 last padrow
+  chain->SetAlias("ky","(kY.fElements)");                   // local inclination angle
+  chain->SetAlias("sec","LTr.fVecSec.fElements");           // sector number
+  chain->SetAlias("ca","cos(LTr.fVecPhi.fElements+0.)");    // cos of global phi position
+  chain->SetAlias("sa","sin(LTr.fVecPhi.fElements+0.)");    // sin of global phi position  
+}
+
+
+TMatrixD * MakeErrVector(TMatrixD & mat){
+  //
+  // get error vector
+  //
+  Int_t    nrows=mat.GetNrows();
+  TMatrixD *err = new TMatrixD(nrows,1);
+  for (Int_t i=0; i<nrows;i++) (*err)(i,0)=TMath::Sqrt(mat(i,i));
+  return err;
+}
+
+
+void MakeFit(Int_t i, TCut cutI, TString  aName){
+  //
+  //
+  //
+  Int_t  ntracks=3000000;
+  TStatToolkit toolkit;
+  Double_t chi2=0;
+  Int_t    npoints=0;
+  //
+  TString  fstringB="";   // magnetic part 
+  TString  fstringT="";    // E file dtilting part
+  TString  fstring0="";   // ROC      part 
+  TString  fstringL="";   // linear part
+  //
+  //  Magnetic field map part
+  // 
+  //                           // 0
+  fstringB+="ablx*dr++";       // 1
+  fstringB+="ablx*ky*dr++";    // 2
+  fstringB+="ably*dr++";       // 3
+  fstringB+="ably*ky*dr++";    // 4
+  //
+  // Electric field tilting part
+  // 
+  //
+  fstringT+="(dr1)++";             // 5     - ex
+  fstringT+="(dr1)*ky++";          // 6     - ex
+  fstringT+="(dr1)*phiL++";        // 7     - ey
+  fstringT+="(dr1)*phiL*ky++";     // 8     - ey
+  //
+  // E field close to the ROC - radius independent part  
+  //                          
+  // 
+  fstring0+="ca++";           // 9
+  fstring0+="sa++";           // 10
+  fstring0+="ky++";           // 11 
+  fstring0+="ky*ca++";        // 12
+  fstring0+="ky*sa++";        // 13
+  //
+  // E field close to the ROC - radius dependent part
+  //
+  fstring0+="r++";              // 14
+  fstring0+="ca*r++";           // 15
+  fstring0+="sa*r++";           // 16
+  fstring0+="ky*r++";           // 17
+  fstring0+="ky*ca*r++";        // 18
+  fstring0+="ky*sa*r++";        // 19
+  //
+  // E field rotation - drift length linear part
+  //  
+  fstringL+="ca*dr1++";         //20
+  fstringL+="sa*dr1++";         //21
+  fstringL+="ky*ca*dr1++";      //22 
+  fstringL+="ky*sa*dr1++";      //23
+  //
+  TString fstringBT   = fstringB +fstringT;
+  TString fstringBT0  = fstringBT+fstring0; 
+  TString fstringA    = fstringBT0+fstringL; 
+  //
+  //
+  TCut cutF=cutI;
+  TString * strFit = 0;
+  //
+  //
+  cutF=cutR+cutI;
+  //
+  strFit = TStatToolkit::FitPlane(chain,"dy:0.1", fstringB.Data(),cutF, chi2,npoints,fitParamB[i],covarB[i],1,0, ntracks);
+  chain->SetAlias(aName+"_FB",strFit->Data());
+  cutF=cutR+cutI + Form("abs(dy-%s)<0.2",(aName+"_FB").Data());
+  //
+  //
+  //
+  strFit = TStatToolkit::FitPlane(chain,"dy:0.1", fstringB.Data(),cutF, chi2,npoints,fitParamB[i],covarB[i],1,0, ntracks);
+  chain->SetAlias(aName+"_FB",strFit->Data()); 
+  chi2B[i]=TMath::Sqrt(chi2/npoints);
+  printf("B fit sqrt(Chi2/npoints)=%f\n",chi2B[i]);
+  //
+  strFit = TStatToolkit::FitPlane(chain,"dy:0.1", fstringBT.Data(),cutF, chi2,npoints,fitParamBT[i],covarBT[i],1,0, ntracks);
+  chain->SetAlias(aName+"_FBT",strFit->Data()); 
+  chi2BT[i]=TMath::Sqrt(chi2/npoints);
+  printf("BT fit sqrt(Chi2/npoints)=%f\n",chi2BT[i]);
+  //
+  strFit = TStatToolkit::FitPlane(chain,"dy:0.1", fstringBT0.Data(),cutF, chi2,npoints,fitParamBT0[i],covarBT0[i],1,0, ntracks);
+  chain->SetAlias(aName+"_FBT0",strFit->Data()); 
+  chi2BT0[i]=TMath::Sqrt(chi2/npoints);
+  printf("BT0 Fit: sqrt(Chi2/npoints)=%f\n",chi2BT0[i]);
+  //
+  //
+  strFit = TStatToolkit::FitPlane(chain,"dy:0.1", fstringA.Data(),cutF, chi2,npoints,fitParamA[i],covarA[i],1, 0, ntracks);
+  chain->SetAlias(aName+"_FA",strFit->Data()); 
+  chi2A[i]=TMath::Sqrt(chi2/npoints);
+  printf("A fit: sqrt(Chi2/npoints)=%f\n",chi2A[i]);
+  errB[i]   =(MakeErrVector(covarB[i]));
+  errBT[i]  =(MakeErrVector(covarBT[i]));
+  errBT0[i] =(MakeErrVector(covarBT0[i]));
+  errA[i]   =(MakeErrVector(covarA[i]));
+}
+
+
+
+
+
+void MakeFitDy(){
+  //
+  //
+  //
+
+
+  MakeFit(0,cutAM5,"dyAM5");
+  MakeFit(1,cutAP2,"dyAP2");
+  MakeFit(2,cutAP5,"dyAP5");
+  //
+  MakeFit(3,cutCM5,"dyCM5");
+  MakeFit(4,cutCP2,"dyCP2");
+  //  MakeFit(5,cutCP5,"dyAP5");
+  DumpFit();
+}
+
+
+void DrawPhi(){
+  //
+  //
+  //
+  chain->Draw("(dy-dyAM5_FA):LTr.fVecPhi.fElements>>hisAM5(60,-3.14,3.14,100,-0.2,0.2)",cutAM5,"");
+  chain->Draw("(dy-dyAP5_FA):LTr.fVecPhi.fElements>>hisAP5(60,-3.14,3.14,100,-0.2,0.2)",cutAP5,"");
+  chain->Draw("(dy-dyAP2_FA):LTr.fVecPhi.fElements>>hisAP2(60,-3.14,3.14,100,-0.2,0.2)",cutAP2,"");
+  hisAM5->FitSlicesY(0,0,-1,20);
+  hisAP5->FitSlicesY(0,0,-1,20);
+  hisAP2->FitSlicesY(0,0,-1,20);
+  hisAM5_1->SetMinimum(-0.2);
+  hisAM5_1->SetMaximum(0.2);
+  hisAM5_1->SetMarkerStyle(20);
+  hisAP5_1->SetMarkerStyle(21);
+  hisAP2_1->SetMarkerStyle(22);
+  hisAM5_1->SetMarkerColor(1);
+  hisAP5_1->SetMarkerColor(2);
+  hisAP2_1->SetMarkerColor(4);
+  hisAM5_1->Draw();
+  hisAP5_1->Draw("same");
+  hisAP2_1->Draw("same");
+}
+
+void MakeGraphs(){
+  Double_t bz[3] ={-5,2,5};
+  Double_t p0A[3]={fitParam[0][0],fitParam[1][0],fitParam[2][0]};
+  Double_t p1A[3]={fitParam[0][1],fitParam[1][1],fitParam[2][1]};
+  Double_t p2A[3]={fitParam[0][2],fitParam[1][2],fitParam[2][2]};
+  Double_t p3A[3]={fitParam[0][3],fitParam[1][3],fitParam[2][3]};
+  Double_t p4A[3]={fitParam[0][4],fitParam[1][4],fitParam[2][4]};
+  Double_t p5A[3]={fitParam[0][5],fitParam[1][5],fitParam[2][5]};
+  Double_t p6A[3]={fitParam[0][6],fitParam[1][6],fitParam[2][6]};
+  Double_t p7A[3]={fitParam[0][7],fitParam[1][7],fitParam[2][7]};
+  Double_t p8A[3]={fitParam[0][8],fitParam[1][8],fitParam[2][8]};
+  Double_t p9A[3]={fitParam[0][9],fitParam[1][9],fitParam[2][9]};
+  TGraph grA0(3,bz,p0A);
+  TGraph grA1(3,bz,p1A);
+  TGraph grA2(3,bz,p2A);
+  TGraph grA3(3,bz,p3A);
+  TGraph grA4(3,bz,p4A);
+  TGraph grA5(3,bz,p5A);
+  TGraph grA6(3,bz,p6A);
+  TGraph grA7(3,bz,p7A);
+  TGraph grA8(3,bz,p8A);
+  TGraph grA9(3,bz,p9A);
+  TF1 f1("f1","[0]*x/(1+([0]*x)^2)");
+  TF1 f2("f2","([0]*x)^2/(1+([0]*x)^2)");
+  TMatrixD matB(5,5);
+  TMatrixD matBT(6,5);
+  TMatrixD matBT0(6,5);
+  for (Int_t i=0; i<5;i++) for (Int_t j=0; j<5;j++) matB[i][j]=fitParamB[j][i];
+  for (Int_t i=0; i<5;i++) for (Int_t j=0; j<5;j++) matBT[i][j]=fitParamBT[j][i];
+  for (Int_t i=0; i<6;i++) for (Int_t j=0; j<5;j++) matB0[i][j]=fitParamBT0[j][i];
+}
+
+
+void DumpFit(){
+  //
+  //
+  //
+  TTreeSRedirector *pcstream = new TTreeSRedirector("exbFits.root");
+
+  for (Int_t i=0; i<5;i++){
+    Double_t bz=0;
+    Double_t side=0;
+    if (i==0) { bz=-5; side=1;}
+    if (i==1) { bz=2;  side=1;}
+    if (i==2) { bz=5;  side=1;}
+    if (i==3) { bz=-5; side=-1;}
+    if (i==4) { bz=2;  side=-1;}
+    (*pcstream)<<"fit"<<
+      "bz="<<bz<<
+      "side="<<side<<
+      "pb.="<<&fitParamB[i]<<
+      "pbT.="<<&fitParamBT[i]<<
+      "pbT0.="<<&fitParamBT0[i]<<
+      "pbA.="<<&fitParamA[i]<<
+      "eb.="<<errB[i]<<
+      "ebT.="<<errBT[i]<<
+      "ebT0.="<<errBT0[i]<<
+      "ebA.="<<errA[i]<<
+      "\n";
+  }
+  delete pcstream;
+}
+
+void MakeFitPic(){
+  TFile f("exbFits.root");
+
+}