AliTPCcalibTimeGain.cxx - Adding the Gamma conversion selected electorns
authormarian <marian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 15 Feb 2011 14:45:32 +0000 (14:45 +0000)
committermarian <marian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 15 Feb 2011 14:45:32 +0000 (14:45 +0000)
AliTPCcalibTracks.cxx AliTPCcalibTracks.h - Use THnSparse for cluster residuals - removed obsolete histograms
AliTPCcalibV0.cxx AliTPCcalibV0.h         - high pt V0  filtering

AliTPCCorrection.cxx AliTPCCorrection.h   - Additional arguemnt in function to create track distortion fit tree
AliTPCExBBShape.cxx AliTPCExBBShape.h     - Adding visualization function for magnetic field
AliTPCExBEffectiveSector.cxx  AliTPCExBEffectiveSector.h - Modified MakeResidualMap function
AliTPCkalmanAlign.cxx AliTPCkalmanAlign.h    - Fit constrain functions
AliTPCReconstructor.cxx AliTPCReconstructor.h - Possibility to use the Altro emulator in reconstruction

Marian

13 files changed:
TPC/AliTPCCorrection.cxx
TPC/AliTPCCorrection.h
TPC/AliTPCExBBShape.cxx
TPC/AliTPCExBBShape.h
TPC/AliTPCExBEffectiveSector.cxx
TPC/AliTPCExBEffectiveSector.h
TPC/AliTPCcalibTimeGain.cxx
TPC/AliTPCcalibTracks.cxx
TPC/AliTPCcalibTracks.h
TPC/AliTPCcalibV0.cxx
TPC/AliTPCcalibV0.h
TPC/AliTPCkalmanAlign.cxx
TPC/AliTPCkalmanAlign.h

index cd2cdf5e8ccd52bae27782a95088f0db4a6caacd..a4317b80e371aa0dc4d46c6bb65bcc2cd2c4014a 100644 (file)
@@ -207,6 +207,41 @@ void AliTPCCorrection::DistortPoint(Float_t x[],const Short_t roc) {
   for (Int_t j=0;j<3;++j) x[j]+=dx[j];
 }
 
+void AliTPCCorrection::DistortPointLocal(Float_t x[],const Short_t roc) {
+  //
+  // Distorts the initial coordinates x (cartesian coordinates)
+  // according to the given effect (inherited classes)
+  // roc represents the TPC read out chamber (offline numbering convention)
+  //
+  Float_t gxyz[3]={0,0,0};
+  Double_t alpha = TMath::Pi()*(roc%18+0.5)/18;
+  Double_t ca=TMath::Cos(alpha), sa= TMath::Sin(alpha);
+  gxyz[0]=  ca*x[0]+sa*x[1];
+  gxyz[1]= -sa*x[0]+ca*x[1];
+  gxyz[2]= x[2];
+  DistortPoint(gxyz,roc);
+  x[0]=  ca*gxyz[0]-sa*gxyz[1];
+  x[1]= +sa*gxyz[0]+ca*gxyz[1];
+  x[2]= gxyz[2];
+}
+void AliTPCCorrection::CorrectPointLocal(Float_t x[],const Short_t roc) {
+  //
+  // Distorts the initial coordinates x (cartesian coordinates)
+  // according to the given effect (inherited classes)
+  // roc represents the TPC read out chamber (offline numbering convention)
+  //
+  Float_t gxyz[3]={0,0,0};
+  Double_t alpha = TMath::Pi()*(roc%18+0.5)/18;
+  Double_t ca=TMath::Cos(alpha), sa= TMath::Sin(alpha);
+  gxyz[0]=  ca*x[0]+sa*x[1];
+  gxyz[1]= -sa*x[0]+ca*x[1];
+  gxyz[2]= x[2];
+  CorrectPoint(gxyz,roc);
+  x[0]=  ca*gxyz[0]-sa*gxyz[1];
+  x[1]=  sa*gxyz[0]+ca*gxyz[1];
+  x[2]=  gxyz[2];
+}
+
 void AliTPCCorrection::DistortPoint(const Float_t x[],const Short_t roc,Float_t xp[]) {
   //
   // Distorts the initial coordinates x (cartesian coordinates) and stores the new 
@@ -1345,6 +1380,7 @@ AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(AliExternalTrackPara
   const Double_t kSigmaZ=0.1;
   const Double_t kMaxR=500;
   const Double_t kMaxZ=500;
+  
   const Double_t kMaxZ0=220;
   const Double_t kZcut=3;
   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
@@ -1368,7 +1404,8 @@ AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(AliExternalTrackPara
     xyz[1]+=gRandom->Gaus(0,0.000005);
     xyz[2]+=gRandom->Gaus(0,0.000005);
     if (TMath::Abs(track.GetZ())>kMaxZ0) continue;
-    if (TMath::Abs(track.GetX())>kMaxR) break;
+    if (TMath::Abs(track.GetX())<kRTPC0) continue;
+    if (TMath::Abs(track.GetX())>kRTPC1) continue;
     AliTrackPoint pIn0;                               // space point          
     AliTrackPoint pIn1;
     Int_t sector= (xyz[2]>0)? 0:18;
@@ -1408,7 +1445,11 @@ AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(AliExternalTrackPara
     pointArray0.GetPoint(point1,npoints-21);
     pointArray0.GetPoint(point2,npoints-11);
     pointArray0.GetPoint(point3,npoints-1);
-  }  
+  } 
+  if ((TMath::Abs(point1.GetX()-point3.GetX())+TMath::Abs(point1.GetY()-point3.GetY()))<10){
+    printf("fit points not properly initialized\n");
+    return 0;
+  }
   track0 = AliTrackerBase::MakeSeed(point1, point2, point3);
   track1 = AliTrackerBase::MakeSeed(point1, point2, point3);
   track0->ResetCovariance(10);
@@ -1426,6 +1467,8 @@ AliExternalTrackParam * AliTPCCorrection::FitDistortedTrack(AliExternalTrackPara
     pointArray1.GetPoint(pIn1,ipoint);
     AliTrackPoint prot0 = pIn0.Rotate(track0->GetAlpha());   // rotate to the local frame - non distoted  point
     AliTrackPoint prot1 = pIn1.Rotate(track1->GetAlpha());   // rotate to the local frame -     distorted point
+    if (TMath::Abs(prot0.GetX())<kRTPC0) continue;
+    if (TMath::Abs(prot0.GetX())>kRTPC1) continue;
     //
     if (!AliTrackerBase::PropagateTrackTo(track0,prot0.GetX(),kMass,5,kFALSE,kMaxSnp)) break;
     if (!AliTrackerBase::PropagateTrackTo(track1,prot0.GetX(),kMass,5,kFALSE,kMaxSnp)) break;
@@ -1545,7 +1588,7 @@ TTree* AliTPCCorrection::CreateDistortionTree(Double_t step){
 
 
 
-void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Bool_t debug ){
+void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Int_t offset, Bool_t debug ){
   //
   // Make a fit tree:
   // For each partial correction (specified in array) and given track topology (phi, theta, snp, refX)
@@ -1571,9 +1614,11 @@ void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t
   //
   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
   //  const Double_t kB2C=-0.299792458e-3;
-  const Int_t kMinEntries=10; 
+  const Int_t kMinEntries=20; 
   Double_t phi,theta, snp, mean,rms, entries,sector,dsec;
-  Float_t refX;
+  Float_t refX;  
+  Int_t run;
+  tinput->SetBranchAddress("run",&run);
   tinput->SetBranchAddress("theta",&theta);
   tinput->SetBranchAddress("phi", &phi);
   tinput->SetBranchAddress("snp",&snp);
@@ -1583,7 +1628,7 @@ void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t
   tinput->SetBranchAddress("sector",&sector);
   tinput->SetBranchAddress("dsec",&dsec);
   tinput->SetBranchAddress("refX",&refX);
-  TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortion%d_%d.root",dtype,ptype));
+  TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortion%d_%d_%d.root",dtype,ptype,offset));
   //
   Int_t nentries=tinput->GetEntries();
   Int_t ncorr=corrArray->GetEntries();
@@ -1598,7 +1643,7 @@ void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t
   if (dtype==3) { dir=1;}
   if (dtype==4) { dir=-1;}
   //
-  for (Int_t ientry=0; ientry<nentries; ientry+=step){
+  for (Int_t ientry=offset; ientry<nentries; ientry+=step){
     tinput->GetEntry(ientry);
     if (TMath::Abs(snp)>kMaxSnp) continue;
     tPar[0]=0;
@@ -1631,9 +1676,12 @@ void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t
     Double_t xyz[3];
     track.GetXYZ(xyz);
     Int_t id=0;
+    Double_t pt=1./tPar[4];
     Double_t dRrec=0; // dummy value - needed for points - e.g for laser
     //if (ptype==4 &&bz<0) mean*=-1;  // interpret as curvature -- COMMENTED out - in lookup signed 1/pt used
+    Double_t refXD=refX;
     (*pcstream)<<"fit"<<
+      "run="<<run<<       // run number
       "bz="<<bz<<         // magnetic filed used
       "dtype="<<dtype<<   // detector match type
       "ptype="<<ptype<<   // parameter type
@@ -1644,15 +1692,18 @@ void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t
       "rms="<<rms<<       // rms
       "sector="<<sector<<
       "dsec="<<dsec<<
-      "refX="<<refX<<         // referece X
+      "refX="<<refXD<<         // referece X as double
       "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
+      "pt="<<pt<<            // pt
       "id="<<id<<             // track id
       "entries="<<entries;// number of entries in bin
     //
     Bool_t isOK=kTRUE;
+    if (entries<kMinEntries) isOK=kFALSE;
+    //
     if (dtype!=4) for (Int_t icorr=0; icorr<ncorr; icorr++) {
       AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
       corrections[icorr]=0;
@@ -1680,11 +1731,8 @@ void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t
        }
        //if (ptype==4 &&bz<0) corrections[icorr]*=-1;  // interpret as curvature - commented out
       }      
-      Double_t dRdummy=0;
       (*pcstream)<<"fit"<<
-       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
+       Form("%s=",corr->GetName())<<corrections[icorr];   // dump correction value
     }
   
     if (dtype==4) for (Int_t icorr=0; icorr<ncorr; icorr++) {
@@ -1694,8 +1742,8 @@ void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t
       AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
       corrections[icorr]=0;
       if (entries>kMinEntries){
-       AliExternalTrackParam trackIn0(refX,phi,tPar,cov);
-       AliExternalTrackParam trackIn1(refX,phi,tPar,cov);
+       AliExternalTrackParam trackIn0(refX,phi,tPar,cov); //Outer - direction to vertex
+       AliExternalTrackParam trackIn1(refX,phi,tPar,cov); //Inner - direction magnet 
        AliExternalTrackParam *trackOut0 = 0;
        AliExternalTrackParam *trackOut1 = 0;
        //
@@ -1718,8 +1766,18 @@ void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t
          //
          corrections[icorr] = (trackOut0->GetParameter()[ptype]-trackIn0.GetParameter()[ptype]);
          corrections[icorr]-= (trackOut1->GetParameter()[ptype]-trackIn1.GetParameter()[ptype]);
+         if (isOK)
+           if ((TMath::Abs(trackOut0->GetX()-trackOut1->GetX())>0.1)||
+               (TMath::Abs(trackOut0->GetX()-trackIn1.GetX())>0.1)||
+               (TMath::Abs(trackOut0->GetAlpha()-trackOut1->GetAlpha())>0.00001)||
+               (TMath::Abs(trackOut0->GetAlpha()-trackIn1.GetAlpha())>0.00001)||
+               (TMath::Abs(trackIn0.GetTgl()-trackIn1.GetTgl())>0.0001)||
+               (TMath::Abs(trackIn0.GetSnp()-trackIn1.GetSnp())>0.0001)
+               ){
+             isOK=kFALSE;
+           }             
          delete trackOut0;      
-         delete trackOut1;      
+         delete trackOut1;       
        }else{
          corrections[icorr]=0;
          isOK=kFALSE;
@@ -1727,11 +1785,8 @@ void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t
        //
        //if (ptype==4 &&bz<0) corrections[icorr]*=-1;  // interpret as curvature - commented out no in lookup
       }      
-      Double_t dRdummy=0;
       (*pcstream)<<"fit"<<
-       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
+       Form("%s=",corr->GetName())<<corrections[icorr];   // dump correction value
     }
     //
     (*pcstream)<<"fit"<<"isOK="<<isOK<<"\n";
@@ -1743,13 +1798,189 @@ void AliTPCCorrection::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t
 
 
 
-void AliTPCCorrection::MakeLaserDistortionTree(TTree* tree, TObjArray *corrArray, Int_t itype){
+void AliTPCCorrection::MakeSectorDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Int_t offset, Bool_t debug ){
+  //
+  // Make a fit tree:
+  // For each partial correction (specified in array) and given track topology (phi, theta, snp, refX)
+  // calculates partial distortions
+  // Partial distortion is stored in the resulting tree
+  // Output is storred in the file distortion_<dettype>_<partype>.root
+  // Partial  distortion is stored with the name given by correction name
+  //
+  //
+  // Parameters of function:
+  // input     - input tree
+  // dtype     - distortion type 10 - IROC-OROC 
+  // ppype     - parameter type
+  // corrArray - array with partial corrections
+  // step      - skipe entries  - if 1 all entries processed - it is slow
+  // debug     0 if debug on also space points dumped - it is slow
+
+  const Double_t kMaxSnp = 0.8;  
+  const Int_t kMinEntries=200; 
+  //  AliTPCROC *tpcRoc =AliTPCROC::Instance();  
+  //
+  const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
+  //  const Double_t kB2C=-0.299792458e-3;
+  Double_t phi,theta, snp, mean,rms, entries,sector,dsec,globalZ;
+  Int_t isec1, isec0;
+  Double_t refXD;
+  Float_t refX;
+  Int_t run;
+  tinput->SetBranchAddress("run",&run);
+  tinput->SetBranchAddress("theta",&theta);
+  tinput->SetBranchAddress("phi", &phi);
+  tinput->SetBranchAddress("snp",&snp);
+  tinput->SetBranchAddress("mean",&mean);
+  tinput->SetBranchAddress("rms",&rms);
+  tinput->SetBranchAddress("entries",&entries);
+  tinput->SetBranchAddress("sector",&sector);
+  tinput->SetBranchAddress("dsec",&dsec);
+  tinput->SetBranchAddress("refX",&refXD);
+  tinput->SetBranchAddress("z",&globalZ);
+  tinput->SetBranchAddress("isec0",&isec0);
+  tinput->SetBranchAddress("isec1",&isec1);
+  TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortionSector%d_%d_%d.root",dtype,ptype,offset));
+  //
+  Int_t nentries=tinput->GetEntries();
+  Int_t ncorr=corrArray->GetEntries();
+  Double_t corrections[100]={0}; //
+  Double_t tPar[5];
+  Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
+  Int_t dir=0;
+  //
+  for (Int_t ientry=offset; ientry<nentries; ientry+=step){
+    tinput->GetEntry(ientry);
+    refX=refXD;
+    Int_t id=-1;
+    if (TMath::Abs(TMath::Abs(isec0%18)-TMath::Abs(isec1%18))==0) id=1;  // IROC-OROC - opposite side
+    if (TMath::Abs(TMath::Abs(isec0%36)-TMath::Abs(isec1%36))==0) id=2;  // IROC-OROC - same side
+    if (dtype==10  && id==-1) continue;
+    //
+    dir=-1;
+    tPar[0]=0;
+    tPar[1]=globalZ;
+    tPar[2]=snp;
+    tPar[3]=theta;
+    tPar[4]=(gRandom->Rndm()-0.1)*0.2;  //
+    Double_t pt=1./tPar[4];
+    //
+    printf("%f\t%f\t%f\t%f\t%f\t%f\n",entries, sector,theta,snp, mean,rms);
+    Double_t bz=AliTrackerBase::GetBz();
+    AliExternalTrackParam track(refX,phi,tPar,cov);    
+    Double_t xyz[3],xyzIn[3],xyzOut[3];
+    track.GetXYZ(xyz);
+    track.GetXYZAt(85,bz,xyzIn);    
+    track.GetXYZAt(245,bz,xyzOut);    
+    Double_t phiIn  = TMath::ATan2(xyzIn[1],xyzIn[0]);
+    Double_t phiOut = TMath::ATan2(xyzOut[1],xyzOut[0]);
+    Double_t phiRef = TMath::ATan2(xyz[1],xyz[0]);
+    Int_t sectorRef = TMath::Nint(9.*phiRef/TMath::Pi()-0.5);
+    Int_t sectorIn  = TMath::Nint(9.*phiIn/TMath::Pi()-0.5);
+    Int_t sectorOut = TMath::Nint(9.*phiOut/TMath::Pi()-0.5);
+    //
+    Bool_t isOK=kTRUE; 
+    if (sectorIn!=sectorOut) isOK=kFALSE;  // requironment - cluster in the same sector
+    if (sectorIn!=sectorRef) isOK=kFALSE;  // requironment - cluster in the same sector
+    if (entries<kMinEntries/(1+TMath::Abs(globalZ/100.))) isOK=kFALSE;  // requironment - minimal amount of tracks in bin
+    // Do downscale
+    if (TMath::Abs(theta)>1) isOK=kFALSE;
+    //
+    Double_t dRrec=0; // dummy value - needed for points - e.g for laser
+    //
+    (*pcstream)<<"fit"<<
+      "run="<<run<<       //run
+      "bz="<<bz<<         // magnetic filed used
+      "dtype="<<dtype<<   // detector match type
+      "ptype="<<ptype<<   // parameter type
+      "theta="<<theta<<   // theta
+      "phi="<<phi<<       // phi 
+      "snp="<<snp<<       // snp
+      "mean="<<mean<<     // mean dist value
+      "rms="<<rms<<       // rms
+      "sector="<<sector<<
+      "dsec="<<dsec<<
+      "refX="<<refXD<<         // referece X
+      "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
+      "pt="<<pt<<      //pt
+      "id="<<id<<             // track id
+      "entries="<<entries;// number of entries in bin
+    //
+    AliExternalTrackParam *trackOut0 = 0;
+    AliExternalTrackParam *trackOut1 = 0;
+    AliExternalTrackParam *ptrackIn0 = 0;
+    AliExternalTrackParam *ptrackIn1 = 0;
+
+    for (Int_t icorr=0; icorr<ncorr; icorr++) {
+      //
+      // special case of the TPC tracks crossing the CE
+      //
+      AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
+      corrections[icorr]=0;
+      if (entries>kMinEntries &&isOK){
+       AliExternalTrackParam trackIn0(refX,phi,tPar,cov);
+       AliExternalTrackParam trackIn1(refX,phi,tPar,cov);
+       ptrackIn1=&trackIn0;
+       ptrackIn0=&trackIn1;
+       //
+       if (debug)  trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,pcstream);
+       if (!debug) trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,0);
+       if (debug)  trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,pcstream);
+       if (!debug) trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,0);
+       //
+       if (trackOut0 && trackOut1){
+         //
+         if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,1,kTRUE,kMaxSnp))  isOK=kFALSE;
+         if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
+         // rotate all tracks to the same frame
+         if (!trackOut0->Rotate(trackIn0.GetAlpha())) isOK=kFALSE;
+         if (!trackIn1.Rotate(trackIn0.GetAlpha()))  isOK=kFALSE;
+         if (!trackOut1->Rotate(trackIn0.GetAlpha())) isOK=kFALSE;       
+         //
+         if (!AliTrackerBase::PropagateTrackTo(trackOut0,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
+         if (!AliTrackerBase::PropagateTrackTo(&trackIn1,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
+         if (!AliTrackerBase::PropagateTrackTo(trackOut1,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
+         //
+         corrections[icorr] = (trackOut0->GetParameter()[ptype]-trackIn0.GetParameter()[ptype]);
+         corrections[icorr]-= (trackOut1->GetParameter()[ptype]-trackIn1.GetParameter()[ptype]);
+         (*pcstream)<<"fitDebug"<< // just to debug the correction
+           "mean="<<mean<<
+           "pIn0.="<<ptrackIn0<<
+           "pIn1.="<<ptrackIn1<<
+           "pOut0.="<<trackOut0<<
+           "pOut1.="<<trackOut1<<
+           "refX="<<refXD<<
+           "\n";
+         delete trackOut0;      
+         delete trackOut1;      
+       }else{
+         corrections[icorr]=0;
+         isOK=kFALSE;
+       }
+      }      
+      (*pcstream)<<"fit"<<
+       Form("%s=",corr->GetName())<<corrections[icorr];   // dump correction value
+    }
+    //
+    (*pcstream)<<"fit"<<"isOK="<<isOK<<"\n";
+  }
+  delete pcstream;
+}
+
+
+
+void AliTPCCorrection::MakeLaserDistortionTreeOld(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;
+  const Double_t kMaxDist=1.;  // max distance - space correction
+  const Double_t kMaxRMS=0.05;  // max distance -between point and local mean
   TVectorD *vecdY=0;
   TVectorD *vecdZ=0;
   TVectorD *veceY=0;
@@ -1773,92 +2004,163 @@ void AliTPCCorrection::MakeLaserDistortionTree(TTree* tree, TObjArray *corrArray
     }
     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=5;
-      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(AliTrackerBase::GetBz());
-      //      const AliTPCRecoParam * recoParam = calib->GetTransform()->GetCurrentRecoParam();
-      //Double_t xyz0[3]={gx,gy,gz};
-      Double_t oldR=TMath::Sqrt(gx*gx+gy*gy);
-      Double_t fphi = TMath::ATan2(gy,gx);      
-      Double_t fsector = 9.*fphi/TMath::Pi();
-      if (fsector<0) fsector+=18;
-      Double_t dsec = fsector-Int_t(fsector)-0.5;
-      Double_t refX=0;
-      //
-      if (1 && oldR>1) {
-       Float_t xyz1[3]={gx,gy,gz};
-       Int_t sector=(gz>0)?0:18;
-       correction->CorrectPoint(xyz1, sector);
-       refX=TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
-       dRrec=oldR-refX;
-      } 
-
-      (*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
-       "sector="<<fsector<<
-       "dsec="<<dsec<<
+    TLinearFitter  fitter(2,"pol1");
+    for (Int_t iter=0; iter<2; iter++){
+      Double_t kfit0=0, kfit1=0;
+      Int_t npoints=fitter.GetNpoints();
+      if (npoints>80){
+       fitter.Eval();
+       kfit0=fitter.GetParameter(0);
+       kfit1=fitter.GetParameter(1);
+      }
+      for (Int_t irow=0; irow<159; irow++){
+       Bool_t isOK=kTRUE;
+       Int_t isOKF=0;
+       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=5;
+       Double_t array[10];
+       Int_t first3=TMath::Max(irow-3,0);
+       Int_t last3 =TMath::Min(irow+3,159);
+       Int_t counter=0;
+       if ((*ltr->GetVecSec())[irow]>=0 && err) {
+         for (Int_t jrow=first3; jrow<=last3; jrow++){
+           if ((*ltr->GetVecSec())[irow]!= (*ltr->GetVecSec())[jrow]) continue;
+           if ((*err)[jrow]<kEpsilon) continue;
+           array[counter]=(*delta)[jrow];
+           counter++;
+         }
+       }    
+       Double_t rms3  = 0;
+       Double_t mean3 = 0;
+       if (counter>2){
+         rms3  = TMath::RMS(counter,array);
+         mean3  = TMath::Mean(counter,array);
+       }else{
+         isOK=kFALSE;
+       }
+       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 dRrec=0;
+       //      Double_t rms = err->GetMatrixArray()[irow];
        //
-       "refX="<<refX<<      // reference radius
-       "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);
+       gx = (*ltr->GetVecGX())[irow];
+       gy = (*ltr->GetVecGY())[irow];
+       gz = (*ltr->GetVecGZ())[irow];
+       //
+       // get delta R used in reconstruction
+       AliTPCcalibDB*  calib=AliTPCcalibDB::Instance();  
+       AliTPCCorrection * correction = calib->GetTPCComposedCorrection(AliTrackerBase::GetBz());
+       //      const AliTPCRecoParam * recoParam = calib->GetTransform()->GetCurrentRecoParam();
+       //Double_t xyz0[3]={gx,gy,gz};
+       Double_t oldR=TMath::Sqrt(gx*gx+gy*gy);
+       Double_t fphi = TMath::ATan2(gy,gx);      
+       Double_t fsector = 9.*fphi/TMath::Pi();
+       if (fsector<0) fsector+=18;
+       Double_t dsec = fsector-Int_t(fsector)-0.5;
+       Double_t refX=0;
+       Int_t id= ltr->GetId();
+       Double_t pt=0;
+       //
+       if (1 && oldR>1) {
+         Float_t xyz1[3]={gx,gy,gz};
+         Int_t sector=(gz>0)?0:18;
+         correction->CorrectPoint(xyz1, sector);
+         refX=TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
+         dRrec=oldR-refX;
+       } 
+       if (TMath::Abs(rms3)>kMaxRMS) isOK=kFALSE;
+       if (TMath::Abs(mean-mean3)>kMaxRMS) isOK=kFALSE;
+       if (counter<4) isOK=kFALSE;     
+       if (npoints<90) isOK=kFALSE;    
+       if (isOK){
+         fitter.AddPoint(&refX,mean);
        }
-       // Double_t value=distPoint[2]-gz;
-       if (itype==0 && r0>1){
-         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;
+       Double_t deltaF=kfit0+kfit1*refX;
+       if (iter==1){
+         (*pcstream)<<"fitFull"<<  // dumpe also intermediate results
+           "bz="<<bz<<         // magnetic filed used
+           "dtype="<<dtype<<   // detector match type
+           "ptype="<<itype<<   // parameter type
+           "theta="<<theta<<   // theta
+           "phi="<<phi<<       // phi 
+           "snp="<<snp<<       // snp
+           "mean="<<mean3<<     // mean dist value
+           "rms="<<rms3<<       // rms
+           "deltaF="<<deltaF<<
+           "npoints="<<npoints<<  //number of points
+           "mean3="<<mean3<<     // mean dist value
+           "rms3="<<rms3<<       // rms
+           "counter="<<counter<<
+           "sector="<<fsector<<
+           "dsec="<<dsec<<
+           //
+           "refX="<<refX<<      // reference radius
+           "gx="<<gx<<         // global position
+           "gy="<<gy<<         // global position
+           "gz="<<gz<<         // global position
+           "dRrec="<<dRrec<<      // delta Radius in reconstruction
+           "id="<<id<<     //bundle    
+           "entries="<<nentries<<// number of entries in bin
+           "\n";
+       }
+       if (iter==1) (*pcstream)<<"fit"<<  // dump valus for 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="<<mean3<<     // mean dist value
+         "rms="<<rms3<<       // rms
+         "sector="<<fsector<<
+         "dsec="<<dsec<<
+         //
+         "refX="<<refX<<      // reference radius
+         "gx="<<gx<<         // global position
+         "gy="<<gy<<         // global position
+         "gz="<<gz<<         // global position
+         "dRrec="<<dRrec<<      // delta Radius in reconstruction
+         "pt="<<pt<<           //pt
+         "id="<<id<<     //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};
+       if (iter==1){
+         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 && r0>1){
+             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;
+           }
+           if (TMath::Abs(distortions[icorr])>kMaxDist) {isOKF=icorr+1; isOK=kFALSE; }
+           if (TMath::Abs(distortionsR[icorr])>kMaxDist) {isOKF=icorr+1; isOK=kFALSE;}
+           (*pcstream)<<"fit"<<
+             Form("%s=",corr->GetName())<<distortions[icorr];    // dump correction value
+         }
+         (*pcstream)<<"fit"<<"isOK="<<isOK<<"\n";
        }
-       (*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;
@@ -2457,3 +2759,184 @@ Double_t AliTPCCorrection::GetCorrXYZ(Double_t gx, Double_t gy, Double_t gz, Int
   if (axisType==2) return distPoint[2]-gz;
   return phi1-phi0;
 }
+
+
+
+
+
+void AliTPCCorrection::MakeLaserDistortionTree(TTree* tree, TObjArray *corrArray, Int_t itype){
+  //
+  // Make a laser fit tree for global minimization
+  //  
+  AliTPCcalibDB*  calib=AliTPCcalibDB::Instance();  
+  AliTPCCorrection * correction = calib->GetTPCComposedCorrection();  
+  if (!correction) correction = calib->GetTPCComposedCorrection(AliTrackerBase::GetBz());  
+  correction->AddVisualCorrection(correction,0);  //register correction
+
+  AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
+  AliTPCParam     *param     = AliTPCcalibDB::Instance()->GetParameters();
+  //
+  const Double_t cutErrY=0.05;
+  const Double_t kSigmaCut=4;
+  //  const Double_t cutErrZ=0.03;
+  const Double_t kEpsilon=0.00000001;
+  const Double_t kMaxDist=1.;  // max distance - space correction
+  TVectorD *vecdY=0;
+  TVectorD *vecdZ=0;
+  TVectorD *veceY=0;
+  TVectorD *veceZ=0;
+  AliTPCLaserTrack *ltr=0;
+  AliTPCLaserTrack::LoadTracks();
+  tree->SetBranchAddress("dY.",&vecdY);
+  tree->SetBranchAddress("dZ.",&vecdZ);
+  tree->SetBranchAddress("eY.",&veceY);
+  tree->SetBranchAddress("eZ.",&veceZ);
+  tree->SetBranchAddress("LTr.",&ltr);
+  Int_t entries= tree->GetEntries();
+  TTreeSRedirector *pcstream= new TTreeSRedirector("distortionLaser_0.root");
+  Double_t bz=AliTrackerBase::GetBz();
+  // 
+  Double_t globalXYZ[3];
+  Double_t globalXYZCorr[3];
+  for (Int_t ientry=0; ientry<entries; ientry++){
+    tree->GetEntry(ientry);
+    if (!ltr->GetVecGX()){
+      ltr->UpdatePoints();
+    }
+    //
+    TVectorD fit10(5);
+    TVectorD fit5(5);
+    printf("Entry\t%d\n",ientry);
+    for (Int_t irow0=0; irow0<158; irow0+=1){
+      //       
+      TLinearFitter fitter10(4,"hyp3");
+      TLinearFitter fitter5(2,"hyp1");
+      Int_t sector= (Int_t)(*ltr->GetVecSec())[irow0];
+      if (sector<0) continue;
+      //if (TMath::Abs(vecdY->GetMatrixArray()[irow0])<kEpsilon) continue;
+
+      Double_t refX= (*ltr->GetVecLX())[irow0];
+      Int_t firstRow1 = TMath::Max(irow0-10,0);
+      Int_t lastRow1  = TMath::Min(irow0+10,158);
+      Double_t padWidth=(irow0<64)?0.4:0.6;
+      // make long range fit
+      for (Int_t irow1=firstRow1; irow1<=lastRow1; irow1++){
+       if (TMath::Abs((*ltr->GetVecSec())[irow1]-sector)>kEpsilon) continue;
+       if (veceY->GetMatrixArray()[irow1]>cutErrY) continue;
+       if (TMath::Abs(vecdY->GetMatrixArray()[irow1])<kEpsilon) continue;
+       Double_t idealX= (*ltr->GetVecLX())[irow1];
+       Double_t idealY= (*ltr->GetVecLY())[irow1];
+       Double_t idealZ= (*ltr->GetVecLZ())[irow1];
+       Double_t gx= (*ltr->GetVecGX())[irow1];
+       Double_t gy= (*ltr->GetVecGY())[irow1];
+       Double_t gz= (*ltr->GetVecGZ())[irow1];
+       Double_t measY=(*vecdY)[irow1]+idealY;
+       Double_t deltaR = GetCorrXYZ(gx, gy, gz, 0,0);
+       // deltaR = R distorted -R ideal
+       Double_t xxx[4]={idealX+deltaR-refX,TMath::Cos(idealY/padWidth), TMath::Sin(idealY/padWidth)};
+       fitter10.AddPoint(xxx,measY,1);
+      }
+      Bool_t isOK=kTRUE;
+      Double_t rms10=0;//TMath::Sqrt(fitter10.GetChisquare()/(fitter10.GetNpoints()-4));
+      Double_t mean10  =0;//   fitter10.GetParameter(0);
+      Double_t slope10  =0;//   fitter10.GetParameter(0);
+      Double_t cosPart10  = 0;//  fitter10.GetParameter(2);
+      Double_t sinPart10   =0;//  fitter10.GetParameter(3); 
+
+      if (fitter10.GetNpoints()>10){
+       fitter10.Eval();
+       rms10=TMath::Sqrt(fitter10.GetChisquare()/(fitter10.GetNpoints()-4));
+       mean10      =   fitter10.GetParameter(0);
+       slope10     =   fitter10.GetParameter(1);
+       cosPart10   =   fitter10.GetParameter(2);
+       sinPart10   =  fitter10.GetParameter(3); 
+       //
+       // make short range fit
+       //
+       for (Int_t irow1=firstRow1+5; irow1<=lastRow1-5; irow1++){
+         if (TMath::Abs((*ltr->GetVecSec())[irow1]-sector)>kEpsilon) continue;
+         if (veceY->GetMatrixArray()[irow1]>cutErrY) continue;
+         if (TMath::Abs(vecdY->GetMatrixArray()[irow1])<kEpsilon) continue;
+         Double_t idealX= (*ltr->GetVecLX())[irow1];
+         Double_t idealY= (*ltr->GetVecLY())[irow1];
+         Double_t idealZ= (*ltr->GetVecLZ())[irow1];
+         Double_t gx= (*ltr->GetVecGX())[irow1];
+         Double_t gy= (*ltr->GetVecGY())[irow1];
+         Double_t gz= (*ltr->GetVecGZ())[irow1];
+         Double_t measY=(*vecdY)[irow1]+idealY;
+         Double_t deltaR = GetCorrXYZ(gx, gy, gz, 0,0);
+         // deltaR = R distorted -R ideal 
+         Double_t expY= mean10+slope10*(idealX+deltaR-refX);
+         if (TMath::Abs(measY-expY)>kSigmaCut*rms10) continue;
+         //
+         Double_t corr=cosPart10*TMath::Cos(idealY/padWidth)+sinPart10*TMath::Sin(idealY/padWidth);
+         Double_t xxx[4]={idealX+deltaR-refX,TMath::Cos(idealY/padWidth), TMath::Sin(idealY/padWidth)};
+         fitter5.AddPoint(xxx,measY-corr,1);
+       }     
+      }else{
+       isOK=kFALSE;
+      }
+      if (fitter5.GetNpoints()<8) isOK=kFALSE;
+
+      Double_t rms5=0;//TMath::Sqrt(fitter5.GetChisquare()/(fitter5.GetNpoints()-4));
+      Double_t offset5  =0;//  fitter5.GetParameter(0);
+      Double_t slope5   =0;//  fitter5.GetParameter(0); 
+      if (isOK){
+       fitter5.Eval();
+       rms5=TMath::Sqrt(fitter5.GetChisquare()/(fitter5.GetNpoints()-4));
+       offset5  =  fitter5.GetParameter(0);
+       slope5   =  fitter5.GetParameter(0); 
+      }
+      //
+      Double_t dtype=5;
+      Double_t ptype=0;
+      Double_t phi   =(*ltr->GetVecPhi())[irow0];
+      Double_t theta =ltr->GetTgl();
+      Double_t mean=(vecdY)->GetMatrixArray()[irow0];
+      Double_t gx=0,gy=0,gz=0;
+      Double_t snp = (*ltr->GetVecP2())[irow0];
+      Int_t bundle= ltr->GetBundle();
+      Int_t id= ltr->GetId();
+      //      Double_t rms = err->GetMatrixArray()[irow];
+      //
+      gx = (*ltr->GetVecGX())[irow0];
+      gy = (*ltr->GetVecGY())[irow0];
+      gz = (*ltr->GetVecGZ())[irow0];
+      Double_t dRrec = GetCorrXYZ(gx, gy, gz, 0,0);
+      fitter10.GetParameters(fit10);
+      fitter5.GetParameters(fit5);      
+      Double_t idealY= (*ltr->GetVecLY())[irow0];
+      Double_t measY=(*vecdY)[irow0]+idealY;
+      Double_t corr=cosPart10*TMath::Cos(idealY/padWidth)+sinPart10*TMath::Sin(idealY/padWidth);
+      if (TMath::Max(rms5,rms10)>0.06) isOK=kFALSE;
+      //
+      (*pcstream)<<"fitFull"<<  // dumpe also intermediate results
+       "bz="<<bz<<         // magnetic filed used
+       "dtype="<<dtype<<   // detector match type
+       "ptype="<<ptype<<   // parameter type
+       "theta="<<theta<<   // theta
+       "phi="<<phi<<       // phi 
+       "snp="<<snp<<       // snp
+       "sector="<<sector<<
+       "bundle="<<bundle<<
+//     //      "dsec="<<dsec<<
+       "refX="<<refX<<      // reference radius
+       "gx="<<gx<<         // global position
+       "gy="<<gy<<         // global position
+       "gz="<<gz<<         // global position
+       "dRrec="<<dRrec<<      // delta Radius in reconstruction
+       "id="<<id<<     //bundle
+       "rms10="<<rms10<<
+       "rms5="<<rms5<<
+       "fit10.="<<&fit10<<
+       "fit5.="<<&fit5<<
+       "measY="<<measY<<
+       "mean="<<mean<<
+       "idealY="<<idealY<<
+       "corr="<<corr<<
+       "isOK="<<isOK<<
+       "\n";
+    }
+  }
+  delete pcstream;
+}
index c8002064bf213738d8009eecb0777288a69bd897..679da2a8007a043a7e568302b80e7e00ea9d6eb4 100644 (file)
@@ -33,11 +33,13 @@ public:
 
   // functions to correct a space point
           void CorrectPoint (      Float_t x[],const Short_t roc);
+          void CorrectPointLocal(Float_t x[],const Short_t roc);
           void CorrectPoint (const Float_t x[],const Short_t roc,Float_t xp[]);
   virtual void GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]);
 
   // functions to distort a space point
           void DistortPoint (      Float_t x[],const Short_t roc);
+          void DistortPointLocal(Float_t x[],const Short_t roc);
           void DistortPoint (const Float_t x[],const Short_t roc,Float_t xp[]);
   virtual void GetDistortion(const Float_t x[],const Short_t roc,Float_t dx[]);
 
@@ -65,7 +67,9 @@ public:
   virtual void SetOmegaTauT1T2(Float_t omegaTau,Float_t t1,Float_t t2);
   AliExternalTrackParam * FitDistortedTrack(AliExternalTrackParam & trackIn, Double_t refX, Int_t dir,TTreeSRedirector *pcstream);
   void StoreInOCDB(Int_t startRun, Int_t endRun, const char *comment=0);
-  static void MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step=1, Bool_t debug=0);
+  static void MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step=1, Int_t offset=0, Bool_t debug=0);
+  static void MakeSectorDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step=1, Int_t offset=0, Bool_t debug=0);
+  static void MakeLaserDistortionTreeOld(TTree* tree, TObjArray *corrArray, Int_t itype);
   static void MakeLaserDistortionTree(TTree* tree, TObjArray *corrArray, Int_t itype);
 
   void FastSimDistortedVertex(Double_t orgVertex[3], Int_t nTracks, AliESDVertex &aV, AliESDVertex &avOrg, AliESDVertex &cV, AliESDVertex &cvOrg, TTreeSRedirector * const pcstream, Double_t etaCuts);
index 51b78b742c4bb199ee41e3d1f54bdcf76c38855f..88c67904dfff9c70fd02ef92a850f7f13a14800e 100644 (file)
@@ -154,3 +154,28 @@ void AliTPCExBBShape::Print(Option_t* option) const {
     printf(" - C1: %1.4f, C2: %1.4f \n",fC1,fC2);
   }    
 }
+
+Double_t AliTPCExBBShape::GetBFieldXYZ(Double_t gx, Double_t gy, Double_t gz, Int_t axisType){
+  //
+  // return B field at given x,y,z
+  // 
+  AliMagF* field = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
+  if (!field) return 0;
+  Double_t xyz[3]={gx,gy,gz};
+  Double_t bxyz[3]={0};
+  field->Field(xyz,bxyz);
+  //
+  Double_t r=TMath::Sqrt(gx*gx+gy*gy);
+  Double_t b=TMath::Sqrt(bxyz[0]*bxyz[0]+bxyz[1]*bxyz[1]);
+  if (axisType==0) {
+    return (xyz[0]*bxyz[1]-xyz[1]*bxyz[0])/(bxyz[2]*r);
+  }
+  if (axisType==1){
+    return (xyz[0]*bxyz[0]+xyz[1]*bxyz[1])/(bxyz[2]*r);
+  }
+  if (axisType==2) return bxyz[2];
+  if (axisType==3) return bxyz[0];
+  if (axisType==4) return bxyz[1];
+  if (axisType==5) return bxyz[2];
+  return bxyz[2];
+}
index defcce7fbb8aabc2a24cada01b3e07b3c82df1be..f2763f8a0376105f51ab6229a963689c5dceb264 100644 (file)
@@ -65,6 +65,7 @@ public:
   void GetBxAndByOverBz(const Float_t x[],const Short_t roc,Float_t BxByOverBz[]);
 
   virtual void Print(Option_t* option="") const;
+  static Double_t GetBFieldXYZ(Double_t gx, Double_t gy, Double_t gz, Int_t axisType);
 
 private:
   Float_t fC1; // coefficient C1          (compare Jim Thomas's notes for definitions)
index c3954a492a3693c9deeb5b8b278cfc798128b07e..fcd1b9cc59f5c3813b28b19e1c9db6236ac94e1c 100644 (file)
@@ -190,7 +190,7 @@ void AliTPCExBEffectiveSector::Print(const Option_t* option) const {
   }    
 }
 
-void  AliTPCExBEffectiveSector::MakeResidualMap(THnSparse * hisInput, const char *sname){
+void  AliTPCExBEffectiveSector::MakeResidualMap(THnSparse * hisInput, const char *sname, Int_t ptype, Int_t dtype){
   //
   // Make cluster residual map from the n-dimensional histogram
   // hisInput supposed to have given format:
@@ -201,7 +201,7 @@ void  AliTPCExBEffectiveSector::MakeResidualMap(THnSparse * hisInput, const char
   Int_t nbins1=hisInput->GetAxis(1)->GetNbins();
   Int_t nbins2=hisInput->GetAxis(2)->GetNbins();
   Int_t nbins3=hisInput->GetAxis(3)->GetNbins();
-
+  TF1 *fgaus=0;
   TH3F * hisResMap3D = 
     new TH3F("his3D","his3D",
             nbins1,hisInput->GetAxis(1)->GetXmin(), hisInput->GetAxis(1)->GetXmax(),
@@ -314,13 +314,36 @@ void  AliTPCExBEffectiveSector::MakeResidualMap(THnSparse * hisInput, const char
        hisResMap3D->SetBinContent(ibin1,ibin2,ibin3, mean);
        Double_t phi=TMath::Pi()*sector/9;
        if (phi>TMath::Pi()) phi+=TMath::Pi();
+       Double_t meanG=0;
+       Double_t rmsG=0;
+       if (entries>50){
+         if (!fgaus) {     
+           his->Fit("gaus","Q","goff");
+           fgaus= (TF1*)((his->GetListOfFunctions()->FindObject("gaus"))->Clone());
+           fgaus->SetName("gausk");
+         }
+         if (fgaus) his->Fit(fgaus,"Q","goff");
+         if (fgaus){
+           meanG=fgaus->GetParameter(1);
+           rmsG=fgaus->GetParameter(2);
+         }
+       }
+       Double_t dsec=sector-Int_t(sector)-0.5;
+       Double_t snp=dsec*TMath::Pi()/9.;
        (*pcstream)<<"delta"<<
+         "ptype="<<ptype<<
+         "dtype="<<dtype<<
          "sector="<<sector<<
+         "dsec="<<dsec<<
+         "snp="<<snp<<
          "phi="<<phi<<
          "localX="<<localX<<
          "kZ="<<kZ<<
+         "theta="<<kZ<<
          "mean="<<mean<<
          "rms="<<rms<<
+         "meanG="<<meanG<<
+         "rmsG="<<rmsG<<
          "entries="<<entries<<
          "meanA="<<meanA<<
          "rmsA="<<rmsA<<
index a019a92721dccc99b5b8f2d842c52cf687f79e93..8bf1d0787c74c0dedf19428be52c8df73effd53b 100644 (file)
@@ -31,7 +31,7 @@ public:
   Float_t GetC1() const {return fC1;}
   Float_t GetC0() const {return fC0;}
   void Print(const Option_t* option) const;
-  static void  MakeResidualMap(THnSparse * hisInput, const char *sname);
+  static void  MakeResidualMap(THnSparse * hisInput, const char *sname, Int_t ptype, Int_t dtype=0);
 public:
   virtual void GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]);
 public:
index 11093d43b9c3533ee7b296b0f2173844f35e9151..d2501e7b1ffbe18a60a021d2d10b8e5679025d82 100644 (file)
@@ -164,6 +164,7 @@ TGaxis *axis = new TGaxis(xmax,ymin,xmax,ymax,ymin,ymax,50510,"+L");
 #include "AliTracker.h"
 #include "AliMagF.h"
 #include "AliTPCCalROC.h"
+#include "AliESDv0.h"
 
 #include "AliLog.h"
 
@@ -234,12 +235,12 @@ AliTPCcalibTimeGain::AliTPCcalibTimeGain(const Text_t *name, const Text_t *title
   Double_t deltaTime = EndTime - StartTime;
   
   
-  // main histogram for time dependence: dE/dx, time, type (1-muon cosmic,2-pion beam data), meanDriftlength, momenta (only filled if enough space is available), run number
+  // main histogram for time dependence: dE/dx, time, type (1-muon cosmic,2-pion beam data, 3&4 - proton points at higher dE/dx), meanDriftlength, momenta (only filled if enough space is available), run number, eta
   Int_t timeBins = TMath::Nint(deltaTime/deltaIntegrationTimeGain);
-  Int_t binsGainTime[6]    = {150,  timeBins,    2,  25, 200, 10000000};
-  Double_t xminGainTime[6] = {0.5, StartTime,  0.5,   0, 0.1,    -0.5};
-  Double_t xmaxGainTime[6] = {  8,   EndTime,  2.5, 250,  50, 9999999.5};
-  fHistGainTime = new THnSparseF("HistGainTime","dEdx time dep.;dEdx,time,type,driftlength,momenta,run number;dEdx",6,binsGainTime,xminGainTime,xmaxGainTime);
+  Int_t binsGainTime[7]    = {150,  timeBins,    4,  25, 200, 10000000, 20};
+  Double_t xminGainTime[7] = {0.5, StartTime,  0.5,   0, 0.1,    -0.5,  -1};
+  Double_t xmaxGainTime[7] = {  8,   EndTime,  4.5, 250,  50, 9999999.5, 1};
+  fHistGainTime = new THnSparseF("HistGainTime","dEdx time dep.;dEdx,time,type,driftlength,momenta,run number, eta;dEdx",7,binsGainTime,xminGainTime,xmaxGainTime);
   BinLogX(fHistGainTime, 4);
   //
   fHistDeDxTotal = new TH2F("DeDx","dEdx; momentum p (GeV); TPC signal (a.u.)",250,0.01,100.,1000,0.,8);
@@ -373,7 +374,7 @@ void AliTPCcalibTimeGain::ProcessBeamEvent(AliESDEvent *event) {
   //
   // track loop
   //
-  for (Int_t i=0;i<ntracks;++i) {
+  for (Int_t i=0;i<ntracks;++i) { // begin track loop
 
     AliESDtrack *track = event->GetTrack(i);
     if (!track) continue;
@@ -392,8 +393,16 @@ void AliTPCcalibTimeGain::ProcessBeamEvent(AliESDEvent *event) {
 
     // exclude tracks which do not look like primaries or are simply too short or on wrong sectors
     if (nclsDeDx < 60) continue;     
-    if (TMath::Abs(trackIn->GetTgl()) > 1) continue;
-    if (TMath::Abs(trackIn->GetSnp()) > 0.6) continue;
+    //if (TMath::Abs(trackIn->GetTgl()) > 1) continue;
+    //if (TMath::Abs(trackIn->GetSnp()) > 0.6) continue;
+    if (TMath::Abs(trackIn->Eta()) > 1) continue;
+    UInt_t status = track->GetStatus();
+    if ((status&AliESDtrack::kTPCrefit)==0) continue;
+    if (track->GetNcls(0) < 3) continue; // ITS clusters
+    Float_t dca[2], cov[3];
+    track->GetImpactParameters(dca,cov);
+    if (dca[0] > 7 || dca[0] < 0.000001) continue; // cut in xy
+    Double_t eta = trackIn->Eta();
     
     // Get seeds
     TObject *calibObject;
@@ -402,19 +411,70 @@ void AliTPCcalibTimeGain::ProcessBeamEvent(AliESDEvent *event) {
       if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
     }    
 
-    if (seed) { 
-      if (fLowMemoryConsumption) {
-       if (meanP > 0.5 || meanP < 0.4) continue;
-       meanP = 0.45; // set all momenta to one in order to save memory
-      }
+    if (seed) {
+      Int_t particleCase = 0;
+      if (meanP > 0.5  || meanP < 0.4)  particleCase = 2; // MIP pions
+      if (meanP > 0.56 || meanP < 0.57) particleCase = 3; // protons 1
+      if (meanP > 0.64 || meanP < 0.66) particleCase = 4; // protons 2
+      //
+      if (fLowMemoryConsumption && particleCase == 0) continue;
+      //
       Double_t tpcSignal = GetTPCdEdx(seed);
       fHistDeDxTotal->Fill(meanP, tpcSignal);
       //
-      //dE/dx, time, type (1-muon cosmic,2-pion beam data), momenta
-      Double_t vec[6] = {tpcSignal,time,2,meanDrift,meanP,runNumber};
+      //dE/dx, time, type (1-muon cosmic,2-pion beam data, 3&4 protons), momenta, runNumner, eta
+      Double_t vec[7] = {tpcSignal,time,particleCase,meanDrift,meanP,runNumber, eta};
       fHistGainTime->Fill(vec);
     }
     
+  } // end track loop
+  //
+  // V0 loop -- in beam events the cosmic part of the histogram is filled with GammaConversions
+  //
+  for(Int_t iv0 = 0; iv0 < event->GetNumberOfV0s(); iv0++) {
+    AliESDv0 * v0 = event->GetV0(iv0);
+    if (!v0->GetOnFlyStatus()) continue;
+    if (v0->GetEffMass(0,0) > 0.02) continue; // select low inv. mass
+    Double_t xyz[3];
+    v0->GetXYZ(xyz[0], xyz[1], xyz[2]);
+    if (TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]) < 3 || TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]) > 30) continue;
+    //
+    // "loop over daughters" 
+    //
+    for(Int_t idaughter = 0; idaughter < 2; idaughter++) { // daughter loop
+      Int_t index = idaughter == 0 ? v0->GetPindex() : v0->GetNindex();
+      AliESDtrack * trackP = event->GetTrack(index);
+      AliESDfriendTrack *friendTrackP = esdFriend->GetTrack(index);
+      if (!friendTrackP) continue;
+      const AliExternalTrackParam * trackPIn = trackP->GetInnerParam();
+      const AliExternalTrackParam * trackPOut = friendTrackP->GetTPCOut();
+      if (!trackPIn) continue;
+      if (!trackPOut) continue;
+      // calculate necessary track parameters
+      Double_t meanP = trackPIn->GetP();
+      Double_t meanDrift = 250 - 0.5*TMath::Abs(trackPIn->GetZ() + trackPOut->GetZ());
+      Int_t nclsDeDx = trackP->GetTPCNcls();
+      // exclude tracks which do not look like primaries or are simply too short or on wrong sectors
+      if (nclsDeDx < 60) continue;     
+      if (TMath::Abs(trackPIn->GetTgl()) > 1) continue;
+      //
+      TObject *calibObject;
+      AliTPCseed *seed = 0;
+      for (Int_t l=0;(calibObject=friendTrackP->GetCalibObject(l));++l) {
+      if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
+      }    
+      if (seed) { 
+       if (fLowMemoryConsumption) {
+         if (meanP > 0.5 || meanP < 0.4) continue;
+         meanP = 0.45; // set all momenta to one in order to save memory
+      }
+       Double_t tpcSignal = GetTPCdEdx(seed);
+       //dE/dx, time, type (1-muon cosmic,2-pion beam data), momenta
+       Double_t vec[6] = {tpcSignal,time,1,meanDrift,meanP,runNumber};
+       fHistGainTime->Fill(vec);
+      }
+    }
+    
   }
 
 }
index 6687d07f1c5b7c081776022cd4b2111259e33127..700bfaef615b99eb545dabdeeeb724b513947ea3 100644 (file)
                Offline/HLT             Offline/HLT                    OCDB entries (AliTPCClusterParam) 
 */            
 
-/*
-
-How to retrive it from file (created using calibration task):
-
-gSystem->Load("libANALYSIS");
-gSystem->Load("libTPCcalib");
-TFile fcalib("CalibObjects.root");
-TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
-AliTPCcalibTracks * calibTracks = ( AliTPCcalibTracks *)array->FindObject("calibTracks");
-
-
-//USAGE of debug stream example
- gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
-  gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
-  AliXRDPROOFtoolkit tool;
-  TChain * chainres = tool.MakeChain("tracks.txt","ResolCl",0,10200);
-  chainres->Lookup();
-*/
-
 
                                                                //
 ///////////////////////////////////////////////////////////////////////////////
@@ -108,7 +89,6 @@ using namespace std;
 #include "AliTPCcalibTracks.h"
 #include "AliTPCClusterParam.h"
 #include "AliTPCcalibTracksCuts.h"
-#include "AliTPCCalPadRegion.h"
 #include "AliTPCCalPad.h"
 #include "AliTPCCalROC.h"
 #include "TText.h"
@@ -117,6 +97,7 @@ using namespace std;
 #include "TStatToolkit.h"
 #include "TCut.h"
 #include "THnSparse.h"
+#include "AliRieman.h"
 
 
 
@@ -133,25 +114,16 @@ AliTPCcalibTracks::AliTPCcalibTracks():
   fHisRMSZ(0),      // THnSparse - rms Z 
   fHisQmax(0),      // THnSparse - qmax 
   fHisQtot(0),      // THnSparse - qtot 
-  fArrayAmpRow(0),
-  fArrayAmp(0), 
   fArrayQDY(0), 
   fArrayQDZ(0), 
   fArrayQRMSY(0),
   fArrayQRMSZ(0),
-  fArrayChargeVsDriftlength(0),
-  fcalPadRegionChargeVsDriftlength(0),
-  fDeltaY(0),
-  fDeltaZ(0),
   fResolY(0),
   fResolZ(0),
   fRMSY(0),
   fRMSZ(0),
   fCuts(0),
-  fHclus(0),
   fRejectedTracksHisto(0),
-  fHclusterPerPadrow(0),
-  fHclusterPerPadrowRaw(0),
   fClusterCutHisto(0),
   fCalPadClusterPerPad(0),
   fCalPadClusterPerPadRaw(0)
@@ -174,25 +146,16 @@ AliTPCcalibTracks::AliTPCcalibTracks(const AliTPCcalibTracks& calibTracks):
   fHisRMSZ(0),      // THnSparse - rms Z 
   fHisQmax(0),      // THnSparse - qmax 
   fHisQtot(0),      // THnSparse - qtot 
-  fArrayAmpRow(0),
-  fArrayAmp(0), 
   fArrayQDY(0), 
   fArrayQDZ(0), 
   fArrayQRMSY(0),
   fArrayQRMSZ(0),
-  fArrayChargeVsDriftlength(0),
-  fcalPadRegionChargeVsDriftlength(0),
-  fDeltaY(0),
-  fDeltaZ(0),
   fResolY(0),
   fResolZ(0),
   fRMSY(0),
   fRMSZ(0),
   fCuts(0),
-  fHclus(0),
   fRejectedTracksHisto(0),
-  fHclusterPerPadrow(0),
-  fHclusterPerPadrowRaw(0),
   fClusterCutHisto(0),
   fCalPadClusterPerPad(0),
   fCalPadClusterPerPadRaw(0)
@@ -206,14 +169,6 @@ AliTPCcalibTracks::AliTPCcalibTracks(const AliTPCcalibTracks& calibTracks):
    TH1::AddDirectory(kFALSE);
    
    Int_t length = -1;
-   // backward compatibility: if the data member doesn't yet exist, it will not be merged
-   (calibTracks.fArrayAmpRow) ? length = calibTracks.fArrayAmpRow->GetEntriesFast() : length = -1;
-   fArrayAmpRow = new TObjArray(length);
-   fArrayAmp = new TObjArray(length);
-   for (Int_t i = 0; i < length; i++) {
-      fArrayAmpRow->AddAt( (TProfile*)calibTracks.fArrayAmpRow->At(i)->Clone(), i);
-      fArrayAmp->AddAt( ((TProfile*)calibTracks.fArrayAmp->At(i)->Clone()), i);
-   }
    
    (calibTracks.fArrayQDY) ? length = calibTracks.fArrayQDY->GetEntriesFast() : length = -1;
    fArrayQDY= new TObjArray(length);
@@ -239,20 +194,9 @@ AliTPCcalibTracks::AliTPCcalibTracks(const AliTPCcalibTracks& calibTracks):
       fRMSZ->AddAt( ((TH1F*)calibTracks.fRMSZ->At(i)->Clone()), i);
    } 
    
-   (calibTracks.fArrayChargeVsDriftlength) ? length = calibTracks.fArrayChargeVsDriftlength->GetEntriesFast() : length = -1;
-   (calibTracks.fArrayChargeVsDriftlength) ? fArrayChargeVsDriftlength = new TObjArray(length) : fArrayChargeVsDriftlength = 0;
-   for (Int_t i = 0; i < length; i++) {
-      fArrayChargeVsDriftlength->AddAt( ((TProfile*)calibTracks.fArrayChargeVsDriftlength->At(i)->Clone()), i);
-   }
    
-   fDeltaY =  (TH1F*)calibTracks.fDeltaY->Clone();
-   fDeltaZ =  (TH1F*)calibTracks.fDeltaZ->Clone();
-   fHclus = (TH1I*)calibTracks.fHclus->Clone();
    fClusterCutHisto = (TH2I*)calibTracks.fClusterCutHisto->Clone();
    fRejectedTracksHisto    = (TH1I*)calibTracks.fRejectedTracksHisto->Clone();
-   fHclusterPerPadrow      = (TH1I*)calibTracks.fHclusterPerPadrow->Clone();
-   fHclusterPerPadrowRaw   = (TH1I*)calibTracks.fHclusterPerPadrowRaw->Clone();
-   fcalPadRegionChargeVsDriftlength = (AliTPCCalPadRegion*)calibTracks.fcalPadRegionChargeVsDriftlength->Clone();
    fCalPadClusterPerPad    = (AliTPCCalPad*)calibTracks.fCalPadClusterPerPad->Clone();
    fCalPadClusterPerPadRaw = (AliTPCCalPad*)calibTracks.fCalPadClusterPerPadRaw->Clone();
 
@@ -286,25 +230,16 @@ AliTPCcalibTracks::AliTPCcalibTracks(const Text_t *name, const Text_t *title, Al
   fHisRMSZ(0),      // THnSparse - rms Z 
   fHisQmax(0),      // THnSparse - qmax 
   fHisQtot(0),      // THnSparse - qtot 
-  fArrayAmpRow(0),
-  fArrayAmp(0), 
   fArrayQDY(0), 
   fArrayQDZ(0), 
   fArrayQRMSY(0),
   fArrayQRMSZ(0),
-  fArrayChargeVsDriftlength(0),
-  fcalPadRegionChargeVsDriftlength(0),
-  fDeltaY(0),
-  fDeltaZ(0),
   fResolY(0),
   fResolZ(0),
   fRMSY(0),
   fRMSZ(0),
   fCuts(0),
-  fHclus(0),
   fRejectedTracksHisto(0),
-  fHclusterPerPadrow(0),
-  fHclusterPerPadrowRaw(0),
   fClusterCutHisto(0),
   fCalPadClusterPerPad(0),
   fCalPadClusterPerPadRaw(0)
@@ -339,61 +274,14 @@ AliTPCcalibTracks::AliTPCcalibTracks(const Text_t *name, const Text_t *title, Al
    
    TH1::AddDirectory(kFALSE);
    
-   char chname[1000];
-   TProfile * prof1=0;
-   TH1F     * his1 =0;
-   fHclus = new TH1I("hclus","Number of clusters per track",160, 0, 160);     // valgrind 3
    fRejectedTracksHisto    = new TH1I("RejectedTracksHisto", "Rejected tracks, sorted by failed cut", 100, -1, 10);
-   fHclusterPerPadrow      = new TH1I("fHclusterPerPadrow", " clusters per padRow, used for the resolution tree", 160, 0, 160);
-   fHclusterPerPadrowRaw   = new TH1I("fHclusterPerPadrowRaw", " clusters per padRow, before cutting clusters", 160, 0, 160);
    fCalPadClusterPerPad    = new AliTPCCalPad("fCalPadClusterPerPad", "clusters per pad");
    fCalPadClusterPerPadRaw = new AliTPCCalPad("fCalPadClusterPerPadRaw", "clusters per pad, before cutting clusters");
    fClusterCutHisto = new TH2I("fClusterCutHisto", "Cutted cluster over padRow; Cut Criterium; PadRow", 5,1,5, 160,0,159);
    
-   // Amplitude  - sector - row histograms 
-   fArrayAmpRow = new TObjArray(72);
-   fArrayAmp    = new TObjArray(72);
-   fArrayChargeVsDriftlength = new TObjArray(72);
-   
-   for (Int_t i = 0; i < 36; i++){   
-     snprintf(chname,100,"Amp_row_Sector%d",i);
-      prof1 = new TProfile(chname,chname,63,0,64);          // valgrind 3   193,536 bytes in 354 blocks are still reachable 
-      prof1->SetXTitle("Pad row");
-      prof1->SetYTitle("Mean Max amplitude");
-      fArrayAmpRow->AddAt(prof1,i);
-      snprintf(chname,100,"Amp_row_Sector%d",i+36);
-      prof1 = new TProfile(chname,chname,96,0,97);       // valgrind 3   3,912 bytes in 6 blocks are possibly lost
-      prof1->SetXTitle("Pad row");
-      prof1->SetYTitle("Mean Max  amplitude");
-      fArrayAmpRow->AddAt(prof1,i+36);
-      
-      // amplitude
-      sprintf(chname,"Amp_Sector%d",i);
-      his1 = new TH1F(chname,chname,100,0,32);         // valgrind 
-      his1->SetXTitle("Max Amplitude (ADC)");
-      fArrayAmp->AddAt(his1,i);
-      sprintf(chname,"Amp_Sector%d",i+36);
-      his1 = new TH1F(chname,chname,100,0,32);         // valgrind 3   13,408,208 bytes in 229 blocks are still reachable
-      his1->SetXTitle("Max Amplitude (ADC)");
-      fArrayAmp->AddAt(his1,i+36);
-      
-      // driftlength
-      snprintf(chname,100, "driftlengt vs. charge, ROC %i", i);
-      prof1 = new TProfile(chname, chname, 25, 0, 250);
-      prof1->SetYTitle("Charge");
-      prof1->SetXTitle("Driftlength");
-      fArrayChargeVsDriftlength->AddAt(prof1,i);
-      snprintf(chname,100, "driftlengt vs. charge, ROC %i", i+36);
-      prof1 = new TProfile(chname, chname, 25, 0, 250);
-      prof1->SetYTitle("Charge");
-      prof1->SetXTitle("Driftlength");
-      fArrayChargeVsDriftlength->AddAt(prof1,i+36);
-   }
    
    TH1::AddDirectory(kFALSE);
    
-   fDeltaY = new TH1F("DeltaY","DeltaY",100,-1,1);
-   fDeltaZ = new TH1F("DeltaZ","DeltaZ",100,-1,1);
    
    fResolY = new TObjArray(3);
    fResolZ = new TObjArray(3);
@@ -456,19 +344,6 @@ AliTPCcalibTracks::AliTPCcalibTracks(const Text_t *name, const Text_t *title, Al
       }
    }
    
-   fcalPadRegionChargeVsDriftlength = new AliTPCCalPadRegion("fcalPadRegionChargeVsDriftlength", "TProfiles with charge vs driftlength for each pad region");
-   TProfile *tempProf;
-   for (UInt_t padSize = 0; padSize < 3; padSize++) {
-      for (UInt_t isector = 0; isector < 36; isector++) {
-         if (padSize == 0) sprintf(chname, "driftlengt vs. charge, sector %i, short pads", isector);
-         if (padSize == 1) sprintf(chname, "driftlengt vs. charge, sector %i, medium  pads", isector);
-         if (padSize == 2) sprintf(chname, "driftlengt vs. charge, sector %i, long  pads", isector);
-         tempProf = new TProfile(chname, chname, 500, 0, 250);
-         tempProf->SetYTitle("Charge");
-         tempProf->SetXTitle("Driftlength");
-         fcalPadRegionChargeVsDriftlength->SetObject(tempProf, isector, padSize);
-      }
-   }
    
 
    if (GetDebugLevel() > 1) cout << "AliTPCcalibTracks object sucessfully constructed: " << GetName() << endl; 
@@ -483,16 +358,7 @@ AliTPCcalibTracks::~AliTPCcalibTracks() {
    
   if (GetDebugLevel() > 0) cout << "AliTPCcalibTracks' destuctor called." << endl;
    Int_t length = 0;
-   if (fArrayAmpRow) length = fArrayAmpRow->GetEntriesFast();
-   for (Int_t i = 0; i < length; i++){
-      delete fArrayAmpRow->At(i);  
-      delete fArrayAmp->At(i);  
-   }
-   delete fArrayAmpRow;
-   delete fArrayAmp;
    
-   delete fDeltaY;
-   delete fDeltaZ;
    
    if (fResolY) {
      length = fResolY->GetEntriesFast();
@@ -518,31 +384,17 @@ AliTPCcalibTracks::~AliTPCcalibTracks() {
      }
    }
    
-   if (fArrayChargeVsDriftlength) {
-     length = fArrayChargeVsDriftlength->GetEntriesFast();
-     for (Int_t i = 0; i < length; i++){
-       delete fArrayChargeVsDriftlength->At(i);
-     }
-   }
    
     
    delete fArrayQDY;
    delete fArrayQDZ;
    delete fArrayQRMSY;
    delete fArrayQRMSZ;
-   delete fArrayChargeVsDriftlength;
    
-  delete fHclus;
   delete fRejectedTracksHisto;
   delete fClusterCutHisto;
-  delete fHclusterPerPadrow;
-  delete fHclusterPerPadrowRaw;
   if (fCalPadClusterPerPad)    delete fCalPadClusterPerPad;
   if (fCalPadClusterPerPadRaw) delete fCalPadClusterPerPadRaw;
-  if(fcalPadRegionChargeVsDriftlength) {
-     fcalPadRegionChargeVsDriftlength->Delete();
-     delete fcalPadRegionChargeVsDriftlength;
-  }
   delete fHisDeltaY;    // THnSparse - delta Y 
   delete fHisDeltaZ;    // THnSparse - delta Z 
   delete fHisRMSY;      // THnSparse - rms Y 
@@ -664,666 +516,272 @@ void  AliTPCcalibTracks::FillResolutionHistoLocal(AliTPCseed * track){
    // if this chi2 is bigger than a given threshold, assume that the current cluster is
    // a kink an goto next padrow
    // if not:
-   // fill fArrayAmpRow, array with amplitudes vs. row for given sector
-   // fill fArrayAmp, array with amplitude histograms for give sector
-   // fill fRMSY, fRMSZ, fArrayQRMSY and fArrayQRMSZ, fDeltaY, fDeltaZ, fResolY, fResolZ, fArrayQDY, fArrayQDY
+   // fill fRMSY, fRMSZ, fArrayQRMSY and fArrayQRMSZ, fResolY, fResolZ, fArrayQDY, fArrayQDY
    // 
    // write debug information to 'TPCSelectorDebug.root'
    // only for every kDeltaWriteDebugStream'th padrow to reduce data volume 
    // and to avoid redundant data
    // 
 
-  static TLinearFitter fFitterLinY1(2,"pol1");   //
-  static TLinearFitter fFitterLinZ1(2,"pol1");   // 
-  static TLinearFitter fFitterLinY2(2,"pol1");   // 
-  static TLinearFitter fFitterLinZ2(2,"pol1");   //
   static TLinearFitter fFitterParY(3,"pol2");    // 
   static TLinearFitter fFitterParZ(3,"pol2");    //
-
-  fFitterLinY1.StoreData(kFALSE);
-  fFitterLinZ1.StoreData(kFALSE);
-  fFitterLinY2.StoreData(kFALSE);
-  fFitterLinZ2.StoreData(kFALSE);
-  fFitterParY.StoreData(kFALSE);
-  fFitterParZ.StoreData(kFALSE);
-
-
+  static TLinearFitter fFitterParYWeight(3,"pol2");    // 
+  static TLinearFitter fFitterParZWeight(3,"pol2");    //
+  fFitterParY.StoreData(kTRUE);
+  fFitterParZ.StoreData(kTRUE);
+  fFitterParYWeight.StoreData(kTRUE);
+  fFitterParZWeight.StoreData(kTRUE);
   if (GetDebugLevel() > 5) Info("FillResolutionHistoLocal"," ***** Start of FillResolutionHistoLocal *****");
-   const Int_t   kDelta    = 10;          // delta rows to fit
-   const Float_t kMinRatio = 0.75;        // minimal ratio
-   //   const Float_t kCutChi2  = 6.;          // cut chi2 - left right  - kink removal
-   const Float_t kErrorFraction = 0.5;    // use only clusters with small interpolation error - for error param
-   const Int_t   kFirstLargePad = 127;    // medium pads -> long pads
-   const Float_t kLargePadSize  = 1.5;    // factor between medium and long pads' area
-   const Int_t   kDeltaWriteDebugStream  = 5;  // only for every kDeltaWriteDebugStream'th padrow debug information is calulated and written to debugstream
-   TVectorD paramY0(2);
-   TVectorD paramZ0(2);
-   TVectorD paramY1(2);
-   TVectorD paramZ1(2);
-   TVectorD paramY2(3);
-   TVectorD paramZ2(3);
-   TMatrixD matrixY0(2,2);
-   TMatrixD matrixZ0(2,2);
-   TMatrixD matrixY1(2,2);
-   TMatrixD matrixZ1(2,2);
-   
-   // estimate mean error
-   Int_t nTrackletsAll = 0;       // number of tracklets for given track
-   Float_t csigmaY     = 0;       // mean sigma for tracklet refit in Y direction
-   Float_t csigmaZ     = 0;       // mean sigma for tracklet refit in Z direction
-   Int_t nClusters     = 0;       // working variable, number of clusters per tracklet
-   Int_t sectorG       = -1;      // working variable, sector of tracklet, has to stay constant for one tracklet
-   
-   fHclus->Fill(track->GetNumberOfClusters());      // for statistics overview
-   // ---------------------------------------------------------------------
-   for (Int_t irow = 0; irow < 159; irow++){
-      // loop over all rows along the track
-      // fit tracklets (length: 13 rows) with pol2 in Y and Z direction
-      // calculate mean chi^2 for this track-fit in Y and Z direction
-      AliTPCclusterMI * cluster0 = track->GetClusterPointer(irow);
-      if (!cluster0) continue;  // no cluster found
-      Int_t sector = cluster0->GetDetector();
-      fHclusterPerPadrowRaw->Fill(irow);
-      
-      Int_t ipad = TMath::Nint(cluster0->GetPad());
-      Float_t value = fCalPadClusterPerPadRaw->GetCalROC(sector)->GetValue((sector<36)?irow:irow-64, TMath::Nint(cluster0->GetPad()));
-      fCalPadClusterPerPadRaw->GetCalROC(sector)->SetValue((sector<36)?irow:irow-64, ipad, value + 1 );
-      
-      if (sector != sectorG){
-         // track leaves sector before it crossed enough rows to fit / initialization
-         nClusters = 0;
-         fFitterParY.ClearPoints();
-         fFitterParZ.ClearPoints();
-         sectorG = sector;
-      }
-      else {
-         nClusters++;
-         Double_t x = cluster0->GetX();
-         fFitterParY.AddPoint(&x, cluster0->GetY(), 1);
-         fFitterParZ.AddPoint(&x, cluster0->GetZ(), 1);
-         //
-         if ( nClusters >= kDelta + 3 ){  
-         // if more than 13 (kDelta+3) clusters were added to the fitters
-         // fit the tracklet, increase trackletCounter
-         fFitterParY.Eval();
-         fFitterParZ.Eval();
-         nTrackletsAll++;
-         csigmaY += fFitterParY.GetChisquare() / (nClusters - 3.);
-         csigmaZ += fFitterParZ.GetChisquare() / (nClusters - 3.);
-         nClusters = -1;
-         fFitterParY.ClearPoints();
-         fFitterParZ.ClearPoints();
-         }
-      }
-   }      // for (Int_t irow = 0; irow < 159; irow++)
-   // mean chi^2 for all tracklet fits in Y and in Z direction: 
-   csigmaY = TMath::Sqrt(TMath::Abs(csigmaY) / (nTrackletsAll+0.1));
-   csigmaZ = TMath::Sqrt(TMath::Abs(csigmaZ) / (nTrackletsAll+0.1));
-   // ---------------------------------------------------------------------
-   for (Int_t irow = 0; irow < 159; irow++){
-      // loop again over all rows along the track
-      // do analysis
-      // 
-      Int_t nclFound = 0;  // number of clusters in the neighborhood
-      Int_t ncl0 = 0;      // number of clusters in rows < rowOfCenterCluster
-      Int_t ncl1 = 0;      // number of clusters in rows > rowOfCenterCluster
-      AliTPCclusterMI * cluster0 = track->GetClusterPointer(irow);
-      if (!cluster0) continue;
-      Int_t sector = cluster0->GetDetector();
-      Float_t xref = cluster0->GetX();
-         
-      // Make Fit
+  const Int_t   kDelta     = 10;          // delta rows to fit
+  const Float_t kMinRatio  = 0.75;        // minimal ratio
+  const Float_t kChi2Cut   =  10.;          // cut chi2 - left right
+  const Float_t kSigmaCut  = 3.;        //sigma cut
+  const Float_t kdEdxCut=300;
+  const Float_t kPtCut=0.300;
+
+  if (track->GetTPCsignal()>kdEdxCut) return;  
+  if (TMath::Abs(AliTracker::GetBz())>0.1 &&TMath::Abs(track->Pt())<kPtCut) return;  
+
+  // estimate mean error
+  Int_t nTrackletsAll = 0;       // number of tracklets for given track
+  Float_t csigmaY     = 0;       // mean sigma for tracklet refit in Y direction
+  Float_t csigmaZ     = 0;       // mean sigma for tracklet refit in Z direction
+  Int_t nClusters     = 0;       // working variable, number of clusters per tracklet
+  Int_t sectorG       = -1;      // working variable, sector of tracklet, has to stay constant for one tracklet
+  Double_t refX=0;
+  // ---------------------------------------------------------------------
+  for (Int_t irow = 0; irow < 159; irow++){
+    // loop over all rows along the track
+    // fit tracklets (length: 13 rows) with pol2 in Y and Z direction
+    // calculate mean chi^2 for this track-fit in Y and Z direction
+    AliTPCclusterMI * cluster0 = track->GetClusterPointer(irow);
+    if (!cluster0) continue;  // no cluster found
+    Int_t sector = cluster0->GetDetector();
+    
+    Int_t ipad = TMath::Nint(cluster0->GetPad());
+    Float_t value = fCalPadClusterPerPadRaw->GetCalROC(sector)->GetValue((sector<36)?irow:irow-64, TMath::Nint(cluster0->GetPad()));
+    fCalPadClusterPerPadRaw->GetCalROC(sector)->SetValue((sector<36)?irow:irow-64, ipad, value + 1 );
+    
+    if (sector != sectorG){
+      // track leaves sector before it crossed enough rows to fit / initialization
+      nClusters = 0;
       fFitterParY.ClearPoints();
       fFitterParZ.ClearPoints();
-      fFitterLinY1.ClearPoints();
-      fFitterLinZ1.ClearPoints();
-      fFitterLinY2.ClearPoints();
-      fFitterLinZ2.ClearPoints();
-      
-      // fit tracklet (clusters in given padrow +- kDelta padrows) 
-      // with polynom of 2nd order and two polynoms of 1st order
-      // take both polynoms of 1st order, calculate difference of their parameters
-      // add covariance matrixes and calculate chi2 of this difference
-      // if this chi2 is bigger than a given threshold, assume that the current cluster is
-      // a kink an goto next padrow
-      
-      for (Int_t idelta = -kDelta; idelta <= kDelta; idelta++){
-         // loop over irow +- kDelta rows (neighboured rows)
-         // 
-         // 
-         if (idelta == 0) continue;                                // don't use center cluster
-         if (idelta + irow < 0 || idelta + irow > 159) continue;   // don't go out of ROC
-         AliTPCclusterMI * currentCluster = track->GetClusterPointer(irow + idelta);
-         if (!currentCluster) continue;
-         if (currentCluster->GetType() < 0) continue;
-         if (currentCluster->GetDetector() != sector) continue;
-         Double_t x = currentCluster->GetX() - xref;  // x = differece: current cluster - cluster @ irow
-         nclFound++;
-         if (idelta < 0){
-         ncl0++;
-         fFitterLinY1.AddPoint(&x, currentCluster->GetY(), csigmaY);
-         fFitterLinZ1.AddPoint(&x, currentCluster->GetZ(), csigmaZ);
-         }
-         if (idelta > 0){
-         ncl1++;
-         fFitterLinY2.AddPoint(&x, currentCluster->GetY(), csigmaY);
-         fFitterLinZ2.AddPoint(&x, currentCluster->GetZ(), csigmaZ);
-         }
-         fFitterParY.AddPoint(&x, currentCluster->GetY(), csigmaY);  
-         fFitterParZ.AddPoint(&x, currentCluster->GetZ(), csigmaZ);  
-      }  // loop over neighbourhood for fitter filling 
-
-
-      
-      if (nclFound < kDelta * kMinRatio) fRejectedTracksHisto->Fill(10);
-      if (nclFound < kDelta * kMinRatio) fClusterCutHisto->Fill(1, irow);
-      if (nclFound < kDelta * kMinRatio) continue;    // if not enough clusters (7.5) found in neighbourhood goto next padrow
-      fFitterParY.Eval();
-      fFitterParZ.Eval();
-      Double_t chi2 = (fFitterParY.GetChisquare() + fFitterParZ.GetChisquare()) / (2. * nclFound - 6.);
-      TTreeSRedirector *cstream = GetDebugStreamer();
-      if (cstream){
-       (*cstream)<<"Cut9"<<
-         "chi2="<<chi2<<
-         "\n";
-      }
-      // REMOVE KINK
-      // only when there are enough clusters (4) in each direction
-      if (ncl0 > 4){
-         fFitterLinY1.Eval();
-         fFitterLinZ1.Eval();
-      }
-      if (ncl1 > 4){
-         fFitterLinY2.Eval();
-         fFitterLinZ2.Eval();
-      }
-      
-      if (ncl0 > 4 && ncl1 > 4){
-         fFitterLinY1.GetCovarianceMatrix(matrixY0);
-         fFitterLinY2.GetCovarianceMatrix(matrixY1);
-         fFitterLinZ1.GetCovarianceMatrix(matrixZ0);
-         fFitterLinZ2.GetCovarianceMatrix(matrixZ1);
-         fFitterLinY2.GetParameters(paramY1);
-         fFitterLinZ2.GetParameters(paramZ1);
-         fFitterLinY1.GetParameters(paramY0);
-         fFitterLinZ1.GetParameters(paramZ0);
-         paramY0 -= paramY1;
-         paramZ0 -= paramZ1;
-         matrixY0 += matrixY1;
-         matrixZ0 += matrixZ1;
-         Double_t cchi2 = 0;
-         
-         TMatrixD difY(2, 1, paramY0.GetMatrixArray());
-         TMatrixD difYT(1, 2, paramY0.GetMatrixArray());
-         matrixY0.Invert();
-         TMatrixD mulY(matrixY0, TMatrixD::kMult, difY);
-         TMatrixD chi2Y(difYT, TMatrixD::kMult, mulY);
-         cchi2 += chi2Y(0, 0);
-         
-         TMatrixD difZ(2, 1, paramZ0.GetMatrixArray());
-         TMatrixD difZT(1, 2, paramZ0.GetMatrixArray());
-         matrixZ0.Invert();
-         TMatrixD mulZ(matrixZ0, TMatrixD::kMult, difZ);
-         TMatrixD chi2Z(difZT, TMatrixD::kMult, mulZ);
-         cchi2 += chi2Z(0, 0);      
-         
-        if (cstream){
-          (*cstream)<<"Cut8"<<
-            "chi2="<<cchi2<<
-            "\n";
-        }       
-      }
-      
-      // current padrow has no kink
-      
-      // get fit parameters from pol2 fit: 
-      Double_t paramY[4], paramZ[4];
-      paramY[0] = fFitterParY.GetParameter(0);
-      paramY[1] = fFitterParY.GetParameter(1);
-      paramY[2] = fFitterParY.GetParameter(2);
-      paramZ[0] = fFitterParZ.GetParameter(0);
-      paramZ[1] = fFitterParZ.GetParameter(1);
-      paramZ[2] = fFitterParZ.GetParameter(2);    
-      
-      Double_t tracky = paramY[0];
-      Double_t trackz = paramZ[0];
-      Float_t  deltay = tracky - cluster0->GetY();
-      Float_t  deltaz = trackz - cluster0->GetZ();
-      Float_t  angley = paramY[1] - paramY[0] / xref;
-      Float_t  anglez = paramZ[1];
-      
-      Float_t max = cluster0->GetMax();
-      UInt_t isegment = cluster0->GetDetector() % 36;
-      Int_t padSize = 0;                          // short pads
-      if (cluster0->GetDetector() >= 36) {
-         padSize = 1;                              // medium pads 
-         if (cluster0->GetRow() > 63) padSize = 2; // long pads
+      sectorG = sector;
+    }
+    else {
+      nClusters++;
+      if (refX<1) refX=cluster0->GetX()+kDelta*0.5;
+      Double_t x = cluster0->GetX()-refX;
+      fFitterParY.AddPoint(&x, cluster0->GetY(), 1);
+      fFitterParZ.AddPoint(&x, cluster0->GetZ(), 1);
+      //
+      if ( nClusters >= kDelta + 3 ){  
+       // if more than 13 (kDelta+3) clusters were added to the fitters
+       // fit the tracklet, increase trackletCounter
+       fFitterParY.Eval();
+       fFitterParZ.Eval();
+       nTrackletsAll++;
+       csigmaY += fFitterParY.GetChisquare() / (nClusters - 3.);
+       csigmaZ += fFitterParZ.GetChisquare() / (nClusters - 3.);
+       nClusters = -1;
+       fFitterParY.ClearPoints();
+       fFitterParZ.ClearPoints();
+       refX=0;
       }
+    }
+  }      // for (Int_t irow = 0; irow < 159; irow++)
+  // mean chi^2 for all tracklet fits in Y and in Z direction: 
+  csigmaY = TMath::Sqrt(TMath::Abs(csigmaY) / (nTrackletsAll+0.1));
+  csigmaZ = TMath::Sqrt(TMath::Abs(csigmaZ) / (nTrackletsAll+0.1));
+  // ---------------------------------------------------------------------
+  //
+  //
 
-      // =========================================
-      // wirte collected information to histograms
-      // =========================================
-      
-      TProfile *profAmpRow =  (TProfile*)fArrayAmpRow->At(sector);
-      if ( irow >= kFirstLargePad) max /= kLargePadSize;
-      Double_t smax = TMath::Sqrt(max);
-      profAmpRow->Fill( (Double_t)cluster0->GetRow(), smax );
-      TH1F *hisAmp =  (TH1F*)fArrayAmp->At(sector);
-      hisAmp->Fill(smax);
-      
-      // remove the following two lines one day:
-      TProfile *profDriftLength = (TProfile*)fArrayChargeVsDriftlength->At(sector);
-      profDriftLength->Fill( 250.-(Double_t)TMath::Abs(cluster0->GetZ()), smax );
-      
-      TProfile *profDriftLengthTmp = (TProfile*)(fcalPadRegionChargeVsDriftlength->GetObject(isegment, padSize));
-      profDriftLengthTmp->Fill( 250.-(Double_t)TMath::Abs(cluster0->GetZ()), smax );
-      
-      fHclusterPerPadrow->Fill(irow);   // fill histogram showing clusters per padrow
-      Int_t ipad = TMath::Nint(cluster0->GetPad());
-      Float_t value = fCalPadClusterPerPad->GetCalROC(sector)->GetValue((sector<36)?irow:irow-64, TMath::Nint(cluster0->GetPad()));
-      fCalPadClusterPerPad->GetCalROC(sector)->SetValue((sector<36)?irow:irow-64, ipad, value + 1 );
-   
-         
-      TH3F * his3 = 0;
-      his3 = (TH3F*)fRMSY->At(padSize);
-      if (his3) his3->Fill(250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), TMath::Sqrt(cluster0->GetSigmaY2()) );
-      his3 = (TH3F*)fRMSZ->At(padSize);
-      if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), TMath::Sqrt(cluster0->GetSigmaZ2()) );
-      
-      his3 = (TH3F*)fArrayQRMSY->At(GetBin(cluster0->GetMax(), padSize));
-      if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), TMath::Sqrt(cluster0->GetSigmaY2()) );
-      his3 = (TH3F*)fArrayQRMSZ->At(GetBin(cluster0->GetMax(), padSize));
-      if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), TMath::Sqrt(cluster0->GetSigmaZ2()) );
-   
-      
-      // Fill resolution histograms
-      Bool_t useForResol = kTRUE;
-      if (fFitterParY.GetParError(0) > kErrorFraction * csigmaY) useForResol = kFALSE;
-   
-      if (cstream){
-       Float_t zdrift = 250 - TMath::Abs(cluster0->GetZ());
-       Float_t sy = cluster0->GetSigmaY2();
-       Float_t sz = cluster0->GetSigmaZ2();
-       (*cstream)<<"Resol0"<<
-         "run="<<fRun<<              //  run number
-         "event="<<fEvent<<          //  event number
-         "time="<<fTime<<            //  time stamp of event
-         "trigger="<<fTrigger<<      //  trigger
-         "mag="<<fMagF<<             //  magnetic field              
-         "padSize="<<padSize<<
-         "angley="<<angley<<
-         "anglez="<<anglez<<
-         "zdr="<<zdrift<<
-         "dy="<<deltay<<
-         "dz="<<deltaz<<
-         "sy="<<sy<<
-         "sz="<<sz<<
-         "\n";
+  for (Int_t irow = kDelta; irow < 159-kDelta; irow++){
+    // loop again over all rows along the track
+    // do analysis
+    // 
+    Int_t nclFound = 0;  // number of clusters in the neighborhood
+    Int_t ncl0 = 0;      // number of clusters in rows < rowOfCenterCluster
+    Int_t ncl1 = 0;      // number of clusters in rows > rowOfCenterCluster
+    AliTPCclusterMI * cluster0 = track->GetClusterPointer(irow);
+    if (!cluster0) continue;
+    Int_t sector = cluster0->GetDetector();
+    Float_t xref = cluster0->GetX();
+    
+    // Make Fit
+    fFitterParY.ClearPoints();
+    fFitterParZ.ClearPoints();    
+    fFitterParYWeight.ClearPoints();
+    fFitterParZWeight.ClearPoints();    
+    // fit tracklet (clusters in given padrow +- kDelta padrows) 
+    // with polynom of 2nd order and two polynoms of 1st order
+    // take both polynoms of 1st order, calculate difference of their parameters
+    // add covariance matrixes and calculate chi2 of this difference
+    // if this chi2 is bigger than a given threshold, assume that the current cluster is
+    // a kink an goto next padrow    
+    AliRieman riemanFit(2*kDelta);
+    AliRieman riemanFitW(2*kDelta);
+    for (Int_t idelta = -kDelta; idelta <= kDelta; idelta++){
+      // loop over irow +- kDelta rows (neighboured rows)
+      // 
+      // 
+      if (idelta == 0) continue;                                // don't use center cluster
+      if (idelta + irow < 0 || idelta + irow > 159) continue;   // don't go out of ROC
+      AliTPCclusterMI * currentCluster = track->GetClusterPointer(irow + idelta);
+      if (!currentCluster) continue;
+      if (currentCluster->GetType() < 0) continue;
+      if (currentCluster->GetDetector() != sector) continue;
+      nclFound++;
+      if (idelta < 0){
+       ncl0++;
       }
-
-      if (useForResol){
-         fDeltaY->Fill(deltay);
-         fDeltaZ->Fill(deltaz);
-         his3 = (TH3F*)fResolY->At(padSize);
-         if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), deltay );
-         his3 = (TH3F*)fResolZ->At(padSize);
-         if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), deltaz );
-         his3 = (TH3F*)fArrayQDY->At(GetBin(cluster0->GetMax(), padSize));
-         if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()),TMath::Abs(angley), deltay );
-         his3 = (TH3F*)fArrayQDZ->At(GetBin(cluster0->GetMax(), padSize));
-         if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()),TMath::Abs(anglez), deltaz );
+      if (idelta > 0){
+       ncl1++;
       }
-      
-      //=============================================================================================
-      
-      if (useForResol && nclFound > 2 * kMinRatio * kDelta 
-         && irow % kDeltaWriteDebugStream == 0 && GetDebugLevel() > 4){
-       if (GetDebugLevel() > 20) Info("FillResolutionHistoLocal","Filling 'TPCSelectorDebug.root', irow = %i", irow);
-         FillResolutionHistoLocalDebugPart(track, cluster0, irow, angley, anglez, nclFound, kDelta);
-      }  // if (useForResol && nclFound > 2 * kMinRatio * kDelta)
-      //
-      // Fill THN histograms
-      //
-      Double_t xvar[9];
-      xvar[1]=padSize;
-      xvar[2]=1.-TMath::Abs(cluster0->GetZ())/250.;
-      xvar[3]=cluster0->GetMax();
-      xvar[5]=angley;
-      xvar[6]=anglez;
-
-      xvar[4]=cluster0->GetPad()-TMath::Nint(cluster0->GetPad());      
-      xvar[0]=deltay;
-      fHisDeltaY->Fill(xvar);
-      xvar[0]=TMath::Sqrt(cluster0->GetSigmaY2());
-      fHisRMSY->Fill(xvar);
-
-      xvar[4]=cluster0->GetTimeBin()-TMath::Nint(cluster0->GetTimeBin());
-      xvar[0]=deltaz;
-      fHisDeltaZ->Fill(xvar);
-      xvar[0]=TMath::Sqrt(cluster0->GetSigmaZ2());
-      fHisRMSZ->Fill(xvar);
-   
-   }    // loop over all padrows along the track: for (Int_t irow = 0; irow < 159; irow++)
-}  // FillResolutionHistoLocal(...)
-
-
-
-void AliTPCcalibTracks::FillResolutionHistoLocalDebugPart(AliTPCseed *track, AliTPCclusterMI *cluster0, Int_t irow, Float_t  angley, Float_t  anglez, Int_t nclFound, Int_t kDelta) {
-   // 
-   //  - debug part of FillResolutionHistoLocal - 
-   // called only for every kDeltaWriteDebugStream'th padrow, to avoid to much redundant data
-   // called only for GetStreamLevel() > 4
-   // fill resolution trees
-   //
-      
-   Int_t sector = cluster0->GetDetector();
-   Float_t xref = cluster0->GetX();
-   Int_t padSize = 0;                          // short pads
-   if (cluster0->GetDetector() >= 36) {
+      riemanFit.AddPoint(currentCluster->GetX(), currentCluster->GetY(),currentCluster->GetZ(), csigmaY,csigmaZ);
+      riemanFitW.AddPoint(currentCluster->GetX(), currentCluster->GetY(),currentCluster->GetZ(), csigmaY*TMath::Sqrt(1+TMath::Abs(idelta)),csigmaZ*TMath::Sqrt(1+TMath::Abs(idelta)));
+    }  // loop over neighbourhood for fitter filling 
+    if (nclFound < kDelta * kMinRatio) fRejectedTracksHisto->Fill(10);
+    if (nclFound < kDelta * kMinRatio) fClusterCutHisto->Fill(1, irow);
+    if (nclFound < kDelta * kMinRatio) continue;    // if not enough clusters (7.5) found in neighbourhood goto next padrow
+    riemanFit.Update();
+    riemanFitW.Update();
+    Double_t chi2R=TMath::Sqrt(riemanFit.CalcChi2()/(2*nclFound-5));
+    Double_t chi2RW=TMath::Sqrt(riemanFitW.CalcChi2()/(2*nclFound-5));
+    Double_t paramYR[3], paramZR[3];
+    Double_t paramYRW[3], paramZRW[3];
+    //
+    paramYR[0]    = riemanFit.GetYat(xref);
+    paramZR[0]    = riemanFit.GetZat(xref);
+    paramYRW[0]   = riemanFitW.GetYat(xref);
+    paramZRW[0]   = riemanFitW.GetZat(xref);
+    //
+    paramYR[1]    = riemanFit.GetDYat(xref);
+    paramZR[1]    = riemanFit.GetDZat(xref);
+    paramYRW[1]   = riemanFitW.GetDYat(xref);
+    paramZRW[1]   = riemanFitW.GetDZat(xref);
+    //
+    Int_t reject=0;
+    if (chi2R>kChi2Cut) reject+=1;
+    if (chi2RW>kChi2Cut) reject+=2;
+    if (TMath::Abs(paramYR[0]-paramYRW[0])>kSigmaCut*csigmaY) reject+=4;
+    if (TMath::Abs(paramZR[0]-paramZRW[0])>kSigmaCut*csigmaZ) reject+=8;
+    if (TMath::Abs(paramYR[1]-paramYRW[1])>csigmaY) reject+=16;
+    if (TMath::Abs(paramZR[1]-paramZRW[1])>csigmaZ) reject+=32;
+    //
+    TTreeSRedirector *cstream = GetDebugStreamer();    
+    // get fit parameters from pol2 fit:     
+    Double_t tracky = paramYR[0];
+    Double_t trackz = paramZR[0];
+    Float_t  deltay = cluster0->GetY()-tracky;
+    Float_t  deltaz = cluster0->GetZ()-trackz;
+    Float_t  angley = paramYR[1];
+    Float_t  anglez = paramZR[1];
+    
+    Int_t padSize = 0;                          // short pads
+    if (cluster0->GetDetector() >= 36) {
       padSize = 1;                              // medium pads 
       if (cluster0->GetRow() > 63) padSize = 2; // long pads
-   }
-      
-   static TLinearFitter fitY0(3, "pol2");
-   static TLinearFitter fitZ0(3, "pol2");
-   static TLinearFitter fitY2(5, "hyp4");
-   static TLinearFitter fitZ2(5, "hyp4");
-   static TLinearFitter fitY2Q(5, "hyp4");
-   static TLinearFitter fitZ2Q(5, "hyp4");
-   static TLinearFitter fitY2S(5, "hyp4");
-   static TLinearFitter fitZ2S(5, "hyp4");
-   fitY0.ClearPoints();
-   fitZ0.ClearPoints();
-   fitY2.ClearPoints();
-   fitZ2.ClearPoints();
-   fitY2Q.ClearPoints();
-   fitZ2Q.ClearPoints();
-   fitY2S.ClearPoints();
-   fitZ2S.ClearPoints();
-   
-   for (Int_t idelta = -kDelta; idelta <= kDelta; idelta++){
-     // loop over irow +- kDelta rows (neighboured rows)
-     // 
-     // 
-     if (idelta == 0) continue;
-     if (idelta + irow < 0 || idelta + irow > 159) continue;   // don't go out of ROC
-     AliTPCclusterMI * cluster = track->GetClusterPointer(irow + idelta);
-     if (!cluster) continue;
-     if (cluster->GetType() < 0) continue;
-     if (cluster->GetDetector() != sector) continue;
-     Double_t x = cluster->GetX() - xref;
-     Double_t sigmaY0 = fClusterParam->GetError0Par( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(angley) );
-     Double_t sigmaZ0 = fClusterParam->GetError0Par( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(anglez) );
-     //
-     Double_t sigmaYQ = fClusterParam->GetErrorQPar( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(angley), TMath::Abs(cluster->GetMax()) );
-     Double_t sigmaZQ = fClusterParam->GetErrorQPar( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(anglez), TMath::Abs(cluster->GetMax()) );
-     Double_t sigmaYS = fClusterParam->GetErrorQParScaled( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())), 
-                                                           TMath::Abs(angley), TMath::Abs(cluster->GetMax()) );
-     Double_t sigmaZS = fClusterParam->GetErrorQParScaled( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())), 
-                                                           TMath::Abs(anglez), TMath::Abs(cluster->GetMax()) );
-     Float_t rmsYFactor = fClusterParam->GetShapeFactor( 0, padSize,(250.0 - TMath::Abs(cluster->GetZ())),
-                                                        TMath::Abs(anglez), TMath::Abs(cluster->GetMax()),
-                                                        TMath::Sqrt(cluster0->GetSigmaY2()), 0 );
-     Float_t rmsZFactor = fClusterParam->GetShapeFactor(0, padSize,(250.0 - TMath::Abs(cluster->GetZ())),
-                                                       TMath::Abs(anglez), TMath::Abs(cluster->GetMax()),
-                                                       TMath::Sqrt(cluster0->GetSigmaZ2()),0 );
-     sigmaYS  = TMath::Sqrt(sigmaYS * sigmaYS + rmsYFactor * rmsYFactor / 12.);
-     sigmaZS  = TMath::Sqrt(sigmaZS * sigmaZS + rmsZFactor * rmsZFactor / 12. + rmsYFactor * rmsYFactor / 24.);
-     //
-     if (kDelta != 0){
-       fitY0.AddPoint(&x, cluster->GetY(), sigmaY0);
-       fitZ0.AddPoint(&x, cluster->GetZ(), sigmaZ0);
-     }
-     Double_t xxx[4];
-     xxx[0] = ( (idelta+irow) % 2 == 0 ) ? 1 : 0;
-     xxx[1] = x;
-     xxx[2] = ( (idelta+irow) % 2 == 0 ) ? x : 0;
-     xxx[3] = x * x;   
-     fitY2.AddPoint(xxx, cluster->GetY(), sigmaY0);
-     fitY2Q.AddPoint(xxx, cluster->GetY(), sigmaYQ);
-     fitY2S.AddPoint(xxx, cluster->GetY(), sigmaYS);
-     fitZ2.AddPoint(xxx, cluster->GetZ(), sigmaZ0);
-     fitZ2Q.AddPoint(xxx, cluster->GetZ(), sigmaZQ);
-     fitZ2S.AddPoint(xxx, cluster->GetZ(), sigmaZS);
-     //
-   }  // neigbouhood-loop
-      //
-   fitY0.Eval();
-   fitZ0.Eval();
-   fitY2.Eval();
-   fitZ2.Eval();
-   fitY2Q.Eval();
-   fitZ2Q.Eval();
-   fitY2S.Eval();
-   fitZ2S.Eval();
-   Float_t chi2Y0 = fitY0.GetChisquare() / (nclFound-3.);
-   Float_t chi2Z0 = fitZ0.GetChisquare() / (nclFound-3.);
-   Float_t chi2Y2 = fitY2.GetChisquare() / (nclFound-5.);
-   Float_t chi2Z2 = fitZ2.GetChisquare() / (nclFound-5.);
-   Float_t chi2Y2Q = fitY2Q.GetChisquare() / (nclFound-5.);
-   Float_t chi2Z2Q = fitZ2Q.GetChisquare() / (nclFound-5.);
-   Float_t chi2Y2S = fitY2S.GetChisquare() / (nclFound-5.);
-   Float_t chi2Z2S = fitZ2S.GetChisquare() / (nclFound-5.);
-   //
-   static  TVectorD    parY0(3);
-   static  TMatrixD    matY0(3, 3);
-   static  TVectorD    parZ0(3);
-   static  TMatrixD    matZ0(3, 3);
-   fitY0.GetParameters(parY0);
-   fitY0.GetCovarianceMatrix(matY0);
-   fitZ0.GetParameters(parZ0);
-   fitZ0.GetCovarianceMatrix(matZ0);
-   //
-   static  TVectorD    parY2(5);
-   static  TMatrixD    matY2(5,5);
-   static  TVectorD    parZ2(5);
-   static  TMatrixD    matZ2(5,5);
-   fitY2.GetParameters(parY2);
-   fitY2.GetCovarianceMatrix(matY2);
-   fitZ2.GetParameters(parZ2);
-   fitZ2.GetCovarianceMatrix(matZ2);
-   //
-   static  TVectorD    parY2Q(5);
-   static  TMatrixD    matY2Q(5,5);
-   static  TVectorD    parZ2Q(5);
-   static  TMatrixD    matZ2Q(5,5);
-   fitY2Q.GetParameters(parY2Q);
-   fitY2Q.GetCovarianceMatrix(matY2Q);
-   fitZ2Q.GetParameters(parZ2Q);
-   fitZ2Q.GetCovarianceMatrix(matZ2Q);
-   static  TVectorD    parY2S(5);
-   static  TMatrixD    matY2S(5,5);
-   static  TVectorD    parZ2S(5);
-   static  TMatrixD    matZ2S(5,5);
-   fitY2S.GetParameters(parY2S);
-   fitY2S.GetCovarianceMatrix(matY2S);
-   fitZ2S.GetParameters(parZ2S);
-   fitZ2S.GetCovarianceMatrix(matZ2S);
-   Float_t sigmaY0   = TMath::Sqrt(matY0(0,0));
-   Float_t sigmaZ0   = TMath::Sqrt(matZ0(0,0));
-   Float_t sigmaDY0  = TMath::Sqrt(matY0(1,1));
-   Float_t sigmaDZ0  = TMath::Sqrt(matZ0(1,1));
-   Float_t sigmaY2   = TMath::Sqrt(matY2(1,1));
-   Float_t sigmaZ2   = TMath::Sqrt(matZ2(1,1));
-   Float_t sigmaDY2  = TMath::Sqrt(matY2(3,3));
-   Float_t sigmaDZ2  = TMath::Sqrt(matZ2(3,3));
-   Float_t sigmaY2Q  = TMath::Sqrt(matY2Q(1,1));
-   Float_t sigmaZ2Q  = TMath::Sqrt(matZ2Q(1,1));
-   Float_t sigmaDY2Q = TMath::Sqrt(matY2Q(3,3));
-   Float_t sigmaDZ2Q = TMath::Sqrt(matZ2Q(3,3));
-   Float_t sigmaY2S  = TMath::Sqrt(matY2S(1,1));
-   Float_t sigmaZ2S  = TMath::Sqrt(matZ2S(1,1));
-   Float_t sigmaDY2S = TMath::Sqrt(matY2S(3,3));
-   Float_t sigmaDZ2S = TMath::Sqrt(matZ2S(3,3));
-   
-   // Error parameters
-   Float_t csigmaY0 = fClusterParam->GetError0Par(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),TMath::Abs(angley));
-   Float_t csigmaZ0 = fClusterParam->GetError0Par(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),TMath::Abs(anglez));
-   //
-   Float_t csigmaYQ = fClusterParam->GetErrorQPar(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
-                                                 TMath::Abs(angley), TMath::Abs(cluster0->GetMax()));
-   Float_t csigmaZQ = fClusterParam->GetErrorQPar(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
-                                                 TMath::Abs(anglez),TMath::Abs(cluster0->GetMax()));
-   Float_t csigmaYS = fClusterParam->GetErrorQParScaled(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
-                                                       TMath::Abs(angley), TMath::Abs(cluster0->GetMax()));
-   Float_t csigmaZS = fClusterParam->GetErrorQParScaled(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
-                                                       TMath::Abs(anglez),TMath::Abs(cluster0->GetMax()));
-   
-   
-   // RMS parameters
-   Float_t meanRMSY = 0;
-   Float_t meanRMSZ = 0;
-   Int_t   nclRMS = 0;
-   for (Int_t idelta = -2; idelta <= 2; idelta++){
-     // loop over neighbourhood
-     if (idelta+irow < 0 || idelta+irow > 159) continue;
-     //        if (idelta+irow>159) continue;
-     AliTPCclusterMI * cluster = track->GetClusterPointer(irow+idelta);
-     if (!cluster) continue;
-     meanRMSY += TMath::Sqrt(cluster->GetSigmaY2());
-     meanRMSZ += TMath::Sqrt(cluster->GetSigmaZ2());
-     nclRMS++;
-   }
-   meanRMSY /= nclRMS; 
-   meanRMSZ /= nclRMS; 
-   
-   Float_t rmsY      = TMath::Sqrt(cluster0->GetSigmaY2());  
-   Float_t rmsZ      = TMath::Sqrt(cluster0->GetSigmaZ2());
-   Float_t rmsYT     = fClusterParam->GetRMSQ(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
-                                             TMath::Abs(angley), TMath::Abs(cluster0->GetMax()));
-   Float_t rmsZT     = fClusterParam->GetRMSQ(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
-                                             TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()));
-   Float_t rmsYT0    = fClusterParam->GetRMS0(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
-                                             TMath::Abs(angley));
-   Float_t rmsZT0    = fClusterParam->GetRMS0(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
-                                                TMath::Abs(anglez));
-   Float_t rmsYSigma = fClusterParam->GetRMSSigma(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
-                                                 TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()));
-   Float_t rmsZSigma = fClusterParam->GetRMSSigma(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
-                                                 TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()));
-   Float_t rmsYFactor = fClusterParam->GetShapeFactor(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
-                                                     TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()),
-                                                     rmsY,meanRMSY);
-   Float_t rmsZFactor = fClusterParam->GetShapeFactor(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
-                                                     TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()),
-                                                     rmsZ,meanRMSZ);
-   //
-   // cluster debug
-   TTreeSRedirector *cstream = GetDebugStreamer();
-   if (cstream){
-     (*cstream)<<"ResolCl"<<   // valgrind 3   40,000 bytes in 1 blocks are possibly 
-       "run="<<fRun<<              //  run number
-       "event="<<fEvent<<          //  event number
-       "time="<<fTime<<            //  time stamp of event
-       "trigger="<<fTrigger<<      //  trigger
-       "mag="<<fMagF<<             //  magnetic field        
-       "Sector="<<sector<<
-       "Cl.="<<cluster0<<
-       "CSigmaY0="<<csigmaY0<<   // cluster errorY
-       "CSigmaYQ="<<csigmaYQ<<   // cluster errorY - q scaled
-       "CSigmaYS="<<csigmaYS<<   // cluster errorY - q scaled
-       "CSigmaZ0="<<csigmaZ0<<   // 
-       "CSigmaZQ="<<csigmaZQ<<
-       "CSigmaZS="<<csigmaZS<<
-       "shapeYF="<<rmsYFactor<<
-       "shapeZF="<<rmsZFactor<<
-       "rmsY="<<rmsY<<
-       "rmsZ="<<rmsZ<<
-       "rmsYM="<<meanRMSY<<
-       "rmsZM="<<meanRMSZ<<
-       "rmsYT="<<rmsYT<<
-       "rmsZT="<<rmsZT<<
-       "rmsYT0="<<rmsYT0<<
-       "rmsZT0="<<rmsZT0<<
-       "rmsYS="<<rmsYSigma<<  
-       "rmsZS="<<rmsZSigma<<
-       "padSize="<<padSize<<
-       "Ncl="<<nclFound<<      
-       "PY0.="<<&parY0<<
-       "PZ0.="<<&parZ0<<
-       "SigmaY0="<<sigmaY0<< 
-       "SigmaZ0="<<sigmaZ0<< 
-       "angley="<<angley<<
-       "anglez="<<anglez<<
-       "\n";
-     
-     //       tracklet dubug
-     (*cstream)<<"ResolTr"<<   
-       "run="<<fRun<<              //  run number
-       "event="<<fEvent<<          //  event number
-       "time="<<fTime<<            //  time stamp of event
-       "trigger="<<fTrigger<<      //  trigger
-       "mag="<<fMagF<<             //  magnetic field        
-       "padSize="<<padSize<<
-       "IPad="<<padSize<<
-       "Sector="<<sector<<
-       "Ncl="<<nclFound<<      
-       "chi2Y0="<<chi2Y0<<
-       "chi2Z0="<<chi2Z0<<
-       "chi2Y2="<<chi2Y2<<
-       "chi2Z2="<<chi2Z2<<
-       "chi2Y2Q="<<chi2Y2Q<<
-       "chi2Z2Q="<<chi2Z2Q<<
-       "chi2Y2S="<<chi2Y2S<<
-       "chi2Z2S="<<chi2Z2S<<
-       "PY0.="<<&parY0<<
-       "PZ0.="<<&parZ0<<
-       "PY2.="<<&parY2<<
-       "PZ2.="<<&parZ2<<
-       "PY2Q.="<<&parY2Q<<
-       "PZ2Q.="<<&parZ2Q<<
-       "PY2S.="<<&parY2S<<
-       "PZ2S.="<<&parZ2S<<
-       "SigmaY0="<<sigmaY0<< 
-       "SigmaZ0="<<sigmaZ0<< 
-       "SigmaDY0="<<sigmaDY0<< 
-       "SigmaDZ0="<<sigmaDZ0<< 
-       "SigmaY2="<<sigmaY2<< 
-       "SigmaZ2="<<sigmaZ2<< 
-       "SigmaDY2="<<sigmaDY2<< 
-       "SigmaDZ2="<<sigmaDZ2<< 
-       "SigmaY2Q="<<sigmaY2Q<< 
-       "SigmaZ2Q="<<sigmaZ2Q<< 
-       "SigmaDY2Q="<<sigmaDY2Q<< 
-       "SigmaDZ2Q="<<sigmaDZ2Q<< 
-       "SigmaY2S="<<sigmaY2S<< 
-       "SigmaZ2S="<<sigmaZ2S<< 
-       "SigmaDY2S="<<sigmaDY2S<< 
-       "SigmaDZ2S="<<sigmaDZ2S<< 
-       "angley="<<angley<<
-       "anglez="<<anglez<<
-       "\n";
-   }  
-}
+    }
+    Int_t ipad = TMath::Nint(cluster0->GetPad());
+    if (cstream){
+      Float_t zdrift = 250 - TMath::Abs(cluster0->GetZ());
+      (*cstream)<<"Resol0"<<   
+       "run="<<fRun<<              //  run number
+       "event="<<fEvent<<          //  event number
+       "time="<<fTime<<            //  time stamp of event
+       "trigger="<<fTrigger<<      //  trigger
+       "mag="<<fMagF<<             //  magnetic field        
+       "padSize="<<padSize<<
+       //
+       "reject="<<reject<<
+       "cl.="<<cluster0<<          // cluster info
+       "xref="<<xref<<             // cluster reference X
+       //rieman fit
+       "yr="<<paramYR[0]<<         // track position y - rieman fit
+       "zr="<<paramZR[0]<<         // track position z - rieman fit 
+       "yrW="<<paramYRW[0]<<         // track position y - rieman fit - weight
+       "zrW="<<paramZRW[0]<<         // track position z - rieman fit - weight
+       "dyr="<<paramYR[1]<<         // track position y - rieman fit
+       "dzr="<<paramZR[1]<<         // track position z - rieman fit 
+       "dyrW="<<paramYRW[1]<<         // track position y - rieman fit - weight
+       "dzrW="<<paramZRW[1]<<         // track position z - rieman fit - weight
+       //
+       "angley="<<angley<<         // angle par fit
+       "anglez="<<anglez<<         // angle par fit
+       "zdr="<<zdrift<<            //
+       "dy="<<deltay<<
+       "dz="<<deltaz<<        
+       "sy="<<csigmaY<<
+       "sz="<<csigmaZ<<
+       "chi2R="<<chi2R<<
+       "chi2RW="<<chi2RW<<
+       "\n";
+    }    
+    if (reject) continue;
+
+    
+    // =========================================
+    // wirte collected information to histograms
+    // =========================================
+        
+    Float_t value = fCalPadClusterPerPad->GetCalROC(sector)->GetValue((sector<36)?irow:irow-64, TMath::Nint(cluster0->GetPad()));
+    fCalPadClusterPerPad->GetCalROC(sector)->SetValue((sector<36)?irow:irow-64, ipad, value + 1 );
+    
+    
+    TH3F * his3 = 0;
+    his3 = (TH3F*)fRMSY->At(padSize);
+    if (his3) his3->Fill(250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), TMath::Sqrt(cluster0->GetSigmaY2()) );
+    his3 = (TH3F*)fRMSZ->At(padSize);
+    if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), TMath::Sqrt(cluster0->GetSigmaZ2()) );
+    
+    his3 = (TH3F*)fArrayQRMSY->At(GetBin(cluster0->GetMax(), padSize));
+    if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), TMath::Sqrt(cluster0->GetSigmaY2()) );
+    his3 = (TH3F*)fArrayQRMSZ->At(GetBin(cluster0->GetMax(), padSize));
+    if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), TMath::Sqrt(cluster0->GetSigmaZ2()) );
+    
+    
+    his3 = (TH3F*)fResolY->At(padSize);
+    if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), deltay );
+    his3 = (TH3F*)fResolZ->At(padSize);
+    if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), deltaz );
+    his3 = (TH3F*)fArrayQDY->At(GetBin(cluster0->GetMax(), padSize));
+    if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()),TMath::Abs(angley), deltay );
+    his3 = (TH3F*)fArrayQDZ->At(GetBin(cluster0->GetMax(), padSize));
+    if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()),TMath::Abs(anglez), deltaz );        
+    //=============================================================================================    
+    //
+    // Fill THN histograms
+    //
+    Double_t xvar[9];
+    xvar[1]=padSize;
+    xvar[2]=cluster0->GetZ();
+    xvar[3]=cluster0->GetMax();
+    
+    xvar[0]=deltay;
+    xvar[4]=cluster0->GetPad()-Int_t(cluster0->GetPad());      
+    xvar[5]=angley;
+    fHisDeltaY->Fill(xvar);
+    xvar[0]=TMath::Sqrt(cluster0->GetSigmaY2());
+    fHisRMSY->Fill(xvar);
+    
+    xvar[0]=deltaz;
+    xvar[4]=cluster0->GetTimeBin()-Int_t(cluster0->GetTimeBin());
+    xvar[5]=anglez;
+    fHisDeltaZ->Fill(xvar);
+    xvar[0]=TMath::Sqrt(cluster0->GetSigmaZ2());
+    fHisRMSZ->Fill(xvar);
+    
+  }    // loop over all padrows along the track: for (Int_t irow = 0; irow < 159; irow++)
+}  // FillResolutionHistoLocal(...)
+
 
 
 
 
 
-TH2D * AliTPCcalibTracks::MakeDiff(TH2D * hfit, TF2 * func){
-   // 
-   // creates a new histogram which contains the difference between
-   // the histogram hfit and the function func
-   // 
-  TH2D * result = (TH2D*)hfit->Clone();      // valgrind 3   40,139 bytes in 11 blocks are still reachable
-  result->SetTitle(Form("%s fit residuals",result->GetTitle()));
-  result->SetName(Form("%s fit residuals",result->GetName()));
-  TAxis *xaxis = hfit->GetXaxis();
-  TAxis *yaxis = hfit->GetYaxis();
-  Double_t x[2];
-  for (Int_t biny = 0; biny <= yaxis->GetNbins(); biny++) {
-    x[1]  = yaxis->GetBinCenter(biny);
-    for (Int_t binx = 0; binx <= xaxis->GetNbins(); binx++) {
-      x[0]  = xaxis->GetBinCenter(binx);
-      Int_t bin = hfit->GetBin(binx, biny);
-      Double_t val = hfit->GetBinContent(bin);
-//      result->SetBinContent( bin, (val - func->Eval(x[0], x[1])) / func->Eval(x[0], x[1]) );
-      result->SetBinContent( bin, (val / func->Eval(x[0], x[1])) - 1 );
-    }    
-  }
-  return result;
-}
 
 
 void  AliTPCcalibTracks::SetStyle() const {
@@ -1340,67 +798,6 @@ void  AliTPCcalibTracks::SetStyle() const {
 }
 
 
-void AliTPCcalibTracks::Draw(Option_t* opt){
-   // 
-   // draw-function of AliTPCcalibTracks
-   // will draws some exemplaric pictures
-   // 
-
-  if (GetDebugLevel() > 6) Info("Draw", "Drawing an exemplaric picture.");
-   SetStyle();
-   Double_t min = 0;
-   Double_t max = 0;
-   TCanvas *c1 = new TCanvas();
-   c1->Divide(0, 3);
-   TVirtualPad *upperThird = c1->GetPad(1);
-   TVirtualPad *middleThird = c1->GetPad(2);
-   TVirtualPad *lowerThird = c1->GetPad(3);
-   upperThird->Divide(2,0);
-   TVirtualPad *upleft  = upperThird->GetPad(1);
-   TVirtualPad *upright = upperThird->GetPad(2);
-   middleThird->Divide(2,0);
-   TVirtualPad *middleLeft  = middleThird->GetPad(1);
-   TVirtualPad *middleRight = middleThird->GetPad(2);
-   lowerThird->Divide(2,0);
-   TVirtualPad *downLeft  = lowerThird->GetPad(1);
-   TVirtualPad *downRight = lowerThird->GetPad(2);
-   
-   
-   upleft->cd(0);
-   min = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())-20;
-   max = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())+20;
-   fDeltaY->SetAxisRange(min, max);
-   fDeltaY->Fit("gaus","q","",min, max);        // valgrind 3  7 block possibly lost   2,400 bytes in 1 blocks are still reachable
-   c1->Update();
-   
-   upright->cd(0);
-   max = fDeltaZ->GetBinCenter(fDeltaZ->GetMaximumBin())+20;
-   min = fDeltaZ->GetBinCenter(fDeltaZ->GetMaximumBin())-20;
-   fDeltaZ->SetAxisRange(min, max);
-   fDeltaZ->Fit("gaus","q","",min, max);
-   c1->Update();
-   
-   middleLeft->cd();
-   fHclus->Draw(opt);
-   
-   middleRight->cd();
-   fRejectedTracksHisto->Draw(opt);
-   TPaveText *pt = new TPaveText(0.6,0.6, 0.8,0.8, "NDC");
-   TText *t1 = pt->AddText("1: kEdgeThetaCutNoise");
-   TText *t2 = pt->AddText("2: kMinClusters");
-   TText *t3 = pt->AddText("3: kMinRatio");
-   TText *t4 = pt->AddText("4: kMax1pt");
-   t1 = t1; t2 = t2; t3 = t3; t4 = t4;    // avoid compiler warnings
-   pt->SetToolTipText("Legend for failed cuts");
-   pt->Draw();
-   
-   downLeft->cd();
-   fHclusterPerPadrowRaw->Draw(opt);
-   
-   downRight->cd();
-   fHclusterPerPadrow->Draw(opt);
-}
-
 
 void AliTPCcalibTracks::MakeReport(Int_t stat, const char* pathName){ 
    // 
@@ -1411,471 +808,10 @@ void AliTPCcalibTracks::MakeReport(Int_t stat, const char* pathName){
    // 
 
   if (GetDebugLevel() > 0) Info("MakeReport","Writing plots and trees to '%s'.", pathName);
-   MakeAmpPlots(stat, pathName);
-   MakeDeltaPlots(pathName);
-   FitResolutionNew(pathName);
-   FitRMSNew(pathName);
-   MakeChargeVsDriftLengthPlots(pathName);
-//    MakeResPlotsQ(1, 1); 
-   MakeResPlotsQTree(stat, pathName);
+  MakeResPlotsQTree(stat, pathName);
 }
    
 
-void AliTPCcalibTracks::MakeAmpPlots(Int_t stat, const char* pathName){ 
-   // 
-   // creates several plots:
-   // fArrayAmp.ps, fArrayAmpRow.ps and DeltaYZ.ps
-   // fArrayAmp.ps: one histogram per sector, the histogram shows the charge per cluster
-   // fArrayAmpRow.ps: one histogram per sector, mean max. amplitude vs. pad row with landau fit
-   // DeltaYZ.ps: DeltaY and DeltaZ histogram with gaus fit
-   // Empty histograms (sectors without data) are not written to file
-   // the ps-files are written to the directory 'pathName', that is created if it does not exist
-   // 'stat': only histograms with more than 'stat' entries are written to file.
-   // 
-   
-   SetStyle();
-   gSystem->MakeDirectory(pathName);
-   gSystem->ChangeDirectory(pathName);
-   
-   TCanvas* c1 = new TCanvas();     // valgrind 3 ???  634 bytes in 28 blocks are still reachable
-   TPostScript *ps; 
-   // histograms with accumulated amplitude for all IROCs and OROCs
-   TH1F *allAmpHisIROC = ((TH1F*)(fArrayAmp->At(0))->Clone());
-   allAmpHisIROC->SetName("Amp all IROCs");
-   allAmpHisIROC->SetTitle("Amp all IROCs");
-   TH1F *allAmpHisOROC = ((TH1F*)(fArrayAmp->At(36))->Clone());
-   allAmpHisOROC->SetName("Amp all OROCs");
-   allAmpHisOROC->SetTitle("Amp all OROCs");
-   
-   ps = new TPostScript("fArrayAmp.ps", 112);
-   if (GetDebugLevel() > -1) cout << "creating fArrayAmp.ps..." << endl;
-   for (Int_t i = 0; i < fArrayAmp->GetEntriesFast(); i++){
-      if ( ((TH1F*)fArrayAmp->At(i))->GetEntries() < stat  ) continue;
-      ps->NewPage();
-      ((TH1F*)fArrayAmp->At(i))->Draw();
-      c1->Update();              // valgrind 3
-      if (i > 0 && i < 36) {
-         allAmpHisIROC->Add(((TH1F*)fArrayAmp->At(i)));
-         allAmpHisOROC->Add(((TH1F*)fArrayAmp->At(i+36)));
-      }
-   }
-   ps->NewPage();
-   allAmpHisIROC->Draw();
-   c1->Update();              // valgrind
-   ps->NewPage();
-   allAmpHisOROC->Draw();
-   c1->Update();
-   ps->Close();
-   delete ps;
-   
-   TH1F *his = 0;
-   Double_t min = 0;
-   Double_t max = 0;
-   ps = new TPostScript("fArrayAmpRow.ps", 112);
-   if (GetDebugLevel() > -1) cout << "creating fArrayAmpRow.ps..." << endl;
-   for (Int_t i = 0; i < fArrayAmpRow->GetEntriesFast(); i++){
-      his = (TH1F*)fArrayAmpRow->At(i);
-      if (his->GetEntries() < stat) continue;
-      ps->NewPage();
-      min = TMath::Max( his->GetBinCenter(his->GetMaximumBin() )-100., 0.);
-      max = his->GetBinCenter(5*his->GetMaximumBin()) + 100;
-      his->SetAxisRange(min, max);
-      his->Fit("pol3", "q", "", min, max);
-      // his->Draw("error");    // don't use this line when you don't want to have empty pages in the ps-file
-      c1->Update();
-   }
-   ps->Close();
-   delete ps;
-   delete c1;
-   gSystem->ChangeDirectory("..");
-}
-
-
-void AliTPCcalibTracks::MakeDeltaPlots(const char* pathName){
-   // 
-   // creates several plots:
-   // DeltaYZ.ps: DeltaY and DeltaZ histogram with gaus fit
-   // the ps-files are written to the directory 'pathName', that is created if it does not exist
-   // 
-   
-   SetStyle();
-   gSystem->MakeDirectory(pathName);
-   gSystem->ChangeDirectory(pathName);
-   
-   TCanvas* c1 = new TCanvas();     // valgrind 3 ???  634 bytes in 28 blocks are still reachable
-   TPostScript *ps; 
-   Double_t min = 0;
-   Double_t max = 0;
-   
-   ps = new TPostScript("DeltaYZ.ps", 112);
-   if (GetDebugLevel() > -1) cout << "creating DeltaYZ.ps..." << endl;
-   min = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())-20;
-   max = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())+20;
-   fDeltaY->SetAxisRange(min, max);
-   ps->NewPage();
-   fDeltaY->Fit("gaus","q","",min, max);        // valgrind 3  7 block possibly lost   2,400 bytes in 1 blocks are still reachable
-   c1->Update();
-   ps->NewPage();
-   max = fDeltaZ->GetBinCenter(fDeltaZ->GetMaximumBin())+20;
-   min = fDeltaZ->GetBinCenter(fDeltaZ->GetMaximumBin())-20;
-   fDeltaZ->SetAxisRange(min, max);
-   fDeltaZ->Fit("gaus","q","",min, max);
-   c1->Update();
-   ps->Close();
-   delete ps;
-   delete c1;
-   gSystem->ChangeDirectory("..");
-}
-
-
-void AliTPCcalibTracks::MakeChargeVsDriftLengthPlotsOld(const char* pathName){
-   // 
-   // creates charge vs. driftlength plots, one TProfile for each ROC
-   // is not correct like this, should be one TProfile for each sector and padsize
-   // 
-   
-   SetStyle();
-   gSystem->MakeDirectory(pathName);
-   gSystem->ChangeDirectory(pathName);
-   
-   TCanvas* c1 = new TCanvas();     // valgrind 3 ???  634 bytes in 28 blocks are still reachable
-   TPostScript *ps; 
-   ps = new TPostScript("chargeVsDriftlengthOld.ps", 112);
-   if (GetDebugLevel() > -1) cout << "creating chargeVsDriftlength.ps..." << endl;
-   TProfile *chargeVsDriftlengthAllIROCs = ((TProfile*)fArrayChargeVsDriftlength->At(0)->Clone());
-   TProfile *chargeVsDriftlengthAllOROCs = ((TProfile*)fArrayChargeVsDriftlength->At(36)->Clone());
-   chargeVsDriftlengthAllIROCs->SetName("allAmpHisIROC");
-   chargeVsDriftlengthAllIROCs->SetTitle("charge vs. driftlength, all IROCs");
-   chargeVsDriftlengthAllOROCs->SetName("allAmpHisOROC");
-   chargeVsDriftlengthAllOROCs->SetTitle("charge vs. driftlength, all OROCs");
-   
-   for (Int_t i = 0; i < fArrayChargeVsDriftlength->GetEntriesFast(); i++) {
-      ((TProfile*)fArrayChargeVsDriftlength->At(i))->Draw();
-      c1->Update();
-      if (i > 0 && i < 36) { 
-         chargeVsDriftlengthAllIROCs->Add(((TProfile*)fArrayChargeVsDriftlength->At(i)));
-         chargeVsDriftlengthAllOROCs->Add(((TProfile*)fArrayChargeVsDriftlength->At(i+36)));
-      }
-      ps->NewPage();
-   }
-   chargeVsDriftlengthAllIROCs->Draw();
-   c1->Update();              // valgrind
-   ps->NewPage();
-   chargeVsDriftlengthAllOROCs->Draw();
-   c1->Update();
-   ps->Close();  
-   delete ps;
-   delete c1;
-   gSystem->ChangeDirectory("..");
-}   
-
-
-void AliTPCcalibTracks::MakeChargeVsDriftLengthPlots(const char* pathName){
-   // 
-   // creates charge vs. driftlength plots, one TProfile for each ROC
-   // under development....
-   // 
-   
-   SetStyle();
-   gSystem->MakeDirectory(pathName);
-   gSystem->ChangeDirectory(pathName);
-   
-   TCanvas* c1 = new TCanvas("c1", "c1", 700,(Int_t)(TMath::Sqrt(2)*700));     // valgrind 3 ???  634 bytes in 28 blocks are still reachable
-//    TCanvas c1("c1", "c1", 500,(sqrt(2)*500))
-   c1->Divide(0,3);
-   TPostScript *ps; 
-   ps = new TPostScript("chargeVsDriftlength.ps", 111);
-   if (GetDebugLevel() > -1) cout << "creating chargeVsDriftlengthNew.ps..." << endl;
-   
-   TProfile *chargeVsDriftlengthAllShortPads  = ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,0)->Clone());
-   TProfile *chargeVsDriftlengthAllMediumPads = ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,1)->Clone());
-   TProfile *chargeVsDriftlengthAllLongPads   = ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,2)->Clone());
-   chargeVsDriftlengthAllShortPads->SetName("allAmpHisShortPads");
-   chargeVsDriftlengthAllShortPads->SetTitle("charge vs. driftlength, all sectors, short pads");
-   chargeVsDriftlengthAllMediumPads->SetName("allAmpHisMediumPads");
-   chargeVsDriftlengthAllMediumPads->SetTitle("charge vs. driftlength, all sectors, medium pads");
-   chargeVsDriftlengthAllLongPads->SetName("allAmpHisLongPads");
-   chargeVsDriftlengthAllLongPads->SetTitle("charge vs. driftlength, all sectors, long pads");
-   
-   for (Int_t i = 0; i < 36; i++) {
-      c1->cd(1)->SetGridx();
-      c1->cd(1)->SetGridy();
-      ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(i,0))->Draw();
-      c1->cd(2)->SetGridx();
-      c1->cd(2)->SetGridy();
-      ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(i,1))->Draw();
-      c1->cd(3)->SetGridx();
-      c1->cd(3)->SetGridy();
-      ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(i,2))->Draw();
-      c1->Update();
-      chargeVsDriftlengthAllShortPads->Add( (TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,0));
-      chargeVsDriftlengthAllMediumPads->Add((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,1));
-      chargeVsDriftlengthAllLongPads->Add(  (TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,2));
-      ps->NewPage();
-   }
-   c1->cd(1)->SetGridx();
-   c1->cd(1)->SetGridy();
-   chargeVsDriftlengthAllShortPads->Draw();
-   c1->cd(2)->SetGridx();
-   c1->cd(2)->SetGridy();
-   chargeVsDriftlengthAllMediumPads->Draw();
-   c1->cd(3)->SetGridx();
-   c1->cd(3)->SetGridy();
-   chargeVsDriftlengthAllLongPads->Draw();
-   c1->Update();              // valgrind
-//    ps->NewPage();
-   ps->Close();  
-   delete ps;
-   delete c1;
-   gSystem->ChangeDirectory("..");
-}   
-
-
-
-void AliTPCcalibTracks::FitResolutionNew(const char* pathName){
-   // 
-   // calculates different resulution fits in Y and Z direction
-   // the histograms are written to 'ResolutionYZ.ps'
-   // writes calculated resolution to 'resol.txt'
-   // all files are stored in the directory pathName
-   // 
-   
-   SetStyle();
-   gSystem->MakeDirectory(pathName);
-   gSystem->ChangeDirectory(pathName);
-   
-   TCanvas c;
-   c.Divide(2,1); 
-   if (GetDebugLevel() > -1) cout << "creating ResolutionYZ.ps..." << endl;
-   TPostScript *ps = new TPostScript("ResolutionYZ.ps", 112); 
-   TF2 *fres = new TF2("fres","TMath::Sqrt([0]*[0]+[1]*[1]*x+[2]*[2]*y*y)",0,250,0,1);
-   fres->SetParameter(0,0.02);
-   fres->SetParameter(1,0.0054);
-   fres->SetParameter(2,0.13);  
-   
-   TH1::AddDirectory(kTRUE);  // TH3F::FitSlicesZ() writes histograms into the current directory
-   
-   // create histogramw for Y-resolution
-   TH3F * hisResY0 = (TH3F*)fResolY->At(0);
-   hisResY0->FitSlicesZ();
-   TH2D * hisResY02 = (TH2D*)gDirectory->Get("Resol Y0_2");
-   TH3F * hisResY1 = (TH3F*)fResolY->At(1); 
-   hisResY1->FitSlicesZ();
-   TH2D * hisResY12 = (TH2D*)gDirectory->Get("Resol Y1_2");
-   TH3F * hisResY2 = (TH3F*)fResolY->At(2);
-   hisResY2->FitSlicesZ();
-   TH2D * hisResY22 = (TH2D*)gDirectory->Get("Resol Y2_2");
-    //
-   ps->NewPage();
-   c.cd(1);
-   hisResY02->Fit(fres, "q");      // valgrind    132,072 bytes in 6 blocks are indirectly lost
-   hisResY02->Draw("surf1"); 
-   c.cd(2);
-   MakeDiff(hisResY02,fres)->Draw("surf1");
-   c.Update();
-   //   c.SaveAs("ResolutionYPad0.eps");
-   ps->NewPage();
-   c.cd(1);
-   hisResY12->Fit(fres, "q");
-   hisResY12->Draw("surf1");
-   c.cd(2);
-   MakeDiff(hisResY12,fres)->Draw("surf1");
-   c.Update();
-   //   c.SaveAs("ResolutionYPad1.eps");
-   ps->NewPage();
-   c.cd(1);
-   hisResY22->Fit(fres, "q");
-   hisResY22->Draw("surf1"); 
-   c.cd(2);
-   MakeDiff(hisResY22,fres)->Draw("surf1");
-   c.Update();
-//    c.SaveAs("ResolutionYPad2.eps");
-   
-   // create histogramw for Z-resolution
-   TH3F * hisResZ0 = (TH3F*)fResolZ->At(0);
-   hisResZ0->FitSlicesZ();
-   TH2D * hisResZ02 = (TH2D*)gDirectory->Get("Resol Z0_2");
-   TH3F * hisResZ1 = (TH3F*)fResolZ->At(1);
-   hisResZ1->FitSlicesZ();
-   TH2D * hisResZ12 = (TH2D*)gDirectory->Get("Resol Z1_2");
-   TH3F * hisResZ2 = (TH3F*)fResolZ->At(2);
-   hisResZ2->FitSlicesZ();
-   TH2D * hisResZ22 = (TH2D*)gDirectory->Get("Resol Z2_2");
-   
-   ps->NewPage();
-   c.cd(1);
-   hisResZ02->Fit(fres, "q");
-   hisResZ02->Draw("surf1");
-   c.cd(2);
-   MakeDiff(hisResZ02,fres)->Draw("surf1");
-   c.Update();
-//    c.SaveAs("ResolutionZPad0.eps");
-   ps->NewPage();
-   c.cd(1);
-   hisResZ12->Fit(fres, "q");
-   hisResZ12->Draw("surf1"); 
-   c.cd(2);
-   MakeDiff(hisResZ12,fres)->Draw("surf1");
-   c.Update();
-//    c.SaveAs("ResolutionZPad1.eps");
-   ps->NewPage();
-   c.cd(1);
-   hisResZ22->Fit(fres, "q");
-   hisResZ22->Draw("surf1");  
-   c.cd(2);
-   MakeDiff(hisResZ22,fres)->Draw("surf1");
-   c.Update();
-//    c.SaveAs("ResolutionZPad2.eps");
-   ps->Close();
-   delete ps;
-   
-   // write calculated resoltuions to 'resol.txt'
-   ofstream fresol("resol.txt");
-   fresol<<"Pad 0.75 cm"<<"\n";
-   hisResY02->Fit(fres, "q");                     // valgrind
-   fresol<<"Y\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
-   hisResZ02->Fit(fres, "q");
-   fresol<<"Z\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
-   //
-   fresol<<"Pad 1.00 cm"<<1<<"\n";
-   hisResY12->Fit(fres, "q");                     // valgrind
-   fresol<<"Y\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
-   hisResZ12->Fit(fres, "q");
-   fresol<<"Z\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
-   //
-   fresol<<"Pad 1.50 cm"<<0<<"\n";
-   hisResY22->Fit(fres, "q");
-   fresol<<"Y\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
-   hisResZ22->Fit(fres, "q");
-   fresol<<"Z\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
-   
-   TH1::AddDirectory(kFALSE);
-   gSystem->ChangeDirectory("..");
-   delete fres;
-}
-
-
-void AliTPCcalibTracks::FitRMSNew(const char* pathName){
-   // 
-   // calculates different resulution-rms fits in Y and Z direction
-   // the histograms are written to 'RMS_YZ.ps'
-   // writes calculated resolution rms to 'rms.txt'
-   // all files are stored in the directory pathName
-   // 
-   
-   SetStyle();
-   gSystem->MakeDirectory(pathName);
-   gSystem->ChangeDirectory(pathName);
-   
-   TCanvas c;        // valgrind 3   42,120 bytes in 405 blocks are still reachable   23,816 bytes in 229 blocks are still reachable
-   c.Divide(2,1); 
-   if (GetDebugLevel() > -1) cout << "creating RMS_YZ.ps..." << endl;
-   TPostScript *ps = new TPostScript("RMS_YZ.ps", 112); 
-   TF2 *frms = new TF2("fres","TMath::Sqrt([0]*[0]+[1]*[1]*x+[2]*[2]*y*y)",0,250,0,1);
-   frms->SetParameter(0,0.02);
-   frms->SetParameter(1,0.0054);
-   frms->SetParameter(2,0.13);  
-   
-   TH1::AddDirectory(kTRUE);  // TH3F::FitSlicesZ() writes histograms into the current directory
-   
-   // create histogramw for Y-RMS   
-   TH3F * hisResY0 = (TH3F*)fRMSY->At(0);
-   hisResY0->FitSlicesZ();
-   TH2D * hisResY02 = (TH2D*)gDirectory->Get("RMS Y0_1");
-   TH3F * hisResY1 = (TH3F*)fRMSY->At(1);
-   hisResY1->FitSlicesZ();
-   TH2D * hisResY12 = (TH2D*)gDirectory->Get("RMS Y1_1");
-   TH3F * hisResY2 = (TH3F*)fRMSY->At(2);
-   hisResY2->FitSlicesZ();
-   TH2D * hisResY22 = (TH2D*)gDirectory->Get("RMS Y2_1");
-   //
-   ps->NewPage();
-   c.cd(1);
-   hisResY02->Fit(frms, "qn0"); 
-   hisResY02->Draw("surf1"); 
-   c.cd(2);
-   MakeDiff(hisResY02,frms)->Draw("surf1");
-   c.Update();
-//    c.SaveAs("RMSYPad0.eps");
-   ps->NewPage();
-   c.cd(1);
-   hisResY12->Fit(frms, "qn0");               // valgrind   several blocks possibly lost
-   hisResY12->Draw("surf1");
-   c.cd(2);
-   MakeDiff(hisResY12,frms)->Draw("surf1");
-   c.Update();
-//    c.SaveAs("RMSYPad1.eps");
-   ps->NewPage();
-   c.cd(1);
-   hisResY22->Fit(frms, "qn0");
-   hisResY22->Draw("surf1"); 
-   c.cd(2);
-   MakeDiff(hisResY22,frms)->Draw("surf1");
-   c.Update();
-//    c.SaveAs("RMSYPad2.eps");
-   
-   // create histogramw for Z-RMS   
-   TH3F * hisResZ0 = (TH3F*)fRMSZ->At(0);
-   hisResZ0->FitSlicesZ();
-   TH2D * hisResZ02 = (TH2D*)gDirectory->Get("RMS Z0_1");
-   TH3F * hisResZ1 = (TH3F*)fRMSZ->At(1); 
-   hisResZ1->FitSlicesZ();
-   TH2D * hisResZ12 = (TH2D*)gDirectory->Get("RMS Z1_1");
-   TH3F * hisResZ2 = (TH3F*)fRMSZ->At(2); 
-   hisResZ2->FitSlicesZ();
-   TH2D * hisResZ22 = (TH2D*)gDirectory->Get("RMS Z2_1");
-   //
-   ps->NewPage();
-   c.cd(1);
-   hisResZ02->Fit(frms, "qn0");         // valgrind
-   hisResZ02->Draw("surf1");
-   c.cd(2);
-   MakeDiff(hisResZ02,frms)->Draw("surf1");
-   c.Update();
-//    c.SaveAs("RMSZPad0.eps");
-   ps->NewPage();
-   c.cd(1);
-   hisResZ12->Fit(frms, "qn0");
-   hisResZ12->Draw("surf1"); 
-   c.cd(2);
-   MakeDiff(hisResZ12,frms)->Draw("surf1");
-   c.Update();
-//    c.SaveAs("RMSZPad1.eps");
-   ps->NewPage();
-   c.cd(1);
-   hisResZ22->Fit(frms, "qn0");         // valgrind  1 block possibly lost
-   hisResZ22->Draw("surf1");  
-   c.cd(2);
-   MakeDiff(hisResZ22,frms)->Draw("surf1");
-   c.Update();
-//    c.SaveAs("RMSZPad2.eps");
-   
-   // write calculated resoltuion rms to 'rms.txt'
-   ofstream filerms("rms.txt");
-   filerms<<"Pad 0.75 cm"<<"\n";
-   hisResY02->Fit(frms, "qn0");         // valgrind   23 blocks indirectly lost
-   filerms<<"Y\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
-   hisResZ02->Fit(frms, "qn0");         // valgrind   23 blocks indirectly lost
-   filerms<<"Z\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
-   //
-   filerms<<"Pad 1.00 cm"<<1<<"\n";
-   hisResY12->Fit(frms, "qn0");         // valgrind      3,256 bytes in 22 blocks are indirectly lost 
-   filerms<<"Y\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
-   hisResZ12->Fit(frms, "qn0");         // valgrind    66,036 bytes in 3 blocks are still reachable
-   filerms<<"Z\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
-   //
-   filerms<<"Pad 1.50 cm"<<0<<"\n";
-   hisResY22->Fit(frms, "qn0");      // valgrind   40,139 bytes in 11 blocks are still reachable   330,180 bytes in 15 blocks are possibly lost
-   filerms<<"Y\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
-   hisResZ22->Fit(frms, "qn0");
-   filerms<<"Z\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
-
-   TH1::AddDirectory(kFALSE);
-   gSystem->ChangeDirectory("..");
-   ps->Close();
-   delete ps;
-   delete frms;
-}
 
 
 void AliTPCcalibTracks::MakeResPlotsQTree(Int_t minEntries, const char* pathName){
@@ -2209,17 +1145,12 @@ Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) {
    TList* deltaZList = new TList;
    TList* arrayAmpRowList = new TList;
    TList* rejectedTracksList = new TList;
-   TList* hclusList = new TList;
    TList* clusterCutHistoList = new TList;
    TList* arrayAmpList = new TList;
    TList* arrayQDYList = new TList;
    TList* arrayQDZList = new TList;
    TList* arrayQRMSYList = new TList;
    TList* arrayQRMSZList = new TList;
-   TList* arrayChargeVsDriftlengthList = new TList;
-   TList* calPadRegionChargeVsDriftlengthList = new TList;
-   TList* hclusterPerPadrowList = new TList;
-   TList* hclusterPerPadrowRawList = new TList;
    TList* resolYList = new TList;
    TList* resolZList = new TList;
    TList* rMSYList = new TList;
@@ -2236,10 +1167,6 @@ Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) {
    while ( (calibTracks = dynamic_cast<AliTPCcalibTracks*> (listIterator->Next())) ){
       // loop over all entries in the collectionList and get dataMembers into lists
       
-      deltaYList->Add( calibTracks->GetfDeltaY() );
-      deltaZList->Add( calibTracks->GetfDeltaZ() );
-      arrayAmpRowList->Add(calibTracks->GetfArrayAmpRow());
-      arrayAmpList->Add(calibTracks->GetfArrayAmp());
       arrayQDYList->Add(calibTracks->GetfArrayQDY());
       arrayQDZList->Add(calibTracks->GetfArrayQDZ());
       arrayQRMSYList->Add(calibTracks->GetfArrayQRMSY());
@@ -2248,13 +1175,8 @@ Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) {
       resolZList->Add(calibTracks->GetfResolZ());
       rMSYList->Add(calibTracks->GetfRMSY());
       rMSZList->Add(calibTracks->GetfRMSZ());
-      arrayChargeVsDriftlengthList->Add(calibTracks->GetfArrayChargeVsDriftlength());
-      calPadRegionChargeVsDriftlengthList->Add(calibTracks->GetCalPadRegionchargeVsDriftlength());
-      hclusList->Add(calibTracks->GetfHclus());
       rejectedTracksList->Add(calibTracks->GetfRejectedTracksHisto());
       clusterCutHistoList->Add(calibTracks->GetfClusterCutHisto());
-      hclusterPerPadrowList->Add(calibTracks->GetfHclusterPerPadrow());
-      hclusterPerPadrowRawList->Add(calibTracks->GetfHclusterPerPadrowRaw());
       //
       if (fCalPadClusterPerPad && calibTracks->GetfCalPadClusterPerPad())
        fCalPadClusterPerPad->Add(calibTracks->GetfCalPadClusterPerPad());      
@@ -2267,49 +1189,15 @@ Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) {
    
    // merge data members
    if (GetDebugLevel() > 0) cout << "histogram's merge-functins are called... " << endl; 
-   fDeltaY->Merge(deltaYList);
-   fDeltaZ->Merge(deltaZList);
-   fHclus->Merge(hclusList);
    fClusterCutHisto->Merge(clusterCutHistoList);
    fRejectedTracksHisto->Merge(rejectedTracksList);
-   fHclusterPerPadrow->Merge(hclusterPerPadrowList);
-   fHclusterPerPadrowRaw->Merge(hclusterPerPadrowRawList);
    
    TObjArray* objarray = 0;
    TH1* hist = 0;
    TList* histList = 0;
    TIterator *objListIterator = 0;
    
-   if (GetDebugLevel() > 0) cout << "merging fArrayAmpRows..." << endl;
-   // merge fArrayAmpRows
-   for (Int_t i = 0; i < fArrayAmpRow->GetEntriesFast(); i++ ) {  // loop over data member, i<72
-      objListIterator = arrayAmpRowList->MakeIterator();
-      histList = new TList;
-      while (( objarray =  (TObjArray*)objListIterator->Next() )) { 
-         // loop over arrayAmpRowList, get TObjArray, get object at position i, cast it into TProfile
-         hist = (TProfile*)(objarray->At(i));
-         histList->Add(hist);
-      }
-      ((TProfile*)(fArrayAmpRow->At(i)))->Merge(histList);
-      delete histList;
-      delete objListIterator;
-   }
-   
-   if (GetDebugLevel() > 0) cout << "merging fArrayAmps..." << endl;
-   // merge fArrayAmps
-   for (Int_t i = 0; i < fArrayAmp->GetEntriesFast(); i++ ) {  // loop over data member, i<72
-      TIterator *cobjListIterator = arrayAmpList->MakeIterator();
-      histList = new TList;
-      while (( objarray =  (TObjArray*)cobjListIterator->Next() )) { 
-         // loop over arrayAmpList, get TObjArray, get object at position i, cast it into TH1F
-         hist = (TH1F*)(objarray->At(i));
-         histList->Add(hist);
-      }
-      ((TH1F*)(fArrayAmp->At(i)))->Merge(histList);
-      delete histList;
-      delete cobjListIterator;
-   }
-   
+      
    if (GetDebugLevel() > 0) cout << "merging fArrayQDY..." << endl;
    // merge fArrayQDY
    for (Int_t i = 0; i < fArrayQDY->GetEntriesFast(); i++) { // loop over data member, i < 300
@@ -2373,52 +1261,7 @@ Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) {
       delete objListIterator;
    } 
    
-   if (GetDebugLevel() > 0) cout << "merging fArrayChargeVsDriftlength..." << endl;
-   // merge fArrayChargeVsDriftlength
-   for (Int_t i = 0; i < fArrayChargeVsDriftlength->GetEntriesFast(); i++) { // loop over data member, i < 300
-      objListIterator = arrayChargeVsDriftlengthList->MakeIterator();
-      histList = new TList;
-      while (( objarray =  (TObjArray*)objListIterator->Next() )) { 
-         // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TProfile
-         if (!objarray) continue;
-         hist = (TProfile*)(objarray->At(i));
-         histList->Add(hist);
-      }
-      ((TProfile*)(fArrayChargeVsDriftlength->At(i)))->Merge(histList);
-      delete histList;
-      delete objListIterator;
-   }    
-   
-   if (GetDebugLevel() > 0) cout << "merging fcalPadRegionChargeVsDriftlength..." << endl;
-   // merge fcalPadRegionChargeVsDriftlength
-   AliTPCCalPadRegion *cpr = 0x0;
-   
-   /*
-   TIterator *regionIterator = fcalPadRegionChargeVsDriftlength->MakeIterator();
-   while (hist = (TProfile*)regionIterator->Next()) {
-      // loop over all calPadRegion's in destination calibTracks object
-         objListIterator = calPadRegionChargeVsDriftlengthList->MakeIterator();
-         while (( cpr =  (AliTPCCalPadRegion*)objListIterator->Next() )) { 
-
-      
-      hist->Merge(...);
-   }
-   */
    
-   for (UInt_t isec = 0; isec < 36; isec++) {
-      for (UInt_t padSize = 0; padSize < 3; padSize++){
-         objListIterator = calPadRegionChargeVsDriftlengthList->MakeIterator();
-         histList = new TList;
-         while (( cpr =  (AliTPCCalPadRegion*)objListIterator->Next() )) { 
-            // loop over calPadRegionChargeVsDriftlengthList, get AliTPCCalPadRegion, get object 
-            hist = (TProfile*)cpr->GetObject(isec, padSize);
-            histList->Add(hist);
-         }
-         ((TProfile*)(fcalPadRegionChargeVsDriftlength->GetObject(isec, padSize)))->Merge(histList);
-         delete histList;
-         delete objListIterator;
-      }
-   }
   
    
         
@@ -2504,179 +1347,6 @@ Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) {
 }
 
 
-AliTPCcalibTracks* AliTPCcalibTracks::TestMerge(AliTPCcalibTracks *ct, AliTPCClusterParam *clusterParam, Int_t nCalTracks){
-   // 
-   // function to test AliTPCcalibTrack::Merge:
-   // in the file 'f' is a AliTPCcalibTrack object with name "calibTracks"
-   // this object is appended 'nCalTracks' times to a TList
-   // A new AliTPCcalibTrack object is created which merges the list
-   // this object is returned
-   // 
-   /*
-   // .L AliTPCcalibTracks.cxx+g
-   .L libTPCcalib.so
-   TFile f("Output.root");
-   AliTPCcalibTracks* calTracks = (AliTPCcalibTracks*)f.Get("calibTracks");
-   //f.Close();
-   TFile clusterParamFile("/u/lbozyk/calibration/workdir/calibTracks/TPCClusterParam.root");
-   AliTPCClusterParam *clusterParam  =  (AliTPCClusterParam *) clusterParamFile.Get("Param"); 
-   clusterParamFile.Close();
-
-   AliTPCcalibTracks::TestMerge(calTracks, clusterParam);
-   */
-   TList *list = new TList();
-   if (ct == 0 || clusterParam == 0) return 0;
-   cout << "making list with " << nCalTracks << " AliTPCcalibTrack objects" << endl;
-   for (Int_t i = 0; i < nCalTracks; i++) {
-      if (i%10==0) cout << "Adding element " << i << " of " << nCalTracks << endl;
-      list->Add(new AliTPCcalibTracks(*ct));
-   }
-   
-   // only for check at the end
-   AliTPCcalibTracks* cal1 = new AliTPCcalibTracks(*ct);
-   Double_t cal1Entries = ((TH1F*)cal1->GetfArrayAmpRow()->At(5))->GetEntries();
-//    Double_t cal1Entries = 5; //((TH1F*)ct->GetfArrayAmpRow()->At(5))->GetEntries();
-
-   cout  << "The list contains " << list->GetEntries() << " entries. " << endl;
-  
-   
-   AliTPCcalibTracksCuts *cuts = new AliTPCcalibTracksCuts(20, 0.4, 0.5, 0.13, 0.018);
-   AliTPCcalibTracks* cal = new AliTPCcalibTracks("calTracksMerged", "calTracksMerged", clusterParam, cuts, 5);
-   cal->Merge(list);
-   
-   cout << "cal->GetfArrayAmpRow()->At(5)->Print():" << endl;
-   cal->GetfArrayAmpRow()->At(5)->Print();
-   Double_t calEntries = ((TH1F*)cal->GetfArrayAmpRow()->At(5))->GetEntries();
-   
-   cout << "cal1->GetfArrayAmpRow()->At(5))->GetEntries() = " << cal1Entries << endl;
-   cout << " cal->GetfArrayAmpRow()->At(5))->GetEntries() = " <<  calEntries << endl;
-   printf("That means there were %f / %f = %f AliTPCcalibTracks-Objects merged. \n", 
-      calEntries, cal1Entries, ((Double_t)calEntries/cal1Entries));
-   
-   return cal;
-
-}
-
-
-void  AliTPCcalibTracks::MakeQPosNormAll(TTree * chainres, AliTPCClusterParam * param, Int_t maxPoints){
-  //
-  // Make position corrections
-  // for the moment Only using debug streamer 
-  // chainres  - debug tree
-  // param     - parameters to be updated
-  // maxPoints - maximal number of points using for fit
-  // verbose   - print info flag
-  //
-  // Current implementation - only using debug streamers
-  // 
-  
-  /*    
-    //Defaults
-    Int_t maxPoints=100000;
-  */
-  /*
-    Usage: 
-    //0. Load libraries
-    gSystem->Load("libANALYSIS");
-    gSystem->Load("libSTAT");
-    gSystem->Load("libTPCcalib");
-    
-
-    //1. Load Parameters to be modified:
-    //e.g:
-    
-    AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
-    AliCDBManager::Instance()->SetRun(0) 
-    AliTPCClusterParam * param = AliTPCcalibDB::Instance()->GetClusterParam();
-
-    //2. Load chain from debug streamers
-    //
-    //e.g
-    gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
-    gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
-    AliXRDPROOFtoolkit tool;  
-    TChain * chainres = tool.MakeChain("tracks.txt","ResolCl",0,10200);
-    chainres->Lookup();
-    //3. Do fits and store results
-    // 
-    AliTPCcalibTracks::MakeQPosNormAll(chainres,param,200000,0) ;
-    TFile f("paramout.root","recreate");
-    param->Write("clusterParam");
-    f.Close();
-    //4. Verification
-    TFile f2("paramout.root");
-    AliTPCClusterParam *param2 = (AliTPCClusterParam*)f2.Get("clusterParam");
-    param2->SetInstance(param2);
-    chainres->Draw("fitZ0:AliTPCClusterParam::SPosCorrection(1,0,Cl.fPad,Cl.fTimeBin,Cl.fZ,Cl.fSigmaY2,Cl.fSigmaZ2,Cl.fMax)","Cl.fDetector<36","",10000); // should be line 
-    
-   */
-
-
-  TStatToolkit toolkit;
-  Double_t chi2;
-  TVectorD fitParamY0;
-  TVectorD fitParamY1;
-  TVectorD fitParamZ0;
-  TVectorD fitParamZ1;
-  TMatrixD covMatrix;
-  Int_t npoints;
-  
-  chainres->SetAlias("dp","(-1+(Cl.fZ>0)*2)*((Cl.fPad-int(Cl.fPad))-0.5)");
-  chainres->SetAlias("dt","(-1+(Cl.fZ>0)*2)*((Cl.fTimeBin-0.66-int(Cl.fTimeBin-0.66))-0.5)");
-  chainres->SetAlias("sp","(sin(dp*pi)-dp*pi)");
-  chainres->SetAlias("st","(sin(dt)-dt)");
-  //
-  chainres->SetAlias("di","sqrt(1.-abs(Cl.fZ/250.))");
-  //
-  //
-  //
-  TCut cutA("1");
-  TString fstringY="";  
-  //
-  fstringY+="(dp)++";            //1
-  fstringY+="(dp)*di++";         //2
-  fstringY+="(sp)++";            //3
-  fstringY+="(sp)*di++";         //4
-  TString fstringZ="";  
-  fstringZ+="(dt)++";            //1
-  fstringZ+="(dt)*di++";         //2
-  fstringZ+="(st)++";            //3
-  fstringZ+="(st)*di++";         //4
-  //
-  // Z corrections
-  //
-  TString *strZ0 = toolkit.FitPlane(chainres,"(Cl.fZ-PZ0.fElements[0]):CSigmaZ0",fstringZ.Data(), "Cl.fDetector<36"+cutA, chi2,npoints,fitParamZ0,covMatrix,-1,0,maxPoints);
-  printf("Z0 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
-  param->PosZcor(0) = (TVectorD*) fitParamZ0.Clone();
-  //
-  TString *strZ1 = toolkit.FitPlane(chainres,"(Cl.fZ-PZ0.fElements[0]):CSigmaZ0",fstringZ.Data(), "Cl.fDetector>36"+cutA, chi2,npoints,fitParamZ1,covMatrix,-1,0,maxPoints);
-  //
-  printf("Z1 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
-  param->PosZcor(1) = (TVectorD*) fitParamZ1.Clone();
-  param->PosZcor(2) = (TVectorD*) fitParamZ1.Clone();
-  //
-  // Y corrections
-  //   
-  TString *strY0 = toolkit.FitPlane(chainres,"(Cl.fY-PY0.fElements[0]):CSigmaY0",fstringY.Data(), "Cl.fDetector<36"+cutA, chi2,npoints,fitParamY0,covMatrix,-1,0,maxPoints);
-  printf("Y0 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
-  param->PosYcor(0) = (TVectorD*) fitParamY0.Clone();
-  
-
-  TString *strY1 = toolkit.FitPlane(chainres,"(Cl.fY-PY0.fElements[0]):CSigmaY0",fstringY.Data(), "Cl.fDetector>36"+cutA, chi2,npoints,fitParamY1,covMatrix,-1,0,maxPoints);
-  //
-  printf("Y1 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
-  param->PosYcor(1) = (TVectorD*) fitParamY1.Clone();
-  param->PosYcor(2) = (TVectorD*) fitParamY1.Clone();
-  //
-  //
-  //
-  chainres->SetAlias("fitZ0",strZ0->Data());
-  chainres->SetAlias("fitZ1",strZ1->Data());
-  chainres->SetAlias("fitY0",strY0->Data());
-  chainres->SetAlias("fitY1",strY1->Data());
-  //  chainres->Draw("Cl.fZ-PZ0.fElements[0]","CSigmaY0<0.7&&CSigmaZ0<0.7"+cutA,"",10000);   
-}
-
 
 
 void AliTPCcalibTracks::MakeHistos(){
@@ -2703,46 +1373,61 @@ void AliTPCcalibTracks::MakeHistos(){
   TString axisName[8];
   
   //
-  binsTrack[0]=100;
+  binsTrack[0]=200;
   axisName[0]  ="var";
 
   binsTrack[1] =3;
   xminTrack[1] =0; xmaxTrack[1]=3;
   axisName[1]  ="pad type";
   //
-  binsTrack[2] =10;
-  xminTrack[2] =0; xmaxTrack[2]=1;
-  axisName[2]  ="drift length";
+  binsTrack[2] =20;
+  xminTrack[2] =-250; xmaxTrack[2]=250;
+  axisName[2]  ="z";
   //
   binsTrack[3] =10;
   xminTrack[3] =1; xmaxTrack[3]=400;
   axisName[3]  ="Qmax";
   //
-  binsTrack[4] =10;
+  binsTrack[4] =20;
   xminTrack[4] =0; xmaxTrack[4]=1;
   axisName[4]  ="cog";
   //
-  binsTrack[5] =10;
-  xminTrack[5] =0; xmaxTrack[5]=2;
-  axisName[5]  ="tan(phi)";
+  binsTrack[5] =15;
+  xminTrack[5] =-1.5; xmaxTrack[5]=1.5;
+  axisName[5]  ="tan(angle)";
   //
-  binsTrack[6] =10;
-  xminTrack[6] =0; xmaxTrack[6]=2;
-  axisName[6]  ="tan(theta)";
   //
-  xminTrack[0] =-0.5; xmaxTrack[0]=0.5;
-  fHisDeltaY=new THnSparseS("#Delta_{y} (cm)","#Delta_{y} (cm)", 7, binsTrack,xminTrack, xmaxTrack);
-  xminTrack[0] =-0.5; xmaxTrack[0]=0.5;
-  fHisDeltaZ=new THnSparseS("#Delta_{z} (cm)","#Delta_{z} (cm)", 7, binsTrack,xminTrack, xmaxTrack);
+  xminTrack[0] =-1.5; xmaxTrack[0]=1.5;
+  fHisDeltaY=new THnSparseF("#Delta_{y} (cm)","#Delta_{y} (cm)", 6, binsTrack,xminTrack, xmaxTrack);
+  xminTrack[0] =-1.5; xmaxTrack[0]=1.5;
+  fHisDeltaZ=new THnSparseF("#Delta_{z} (cm)","#Delta_{z} (cm)", 6, binsTrack,xminTrack, xmaxTrack);
   xminTrack[0] =0.; xmaxTrack[0]=0.5;
-  fHisRMSY=new THnSparseS("#RMS_{y} (cm)","#RMS_{y} (cm)", 7, binsTrack,xminTrack, xmaxTrack);
+  fHisRMSY=new THnSparseF("#RMS_{y} (cm)","#RMS_{y} (cm)", 6, binsTrack,xminTrack, xmaxTrack);
   xminTrack[0] =0.; xmaxTrack[0]=0.5;
-  fHisRMSZ=new THnSparseS("#RMS_{z} (cm)","#RMS_{z} (cm)", 7, binsTrack,xminTrack, xmaxTrack);
+  fHisRMSZ=new THnSparseF("#RMS_{z} (cm)","#RMS_{z} (cm)", 6, binsTrack,xminTrack, xmaxTrack);
   xminTrack[0] =0.; xmaxTrack[0]=100;
-  fHisQmax=new THnSparseS("Qmax (ADC)","Qmax (ADC)", 7, binsTrack,xminTrack, xmaxTrack);
+  fHisQmax=new THnSparseF("Qmax (ADC)","Qmax (ADC)", 6, binsTrack,xminTrack, xmaxTrack);
 
   xminTrack[0] =0.; xmaxTrack[0]=250;
-  fHisQtot=new THnSparseS("Qtot (ADC)","Qtot (ADC)", 7, binsTrack,xminTrack, xmaxTrack);
+  fHisQtot=new THnSparseF("Qtot (ADC)","Qtot (ADC)", 6, binsTrack,xminTrack, xmaxTrack);
+
+
+  for (Int_t ivar=0;ivar<6;ivar++){
+    fHisDeltaY->GetAxis(ivar)->SetName(axisName[ivar].Data());
+    fHisDeltaZ->GetAxis(ivar)->SetName(axisName[ivar].Data());
+    fHisRMSY->GetAxis(ivar)->SetName(axisName[ivar].Data());
+    fHisRMSZ->GetAxis(ivar)->SetName(axisName[ivar].Data());
+    fHisQmax->GetAxis(ivar)->SetName(axisName[ivar].Data());
+    fHisQtot->GetAxis(ivar)->SetName(axisName[ivar].Data());
+    fHisDeltaY->GetAxis(ivar)->SetTitle(axisName[ivar].Data());
+    fHisDeltaZ->GetAxis(ivar)->SetName(axisName[ivar].Data());
+    fHisRMSY->GetAxis(ivar)->SetName(axisName[ivar].Data());
+    fHisRMSZ->GetAxis(ivar)->SetName(axisName[ivar].Data());
+    fHisQmax->GetAxis(ivar)->SetName(axisName[ivar].Data());
+    fHisQtot->GetAxis(ivar)->SetName(axisName[ivar].Data());
+  }
+
+
   BinLogX(fHisDeltaY,3);
   BinLogX(fHisDeltaZ,3);
   BinLogX(fHisRMSY,3);
@@ -2761,3 +1446,128 @@ void    AliTPCcalibTracks::AddHistos(AliTPCcalibTracks* calib){
   if (calib->fHisRMSY)   fHisRMSY->Add(calib->fHisRMSY);
   if (calib->fHisRMSZ)   fHisRMSZ->Add(calib->fHisRMSZ);
 }
+
+
+
+void AliTPCcalibTracks::MakeSummaryTree(THnSparse *hisInput, TTreeSRedirector *pcstream, Int_t ptype){
+  //
+  // Dump summary info
+  //
+  //  0.OBJ: TAxis     var
+  //  1.OBJ: TAxis     pad type
+  //  2.OBJ: TAxis     z
+  //  3.OBJ: TAxis     Qmax
+  //  4.OBJ: TAxis     cog
+  //  5.OBJ: TAxis     tan(angle)
+  //
+  if (ptype>3) return;
+  Int_t idim[6]={0,1,2,3,4,5};
+  TString hname[4]={"dy","dz","rmsy","rmsz"};
+  //
+  Int_t nbins5=hisInput->GetAxis(5)->GetNbins();
+  Int_t first5=hisInput->GetAxis(5)->GetFirst();
+  Int_t last5 =hisInput->GetAxis(5)->GetLast();
+  //
+  for (Int_t ibin5=first5-1; ibin5<=last5; ibin5+=1){   // axis 5 - local angle
+    Bool_t bin5All=kTRUE;
+    if (ibin5>=first5){
+      hisInput->GetAxis(5)->SetRange(TMath::Max(ibin5-1,first5),TMath::Min(ibin5+1,nbins5));
+      if (ibin5==first5) hisInput->GetAxis(5)->SetRange(TMath::Max(ibin5,first5),TMath::Min(ibin5,nbins5));
+      bin5All=kFALSE;
+    }
+    Double_t      x5= hisInput->GetAxis(5)->GetBinCenter(ibin5);
+    THnSparse * his5= hisInput->Projection(5,idim);         //projected histogram according selection 5    
+    Int_t nbins4=his5->GetAxis(4)->GetNbins();
+    Int_t first4=his5->GetAxis(4)->GetFirst();
+    Int_t last4 =his5->GetAxis(4)->GetLast();
+    
+    for (Int_t ibin4=first4-1; ibin4<=last4; ibin4+=1){   // axis 4 - cog
+      Bool_t bin4All=kTRUE;
+      if (ibin4>=first4){
+       his5->GetAxis(4)->SetRange(TMath::Max(ibin4+1,first4),TMath::Min(ibin4-1,nbins4));
+       if (ibin4==first4||ibin4==last4) his5->GetAxis(4)->SetRange(TMath::Max(ibin4,first4),TMath::Min(ibin4,nbins4));
+       bin4All=kFALSE;
+      }
+      Double_t      x4= his5->GetAxis(4)->GetBinCenter(ibin4);
+      THnSparse * his4= his5->Projection(4,idim);         //projected histogram according selection 4
+      //
+      Int_t nbins3=his4->GetAxis(3)->GetNbins();
+      Int_t first3=his4->GetAxis(3)->GetFirst();
+      Int_t last3 =his4->GetAxis(3)->GetLast();
+      //
+      for (Int_t ibin3=first3-1; ibin3<=last3; ibin3+=1){   // axis 3 - Qmax
+       Bool_t bin3All=kTRUE;
+       if (ibin3>=first3){
+         his4->GetAxis(3)->SetRange(TMath::Max(ibin3,first3),TMath::Min(ibin3,nbins3));
+         bin3All=kFALSE;
+       }
+       Double_t      x3= his4->GetAxis(3)->GetBinCenter(ibin3);
+       THnSparse * his3= his4->Projection(3,idim);         //projected histogram according selection 3
+       //
+       Int_t nbins2    = his3->GetAxis(2)->GetNbins();
+       Int_t first2    = his3->GetAxis(2)->GetFirst();
+       Int_t last2     = his3->GetAxis(2)->GetLast();
+       //
+       for (Int_t ibin2=first2-1; ibin2<=last2; ibin2+=1){   // axis 2 - z     
+         Bool_t bin2All=kTRUE;
+         Double_t      x2= his3->GetAxis(2)->GetBinCenter(ibin2);
+         if (ibin2>=first2){
+           his3->GetAxis(2)->SetRange(TMath::Max(ibin2-1,first2),TMath::Min(ibin2+1,nbins2));
+           if (ibin2==first2||ibin2==last2||TMath::Abs(x2)<20) his3->GetAxis(2)->SetRange(TMath::Max(ibin2,first2),TMath::Min(ibin2,nbins2));
+           bin2All=kFALSE;
+         }
+         THnSparse * his2= his3->Projection(2,idim);         //projected histogram according selection 2
+         //
+         Int_t nbins1     = his2->GetAxis(1)->GetNbins();
+         Int_t first1     = his2->GetAxis(1)->GetFirst();
+         Int_t last1      = his2->GetAxis(1)->GetLast();
+         for (Int_t ibin1=first1-1; ibin1<=last1; ibin1++){   //axis 1 - pad type 
+           Bool_t bin1All=kTRUE;
+           if (ibin1>=first1){
+             his2->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1,nbins1));           
+             bin1All=kFALSE;
+           }
+           Double_t       x1= TMath::Nint(his2->GetAxis(1)->GetBinCenter(ibin1)-0.5);
+           TH1 * hisDelta = his2->Projection(0);
+           Double_t entries = hisDelta->GetEntries();
+           Double_t mean=0, rms=0;
+           if (entries>10){
+             mean    = hisDelta->GetMean();
+             rms = hisDelta->GetRMS();   
+             hisDelta->GetXaxis()->SetRangeUser(mean-4*rms,mean+4*rms);
+             mean    = hisDelta->GetMean();
+             rms = hisDelta->GetRMS();
+           }
+           //
+           (*pcstream)<<hname[ptype].Data()<<
+             // flag - intgrated values for given bin
+             "angleA="<<bin5All<<
+             "cogA="<<bin4All<<
+             "qmaxA="<<bin3All<<
+             "zA="<<bin2All<<
+             "ipadA="<<bin1All<<
+             // center of bin value
+             "angle="<<x5<<
+             "cog="<<x4<<
+             "qmax="<<x3<<
+             "z="<<x2<<
+             "ipad="<<x1<<
+             // mean values
+             //
+             "entries="<<entries<<
+             "mean="<<mean<<
+             "rms="<<rms<<
+             "ptype="<<ptype<<   //
+             "\n";
+           delete hisDelta;
+           printf("%f\t%f\t%f\t%f\t%f\t%f\t%f\n",x5,x4,x3,x2,x1, entries,mean);          
+         }//loop z
+         delete his2;
+       } // loop Q max
+       delete his3;
+      } // loop COG
+      delete his4;
+    }//loop local angle
+    delete his5;
+  }
+}
index 7c7cab4a8ee4b10f59c1174f7a9241382cc35ce2..d940c7185d8bff7b92e5cdffa509e853b2758244 100644 (file)
@@ -5,7 +5,7 @@
 ///////////////////////////////////////////////////////////////////////////////
 //                                                                           //
 //     Class to analyse tracks for calibration                               //
-//     to be used as a component in AliTPCSelectorTracks                     //
+//    to be used as a component in AliTPCSelectorTracks                     //
 //     In the constructor you have to specify name and title                 //
 //     to get the Object out of a file.                                      //
 //     The parameter 'clusterParam', a AliTPCClusterParam object             //
@@ -36,7 +36,6 @@ class AliTPCseed;
 class AliESDtrack;
 class AliTPCclusterMI;
 class AliTPCcalibTracksCuts;
-class AliTPCCalPadRegion;
 class AliTPCCalPad;
 class TChain;
 class TTree;
@@ -60,49 +59,29 @@ public :
   virtual Long64_t Merge(TCollection *li);
   void    AddHistos(AliTPCcalibTracks* calib);
   void     MakeResPlotsQTree(Int_t minEntries = 100, const char* pathName = "plots");
-  static void MakeQPosNormAll(TTree * chain, AliTPCClusterParam * param, Int_t maxPoints=1000000);
    void     MakeReport(Int_t stat, const char* pathName = "plots");     // calls all functions that procude pictures, results are written to pathName, stat is the minimal statistic threshold
-  //
-
-   
+  //   
   Int_t           AcceptTrack(AliTPCseed * track);
   void            FillResolutionHistoLocal(AliTPCseed * track);  // the MAIN-FUNCTION, called for each track to fill the histograms, called by Process(...)
-  static  TH2D*   MakeDiff(TH2D * hfit, TF2 * func);
   
-  static AliTPCcalibTracks* TestMerge(AliTPCcalibTracks *ct, AliTPCClusterParam *clusterParam, Int_t nCalTracks = 50);
   
   void     SetStyle() const;
-  void     Draw(Option_t* opt);                                  // draw some exemplaric histograms for fast result-check
-  void     MakeAmpPlots(Int_t stat, const char* pathName = "plots");
-  void     MakeDeltaPlots(const char* pathName = "plots");
-  void     MakeChargeVsDriftLengthPlotsOld(const char* pathName = "plots");
-  void     MakeChargeVsDriftLengthPlots(const char* pathName = "plots");
-  void     FitResolutionNew(const char* pathName = "plots");
-  void     FitRMSNew(const char* pathName = "plots");
   
-  TObjArray* GetfArrayAmpRow() const {return fArrayAmpRow;}
-  TObjArray* GetfArrayAmp() const {return fArrayAmp;}
   TObjArray* GetfArrayQDY() const {return fArrayQDY;}
   TObjArray* GetfArrayQDZ() const {return fArrayQDZ;}
   TObjArray* GetfArrayQRMSY() const {return fArrayQRMSY;}
   TObjArray* GetfArrayQRMSZ() const {return fArrayQRMSZ;}
-  TObjArray* GetfArrayChargeVsDriftlength() const {return fArrayChargeVsDriftlength;}
-  TH1F*      GetfDeltaY() const {return fDeltaY;}
-  TH1F*      GetfDeltaZ() const {return fDeltaZ;}
   TObjArray* GetfResolY() const {return fResolY;}
   TObjArray* GetfResolZ() const {return fResolZ;}
   TObjArray* GetfRMSY() const {return fRMSY;}
   TObjArray* GetfRMSZ() const {return fRMSZ;}
-  TH1I*      GetfHclus() const {return fHclus;}
   TH1I*      GetfRejectedTracksHisto() const {return fRejectedTracksHisto;}
-  TH1I*      GetfHclusterPerPadrow() const {return fHclusterPerPadrow;}
-  TH1I*      GetfHclusterPerPadrowRaw() const {return fHclusterPerPadrowRaw;}
   TH2I*      GetfClusterCutHisto() const  {return fClusterCutHisto;}
   AliTPCCalPad*          GetfCalPadClusterPerPad() const {return fCalPadClusterPerPad; }
   AliTPCCalPad*          GetfCalPadClusterPerPadRaw() const {return fCalPadClusterPerPadRaw;}
-  AliTPCCalPadRegion*    GetCalPadRegionchargeVsDriftlength() const {return fcalPadRegionChargeVsDriftlength;}
   AliTPCcalibTracksCuts* GetCuts() {return fCuts;}
   void MakeHistos();  //make THnSparse
+  static void MakeSummaryTree(THnSparse *hisInput, TTreeSRedirector *pcstream, Int_t ptype);
 protected:         
   
 private:
@@ -111,9 +90,9 @@ private:
    static Int_t   GetBin(Int_t  iq, Int_t pad);
    static Float_t GetQ(Int_t bin);
    static Float_t GetPad(Int_t bin);
-   void FillResolutionHistoLocalDebugPart(AliTPCseed *track, AliTPCclusterMI *cluster0, Int_t irow, Float_t  angley, Float_t  anglez, Int_t nclFound, Int_t kDelta);
    AliTPCClusterParam *fClusterParam; // pointer to cluster parameterization
    AliTPCROC *fROC;          //!
+public:
   THnSparse  *fHisDeltaY;    // THnSparse - delta Y 
   THnSparse  *fHisDeltaZ;    // THnSparse - delta Z 
   THnSparse  *fHisRMSY;      // THnSparse - rms Y 
@@ -121,31 +100,22 @@ private:
   THnSparse  *fHisQmax;      // THnSparse - qmax 
   THnSparse  *fHisQtot;      // THnSparse - qtot 
 
-
-   TObjArray *fArrayAmpRow; // array with amplitudes versus row for given sector 
-   TObjArray *fArrayAmp;    // array with amplitude for sectors
+private:
    TObjArray *fArrayQDY;    // q binned delta Y histograms
    TObjArray *fArrayQDZ;    // q binned delta Z histograms 
    TObjArray *fArrayQRMSY;  // q binned delta Y histograms
    TObjArray *fArrayQRMSZ;  // q binned delta Z histograms 
-   TObjArray *fArrayChargeVsDriftlength; // array of arrays of TProfiles with charge vs. driftlength for each padsize and sector
-   AliTPCCalPadRegion *fcalPadRegionChargeVsDriftlength;  // CalPadRegion, one TProfile for charge vs. driftlength for each padsize and sector
 
-   TH1F      *fDeltaY;      // integrated delta y histo
-   TH1F      *fDeltaZ;      // integrated delta z histo
    TObjArray *fResolY;      // array of resolution histograms Y
    TObjArray *fResolZ;      // array of resolution histograms Z
    TObjArray *fRMSY;        // array of RMS histograms Y
    TObjArray *fRMSZ;        // array of RMS histograms Z
    AliTPCcalibTracksCuts *fCuts; // object with cuts, that is passed to the constructor
-   TH1I      *fHclus;       // number of clusters per track
    TH1I      *fRejectedTracksHisto; // histogram of rejecteced tracks, the number coresponds to the failed cut
-   TH1I      *fHclusterPerPadrow;   // histogram showing the number of clusters per padRow
-   TH1I      *fHclusterPerPadrowRaw;// histogram showing the number of clusters per padRow before cuts on clusters are applied
    TH2I      *fClusterCutHisto;     // histogram showing in which padRow the clusters were cutted by which criterium
    AliTPCCalPad *fCalPadClusterPerPad;    // AliTPCCalPad showing the number of clusters per Pad
    AliTPCCalPad *fCalPadClusterPerPadRaw; // AliTPCCalPad showing the number of clusters per Pad before cuts on clusters are applied
-   ClassDef(AliTPCcalibTracks,1)
+   ClassDef(AliTPCcalibTracks,2)
    
 };
 
index 96c0ca4f9c5a16a98cc3448c6c428357948e99a5..4808ea45ba651ef734db44444a4e9e62a627271b 100644 (file)
 #include "TCint.h"
 #include "AliTPCcalibV0.h"
 #include "AliV0.h"
+#include "TRandom.h"
+#include "TTreeStream.h"
+#include "AliTPCcalibDB.h"
+#include "AliTPCCorrection.h"
+#include "AliGRPObject.h"
+#include "AliTPCTransform.h"
 
 
 
@@ -48,33 +54,40 @@ ClassImp(AliTPCcalibV0)
 
 AliTPCcalibV0::AliTPCcalibV0() : 
    AliTPCcalibBase(),
+   fV0Tree(0),
+   fHPTTree(0),
    fStack(0),
    fESD(0),
    fPdg(0),
    fParticles(0),
    fV0s(0),
-   fGammas(0),
-   fV0type(0),
-   fTPCdEdx(0),
-   fTPCdEdxPi(0),
-   fTPCdEdxEl(0),
-   fTPCdEdxP(0)
+   fGammas(0)
+{
+  
+}   
+AliTPCcalibV0::AliTPCcalibV0(const Text_t *name, const Text_t *title):
+   AliTPCcalibBase(),
+   fV0Tree(0),
+   fHPTTree(0),
+   fStack(0),
+   fESD(0),
+   fPdg(0),
+   fParticles(0),
+   fV0s(0),
+   fGammas(0)
 {
   fPdg = new TDatabasePDG;       
   // create output histograms
-  fTPCdEdx   = new TH2F("TPCdEdX",  "dE/dX; BetaGamma; TPC signal (a.u.)", 1000, 0.1, 10000, 300, 0, 300);
-  BinLogX(fTPCdEdx); 
-  fTPCdEdxPi = new TH2F("TPCdEdXPi","dE/dX; BetaGamma; TPC signal (a.u.)", 1000, 0.1, 10000, 300, 0, 300);
-  fTPCdEdxEl = new TH2F("TPCdEdXEl","dE/dX; BetaGamma; TPC signal (a.u.)", 1000, 0.1, 10000, 300, 0, 300);  
-  fTPCdEdxP  = new TH2F("TPCdEdXP", "dE/dX; BetaGamma; TPC signal (a.u.)", 1000, 0.1, 10000, 300, 0, 300);        
-  
-  
+  SetName(name);
+  SetTitle(title);
 }   
 
 AliTPCcalibV0::~AliTPCcalibV0(){
   //
   //
   //
+  delete fV0Tree;
+  delete fHPTTree;
 }
 
 
@@ -87,6 +100,11 @@ void  AliTPCcalibV0::ProcessESD(AliESDEvent *esd, AliStack *stack){
   //
   fESD = esd;
   AliKFParticle::SetField(esd->GetMagneticField());
+  if (TMath::Abs(AliTracker::GetBz())<1) return;  
+  DumpToTree(esd);
+  DumpToTreeHPT(esd);
+  return;
+  //
   MakeV0s();
   if (stack) {
     fStack = stack;
@@ -96,6 +114,304 @@ void  AliTPCcalibV0::ProcessESD(AliESDEvent *esd, AliStack *stack){
   }
 }
 
+void  AliTPCcalibV0::DumpToTreeHPT(AliESDEvent *esd){
+  //
+  // Dump V0s fith full firend information to the 
+  // 
+  if (TMath::Abs(AliTracker::GetBz())<1) return;
+  const Int_t kMinCluster=110;
+  const Float_t kMinPt   =3.;
+  AliESDfriend *esdFriend=static_cast<AliESDfriend*>(esd->FindListObject("AliESDfriend"));
+  if (!esdFriend) {
+    Printf("ERROR: esdFriend not available");
+    return;
+  }
+  //
+  Int_t ntracks=esd->GetNumberOfTracks();
+  for (Int_t i=0;i<ntracks;++i) {
+    Bool_t isOK=kFALSE;
+    AliESDtrack *track = esd->GetTrack(i);
+    if (track->GetTPCncls()<kMinCluster) continue;
+    AliESDfriendTrack *friendTrack = esdFriend->GetTrack(i);
+    if (!friendTrack) continue;
+    if (TMath::Abs(AliTracker::GetBz())>1){ // cut on momenta if measured
+      if (track->Pt()>kMinPt) isOK=kTRUE;
+    }
+    if (TMath::Abs(AliTracker::GetBz())<1){  // require primary track for the B field OFF data
+      Bool_t isAccepted=kTRUE;
+      if (!track->IsOn(AliESDtrack::kITSrefit)) isAccepted=kFALSE;
+      if (!track->IsOn(AliESDtrack::kTPCrefit)) isAccepted=kFALSE;
+      if (!track->IsOn(AliESDtrack::kTOFout))   isAccepted=kFALSE;
+      Float_t dvertex[2],cvertex[3]; 
+      track->GetImpactParametersTPC(dvertex,cvertex);
+      if (TMath::Abs(dvertex[0]/TMath::Sqrt(cvertex[0]+0.01))>20) isAccepted=kFALSE;
+      if (TMath::Abs(dvertex[1]/TMath::Sqrt(TMath::Abs(cvertex[2]+0.01)))>20) isAccepted=kFALSE;
+      track->GetImpactParameters(dvertex,cvertex);
+      if (TMath::Abs(dvertex[0]/TMath::Sqrt(cvertex[0]+0.01))>10) isAccepted=kFALSE;
+      if (TMath::Abs(dvertex[1]/TMath::Sqrt(TMath::Abs(cvertex[2]+0.01)))>10) isAccepted=kFALSE;
+      if (!isAccepted) isOK=kFALSE;
+    }
+    
+    if (!isOK) continue;
+    //
+
+    TObject *calibObject;
+    AliTPCseed *seed = 0;
+    for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
+      if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
+    }
+    if (!seed) continue;
+      if (!fHPTTree) {
+      fHPTTree = new TTree("HPT","HPT");
+      fHPTTree->SetDirectory(0);
+    }
+    if (fHPTTree->GetEntries()==0){
+      //
+      fHPTTree->SetDirectory(0);
+      fHPTTree->Branch("t.",&track);
+      fHPTTree->Branch("ft.",&friendTrack);
+      fHPTTree->Branch("s.",&seed);
+    }else{
+      fHPTTree->SetBranchAddress("t.",&track);
+      fHPTTree->SetBranchAddress("ft.",&friendTrack);
+      fHPTTree->SetBranchAddress("s.",&seed);
+    }
+    fHPTTree->Fill();
+    //
+  }
+}
+
+
+
+void  AliTPCcalibV0::DumpToTree(AliESDEvent *esd){
+  //
+  // Dump V0s fith full firend information to the 
+  // 
+  Int_t nV0s  = fESD->GetNumberOfV0s();
+  const Int_t kMinCluster=110;
+  const Double_t kDownscale=0.01;
+  const Float_t kMinR    =0;
+  const Float_t kMinPt   =1.0;
+  const Float_t kMinMinPt   =0.7;
+  AliESDfriend *esdFriend=static_cast<AliESDfriend*>(esd->FindListObject("AliESDfriend"));
+  if (!esdFriend) {
+    Printf("ERROR: esdFriend not available");
+    return;
+  }
+  //
+  
+  for (Int_t ivertex=0; ivertex<nV0s; ivertex++){
+    Bool_t isOK=kFALSE;
+    AliESDv0 * v0 = (AliESDv0*) esd->GetV0(ivertex);
+    AliESDtrack * track0 = fESD->GetTrack(v0->GetIndex(0)); // negative track
+    AliESDtrack * track1 = fESD->GetTrack(v0->GetIndex(1)); // positive track 
+    if (track0->GetTPCNcls()<kMinCluster) continue;
+    if (track0->GetKinkIndex(0)>0) continue;    
+    if (track1->GetTPCNcls()<kMinCluster) continue;
+    if (track1->GetKinkIndex(0)>0) continue;
+    if (v0->GetOnFlyStatus()==kFALSE) continue;
+    //
+    if (TMath::Min(track0->Pt(),track1->Pt())<kMinMinPt) continue;
+    //
+    //
+    if (TMath::Max(track0->Pt(),track1->Pt())>kMinPt) isOK=kTRUE;
+    if (gRandom->Rndm()<kDownscale) isOK=kTRUE;  
+    if (!isOK) continue;
+    //
+    AliESDfriendTrack *ftrack0 = esdFriend->GetTrack(v0->GetIndex(0));
+    if (!ftrack0) continue;
+    AliESDfriendTrack *ftrack1 = esdFriend->GetTrack(v0->GetIndex(1));
+    if (!ftrack1) continue;
+    //
+    TObject *calibObject;
+    AliTPCseed *seed0 = 0;
+    AliTPCseed *seed1 = 0;
+    for (Int_t l=0;(calibObject=ftrack0->GetCalibObject(l));++l) {
+      if ((seed0=dynamic_cast<AliTPCseed*>(calibObject))) break;
+    }
+    for (Int_t l=0;(calibObject=ftrack1->GetCalibObject(l));++l) {
+      if ((seed1=dynamic_cast<AliTPCseed*>(calibObject))) break;
+    }
+    if (!seed0) continue;
+    if (!seed1) continue;
+    AliExternalTrackParam * paramIn0 = (AliExternalTrackParam *)track0->GetInnerParam();
+    AliExternalTrackParam * paramIn1 = (AliExternalTrackParam *)track1->GetInnerParam();
+    if (!paramIn0) continue;
+    if (!paramIn1) continue;
+    //
+    //
+    if (!fV0Tree) {
+      fV0Tree = new TTree("V0s","V0s");
+      fV0Tree->SetDirectory(0);
+    }
+    if (fV0Tree->GetEntries()==0){
+      //
+      fV0Tree->SetDirectory(0);
+      fV0Tree->Branch("v0.",&v0);
+      fV0Tree->Branch("t0.",&track0);
+      fV0Tree->Branch("t1.",&track1);
+      fV0Tree->Branch("ft0.",&ftrack0);
+      fV0Tree->Branch("ft1.",&ftrack1);
+      fV0Tree->Branch("s0.",&seed0);
+      fV0Tree->Branch("s1.",&seed1);
+    }else{
+      fV0Tree->SetBranchAddress("v0.",&v0);
+      fV0Tree->SetBranchAddress("t0.",&track0);
+      fV0Tree->SetBranchAddress("t1.",&track1);
+      fV0Tree->SetBranchAddress("ft0.",&ftrack0);
+      fV0Tree->SetBranchAddress("ft1.",&ftrack1);
+      fV0Tree->SetBranchAddress("s0.",&seed0);
+      fV0Tree->SetBranchAddress("s1.",&seed1);
+    }
+    fV0Tree->Fill();
+  }
+}
+
+
+Long64_t AliTPCcalibV0::Merge(TCollection *const li) {
+
+  TIterator* iter = li->MakeIterator();
+  AliTPCcalibV0* cal = 0;
+
+  while ((cal = (AliTPCcalibV0*)iter->Next())) {
+    if (cal->fV0Tree){
+      if (!fV0Tree) {
+       fV0Tree = new TTree("V0s","V0s");
+       fV0Tree->SetDirectory(0);
+      }
+      if (cal->fV0Tree->GetEntries()>0) AliTPCcalibV0::AddTree(cal->fV0Tree);
+      if (cal->fHPTTree->GetEntries()>0) AliTPCcalibV0::AddTreeHPT(cal->fHPTTree);
+    }    
+  }
+  return 0;
+}
+
+
+void AliTPCcalibV0::AddTree(TTree * treeInput){
+  //
+  // Add the content of tree: 
+  // Notice automatic copy of tree in ROOT does not work for such complicated tree
+  //  
+  AliESDv0 * v0 = new AliESDv0;
+  Double_t kMinPt=0.8;
+  AliESDtrack * track0 = 0; // negative track
+  AliESDtrack * track1 = 0; // positive track 
+  AliESDfriendTrack *ftrack0 = 0;
+  AliESDfriendTrack *ftrack1 = 0;
+  AliTPCseed *seed0 = 0;
+  AliTPCseed *seed1 = 0;
+  treeInput->SetBranchStatus("ft0.",kFALSE);
+  treeInput->SetBranchStatus("ft1.",kFALSE);
+  TDatabasePDG pdg;
+  Double_t massK0= pdg.GetParticle("K0")->Mass();
+  Double_t massLambda= pdg.GetParticle("Lambda0")->Mass();
+
+  Int_t entries= treeInput->GetEntries();
+  for (Int_t i=0; i<entries; i++){
+    treeInput->SetBranchAddress("v0.",&v0);
+    treeInput->SetBranchAddress("t0.",&track0);
+    treeInput->SetBranchAddress("t1.",&track1);
+    treeInput->SetBranchAddress("ft0.",&ftrack0);
+    treeInput->SetBranchAddress("ft1.",&ftrack1);
+    treeInput->SetBranchAddress("s0.",&seed0);
+    treeInput->SetBranchAddress("s1.",&seed1);
+    if (fV0Tree->GetEntries()==0){
+      fV0Tree->SetDirectory(0);
+      fV0Tree->Branch("v0.",&v0);
+      fV0Tree->Branch("t0.",&track0);
+      fV0Tree->Branch("t1.",&track1);
+      fV0Tree->Branch("ft0.",&ftrack0);
+      fV0Tree->Branch("ft1.",&ftrack1);
+      fV0Tree->Branch("s0.",&seed0);
+      fV0Tree->Branch("s1.",&seed1);
+    }else{
+      fV0Tree->SetBranchAddress("v0.",&v0);
+      fV0Tree->SetBranchAddress("t0.",&track0);
+      fV0Tree->SetBranchAddress("t1.",&track1);
+      fV0Tree->SetBranchAddress("ft0.",&ftrack0);
+      fV0Tree->SetBranchAddress("ft1.",&ftrack1);
+      fV0Tree->SetBranchAddress("s0.",&seed0);
+      fV0Tree->SetBranchAddress("s1.",&seed1);
+    }
+    //
+    treeInput->GetEntry(i);
+    ftrack0->GetCalibContainer()->SetOwner(kTRUE);
+    ftrack1->GetCalibContainer()->SetOwner(kTRUE);
+    Bool_t isOK=kTRUE;
+    if (v0->GetOnFlyStatus()==kFALSE) isOK=kFALSE;
+    if (track0->GetTPCncls()<100) isOK=kFALSE;
+    if (track1->GetTPCncls()<100) isOK=kFALSE;    
+    if (TMath::Min(seed0->Pt(),seed1->Pt())<kMinPt) isOK=kFALSE;
+    if (TMath::Min(track0->Pt(),track1->Pt())<kMinPt) isOK=kFALSE;
+    Bool_t isV0=kFALSE;    
+    if (TMath::Abs(v0->GetEffMass(2,2)-massK0)<0.05)     isV0=kTRUE;
+    if (TMath::Abs(v0->GetEffMass(4,2)-massLambda)<0.05) isV0=kTRUE; 
+    if (TMath::Abs(v0->GetEffMass(2,4)-massLambda)<0.05) isV0=kTRUE;
+    if (TMath::Abs(v0->GetEffMass(0,0))<0.02) isV0=kFALSE; //reject electrons
+    if (!isV0) isOK=kFALSE;
+    if (isOK) fV0Tree->Fill();
+    delete v0;
+    delete track0;
+    delete track1;
+    delete ftrack0;
+    delete ftrack1;
+    delete seed0;
+    delete seed1;
+    v0=0;
+    track0=0;
+    track1=0;
+    ftrack0=0;
+    ftrack1=0;
+    seed0=0;
+    seed1=0;
+  }
+}
+
+void AliTPCcalibV0::AddTreeHPT(TTree * treeInput){
+  //
+  // Add the content of tree: 
+  // Notice automatic copy of tree in ROOT does not work for such complicated tree
+  //  
+  AliESDtrack *track = 0;
+  AliESDfriendTrack *friendTrack = 0;
+  AliTPCseed *seed = 0;
+  if (!treeInput) return;
+  if (treeInput->GetEntries()==0) return;
+  //
+  Int_t entries= treeInput->GetEntries();  
+  //
+  for (Int_t i=0; i<entries; i++){
+    track=0;
+    friendTrack=0;
+    seed=0;
+    //
+    treeInput->SetBranchAddress("t.",&track);
+    treeInput->SetBranchAddress("ft.",&friendTrack);
+    treeInput->SetBranchAddress("s.",&seed);
+    treeInput->GetEntry(i);
+    //
+    if (fHPTTree->GetEntries()==0){
+      fHPTTree->SetDirectory(0);
+      fHPTTree->Branch("t.",&track);
+      fHPTTree->Branch("ft.",&friendTrack);
+      fHPTTree->Branch("s.",&seed);
+    }else{
+      fHPTTree->SetBranchAddress("t.",&track);
+      fHPTTree->SetBranchAddress("ft.",&friendTrack);
+      fHPTTree->SetBranchAddress("s.",&seed);
+    }    
+    Bool_t isOK=kTRUE;
+    if (!track->IsOn(AliESDtrack::kITSrefit)) isOK=kFALSE;
+    if (!track->IsOn(AliESDtrack::kTOFout)) isOK=kFALSE;
+    if (isOK) fHPTTree->Fill();
+    //
+    delete track;
+    delete friendTrack;
+    delete seed;
+  }
+}
+
+
 void AliTPCcalibV0::MakeMC(){
   //
   // MC comparison 
@@ -257,11 +573,430 @@ void AliTPCcalibV0::MakeMC(){
   }
 }
 
+void AliTPCcalibV0::MakeFitTreeTrack(const TObjArray * corrArray, Double_t ptCut, Int_t run){
+  //
+  // Make a fit tree
+  //
+  // 0. Loop over selected tracks
+  // 1. Loop over all transformation - refit the track with and without the
+  //    transformtation
+  // 2. Dump the matching paramaeters to the debugStremer
+  //
+  
+  //Connect input
+  const Int_t kMinNcl=120;
+  TFile f("TPCV0Objects.root");
+  AliTPCcalibV0 *v0TPC = (AliTPCcalibV0*) f.Get("v0TPC");
+  TTree * treeInput = v0TPC->GetHPTTree();
+  TTreeSRedirector *pcstream = new TTreeSRedirector("fitHPT.root");
+  AliESDtrack *track = 0;
+  AliESDfriendTrack *friendTrack = 0;
+  AliTPCseed *seed = 0;
+  if (!treeInput) return;
+  if (treeInput->GetEntries()==0) return;
+  //
+  treeInput->SetBranchAddress("t.",&track);
+  treeInput->SetBranchAddress("ft.",&friendTrack);
+  treeInput->SetBranchAddress("s.",&seed);
+  //
+  Int_t ncorr=0;
+  if (corrArray) ncorr = corrArray->GetEntries();
+  AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
+  AliTPCParam     *param     = AliTPCcalibDB::Instance()->GetParameters();
+  AliGRPObject*  grp = AliTPCcalibDB::Instance()->GetGRP(run);
+  Double_t time=0.5*(grp->GetTimeStart() +grp->GetTimeEnd());
+  //
+  //
+  //  
+  Int_t ntracks= treeInput->GetEntries();
+  for (Int_t itrack=0; itrack<ntracks; itrack++){
+    treeInput->GetEntry(itrack);
+    if (!track) continue;
+    if (seed->Pt()<ptCut) continue;
+    if (track->Pt()<ptCut) continue;
+    if (track->GetTPCncls()<kMinNcl) continue;
+    //
+    // Reapply transformation
+    //
+    for (Int_t irow=0; irow<159; irow++){
+      AliTPCclusterMI *cluster=seed->GetClusterPointer(irow);
+      if (cluster &&cluster->GetX()>10){
+        Double_t x0[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
+        Int_t index0[1]={cluster->GetDetector()};
+        transform->Transform(x0,index0,0,1);
+        cluster->SetX(x0[0]);
+        cluster->SetY(x0[1]);
+        cluster->SetZ(x0[2]);
+        //
+      }
+    }    
+    //
+    Double_t xref=134;
+    AliExternalTrackParam* paramInner=0;
+    AliExternalTrackParam* paramOuter=0;
+    AliExternalTrackParam* paramIO=0;
+    Bool_t isOK=kTRUE;
+    for (Int_t icorr=-1; icorr<ncorr; icorr++){
+      //
+      AliTPCCorrection *corr = 0;
+      if (icorr>=0) corr = (AliTPCCorrection*)corrArray->At(icorr);
+      AliExternalTrackParam * trackInner = RefitTrack(seed, corr,85,134,0.1);      
+      AliExternalTrackParam * trackIO = RefitTrack(seed, corr,245,85,0.1);      
+      AliExternalTrackParam * trackOuter = RefitTrack(seed, corr,245,134,0.1 ); 
+      if (trackInner&&trackOuter&&trackIO){
+       trackOuter->Rotate(trackInner->GetAlpha());
+       trackOuter->PropagateTo(trackInner->GetX(),AliTracker::GetBz());
+       if (icorr<0) {
+         paramInner=trackInner;
+         paramOuter=trackOuter;
+         paramIO=trackIO;
+         paramIO->Rotate(seed->GetAlpha());
+         paramIO->PropagateTo(seed->GetX(),AliTracker::GetBz());
+       }
+      }else{
+       isOK=kFALSE;
+      }
+      
+    }
+    if (paramOuter&& paramInner) {
+      //      Bool_t isOK=kTRUE;
+      if (paramInner->GetSigmaY2()>0.01) isOK&=kFALSE;
+      if (paramOuter->GetSigmaY2()>0.01) isOK&=kFALSE;
+      if (paramInner->GetSigmaZ2()>0.01) isOK&=kFALSE;
+      if (paramOuter->GetSigmaZ2()>0.01) isOK&=kFALSE;      
+      (*pcstream)<<"fit"<<
+       "s.="<<seed<<
+       "io.="<<paramIO<<
+       "pIn.="<<paramInner<<
+       "pOut.="<<paramOuter;      
+    }
+    //
+  }
+  delete pcstream;
+  /*
+    .x ~/rootlogon.C
+    Int_t run=117112;
+    .x ../ConfigCalibTrain.C(run)
+    .L ../AddTaskTPCCalib.C
+    ConfigOCDB(run)
+    TFile fexb("../../RegisterCorrectionExB.root");
+    AliTPCComposedCorrection *cc=  (AliTPCComposedCorrection*) fexb.Get("ComposedExB");
+    cc->Init();
+    cc->Print("DA"); // Print used correction classes
+    TObjArray *array = cc->GetCorrections()
+    AliTPCcalibV0::MakeFitTreeTrack(array,2,run);
+   
+   */
+}
 
+void AliTPCcalibV0::MakeFitTreeV0(const TObjArray * corrArray, Double_t ptCut, Int_t run){
+  //
+  // Make a fit tree
+  //
+  // 0. Loop over selected tracks
+  // 1. Loop over all transformation - refit the track with and without the
+  //    transformtation
+  // 2. Dump the matching paramaeters to the debugStremer
+  //
+  
+  //Connect input
+  const Int_t kMinNcl=120;
+  TFile f("TPCV0Objects.root");
+  AliTPCcalibV0 *v0TPC = (AliTPCcalibV0*) f.Get("v0TPC");
+  TTree * treeInput = v0TPC->GetV0Tree();
+  TTreeSRedirector *pcstream = new TTreeSRedirector("fitV0.root");
+  AliESDv0 *v0 = 0;
+  AliESDtrack *track0 = 0;
+  AliESDfriendTrack *friendTrack0 = 0;
+  AliTPCseed *seed0 = 0;
+  AliTPCseed *s0 = 0;
+  AliESDtrack *track1 = 0;
+  AliESDfriendTrack *friendTrack1 = 0;
+  AliTPCseed *s1 = 0;
+  AliTPCseed *seed1 = 0;
+  if (!treeInput) return;
+  if (treeInput->GetEntries()==0) return;
+  //
+  treeInput->SetBranchAddress("v0.",&v0);
+  treeInput->SetBranchAddress("t0.",&track0);
+  treeInput->SetBranchAddress("ft0.",&friendTrack0);
+  treeInput->SetBranchAddress("s0.",&s0);
+  treeInput->SetBranchAddress("t1.",&track1);
+  treeInput->SetBranchAddress("ft1.",&friendTrack1);
+  treeInput->SetBranchAddress("s1.",&s1);
+  //
+  TDatabasePDG pdg;
+  Int_t ncorr=0;
+  if (corrArray) ncorr = corrArray->GetEntries();
+  AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
+  AliTPCParam     *param     = AliTPCcalibDB::Instance()->GetParameters();
+  AliGRPObject*  grp = AliTPCcalibDB::Instance()->GetGRP(run);
+  Double_t time=0.5*(grp->GetTimeStart() +grp->GetTimeEnd());
+  Double_t massK0= pdg.GetParticle("K0")->Mass();
+  Double_t massLambda= pdg.GetParticle("Lambda0")->Mass();
+  Double_t massPion=pdg.GetParticle("pi+")->Mass();
+  Double_t massProton=pdg.GetParticle("proton")->Mass();
+  Double_t pdgPion=pdg.GetParticle("pi+")->PdgCode();
+  Double_t pdgProton=pdg.GetParticle("proton")->PdgCode();
+  Double_t mass0=0;
+  Double_t mass1=0;
+  Double_t massV0=0;
+  Int_t    pdg0=0;
+  Int_t    pdg1=0;
+  //
+  //
+  //  
+  Int_t nv0s= treeInput->GetEntries();
+  for (Int_t iv0=0; iv0<nv0s; iv0++){
+    Int_t  v0Type=0;
+    Int_t isK0=0;
+    Int_t isLambda=0;
+    Int_t isAntiLambda=0;
+    treeInput->GetEntry(iv0);
+    if (TMath::Abs(v0->GetEffMass(2,2)-massK0)<0.03) {isK0=1; v0Type=1;} //select K0s    
+    if (TMath::Abs(v0->GetEffMass(4,2)-massLambda)<0.01) {isLambda=1; v0Type=2;} //select Lambda   
+    if (TMath::Abs(v0->GetEffMass(2,4)-massLambda)<0.01) {isAntiLambda=1;v0Type=3;} //select Anti Lambda
+    if (isK0+isLambda+isAntiLambda!=1) continue;
+    mass0=massPion;
+    mass1=massPion;
+    pdg0=pdgPion;
+    pdg1=pdgPion;
+    if (isLambda) {mass0=massProton; pdg0=pdgProton;}
+    if (isAntiLambda) {mass1=massProton; pdg1=pdgProton;}
+    massV0=massK0;
+    if (isK0==0) massV0=massLambda;
+    //
+    if (!s0) continue;
+    seed0=(s0->GetSign()>0)?s0:s1;
+    seed1=(s0->GetSign()>0)?s1:s0;
+    if (seed0->GetZ()*seed1->GetZ()<0) continue; //remove membrane crossed tracks
+    if (seed0->Pt()<ptCut) continue;
+    if (seed1->Pt()<ptCut) continue;
+    //
+    // Reapply transformation
+    //
+    for  (Int_t itype=0; itype<2; itype++){
+      AliTPCseed * seed = (itype==0) ? seed0: seed1;      
+      for (Int_t irow=0; irow<159; irow++){
+       AliTPCclusterMI *cluster=seed->GetClusterPointer(irow);
+       if (cluster &&cluster->GetX()>10){
+         Double_t x0[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
+         Int_t index0[1]={cluster->GetDetector()};
+         transform->Transform(x0,index0,0,1);
+         cluster->SetX(x0[0]);
+         cluster->SetY(x0[1]);
+         cluster->SetZ(x0[2]);
+         //
+       }
+      }
+    }   
+    Bool_t isOK=kTRUE;
+    Double_t radius = v0->GetRr();
+    Double_t xyz[3];
+    v0->GetXYZ(xyz[0],xyz[1],xyz[2]);
+    Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
+    TObjArray arrayV0in(ncorr+1);
+    TObjArray arrayV0io(ncorr+1);
+    TObjArray arrayT0(ncorr+1);
+    TObjArray arrayT1(ncorr+1);
+    arrayV0in.SetOwner(kTRUE);
+    arrayV0io.SetOwner(kTRUE);
+    //
+    for (Int_t icorr=-1; icorr<ncorr; icorr++){
+      AliTPCCorrection *corr =0;
+      if (icorr>=0) corr = (AliTPCCorrection*)corrArray->At(icorr);
+      //
+      AliExternalTrackParam * trackInner0 = RefitTrack(seed0, corr,160,85,mass0);      
+      AliExternalTrackParam * trackIO0    = RefitTrack(seed0, corr,245,85,mass0);      
+      AliExternalTrackParam * trackInner1 = RefitTrack(seed1, corr,160,85,mass1);      
+      AliExternalTrackParam * trackIO1    = RefitTrack(seed1, corr,245,85,mass1);      
+      if (!trackInner0) isOK=kFALSE;
+      if (!trackInner1) isOK=kFALSE;
+      if (!trackIO0)    isOK=kFALSE;
+      if (!trackIO1)    isOK=kFALSE;
+      if (isOK){
+       if (!trackInner0->Rotate(alpha)) isOK=kFALSE;
+       if (!trackInner1->Rotate(alpha)) isOK=kFALSE;
+       if (!trackIO0->Rotate(alpha)) isOK=kFALSE;
+       if (!trackIO1->Rotate(alpha)) isOK=kFALSE;
+       //
+       if (!AliTracker::PropagateTrackToBxByBz(trackInner0, radius, mass0, 1, kFALSE)) isOK=kFALSE; 
+       if (!AliTracker::PropagateTrackToBxByBz(trackInner1, radius, mass1, 1, kFALSE)) isOK=kFALSE; 
+       if (!AliTracker::PropagateTrackToBxByBz(trackIO0, radius, mass0, 1, kFALSE)) isOK=kFALSE; 
+       if (!AliTracker::PropagateTrackToBxByBz(trackIO1, radius, mass1, 1, kFALSE)) isOK=kFALSE; 
+       if (!isOK) continue;
+       arrayT0.AddAt(trackIO0->Clone(),icorr+1);
+       arrayT1.AddAt(trackIO1->Clone(),icorr+1);
+       Int_t charge=(trackIO0->GetSign());
+       AliKFParticle pin0( *trackInner0,  pdg0*charge);
+       AliKFParticle pin1( *trackInner1, -pdg1*charge);
+       AliKFParticle pio0( *trackIO0,  pdg0*charge);
+       AliKFParticle pio1( *trackIO1, -pdg1*charge);
+       AliKFParticle v0in;
+       AliKFParticle v0io;
+       v0in+=pin0;
+       v0in+=pin1;
+       v0io+=pio0;
+       v0io+=pio1;
+       arrayV0in.AddAt(v0in.Clone(),icorr+1);
+       arrayV0io.AddAt(v0io.Clone(),icorr+1);
+      }
+    }
+    if (!isOK) continue;
+    //
+    //AliKFParticle* pin0= (AliKFParticle*)arrayV0in.At(0);
+    AliKFParticle* pio0= (AliKFParticle*)arrayV0io.At(0);
+    AliExternalTrackParam *param0=(AliExternalTrackParam *)arrayT0.At(0);
+    AliExternalTrackParam *param1=(AliExternalTrackParam *)arrayT1.At(0);
+    Double_t mass0=0, mass0E=0; 
+    pio0->GetMass( mass0,mass0E);
+    //
+    Double_t mean=mass0-massV0;
+    if (TMath::Abs(mean)>0.05) continue;
+    Double_t mass[10000];
+    //
+    Int_t dtype=30;  // id for V0
+    Int_t ptype=5;   // id for invariant mass
+    //    Int_t id=TMath::Nint(100.*(param0->Pt()-param1->Pt())/(param0->Pt()+param1->Pt()));      // K0s V0 asymetry
+    Int_t id=1000.*(param0->Pt()-param1->Pt());      // K0s V0 asymetry
+    Double_t gx,gy,gz, px,py,pz;
+    Double_t pt = v0->Pt();
+    v0->GetXYZ(gx,gy,gz);
+    v0->GetPxPyPz(px,py,pz);
+    Double_t theta=pz/TMath::Sqrt(px*px+py*py);
+    Double_t phi=TMath::ATan2(py,px);
+    Double_t snp=0.5*(seed0->GetSnp()+seed1->GetSnp());
+    Double_t sector=9*phi/TMath::Pi();
+    Double_t dsec=sector-TMath::Nint(sector);
+    Double_t refX=TMath::Sqrt(gx*gx+gy*gy);
+    //Int_t nentries=v0Type;
+    Double_t bz=AliTracker::GetBz();
+    Double_t dRrec=0;
+    (*pcstream)<<"fitDebug"<< 
+      "id="<<id<<
+      "v0.="<<v0<<
+      "mean="<<mean<<
+      "rms="<<mass0E<<
+      "pio0.="<<pio0<<
+      "p0.="<<param0<<
+      "p1.="<<param1;
+    (*pcstream)<<"fit"<<  // dump valus for fit
+      "run="<<run<<       //run number
+      "bz="<<bz<<         // magnetic filed used
+      "dtype="<<dtype<<   // detector match type 30
+      "ptype="<<ptype<<   // parameter type
+      "theta="<<theta<<   // theta
+      "phi="<<phi<<       // phi 
+      "snp="<<snp<<       // snp
+      "mean="<<mean<<     // mean dist value
+      "rms="<<mass0E<<       // rms
+      "sector="<<sector<<
+      "dsec="<<dsec<<
+      //
+      "refX="<<refX<<      // reference radius
+      "gx="<<gx<<         // global position
+      "gy="<<gy<<         // global position
+      "gz="<<gz<<         // global position
+      "dRrec="<<dRrec<<      // delta Radius in reconstruction
+      "pt="<<pt<<         // pt of the particle
+      "id="<<id<<     //delta of the momenta      
+      "entries="<<v0Type;//  type of the V0
+    for (Int_t icorr=0; icorr<ncorr; icorr++){
+      AliTPCCorrection *corr =0;
+      if (icorr>=0) corr = (AliTPCCorrection*)corrArray->At(icorr);
+      AliKFParticle* pin= (AliKFParticle*)arrayV0in.At(icorr+1);
+      AliKFParticle* pio= (AliKFParticle*)arrayV0io.At(icorr+1);
+      AliExternalTrackParam *par0=(AliExternalTrackParam *)arrayT0.At(icorr+1);
+      AliExternalTrackParam *par1=(AliExternalTrackParam *)arrayT1.At(icorr+1);
+      Double_t massE=0; 
+      pio->GetMass( mass[icorr],massE);
+      mass[icorr]-=mass0;
+      (*pcstream)<<"fit"<< 
+       Form("%s=",corr->GetName())<<mass[icorr];
+      (*pcstream)<<"fitDebug"<< 
+       Form("%s=",corr->GetName())<<mass[icorr]<<
+       Form("%sp0.=",corr->GetName())<<par0<<
+       Form("%sp1=",corr->GetName())<<par1;
+    }
+    (*pcstream)<<"fit"<< "isOK="<<isOK<<"\n";        
+    (*pcstream)<<"fitDebug"<< "isOK="<<isOK<<"\n";        
+  }
+  delete pcstream;
+}
 
-void AliTPCcalibV0::MakeV0s(){
+AliExternalTrackParam * AliTPCcalibV0::RefitTrack(AliTPCseed *seed, AliTPCCorrection * corr, Double_t xstart, Double_t xstop, Double_t mass){
   //
+  // Refit the track:
+  //    seed   - tpc track with cluster
+  //    corr   - distrotion/correction class  - apllied to the points
+  //    xstart - radius to start propagate/update
+  //    xstop  - radius to stop propagate/update
+  // 
+  const Double_t kResetCov=20.;
+  const Double_t kSigma=5.;
+  Double_t covar[15];
+  for (Int_t i=0;i<15;i++) covar[i]=0;
+  covar[0]=kSigma*kSigma;
+  covar[2]=kSigma*kSigma;
+  covar[5]=kSigma*kSigma/Float_t(150*150);
+  covar[9]=kSigma*kSigma/Float_t(150*150);
+  covar[14]=1*1;
+  // 
+  AliExternalTrackParam *refit  = new AliExternalTrackParam(*seed);
+  refit->PropagateTo(xstart, AliTracker::GetBz());
+  refit->AddCovariance(covar);
+  refit->ResetCovariance(kResetCov);
+  Double_t xmin = TMath::Min(xstart,xstop);
+  Double_t xmax = TMath::Max(xstart,xstop);
+  Int_t ncl=0;
   //
+  Bool_t isOK=kTRUE;
+  for (Int_t index0=0; index0<159; index0++){
+    Int_t irow= (xstart<xstop)? index0:159-index0;
+    AliTPCclusterMI *cluster=seed->GetClusterPointer(irow);  //cluster in local system
+    if (!cluster) continue;
+    if (cluster->GetX()<xmin) continue;
+    if (cluster->GetX()>xmax) continue;
+    Double_t alpha = TMath::Pi()*(cluster->GetDetector()+0.5)/9.;
+    if (!refit->Rotate(alpha)) isOK=kFALSE;
+    Double_t x     = cluster->GetX();
+    Double_t y     = cluster->GetY();
+    Double_t z     = cluster->GetZ();
+    if (corr){      
+      Float_t xyz[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};  // original position
+      corr->DistortPointLocal(xyz,cluster->GetDetector());
+      x=xyz[0];
+      y=xyz[1];
+      z=xyz[2];
+    }
+    if (!AliTracker::PropagateTrackToBxByBz(refit, x,mass,1.,kFALSE)) isOK=kFALSE;
+    if (!isOK) continue;
+    Double_t cov[3]={0.01,0.,0.01};
+    Double_t yz[2]={y,z};
+    if (!refit->Update(yz,cov)) isOK=kFALSE;    
+    ncl++;
+  }
+  if (!AliTracker::PropagateTrackToBxByBz(refit, xstop, mass,1.,kTRUE)) isOK=kFALSE;
+  //  
+  if (!isOK) {
+    delete refit;
+    return 0;
+  }  
+  return refit;
+}
+
+
+
+
+
+//
+// Obsolete part
+//
+
+
+void AliTPCcalibV0::MakeV0s(){
+  //
+  // 
   //
   const Int_t kMinCluster=40;
   const Float_t kMinR    =0;
@@ -333,166 +1068,12 @@ void AliTPCcalibV0::MakeV0s(){
 
 
 
-// void AliTPCcalibV0::ProcessV0(Int_t ftype){
-//   //
-//   //
-//   const Double_t ktimeK0     = 2.684;
-//   const Double_t ktimeLambda = 7.89; 
-  
-  
-//   if (! fGammas) fGammas = new TObjArray(10);
-//   fGammas->Clear();
-//   Int_t nV0s  = fV0s->GetEntries();
-//   if (nV0s==0) return;
-//   AliKFVertex primVtx(*(fESD->GetPrimaryVertex()));
-//   //
-//   for (Int_t ivertex=0; ivertex<nV0s; ivertex++){
-//     AliESDv0 * v0 = (AliESDv0*)fV0s->At(ivertex);
-//     AliESDtrack * trackN = fESD->GetTrack(v0->GetIndex(0));
-//     AliESDtrack * trackP = fESD->GetTrack(v0->GetIndex(1));
-//     // 
-//     // 
-//     //
-//     AliKFParticle *v0K0       = Fit(primVtx,v0,211,211);
-//     AliKFParticle *v0Gamma    = Fit(primVtx,v0,11,-11);
-//     AliKFParticle *v0Lambda42 = Fit(primVtx,v0,2212,211);
-//     AliKFParticle *v0Lambda24 = Fit(primVtx,v0,211,2212);
-//     //Set production vertex
-//     v0K0->SetProductionVertex( primVtx );
-//     v0Gamma->SetProductionVertex( primVtx );
-//     v0Lambda42->SetProductionVertex( primVtx );
-//     v0Lambda24->SetProductionVertex( primVtx );
-//     Double_t massK0, massGamma, massLambda42,massLambda24, massSigma;
-//     v0K0->GetMass( massK0,massSigma);
-//     v0Gamma->GetMass( massGamma,massSigma);
-//     v0Lambda42->GetMass( massLambda42,massSigma);
-//     v0Lambda24->GetMass( massLambda24,massSigma);
-//     Float_t chi2K0       = v0K0->GetChi2()/v0K0->GetNDF();
-//     Float_t chi2Gamma    = v0Gamma->GetChi2()/v0Gamma->GetNDF();
-//     Float_t chi2Lambda42 = v0Lambda42->GetChi2()/v0Lambda42->GetNDF();
-//     Float_t chi2Lambda24 = v0Lambda24->GetChi2()/v0Lambda24->GetNDF();
-//     //
-//     // Mass Contrained params
-//     //
-//     AliKFParticle *v0K0C       = Fit(primVtx,v0,211,211);
-//     AliKFParticle *v0GammaC    = Fit(primVtx,v0,11,-11);
-//     AliKFParticle *v0Lambda42C = Fit(primVtx,v0,2212,211);
-//     AliKFParticle *v0Lambda24C = Fit(primVtx,v0,211,2212);
-//     //   
-//     v0K0C->SetProductionVertex( primVtx );
-//     v0GammaC->SetProductionVertex( primVtx );
-//     v0Lambda42C->SetProductionVertex( primVtx );
-//     v0Lambda24C->SetProductionVertex( primVtx );
-
-//     v0K0C->SetMassConstraint(fPdg->GetParticle(310)->Mass());
-//     v0GammaC->SetMassConstraint(0);
-//     v0Lambda42C->SetMassConstraint(fPdg->GetParticle(3122)->Mass());
-//     v0Lambda24C->SetMassConstraint(fPdg->GetParticle(-3122)->Mass());
-//     //    
-//     Double_t timeK0, sigmaTimeK0;  
-//     Double_t timeLambda42, sigmaTimeLambda42;  
-//     Double_t timeLambda24, sigmaTimeLambda24;  
-//     v0K0C->GetLifeTime(timeK0, sigmaTimeK0);
-//     //v0K0Gamma->GetLifeTime(timeK0, sigmaTimeK0);
-//     v0Lambda42C->GetLifeTime(timeLambda42, sigmaTimeLambda42);
-//     v0Lambda24C->GetLifeTime(timeLambda24, sigmaTimeLambda24);
-    
-
-//     //
-//     Float_t chi2K0C       = v0K0C->GetChi2()/v0K0C->GetNDF();
-//     if (chi2K0C<0) chi2K0C=100;
-//     Float_t chi2GammaC    = v0GammaC->GetChi2()/v0GammaC->GetNDF();
-//     if (chi2GammaC<0) chi2GammaC=100;
-//     Float_t chi2Lambda42C = v0Lambda42C->GetChi2()/v0Lambda42C->GetNDF();
-//     if (chi2Lambda42C<0) chi2Lambda42C=100;
-//     Float_t chi2Lambda24C = v0Lambda24C->GetChi2()/v0Lambda24C->GetNDF();
-//     if (chi2Lambda24C<0) chi2Lambda24C=100;
-//     //
-//     Float_t  minChi2C=99;
-//     Int_t   type   =-1;
-//     if (chi2K0C<minChi2C) { minChi2C= chi2K0C; type=0;}
-//     if (chi2GammaC<minChi2C) { minChi2C= chi2GammaC; type=1;}
-//     if (chi2Lambda42C<minChi2C) { minChi2C= chi2Lambda42C; type=2;}
-//     if (chi2Lambda24C<minChi2C) { minChi2C= chi2Lambda24C; type=3;}
-//     Float_t  minChi2=99;
-//     Int_t   type0   =-1;
-//     if (chi2K0<minChi2) { minChi2= chi2K0; type0=0;}
-//     if (chi2Gamma<minChi2) { minChi2= chi2Gamma; type0=1;}
-//     if (chi2Lambda42<minChi2) { minChi2= chi2Lambda42; type0=2;}
-//     if (chi2Lambda24<minChi2) { minChi2= chi2Lambda24; type0=3;}
-//     Float_t betaGammaP  = trackN->GetP()/fPdg->GetParticle(-2212)->Mass(); 
-//     Float_t betaGammaPi = trackN->GetP()/fPdg->GetParticle(-211)->Mass();
-//     Float_t betaGammaEl = trackN->GetP()/fPdg->GetParticle(11)->Mass();
-//     Float_t dedxTeorP = BetheBlochAleph(betaGammaP);
-//     Float_t dedxTeorPi = BetheBlochAleph(betaGammaPi);;
-//     Float_t dedxTeorEl = BetheBlochAleph(betaGammaEl);;
-//     //
-//     //
-//     if (minChi2>50) continue;
-//     (*fDebugStream)<<"V0"<<
-//       "ftype="<<ftype<<
-//       "v0.="<<v0<<
-//       "trackN.="<<trackN<<
-//       "trackP.="<<trackP<<
-//       //
-//       "dedxTeorP="<<dedxTeorP<<
-//       "dedxTeorPi="<<dedxTeorPi<<
-//       "dedxTeorEl="<<dedxTeorEl<<
-//       //
-//       "type="<<type<<
-//       "chi2C="<<minChi2C<<
-//       "v0K0.="<<v0K0<<
-//       "v0Gamma.="<<v0Gamma<<
-//       "v0Lambda42.="<<v0Lambda42<<
-//       "v0Lambda24.="<<v0Lambda24<<
-//       //
-//       "chi20K0.="<<chi2K0<<
-//       "chi2Gamma.="<<chi2Gamma<<
-//       "chi2Lambda42.="<<chi2Lambda42<<
-//       "chi2Lambda24.="<<chi2Lambda24<<
-//       //
-//       "chi20K0c.="<<chi2K0C<<
-//       "chi2Gammac.="<<chi2GammaC<<
-//       "chi2Lambda42c.="<<chi2Lambda42C<<
-//       "chi2Lambda24c.="<<chi2Lambda24C<<
-//       //
-//       "v0K0C.="<<v0K0C<<
-//       "v0GammaC.="<<v0GammaC<<
-//       "v0Lambda42C.="<<v0Lambda42C<<
-//       "v0Lambda24C.="<<v0Lambda24C<<
-//       //
-//       "massK0="<<massK0<<
-//       "massGamma="<<massGamma<<
-//       "massLambda42="<<massLambda42<<
-//       "massLambda24="<<massLambda24<<
-//       //
-//       "timeK0="<<timeK0<<
-//       "timeLambda42="<<timeLambda42<<
-//       "timeLambda24="<<timeLambda24<<
-//       "\n";
-//     if (type==1) fGammas->AddLast(v0); 
-//     //
-//     //
-//     //
-//     delete v0K0;
-//     delete v0Gamma;
-//     delete v0Lambda42;
-//     delete v0Lambda24;    
-//     delete v0K0C;
-//     delete v0GammaC;
-//     delete v0Lambda42C;
-//     delete v0Lambda24C;    
-//   }
-//   ProcessPI0(); 
-// }
-
-
 
 
 void AliTPCcalibV0::ProcessV0(Int_t ftype){
   //
-  //
-  
+  // Obsolete
+  //  
   if (! fGammas) fGammas = new TObjArray(10);
   fGammas->Clear();
   Int_t nV0s  = fV0s->GetEntries();
@@ -613,10 +1194,6 @@ void AliTPCcalibV0::ProcessV0(Int_t ftype){
      if (chi2Lambda42C < 10*minChi2C) numCand++;
      if (chi2Lambda24C < 10*minChi2C) numCand++;
      //
-     if (numCand < 2) {
-      if (paramInNeg->GetP() > 0.4) fTPCdEdx->Fill(betaGamma0, trackN->GetTPCsignal());
-      if (paramInPos->GetP() > 0.4) fTPCdEdx->Fill(betaGamma1, trackP->GetTPCsignal());
-     }    
     }
     
     //
@@ -750,13 +1327,6 @@ AliKFParticle * AliTPCcalibV0::Fit(AliKFVertex & /*primVtx*/, AliESDv0 *v0, Int_
   return V0;  
 }
 
-TH2F * AliTPCcalibV0::GetHistograms() {
-  //
-  //
-  //
- return fTPCdEdx;
-}
-
 
 
 
index 0bbba452dfa00ae0ac61f251dd0a781f50228daa..c4628c179a97d234b20a62eda61eb5fbdb2d9f39 100644 (file)
@@ -21,6 +21,8 @@ class AliESDv0;
 class TArrayI;
 class TTree;
 class AliStack;
+class  AliExternalTrackParam;
+class  AliTPCCorrection;
 
 class AliTPCcalibV0 : public AliTPCcalibBase {
 public :
@@ -28,47 +30,53 @@ public :
    // List of branches
 
   AliTPCcalibV0();
+  AliTPCcalibV0(const Text_t *name, const Text_t *title);
   virtual ~AliTPCcalibV0();
   virtual void     Process(AliESDEvent *event) {return ProcessESD(event,0);}
-
+  Long64_t Merge(TCollection *const li);
+  void AddTree(TTree * treeInput);
+  void AddTreeHPT(TTree * treeInput);
+  static AliExternalTrackParam * RefitTrack(AliTPCseed *seed, AliTPCCorrection * corr, Double_t xstart, Double_t xstop, Double_t mass);
+  static void MakeFitTreeTrack(const TObjArray * corrArray, Double_t ptCut, Int_t run);
+  static void MakeFitTreeV0(const TObjArray * corrArray, Double_t ptCut, Int_t run);
   //
   //
   //
   void ProcessESD(AliESDEvent *esd, AliStack *stack=0);
+  void DumpToTree(AliESDEvent *esd);
+  void DumpToTreeHPT(AliESDEvent *esd);
+  TTree * GetV0Tree(){return fV0Tree;}
+  TTree * GetHPTTree(){return fHPTTree;}
+  //  
   void MakeMC();
   void MakeV0s();
   void ProcessV0(Int_t ftype);
   void ProcessPI0();
-  TH2F * GetHistograms();
   void BinLogX(TH2F * h);
   //
   //
   //  
   static AliKFParticle * Fit(AliKFVertex & primVtx, AliESDv0 *v0, Int_t PDG1, Int_t PDG2);
-  void     Process(AliESDtrack *track, Int_t runNo=-1){AliTPCcalibBase::Process(track,runNo);};
-  void     Process(AliTPCseed *track){return AliTPCcalibBase::Process(track);}
 
 protected:
 private:
 
-   AliTPCcalibV0(const AliTPCcalibV0&); // Not implemented
-   AliTPCcalibV0& operator=(const AliTPCcalibV0&); // Not implemented
-
-
-   AliStack       *fStack;        // pointer to kinematic tree        
-   AliESDEvent    *fESD;              //! current ED to proccess - NOT OWNER
-   TDatabasePDG   *fPdg;              // particle database
-   TObjArray      *fParticles;         // array of selected MC particles
-   TObjArray      *fV0s;               // array of V0s
-   TObjArray      *fGammas;           // gamma conversion candidates
-   //
-   TArrayI        *fV0type;            // array of types for V0s       
-   TH2F           *fTPCdEdx;              // dEdx spectra
-   TH2F           *fTPCdEdxPi;            // dEdx spectra - pion anti-pion
-   TH2F           *fTPCdEdxEl;            // dEdx spectra - electroen -positrons from gamma
-   TH2F           *fTPCdEdxP;             // dEdx spectra - proton antiproton - lambda -  antilambda
-   //       
-   ClassDef(AliTPCcalibV0,1);
+  AliTPCcalibV0(const AliTPCcalibV0&); // Not implemented
+  AliTPCcalibV0& operator=(const AliTPCcalibV0&); // Not implemented
+  TTree          *fV0Tree;      // tree with full V0 information
+  TTree          *fHPTTree;      // tree with high mometa tracks - full calib info
+  //
+  AliStack       *fStack;        // pointer to kinematic tree        
+  AliESDEvent    *fESD;              //! current ED to proccess - NOT OWNER
+  TDatabasePDG   *fPdg;              //! particle database
+  TObjArray      *fParticles;         // array of selected MC particles
+  TObjArray      *fV0s;               // array of V0s
+  TObjArray      *fGammas;           // gamma conversion candidates
+  //
+  void     Process(AliESDtrack *track, Int_t runNo=-1){AliTPCcalibBase::Process(track,runNo);};
+  void     Process(AliTPCseed *track){return AliTPCcalibBase::Process(track);}  
+  //       
+  ClassDef(AliTPCcalibV0,3);
 };
 
 
index 513afa63e7a3dccd5ac7732fa107d774f48a0831..e6049f05b6b90a5f1ba375756247858b4c52adc9 100644 (file)
@@ -969,3 +969,111 @@ void AliTPCkalmanAlign::MakeAliasCE(TTree * chain){
   chain->SetAlias("La","(ly.fElements/lx.fElements/0.155)");
   chain->SetAlias("deltaT","(CETmean.fElements-PulserTmean.fElements)");
 }
+
+
+
+
+void AliTPCkalmanAlign::Update1D(Double_t delta, Double_t sigma, Int_t s1, TMatrixD &vecXk, TMatrixD &covXk){
+  //
+  // Update parameters and covariance - with one measurement
+  //
+  const Int_t knMeas=1;
+  Int_t knElem=vecXk.GetNrows();
+  TMatrixD mat1(knElem,knElem);            // update covariance matrix
+  TMatrixD matHk(1,knElem);        // vector to mesurement
+  TMatrixD vecYk(knMeas,1);        // Innovation or measurement residual
+  TMatrixD matHkT(knElem,knMeas);  // helper matrix Hk transpose
+  TMatrixD matSk(knMeas,knMeas);   // Innovation (or residual) covariance
+  TMatrixD matKk(knElem,knMeas);   // Optimal Kalman gain
+  TMatrixD covXk2(knElem,knElem);  // helper matrix
+  TMatrixD covXk3(knElem,knElem);  // helper matrix
+  TMatrixD vecZk(1,1);
+  TMatrixD measR(1,1);
+  vecZk(0,0)=delta;
+  measR(0,0)=sigma*sigma;
+  //
+  // reset matHk
+  for (Int_t iel=0;iel<knElem;iel++) 
+    for (Int_t ip=0;ip<knMeas;ip++) matHk(ip,iel)=0; 
+  //mat1
+  for (Int_t iel=0;iel<knElem;iel++) {
+    for (Int_t jel=0;jel<knElem;jel++) mat1(iel,jel)=0;
+    mat1(iel,iel)=1;
+  }
+  //
+  matHk(0, s1)=1;
+  vecYk = vecZk-matHk*vecXk;               // Innovation or measurement residual
+  matHkT=matHk.T(); matHk.T();
+  matSk = (matHk*(covXk*matHkT))+measR;    // Innovation (or residual) covariance
+  matSk.Invert();
+  matKk = (covXk*matHkT)*matSk;            //  Optimal Kalman gain
+  vecXk += matKk*vecYk;                    //  updated vector 
+  covXk2= (mat1-(matKk*matHk));
+  covXk3 =  covXk2*covXk;          
+  covXk = covXk3;  
+  Int_t nrows=covXk3.GetNrows();
+  
+  for (Int_t irow=0; irow<nrows; irow++)
+    for (Int_t icol=0; icol<nrows; icol++){
+      // rounding problems - make matrix again symteric
+      covXk(irow,icol)=(covXk3(irow,icol)+covXk3(icol,irow))*0.5; 
+    }
+}
+
+
+void   AliTPCkalmanAlign::Update1D(TString &input, TString filter, TVectorD &param, TMatrixD & covar, Double_t mean, Double_t sigma){
+  //
+  //  Update Parameters
+  //  input  - variable name description
+  //  filter - filter string 
+  //  param  - parameter vector
+  //  covar  - covariance
+  //  mean   - value to set
+  //  sigma  - sigma of measurement 
+  TObjArray *array0= input.Tokenize("++");
+  TObjArray *array1= filter.Tokenize("++");
+  TMatrixD paramM(param.GetNrows(),1);
+  for (Int_t i=0; i<array0->GetEntries(); i++){paramM(i,0)=param(i);}
+
+  for (Int_t i=0; i<array0->GetEntries(); i++){
+    Bool_t isOK=kTRUE;
+    TString str(array0->At(i)->GetName());
+    for (Int_t j=0; j<array1->GetEntries(); j++){
+      if (str.Contains(array1->At(j)->GetName())==0) isOK=kFALSE;
+    }
+    if (isOK) {
+      AliTPCkalmanAlign::Update1D(mean, sigma, i, paramM, covar);
+    }
+  }
+  for (Int_t i=0; i<array0->GetEntries(); i++){
+    param(i)=paramM(i,0);
+  }
+}
+
+
+TString  AliTPCkalmanAlign::FilterFit(TString &input, TString filter, TVectorD &param, TMatrixD & covar){
+  //
+  //
+  //
+  TObjArray *array0= input.Tokenize("++");
+  TObjArray *array1= filter.Tokenize("++");
+  //TString *presult=new TString("(0");
+  TString result="(0.0";
+  for (Int_t i=0; i<array0->GetEntries(); i++){
+    Bool_t isOK=kTRUE;
+    TString str(array0->At(i)->GetName());
+    for (Int_t j=0; j<array1->GetEntries(); j++){
+      if (str.Contains(array1->At(j)->GetName())==0) isOK=kFALSE;
+    }
+    if (isOK) {
+      result+="+"+str;
+      result+=Form("*(%f)",param[i+1]);
+      printf("%f\t%f\t%s\n",param[i+1], TMath::Sqrt(covar(i+1,i+1)),str.Data());
+    }
+  }
+  result+="-0.)";
+  return result;
+}
+
+
index 09e54fcba321fb97228997b9c689a8abfc483905..50a1d1564c65ae771e1d2934f62f55ff656e8935 100644 (file)
@@ -6,6 +6,7 @@
 
 #include "TNamed.h"
 #include "TMatrixD.h"
+#include "TString.h"
 class TTreeSRedirector;
 class TObjArray;
 class AliTPCcalibAlign;
@@ -24,6 +25,11 @@ public:
   void UpdateOCDBTime0( AliTPCCalPad  *pad, Int_t ustartRun, Int_t uendRun,  const char* storagePath );
   static void UpdateAlign1D(Double_t delta, Double_t sigma, Int_t s1, Int_t s2, TMatrixD &param, TMatrixD &covar);
   static void UpdateAlign1D(Double_t delta, Double_t sigma, Int_t s1, TMatrixD &param, TMatrixD &covar);
+  //
+  static void Update1D(Double_t delta, Double_t sigma, Int_t s1, TMatrixD &param, TMatrixD &covar);
+  static void Update1D(TString &input, TString filter, TVectorD &param, TMatrixD & covar, Double_t mean, Double_t sigma);
+  static TString  FilterFit(TString &input, TString filter, TVectorD &param, TMatrixD & covar);
+  //
   static void BookAlign1D(TMatrixD &param, TMatrixD &covar, Double_t sigma, Double_t mean);
   void DumpOldAlignment(TTreeSRedirector *pcstream);
   void MakeNewAlignment(Bool_t add,TTreeSRedirector *pcstream=0);