]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TPC/AliTPCcalibLaser.cxx
Using masses form the PDG
[u/mrichter/AliRoot.git] / TPC / AliTPCcalibLaser.cxx
index a13bb7b6544df8e6bf249640df59f031e51f1696..4a91299d4ae0136c2b7e5c07483e6e7828f966a8 100644 (file)
   .x ~/NimStyle.C
   gSystem->Load("libANALYSIS");
   gSystem->Load("libTPCcalib");
-  TFile fcalib("CalibObjects.root");
-  TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
-  AliTPCcalibLaser * laser = ( AliTPCcalibLaser *)array->FindObject("laserTPC");
-  laser->DumpMeanInfo(-0,0)
+  TFile fcalib("CalibObjectsTrain2.root");
+  AliTPCcalibLaser * laser = ( AliTPCcalibLaser *)fcalib->Get("laserTPC");
+  laser->DumpMeanInfo(run)
   TFile fmean("laserMean.root")
   //
   //  laser track clasification;
   //
   TCut cutT("cutT","abs(Tr.fP[3])<0.06");
-  TCut cutPt("cutPt","abs(Tr.fP[4])<0.1");
   TCut cutN("cutN","fTPCncls>70");
   TCut cutP("cutP","abs(atan2(x1,x0)-atan2(lx1,lx0))<0.03")
   TCut cutA = cutT+cutPt+cutP;
   gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
   gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
   AliXRDPROOFtoolkit tool;
-  TChain * chainDrift = tool.MakeChain("laser.txt","driftv",0,10200);
+  AliXRDPROOFtoolkit::FilterList("laserDebug.list","* driftvN",1) 
+  TChain * chainDrift = tool.MakeChainRandom("laser.txt.Good","driftv",0,50);
   chainDrift->Lookup();
-  TChain * chainDriftN = tool.MakeChain("laser.txt","driftvN",0,10200);
+  TChain * chainDriftN = tool.MakeChainRandom("laserDebug.list.Good","driftvN",0,300);
   chainDriftN->Lookup();
 
 
-
   TChain * chain = tool.MakeChain("laser.txt","Residuals",0,10200);
   chain->Lookup();
   TChain * chainFit = tool.MakeChain("laser.txt","FitModels",0,10200);
   chainFit->Lookup();
+  TChain * chainTrack = tool.MakeChainRandom("laser.txt","Track",0,30);
+  chainTrack->Lookup();
 
 */
 
 #include "TH2F.h"
 #include "TStatToolkit.h"
 #include "TROOT.h"
+#include "TDatabasePDG.h"
 
 
 #include "TTreeStream.h"
 #include "TTimeStamp.h"
 #include "AliDCSSensorArray.h"
 #include "AliDCSSensor.h"
+#include "AliGRPObject.h"
 
 using namespace std;
 
@@ -196,8 +198,10 @@ AliTPCcalibLaser::AliTPCcalibLaser():
   fDeltaZres(336),   //->2D histo fo residuals
   fDeltaYres2(336),   //->2D histo of residuals
   fDeltaZres2(336),   //->2D histo fo residuals
-  fDeltaYresAbs(336),   //->2D histo of residuals
-  fDeltaZresAbs(336),   //->2D histo of residuals
+  fDeltaYresAbs(336), //->2D histo of residuals
+  fHisYAbsErrors(0),  //-> total errors per beam in the abs res y analysis
+  fDeltaZresAbs(336), //->2D histo of residuals
+  fHisZAbsErrors(0),  //-> total errors per beam in the abs res z analysis
   //fDeltaYres3(336),   //->2D histo of residuals
   //fDeltaZres3(336),   //->2D histo fo residuals
   fFitAside(new TVectorD(5)),
@@ -217,7 +221,12 @@ AliTPCcalibLaser::AliTPCcalibLaser():
   fBeamSlopeZOuter(336),
   fBeamOffsetZInner(336),
   fBeamSlopeZInner(336),
-  fInverseSlopeZ(kTRUE)
+  fInverseSlopeZ(kTRUE),
+  fUseFixedDriftV(0),
+  fFixedFitAside0(0.0),
+  fFixedFitAside1(1.0),
+  fFixedFitCside0(0.0),
+  fFixedFitCside1(1.0)
 {
   //
   // Constructor
@@ -295,7 +304,9 @@ AliTPCcalibLaser::AliTPCcalibLaser(const Text_t *name, const Text_t *title, Bool
   fDeltaYres2(336),
   fDeltaZres2(336),  
   fDeltaYresAbs(336),
+  fHisYAbsErrors(0),
   fDeltaZresAbs(336),
+  fHisZAbsErrors(0),
   //  fDeltaYres3(336),
   //fDeltaZres3(336),  
   fFitAside(new TVectorD(5)),        // drift fit - A side
@@ -315,7 +326,12 @@ AliTPCcalibLaser::AliTPCcalibLaser(const Text_t *name, const Text_t *title, Bool
   fBeamSlopeZOuter(336),
   fBeamOffsetZInner(336),
   fBeamSlopeZInner(336),
-  fInverseSlopeZ(kTRUE)
+  fInverseSlopeZ(kTRUE),
+  fUseFixedDriftV(0),
+  fFixedFitAside0(0.0),
+  fFixedFitAside1(1.0),
+  fFixedFitCside0(0.0),
+  fFixedFitCside1(1.0)
 {
   SetName(name);
   SetTitle(title);
@@ -394,7 +410,9 @@ AliTPCcalibLaser::AliTPCcalibLaser(const AliTPCcalibLaser& calibLaser):
   fDeltaYres2(((calibLaser.fDeltaYres))),
   fDeltaZres2(((calibLaser.fDeltaZres))),  
   fDeltaYresAbs(((calibLaser.fDeltaYresAbs))),
+  fHisYAbsErrors(new TH1F(*(calibLaser.fHisYAbsErrors))),
   fDeltaZresAbs(((calibLaser.fDeltaZresAbs))),
+  fHisZAbsErrors(new TH1F(*(calibLaser.fHisZAbsErrors))),
   //  fDeltaYres3(((calibLaser.fDeltaYres))),
   //fDeltaZres3(((calibLaser.fDeltaZres))),  
   fFitAside(new TVectorD(5)),        // drift fit - A side
@@ -414,8 +432,12 @@ AliTPCcalibLaser::AliTPCcalibLaser(const AliTPCcalibLaser& calibLaser):
   fBeamSlopeZOuter(336),
   fBeamOffsetZInner(336),
   fBeamSlopeZInner(336),
-  fInverseSlopeZ(calibLaser.fInverseSlopeZ)
-
+  fInverseSlopeZ(calibLaser.fInverseSlopeZ),
+  fUseFixedDriftV(calibLaser.fUseFixedDriftV),
+  fFixedFitAside0(calibLaser.fFixedFitAside0),
+  fFixedFitAside1(calibLaser.fFixedFitAside1),
+  fFixedFitCside0(calibLaser.fFixedFitCside0),
+  fFixedFitCside1(calibLaser.fFixedFitCside1)
 {
   //
   // copy constructor
@@ -506,8 +528,10 @@ AliTPCcalibLaser::~AliTPCcalibLaser() {
 
   fDeltaYres.SetOwner();
   fDeltaYres.Delete();
+  delete fHisYAbsErrors;
   fDeltaZres.SetOwner();
   fDeltaZres.Delete();
+  delete fHisZAbsErrors;
   fDeltaYres2.SetOwner();
   fDeltaYres2.Delete();
   fDeltaZres2.SetOwner();
@@ -526,7 +550,9 @@ void AliTPCcalibLaser::Process(AliESDEvent * event) {
   //
   // Loop over tracks and call  Process function
   //
-  Int_t kMinTracks=20;
+  const Int_t  kMinTracks=20;
+  const Int_t  kMinClusters=40;
+
   fESD = event;
   if (!fESD) {
     return;
@@ -537,6 +563,22 @@ void AliTPCcalibLaser::Process(AliESDEvent * event) {
   }
   if (fESD->GetNumberOfTracks()<kMinTracks) return; //not enough tracks
   AliDebug(4,Form("Event number in current file: %d",event->GetEventNumberInFile()));
+  //
+  // find CE background if present
+  //
+  if (AliTPCLaserTrack::GetTracks()==0) AliTPCLaserTrack::LoadTracks();
+  TH1D hisCE("hhisCE","hhisCE",100,-100,100);  
+  for (Int_t i=0;i<fESD->GetNumberOfTracks();++i) {
+    AliESDtrack *track=fESD->GetTrack(i);
+    if (!track) continue;
+    hisCE.Fill(track->GetZ());
+    hisCE.Fill(track->GetZ()+2);
+    hisCE.Fill(track->GetZ()-2);
+  }
+  //
+  //
+
+
   fTracksTPC.Clear();
   fTracksEsd.Clear();
   fTracksEsdParam.Delete();
@@ -552,23 +594,25 @@ void AliTPCcalibLaser::Process(AliESDEvent * event) {
     AliESDfriendTrack *friendTrack=fESDfriend->GetTrack(i);
     if (!friendTrack) continue;
     AliESDtrack *track=fESD->GetTrack(i);
+    Double_t binC = hisCE.GetBinContent(hisCE.FindBin(track->GetZ()));
+    if (binC>336) continue; //remove CE background
     TObject *calibObject=0;
     AliTPCseed *seed=0;
     for (Int_t j=0;(calibObject=friendTrack->GetCalibObject(j));++j)
       if ((seed=dynamic_cast<AliTPCseed*>(calibObject)))
        break;
-    if (track&&seed &&TMath::Abs(track->Pt()) >1 ) {
+    if (track&&seed&&track->GetTPCNcls()>kMinClusters && seed->GetNumberOfClusters() >kMinClusters) {
       //filter CE tracks
       Int_t id = FindMirror(track,seed);
-      if (id>0) counter++;
+      if (id>=0) counter++;
     }
     //
   } 
   fNtracks=counter;
   if (counter<kMinTracks) return;
 
-  FitDriftV();
-  FitDriftV(0.5);
+  //FitDriftV();
+  FitDriftV(0.2);
   if (!fFullCalib) return;
   static Bool_t init=kFALSE;
   if (!init){
@@ -879,6 +923,14 @@ Bool_t  AliTPCcalibLaser::FitDriftV(Float_t minFraction){
   // Fit corrections to the drift velocity - linear approximation in the z and global y
   //The transfromatiom from the drift time to the z position done in AliTPCTracnsform class
   // 
+  // Source of outlyers : 
+  // 0. Track in the saturation - postpeak
+  // 1. gating grid close the part of the signal for first bundle
+
+  // The robust fit is performed in 2 itterations /robust fraction controlled by kFraction/
+  // 1. Robust fit is used in the itteration number 0
+  // only fraction of laser uted
+  // 2. Only the tracks close to the fit used in the second itteration
   /*
     Formulas:
     
@@ -904,14 +956,15 @@ Bool_t  AliTPCcalibLaser::FitDriftV(Float_t minFraction){
     dzs/dl = dz/dl +s*s*vr*dz/dl 
     d(dz/dl) = vr*dz/dl     
   */
-  const Int_t   knLaser      = 336;    //n laser tracks
-  const Float_t kFraction    = 0.9;   // robust fit fraction
-  const Float_t kSaturCut    = 0.05;   // remove saturated lasers - cut on fraction of saturated 
-  const Float_t kDistCut     = 6;      // distance sigma cut
-  const Float_t kDistCutAbs  = 0.25;  
-  const Float_t kMinClusters = 60;     // minimal amount of the clusters
-  const Float_t kMinSignal   = 16;     // minimal mean height of the signal
-  const Float_t kChi2Cut     = 0.1;    // chi2 cut to accept drift fit
+  const Int_t   knLaser      = 336;     //n laser tracks
+  const Float_t kFraction[2] = {0.70,0.95};    // robust fit fraction
+
+  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 = 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");
   static TLinearFitter fdriftC(3,"hyp2");
@@ -940,19 +993,45 @@ Bool_t  AliTPCcalibLaser::FitDriftV(Float_t minFraction){
   Int_t npointsAC=0;
   Int_t nbA[4]={0,0,0,0};
   Int_t nbC[4]={0,0,0,0};
+  TVectorD vecZM(336);      // measured z potion of laser
+  TVectorD vecA(336);       // accepted laser
+  TVectorD vecZF(336);      // fitted position  
+  TVectorD vecDz(336);      // deltaZ
+  TVectorD vecZS(336);      // surveyed position of laser
+  // additional variable to cut
+  TVectorD vecdEdx(336);    // dEdx
+  TVectorD vecSy(336);      // shape y
+  TVectorD vecSz(336);      // shape z
+  //
   //
   for (Int_t id=0; id<336; id++){
-    if (!fTracksEsdParam.At(id)) continue;
-    if (!AcceptLaser(id)) continue;
-    if ( fClusterSatur[id]>kSaturCut)  continue;
-    if ( fClusterCounter[id]<kMinClusters)  continue;
-    AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
+    Int_t reject=0;
+    AliTPCLaserTrack *ltrp =
+      (AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
+    AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
+    AliTPCseed * seed = (AliTPCseed*)fTracksTPC.At(id);
+    vecZM(id)= (param==0) ? 0:param->GetZ();
+    vecZS(id)= ltrp->GetZ();
+    vecDz(id)= 0;
+    vecA(id)=1;
+    vecdEdx(id)=(seed)?seed->GetdEdx():0;
+    vecSy(id) =(seed)?seed->CookShape(1):0;
+    vecSz(id) =(seed)?seed->CookShape(2):0;
+    //
+    fFitZ[id]=0;
+    if (!fTracksEsdParam.At(id)) reject|=1;
+    if (!param) continue; 
+    if (!AcceptLaser(id))        reject|=2;
+    if ( fClusterSatur[id]>kSaturCut)  reject|=4;
+    if ( fClusterCounter[id]<kMinClusters)  reject|=8;
+    vecA(id)=reject;
+    if (reject>0) continue;
     if (ltrp->GetSide()==0){
       npointsA++;
       nbA[ltrp->GetBundle()]++;
     }
     if (ltrp->GetSide()==1){
-      npointsA++;
+      npointsC++;
       nbC[ltrp->GetBundle()]++;
     }
   }
@@ -978,7 +1057,8 @@ Bool_t  AliTPCcalibLaser::FitDriftV(Float_t minFraction){
   for (Int_t iter=0; iter<2; iter++){
     fdriftA.ClearPoints();
     fdriftC.ClearPoints();
-    fdriftAC.ClearPoints();
+    fdriftAC.ClearPoints(); 
+    npointsA=0;  npointsC=0;  npointsAC=0;
     //
     for (Int_t id=0; id<336; id++){
       if (!fTracksEsdParam.At(id)) continue;
@@ -989,7 +1069,6 @@ Bool_t  AliTPCcalibLaser::FitDriftV(Float_t minFraction){
       if (track->GetTPCsignal()<kMinSignal) continue;
       AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
       AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
-
       Double_t xyz[3];
       Double_t pxyz[3];
       Double_t lxyz[3];
@@ -1000,14 +1079,22 @@ Bool_t  AliTPCcalibLaser::FitDriftV(Float_t minFraction){
       ltrp->GetPxPyPz(lpxyz);
       Float_t sz = (ltrp->GetSide()==0) ? TMath::Sqrt(chi2A): TMath::Sqrt(chi2C);
       if (npointsAC>0) sz =TMath::Sqrt(chi2AC); 
-      if (TMath::Abs(fFitZ[id])>sz*kDistCut) continue;
+      if (iter>0 && TMath::Abs(fFitZ[id])>sz*kDistCut) continue;
       if (iter>0 && TMath::Abs(fFitZ[id])>kDistCutAbs) continue;
 
-      // drift distance
-      Double_t zlength =  tpcparam->GetZLength(0);
-      Double_t ldrift = zlength-TMath::Abs(lxyz[2]);
-      Double_t mdrift = zlength-TMath::Abs(xyz[2]);
+//       // drift distance
+//       Double_t zlength =  tpcparam->GetZLength(0);
+//       Double_t ldrift  = zlength-TMath::Abs(lxyz[2]);
+//       Double_t mdrift  = zlength-TMath::Abs(xyz[2]);
+      //
+      Double_t zlength = (ltrp->GetSide()==0)? tpcparam->GetZLength(36): tpcparam->GetZLength(71);
+      Double_t ldrift  = (ltrp->GetSide()==0)?zlength-lxyz[2]:lxyz[2]+zlength;
+      Double_t mdrift  = (ltrp->GetSide()==0)?zlength-xyz[2]:xyz[2]+zlength;
+
+
       Double_t xxx[2] = {ldrift,lxyz[1]*ldrift/(zlength*250.)};
+      if (iter==0 &&ltrp->GetBundle()==0) continue;  
+      // skip bundle 0 - can be wrong in the 0 iteration
       if (ltrp->GetSide()==0){
        fdriftA.AddPoint(xxx,mdrift,1);
       }else{
@@ -1022,7 +1109,7 @@ Bool_t  AliTPCcalibLaser::FitDriftV(Float_t minFraction){
       fdriftA.Eval();
       npointsA= fdriftA.GetNpoints();
       chi2A = fdriftA.GetChisquare()/fdriftA.GetNpoints();
-      fdriftA.EvalRobust(kFraction);
+      fdriftA.EvalRobust(kFraction[iter]);
       fdriftA.GetParameters(fitA);
       if (chi2A<kChi2Cut ||(*fFitAside)[0]==0 ) {
        if (fFitAside->GetNoElements()<5) fFitAside->ResizeTo(5);
@@ -1037,7 +1124,7 @@ Bool_t  AliTPCcalibLaser::FitDriftV(Float_t minFraction){
       fdriftC.Eval();
       npointsC= fdriftC.GetNpoints();
       chi2C = fdriftC.GetChisquare()/fdriftC.GetNpoints();
-      fdriftC.EvalRobust(kFraction);
+      fdriftC.EvalRobust(kFraction[iter]);
       fdriftC.GetParameters(fitC);
       if (chi2C<kChi2Cut||(*fFitCside)[0]==0) {        
        if (fFitCside->GetNoElements()<5) fFitCside->ResizeTo(5);
@@ -1049,11 +1136,11 @@ Bool_t  AliTPCcalibLaser::FitDriftV(Float_t minFraction){
       }
     }
 
-    if (fdriftAC.GetNpoints()>minFraction*knLaser){
+    if (fdriftAC.GetNpoints()>minFraction*knLaser &&npointsA>0.5*minFraction*knLaser&&npointsC>0.5*minFraction*knLaser){
       fdriftAC.Eval();
       npointsAC= fdriftAC.GetNpoints();
       chi2AC = fdriftAC.GetChisquare()/fdriftAC.GetNpoints();
-      fdriftAC.EvalRobust(kFraction);
+      fdriftAC.EvalRobust(kFraction[iter]);
       fdriftAC.GetParameters(fitAC);
       if (chi2AC<kChi2Cut||(*fFitACside)[0]==0) (*fFitACside) = fitAC;
       (*fFitACside)[0] = fitAC[0];
@@ -1076,9 +1163,13 @@ Bool_t  AliTPCcalibLaser::FitDriftV(Float_t minFraction){
       param->GetPxPyPz(pxyz);
       ltrp->GetXYZ(lxyz);
       ltrp->GetPxPyPz(lpxyz); 
-      Double_t zlength =  tpcparam->GetZLength(0);
-      Double_t ldrift = zlength-TMath::Abs(lxyz[2]);
-      Double_t mdrift = zlength-TMath::Abs(xyz[2]);
+      //Double_t zlength =  tpcparam->GetZLength(0);
+      //Double_t ldrift = zlength-TMath::Abs(lxyz[2]);
+      //Double_t mdrift = zlength-TMath::Abs(xyz[2]);
+      Double_t zlength = (ltrp->GetSide()==0)? tpcparam->GetZLength(36): tpcparam->GetZLength(71);
+      Double_t ldrift  = (ltrp->GetSide()==0)?zlength-lxyz[2]:lxyz[2]+zlength;
+      Double_t mdrift  = (ltrp->GetSide()==0)?zlength-xyz[2]:xyz[2]+zlength;
+
 
       Float_t fz =0;
       if (ltrp->GetSide()==0){
@@ -1087,9 +1178,11 @@ Bool_t  AliTPCcalibLaser::FitDriftV(Float_t minFraction){
        fz = (fitC)[0]+(fitC)[1]*ldrift+(fitC)[2]*lxyz[1]*ldrift/(zlength*250.);        
       }
       if (npointsAC>10){
-       fz = (fitAC)[0]+(fitAC)[1]*ltrp->GetSide()+(fitAC)[2]*ldrift+(fitAC)[3]*lxyz[1]*ldrift/(zlength*250.);
+       //fz = (fitAC)[0]+(fitAC)[1]*ltrp->GetSide()+(fitAC)[2]*ldrift+(fitAC)[3]*lxyz[1]*ldrift/(zlength*250.);
       }
       fFitZ[id]=mdrift-fz;
+      vecZF[id]=fz;
+      vecDz[id]=mdrift-fz;
     }
     if (fStreamLevel>0){
       TTreeSRedirector *cstream = GetDebugStreamer();
@@ -1125,6 +1218,15 @@ Bool_t  AliTPCcalibLaser::FitDriftV(Float_t minFraction){
          "vecGoofie.="<<&vecGoofie<<
          //
          //
+         "vecZM.="<<&vecZM<<   // measured z position
+         "vecZS.="<<&vecZS<<   // surveyed z position
+         "vecZF.="<<&vecZF<<   // fitted   z position
+         "vecDz.="<<&vecDz<<   // fitted   z position
+         "vecA.="<<&vecA<<     // accept laser flag
+         "vecdEdx.="<<&vecdEdx<<   // dEdx  - to cut on
+         "vecSy.="<<&vecSy<<       // shape y - to cut on
+         "vecSz.="<<&vecSz<<       // shape z - to cut on
+         //
          "iter="<<iter<<
          "driftA.="<<fFitAside<<
          "driftC.="<<fFitCside<<
@@ -1136,6 +1238,18 @@ Bool_t  AliTPCcalibLaser::FitDriftV(Float_t minFraction){
          "nC="<<npointsC<<
          "nAC="<<npointsAC<<
          "\n";
+       /*
+         //
+         variables to check in debug mode:
+         //
+         chainDriftN->SetAlias("driftS","250-abs(vecZS.fElements)");
+         chainDriftN->SetAlias("driftM","250-abs(vecZM.fElements)");
+         chainDriftN->SetAlias("driftF","vecZF.fElements");      
+         chainDriftN->SetAlias("deltaZ","driftF-driftM");  //deltaZ
+         TCut cutA="iter==1&&sqrt(chi2A)<0.1&&vecZS.fElements>0&&vecA.fElements==0";
+         TCut cutC="iter==1&&sqrt(chi2C)<0.1&&vecZS.fElements<0&&vecA.fElements==0";
+         
+        */
       }
     }
   }
@@ -1185,7 +1299,6 @@ Bool_t  AliTPCcalibLaser::AcceptLaser(Int_t id){
   TCut cutP1("cutP1","abs(LTr.fP[1]-Tr.fP[1])<30");
   TCut cutP2("cutP2","abs(LTr.fP[2]-Tr.fP[2])<0.03");
   TCut cutP3("cutP3","abs(Tr.fP[3])<0.05");
-  TCut cutP4("cutPt","abs(Tr.fP[4])<0.1");
 
   TCut cutA = cutP0+cutP1+cutP2+cutP3+cutP4;
   */
@@ -1200,7 +1313,6 @@ Bool_t  AliTPCcalibLaser::AcceptLaser(Int_t id){
   if (TMath::Abs(param->GetParameter()[1]-ltrp->GetParameter()[1])>30) return kFALSE;    // cutZ -P1
   if (TMath::Abs(param->GetParameter()[2]-ltrp->GetParameter()[2])>0.03) return kFALSE;  // cut -P2
   if (TMath::Abs(param->GetParameter()[3])>0.05) return kFALSE;   // cut Tl -P3
-  if (TMath::Abs(param->GetParameter()[4])>0.1) return kFALSE;   // cut Pt  -P4
   //
   //
 
@@ -1212,32 +1324,10 @@ Int_t  AliTPCcalibLaser::FindMirror(AliESDtrack *track, AliTPCseed *seed){
   // Find corresponding mirror
   // add the corresponding tracks
  
-
-  Float_t kRadius0  = 252;
-  Float_t kRadius   = 253.4;
-
   if (!track->GetOuterParam()) return -1;
-  AliExternalTrackParam param(*(track->GetOuterParam()));
-  AliTracker::PropagateTrackTo(&param,kRadius0,0.10566,3,kTRUE);
-  AliTracker::PropagateTrackTo(&param,kRadius,0.10566,0.1,kTRUE);
-  AliTPCLaserTrack ltr;
-  AliTPCLaserTrack *ltrp=0x0;
-  //  AliTPCLaserTrack *ltrpjw=0x0;
-  //
-  Int_t id   = AliTPCLaserTrack::IdentifyTrack(&param);
- // Int_t idjw = AliTPCLaserTrack::IdentifyTrackJW(&param);
-  //AliDebug(4,Form("Identified Track: %03d (%03d)",id,idjw));
 
-  if (id!=-1 && (AliTPCLaserTrack::GetTracks()->UncheckedAt(id)))
-    ltrp=(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
-  else
-    ltrp=&ltr;
-
-  if (id<0) return id;
-  fCounter[id]++;
-  //
-  //
-  //
+  Float_t kRadius0  = 252;
+  Float_t kRadius   = 254.2;
   Int_t countercl=0;
   Float_t counterSatur=0;
   Int_t csideA =0;
@@ -1258,15 +1348,41 @@ Int_t  AliTPCcalibLaser::FindMirror(AliESDtrack *track, AliTPCseed *seed){
   //
   //
   //
-  if (csideA>0.1*seed->GetNumberOfClusters() && ltrp->GetSide()==1) return -1;
-  if (csideC>0.1*seed->GetNumberOfClusters() && ltrp->GetSide()==0) return -1;
+  if (csideA<0.9*seed->GetNumberOfClusters() && csideC<0.9*seed->GetNumberOfClusters()) return 0;  // cross laser track can not happen
+
+  Int_t side= 0;
+  if (csideC>0.5*seed->GetNumberOfClusters()) side=1;
+
+
+  AliExternalTrackParam param(*(track->GetOuterParam()));
+  AliTracker::PropagateTrackTo(&param,kRadius0,TDatabasePDG::Instance()->GetParticle("mu+")->Mass(),3,kTRUE);
+  AliTracker::PropagateTrackTo(&param,kRadius,TDatabasePDG::Instance()->GetParticle("mu+")->Mass(),0.1,kTRUE);
+  AliTPCLaserTrack ltr;
+  AliTPCLaserTrack *ltrp=0x0;
+  //  AliTPCLaserTrack *ltrpjw=0x0;
+  //
+  Int_t id   = AliTPCLaserTrack::IdentifyTrack(&param,side);
+ // Int_t idjw = AliTPCLaserTrack::IdentifyTrackJW(&param);
+  //AliDebug(4,Form("Identified Track: %03d (%03d)",id,idjw));
+
+  if (id!=-1 && (AliTPCLaserTrack::GetTracks()->UncheckedAt(id)))
+    ltrp=(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
+  else
+    ltrp=&ltr;
+
+  if (id<0) return -1;
+  if (ltrp->GetSide()!=side) return -1;
+  fCounter[id]++;
+  //
+  //
+  //
   //
   if (counterSatur>fClusterSatur[id]) fClusterSatur[id]=counterSatur;
   //
   //
   Float_t radius=TMath::Abs(ltrp->GetX());
-  AliTracker::PropagateTrackTo(&param,radius,0.10566,0.01,kTRUE);
   param.Rotate(ltrp->GetAlpha());
+  AliTracker::PropagateTrackTo(&param,radius,0.10566,0.01,kFALSE);
   //
   if (!fTracksMirror.At(id)) fTracksMirror.AddAt(ltrp,id);
   Bool_t accept=kTRUE;  
@@ -1994,25 +2110,37 @@ void AliTPCcalibLaser::RefitLaserJW(Int_t id){
        // zreal = (zmeasured - [0] + (1 - [1])*250.0)/[1]
 
        Double_t dzabs = -1000;
+       Double_t zcorrected = -1000;
        if (ltrp->GetSide()==0){
-         if ((*fFitAside)[1]>0.) { 
+         if ((*fFitAside)[1]>0. || fUseFixedDriftV) { 
            // ignore global y dependence for now
-           Double_t zcorrected = (zcl + (*fFitAside)[0] - 
-                                  (1.0-(*fFitAside)[1])*250.0)/(*fFitAside)[1];
+           zcorrected = 0;         
+           if(!fUseFixedDriftV) 
+             zcorrected = (zcl + (*fFitAside)[0] - 
+                           (1.0-(*fFitAside)[1])*250.0)/(*fFitAside)[1];
+           else
+             zcorrected = (zcl + fFixedFitAside0 - 
+                           (1.0-fFixedFitAside1)*250.0)/fFixedFitAside1;
            //      zcorrected = zcl;
-           if(vecSec[irow]==outerSector && (outerSector%36)==fBeamSectorOuter[id])
+           if(vecSec[irow]==outerSector && outerSector==fBeamSectorOuter[id])
              dzabs = zcorrected -fBeamSlopeZOuter[id]*x -fBeamOffsetZOuter[id];
-           else if((innerSector%36)==fBeamSectorInner[id])
+           else if(innerSector==fBeamSectorInner[id])
              dzabs = zcorrected -fBeamSlopeZInner[id]*x -fBeamOffsetZInner[id];
          }
        } else {
-         if ((*fFitCside)[1]>0.) {
-           Double_t zcorrected = (zcl - (*fFitCside)[0] + 
-                                  (1.0-(*fFitCside)[1])*250.0)/(*fFitCside)[1];
+         if ((*fFitCside)[1]>0. || fUseFixedDriftV) {
+           
+           if(!fUseFixedDriftV) 
+             zcorrected = (zcl - (*fFitCside)[0] + 
+                           (1.0-(*fFitCside)[1])*250.0)/(*fFitCside)[1];
+           else
+             zcorrected = (zcl - fFixedFitCside0 + 
+                           (1.0-fFixedFitCside1)*250.0)/fFixedFitCside1;
+           
            //      zcorrected = zcl;
-           if(vecSec[irow]==outerSector && (outerSector%36)==fBeamSectorOuter[id])
+           if(vecSec[irow]==outerSector && outerSector==fBeamSectorOuter[id])
              dzabs = zcorrected -fBeamSlopeZOuter[id]*x -fBeamOffsetZOuter[id];
-           else if((innerSector%36)==fBeamSectorInner[id])
+           else if(innerSector==fBeamSectorInner[id])
              dzabs = zcorrected -fBeamSlopeZInner[id]*x -fBeamOffsetZInner[id];
          }
        }
@@ -2021,9 +2149,10 @@ void AliTPCcalibLaser::RefitLaserJW(Int_t id){
          if (profy){       
              profy->Fill(irow,ycl-yfit);
              profy2->Fill(irow,ycl-yfit2);
-             if(yabsbeam<-100)
-               profyabs->Fill(irow,-0.99);
-             else
+             if(yabsbeam<-100) {
+               fHisYAbsErrors->Fill(id);
+               //              profyabs->Fill(irow,-0.99);
+             } else
                profyabs->Fill(irow,ycl-yabsbeam);
 
              //              profy3->Fill(irow,ycl-yfit3);
@@ -2032,9 +2161,10 @@ void AliTPCcalibLaser::RefitLaserJW(Int_t id){
              profz->Fill(irow,zcl-zfit);
              profz2->Fill(irow,zcl-zfit2);
              //profz3->Fill(irow,zcl-zfit3);
-             if(dzabs<-100)
-               profzabs->Fill(irow,-0.99);
-             else
+             if(dzabs<-100) {
+
+               fHisZAbsErrors->Fill(id);
+             }else
                profzabs->Fill(irow,dzabs);
          }
        }
@@ -2104,7 +2234,7 @@ void AliTPCcalibLaser::RefitLaserJW(Int_t id){
 
 
 
-void AliTPCcalibLaser::DumpMeanInfo(Float_t bfield, Int_t run){
+void AliTPCcalibLaser::DumpMeanInfo(Int_t run){
   //
   //  Dump information about laser beams
   //  isOK variable indicates usability of the beam  
@@ -2119,17 +2249,33 @@ void AliTPCcalibLaser::DumpMeanInfo(Float_t bfield, Int_t run){
   const Float_t krmsCut1=0.16;
   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};         // 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");
   TF1 fg("fg","gaus");
-
+  AliTPCParam  * tpcparam    = 0;   
   // start set up for absolute residuals analysis
-  //AliTPCcalibDB*  calib=AliTPCcalibDB::Instance();
-  //  AliTPCParam  * tpcparam    = calib->GetParameters(); 
-  AliTPCParam  * tpcparam    = 0
+  //
+  AliTPCcalibDB*  calib=AliTPCcalibDB::Instance();
+  tpcparam    = calib->GetParameters()
   if (!tpcparam) tpcparam    = new AliTPCParamSR;
   tpcparam->Update();
+  AliGRPObject *grp = AliTPCcalibDB::GetGRP(run);
+  Float_t current=0;
+  Float_t bfield      = 0, bz=0;
+
+  if (grp){
+    Float_t polarity = (grp->GetL3Polarity()>0) ? -1.:1;
+    current = grp->GetL3Current((AliGRPObject::Stats)0);
+    bfield = polarity*5*current/30000.;
+    bz = polarity*5*current/30000.;
+    printf("Run%d\tL3 current%f\tBz\t%f\n",run,current,bz);
+  }
+
   SetBeamParameters(fBeamOffsetYOuter, fBeamSlopeYOuter, fBeamSectorOuter,0);
   SetBeamParameters(fBeamOffsetYInner, fBeamSlopeYInner, fBeamSectorInner,1);
   TLinearFitter lfabsyInner(2);
@@ -2162,6 +2308,7 @@ void AliTPCcalibLaser::DumpMeanInfo(Float_t bfield, Int_t run){
      AliTPCLaserTrack::LoadTracks();
       ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
     }
+    ltrp->UpdatePoints();
     pcstream->GetFile()->cd();
     if (hisphi)  hisphi->Write();
     if (hisphiP) hisphiP->Write();
@@ -2408,6 +2555,49 @@ void AliTPCcalibLaser::DumpMeanInfo(Float_t bfield, Int_t run){
     TVectorD vecEz(159);       //error z
     TVectorD vecPhi(159);      // local tangent
     TVectorD vecPhiR(159);     // local tangent
+    // 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];
+      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];
+      }
+    }
+
 
     lfabsyInner.ClearPoints();    
     lfabszInner.ClearPoints();    
@@ -2547,7 +2737,8 @@ void AliTPCcalibLaser::DumpMeanInfo(Float_t bfield, Int_t run){
                                  TMath::Max(zprof->GetBinError(bin), 0.001));
          }
        }
-
+       // global position
+       
       }
        
       delete yprof; delete zprof;
@@ -2744,6 +2935,16 @@ void AliTPCcalibLaser::DumpMeanInfo(Float_t bfield, Int_t run){
       "SecOut="<<secOuter<<             // outer sector
       "lasTanPhiLocIn="<<lasTanPhiLocIn<< // laser tan phi in local frame (inner)
       "lasTanPhiLocOut="<<lasTanPhiLocOut<<// laser tan phi in local frame (outer)
+      "ibr.="<<&vecIBR<<   // radial filed integral
+      "ibrphi.="<<&vecIBRPhi<<   // r=phifiled 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 
@@ -3174,6 +3375,9 @@ void   AliTPCcalibLaser::MakeFitHistos(){
   fHisPz2vP2Out= new TH2F("HisPz2vP2Out","HisPz2vP2Out",336,0,336,500,-0.0002,0.0002);
   fHisPz3vP2IO = new TH2F("HisPz3vP2IO", "HisPz3vP2IO",336,0,336,500,-0.0002,0.0002);
   
+  fHisYAbsErrors = new TH1F("HisYabsErrors", "Errors per beam (y)", 336, 0, 336);
+  fHisZAbsErrors = new TH1F("HisZabsErrors", "Errors per beam (z)", 336, 0, 336);
+
   fHisNclIn->SetDirectory(0);      //->Number of clusters inner
   fHisNclOut->SetDirectory(0);     //->Number of clusters outer
   fHisNclIO->SetDirectory(0);      //->Number of cluster inner outer
@@ -3217,6 +3421,10 @@ void   AliTPCcalibLaser::MakeFitHistos(){
   fHisPz2vP2Out->SetDirectory(0);  //-> Curv  P2outer - parabola
   fHisPz3vP2IO->SetDirectory(0);   //-> Curv  P2outerinner - common parabola
 
+  fHisYAbsErrors->SetDirectory(0); //-> total errors per beam in the abs res y analysis
+  fHisZAbsErrors->SetDirectory(0); //-> total errors per beam in the abs res z analysis
+
+  
 
   //
   //
@@ -3237,7 +3445,10 @@ void   AliTPCcalibLaser::MakeFitHistos(){
       profy2=new TH2F(Form("pry%03d",id),Form("Y Residuals for Laser Beam %03d -Parabolic",id),160,0,160,50,-0.5,0.5);
       profy2->SetDirectory(0);
       fDeltaYres2.AddAt(profy2,id);
-      profyabs=new TH2F(Form("pryabs%03d",id),Form("Y Residuals for Laser Beam %03d -Absolute",id),160,0,160,100,-1.0,1.0); // has to be bigger based on earlier studies
+      if(!fUseFixedDriftV)
+       profyabs=new TH2F(Form("pryabs%03d",id),Form("Y Residuals for Laser Beam %03d -Absolute",id),160,0,160,100,-1.0,1.0); // has to be bigger based on earlier studies
+      else
+       profyabs=new TH2F(Form("pryabs%03d",id),Form("Y Residuals for Laser Beam %03d -Absolute",id),160,0,160,200,-2.0,2.0); // has to be bigger based on earlier studies
       profyabs->SetDirectory(0);
       fDeltaYresAbs.AddAt(profyabs,id);
       //profy3=new TH2F(Form("pry%03d",id),Form("Y Residuals for Laser Beam %03d- Parabolic2",id),160,0,160,100,-0.5,0.5);
@@ -3251,7 +3462,10 @@ void   AliTPCcalibLaser::MakeFitHistos(){
       profz2=new TH2F(Form("prz%03d",id),Form("Z Residuals for Laser Beam %03d - Parabolic",id),160,0,160,50,-0.5,0.5);
       profz2->SetDirectory(0);
       fDeltaZres2.AddAt(profz2,id);
-      profzabs=new TH2F(Form("przabs%03d",id),Form("Z Residuals for Laser Beam %03d -Absolute",id),160,0,160,100,-1.0,1.0); // has to be bigger based on earlier studies
+      if(!fUseFixedDriftV)
+       profzabs=new TH2F(Form("przabs%03d",id),Form("Z Residuals for Laser Beam %03d -Absolute",id),160,0,160,100,-1.0,1.0); // has to be bigger based on earlier studies
+      else
+       profzabs=new TH2F(Form("przabs%03d",id),Form("Z Residuals for Laser Beam %03d -Absolute",id),160,0,160,200,-2.0,2.0); // has to be bigger based on earlier studies
       profzabs->SetDirectory(0);
       fDeltaZresAbs.AddAt(profzabs,id);
       //profz3=new TH2F(Form("prz%03d",id),Form("Z Residuals for Laser Beam %03d- Parabolic2",id),160,0,160,100,-0.5,0.5);
@@ -3425,6 +3639,8 @@ void AliTPCcalibLaser::MergeFitHistos(AliTPCcalibLaser * laser){
   fHisPz2vP2In->Add(laser->fHisPz2vP2In  );   //-> Curv  P2inner - parabola
   fHisPz2vP2Out->Add(laser->fHisPz2vP2Out  );  //-> Curv  P2outer - parabola
   fHisPz3vP2IO->Add(laser->fHisPz3vP2IO  );   //-> Curv  P2outerinner - common parabola
+  fHisYAbsErrors->Add(laser->fHisYAbsErrors); //-> total errors per beam in the abs res y analysis
+  fHisZAbsErrors->Add(laser->fHisZAbsErrors); //-> total errors per beam in the abs res z analysis
   //
   //
   //