// -*- mode: c++; mode: auto-fill; mode: flyspell-prog; -*-

// Author: A.F.Zarnecki, University of Warsaw <mailto:zarnecki@fuw.edu.pl>
// Version: $Id: EUTelTestFitter.cc,v 1.39 2009-07-30 14:00:19 zarnecki Exp $
// Date 2007.06.04

/*
 *   This source code is part of the Eutelescope package of Marlin.
 *   You are free to use this source files for your own development as
 *   long as it stays in a public research context. You are not
 *   allowed to use it for commercial purpose. You must put this
 *   header with author names in all development based on this file.
 *
 */

#if defined USE_GEAR

// eutelescope includes:

#include "EUTelTestFitter.h"
#include "EUTELESCOPE.h"
#include "EUTelEventImpl.h"
#include "EUTelRunHeaderImpl.h"
#include "EUTelHistogramManager.h"
#include "EUTelExceptions.h"
#include "EUTelPStream.h" // process streams redi::ipstream

// GBL:

#include "../GeneralBrokenLines/trunk/cpp/GblTrajectory.h"

// AIDA histogram package (on top of ROOT):

#if defined(USE_AIDA) || defined(MARLIN_USE_AIDA)
#include <marlin/AIDAProcessor.h>
#include <AIDA/IHistogramFactory.h>
#include <AIDA/IHistogram1D.h>
#include <AIDA/IHistogram2D.h>
#include <AIDA/IProfile1D.h>
#include <AIDA/IProfile2D.h>
#include <AIDA/ITree.h>
#endif

// marlin includes ".h"
#include "marlin/Processor.h"
#include "marlin/Exceptions.h"
#include "marlin/Global.h"

// gear includes <.h>
#include <gear/GearMgr.h>
#include <gear/SiPlanesParameters.h>

// LCIO includes <.h>
#include <EVENT/LCCollection.h>
#include <EVENT/LCEvent.h>
#include <IMPL/LCCollectionVec.h>
#include <IMPL/TrackerHitImpl.h>
#include <IMPL/TrackImpl.h>
#include <IMPL/LCFlagImpl.h>
#include <Exceptions.h>

// system includes <>
#include <iostream>
#include <fstream>
#include <vector>
#include <cmath>
#include <sstream>
#include <memory>
#include <string>
#include <map>
#include <cstdlib>
#include <limits>

// ROOT includes ".h"
#include <TMath.h>
#include <TVectorD.h>
#include <TMatrixD.h>
#include <TVector3.h>
#include <TRotation.h>
#include "TH1D.h"

using namespace std; //must be before pixelForReadout.h
//AB#include "pixelForReadout.h" //struct pixel, struct cluster
#include "CMSpixelClust.h"

using namespace lcio;
using namespace marlin;
using namespace eutelescope;

// definition of static members mainly used to name histograms
#if defined(USE_AIDA) || defined(MARLIN_USE_AIDA)

AIDA::IHistogram1D * t100Histo;
AIDA::IHistogram1D * t300Histo;
AIDA::IHistogram1D * t1000Histo;
AIDA::IHistogram1D * t3600Histo;
AIDA::IHistogram1D * dtHisto;
AIDA::IHistogram1D * dtmsHisto;
AIDA::IHistogram1D * logdtHisto;
AIDA::IHistogram1D * logdtcmsHisto;

AIDA::IHistogram1D * cmscolHisto;
AIDA::IHistogram1D * cmsrowHisto;
AIDA::IHistogram1D * cmsnpxHisto;
AIDA::IHistogram1D * cmsadcHisto;

AIDA::IHistogram1D * refcolHisto;
AIDA::IHistogram1D * refrowHisto;
AIDA::IHistogram1D * refnpxHisto;
AIDA::IHistogram1D * refadcHisto;

AIDA::IHistogram1D * cmsncolHisto;
AIDA::IHistogram1D * cmsnrowHisto;
AIDA::IHistogram1D * cmsetaHisto;

AIDA::IHistogram1D * cmsdtHisto;
AIDA::IHistogram1D * sysrtHisto;
AIDA::IHistogram1D * sysrdtHisto;
AIDA::IHistogram1D * sysdtHisto;
AIDA::IProfile1D * sysdtvsdt;

AIDA::IHistogram1D * hits0Histo;
AIDA::IHistogram1D * hits1Histo;
AIDA::IHistogram1D * hits2Histo;
AIDA::IHistogram1D * hits3Histo;
AIDA::IHistogram1D * hits4Histo;
AIDA::IHistogram1D * hits5Histo;

AIDA::IHistogram1D * dx01Histo;
AIDA::IHistogram1D * dy01Histo;
AIDA::IHistogram1D * du01Histo;
AIDA::IHistogram1D * dx02Histo;
AIDA::IHistogram1D * dx03Histo;
AIDA::IHistogram1D * dx04Histo;
AIDA::IHistogram1D * dx05Histo;

AIDA::IHistogram1D * dx12Histo;
AIDA::IHistogram1D * dy12Histo;
AIDA::IHistogram1D * du12Histo;

AIDA::IHistogram1D * dx23Histo;
AIDA::IHistogram1D * dy23Histo;
AIDA::IHistogram1D * du23Histo;

AIDA::IHistogram1D * dx34Histo;
AIDA::IHistogram1D * dy34Histo;
AIDA::IHistogram1D * du34Histo;

AIDA::IHistogram1D * dx45Histo;
AIDA::IHistogram1D * dy45Histo;
AIDA::IHistogram1D * du45Histo;

// triplets 0-1-2:

AIDA::IHistogram1D * da02Histo;
AIDA::IHistogram1D * db02Histo;

AIDA::IHistogram1D * tridxHisto;
AIDA::IHistogram1D * tridyHisto;

AIDA::IProfile1D * tridxvsx;
AIDA::IProfile1D * tridxvsy;
AIDA::IProfile1D * tridxvstx;
AIDA::IProfile1D * tridxvsty;
AIDA::IProfile1D * tridyvsx;
AIDA::IProfile1D * tridyvsy;
AIDA::IProfile1D * tridyvstx;
AIDA::IProfile1D * tridyvsty;

AIDA::IHistogram1D * tridx1Histo;
AIDA::IHistogram1D * tridy1Histo;
AIDA::IHistogram1D * tridx3Histo;
AIDA::IHistogram1D * tridy3Histo;
AIDA::IHistogram1D * tridx3bHisto;
AIDA::IHistogram1D * tridy3bHisto;
AIDA::IHistogram1D * tridx4Histo;
AIDA::IHistogram1D * tridy4Histo;
AIDA::IHistogram1D * tridx4bHisto;
AIDA::IHistogram1D * tridy4bHisto;
AIDA::IHistogram1D * tridx5Histo;
AIDA::IHistogram1D * tridy5Histo;
AIDA::IHistogram1D * tridx5bHisto;
AIDA::IHistogram1D * tridy5bHisto;

AIDA::IHistogram1D * trixHisto;
AIDA::IHistogram1D * triyHisto;
AIDA::IHistogram2D * trixyHisto;
AIDA::IHistogram1D * tritxHisto;
AIDA::IHistogram1D * trityHisto;

AIDA::IHistogram1D * dutxHisto; // at DUT
AIDA::IHistogram1D * dutyHisto;
AIDA::IHistogram2D * dutxyHisto;

AIDA::IHistogram2D * cmsxxHisto;
AIDA::IHistogram2D * cmsyyHisto;

AIDA::IHistogram1D * cmssxaHisto;
AIDA::IHistogram1D * cmsdyaHisto;

AIDA::IHistogram1D * cmssxHisto;
AIDA::IHistogram1D * cmsdyHisto;

AIDA::IHistogram1D * cmssxfHisto;
AIDA::IHistogram1D * cmsdyfHisto;
AIDA::IHistogram1D * cmssxfcHisto;
AIDA::IHistogram1D * cmsdyfcHisto;

AIDA::IHistogram1D * cmsdy1Histo;
AIDA::IHistogram1D * cmsdy2Histo;
AIDA::IHistogram1D * cmsdy3Histo;

AIDA::IHistogram1D * cmssxfcnHisto;
AIDA::IHistogram1D * cmsdyfcnHisto;
AIDA::IHistogram1D * cmsdyfcntHisto;
AIDA::IHistogram1D * cmsdyfcntqHisto;
AIDA::IHistogram1D * cmsdyfcntq1Histo;

AIDA::IHistogram1D * cmsdyq0Histo;
AIDA::IHistogram1D * cmsdyq1Histo;
AIDA::IHistogram1D * cmsdyq2Histo;

AIDA::IHistogram1D * cmsqHisto;
AIDA::IHistogram1D * cmsqfHisto;
AIDA::IHistogram1D * cmsqdt0Histo;
AIDA::IHistogram1D * cmsqdt1Histo;
AIDA::IHistogram1D * cmsqdt2Histo;
AIDA::IHistogram1D * cmsqdt3Histo;
AIDA::IHistogram1D * cmsqdt4Histo;
AIDA::IHistogram1D * cmsqdt5Histo;
AIDA::IHistogram1D * cmsxHisto;
AIDA::IHistogram1D * cmsyHisto;
AIDA::IHistogram2D * cmsxyHisto;
AIDA::IProfile1D * cmsdyvsx;
AIDA::IProfile1D * cmsdyvsy;
AIDA::IProfile1D * cmsdyvsty;

AIDA::IHistogram2D * cmspixvsxy;
AIDA::IProfile2D * cmsqvsxy;
AIDA::IProfile1D * cmsqvsdt;

AIDA::IProfile1D * cmsrmsxvsxm;
AIDA::IProfile1D * cmsrmsyvsym;
AIDA::IProfile1D * cmsncolvsxm;
AIDA::IProfile1D * cmsnrowvsym;
AIDA::IProfile1D * cmsetavsym;
AIDA::IProfile1D * cmsetavsym2;
AIDA::IHistogram1D * cmsym1Histo;
AIDA::IHistogram2D * cmsdyvsxHisto;

AIDA::IHistogram2D * effxyHisto;
AIDA::IProfile2D * effvsxy;
AIDA::IProfile1D * effvsx;
AIDA::IProfile1D * effvsxg;
AIDA::IProfile1D * effvsy;
AIDA::IProfile1D * effvst1;
AIDA::IProfile1D * effvst3;
AIDA::IProfile1D * effvsdt;
AIDA::IProfile2D * effvsmm;

AIDA::IProfile2D * rffvsxy;
AIDA::IProfile1D * rffvsx;

AIDA::IHistogram1D * ntriHisto;
AIDA::IProfile1D * lkAvst;

// triplet eff w.r.t. CMS:

AIDA::IHistogram1D * cmsxeHisto;
AIDA::IHistogram1D * cmsyeHisto;

AIDA::IHistogram1D * cmssxeHisto;
AIDA::IHistogram1D * cmsdyeHisto;

AIDA::IHistogram1D * cmsnmHisto;
AIDA::IProfile2D * trieffvsxy;

// driplets 3-4-5:

AIDA::IHistogram1D * dx35Histo;
AIDA::IHistogram1D * dy35Histo;

AIDA::IHistogram1D * dridxHisto;
AIDA::IHistogram1D * dridyHisto;
AIDA::IHistogram1D * drixHisto;
AIDA::IHistogram1D * driyHisto;
AIDA::IHistogram2D * drixyHisto;
AIDA::IHistogram1D * dritxHisto;
AIDA::IHistogram1D * drityHisto;

AIDA::IProfile1D * dridxvsx;
AIDA::IProfile1D * dridxvsy;
AIDA::IProfile1D * dridxvstx;
AIDA::IProfile1D * dridxvsty;
AIDA::IProfile1D * dridyvsx;
AIDA::IProfile1D * dridyvsy;
AIDA::IProfile1D * dridyvstx;
AIDA::IProfile1D * dridyvsty;

AIDA::IHistogram1D * drixcmsHisto;
AIDA::IHistogram1D * driycmsHisto;

AIDA::IHistogram1D * refxHisto; // at REF
AIDA::IHistogram1D * refyHisto;
AIDA::IHistogram2D * refxyHisto;

AIDA::IHistogram1D * refxcmsHisto;
AIDA::IHistogram1D * refycmsHisto;
AIDA::IHistogram2D * refxycmsHisto;
AIDA::IHistogram2D * refpixvsxy;

AIDA::IHistogram1D * refqHisto;
AIDA::IProfile2D * refqvsxy;

AIDA::IHistogram2D * refxxHisto; //REF vs driplet
AIDA::IHistogram2D * refyyHisto;

AIDA::IHistogram1D * refsxaHisto;
AIDA::IHistogram1D * refdyaHisto;
AIDA::IHistogram1D * refsxHisto;
AIDA::IHistogram1D * refdyHisto;
AIDA::IHistogram1D * refsxcHisto;
AIDA::IHistogram1D * refdycHisto;

AIDA::IProfile1D * refdyvsx;
AIDA::IProfile1D * refdyvsy;
AIDA::IProfile1D * refdyvsty;

AIDA::IHistogram1D * bacsxaHisto;
AIDA::IHistogram1D * bacdyaHisto;
AIDA::IHistogram1D * bacsxcHisto;
AIDA::IHistogram1D * bacdycHisto;
AIDA::IHistogram1D * bacsxcqHisto;
AIDA::IHistogram1D * bacdycqHisto;

AIDA::IHistogram1D * ndriHisto;
AIDA::IHistogram1D * ndricmsHisto;
AIDA::IProfile1D * lkBvst;

AIDA::IHistogram1D * sixkxHisto; //driplet-triplet
AIDA::IHistogram1D * sixkyHisto;
AIDA::IHistogram2D * sixxyHisto;
AIDA::IHistogram2D * sixxycHisto;
AIDA::IHistogram2D * sixxylkHisto;
AIDA::IHistogram1D * sixdxHisto;
AIDA::IHistogram1D * sixdyHisto;
AIDA::IProfile2D * kinkvsxy;

AIDA::IHistogram1D * selxHisto;
AIDA::IHistogram1D * selyHisto;
AIDA::IHistogram1D * selaxHisto;
AIDA::IHistogram1D * selayHisto;
AIDA::IHistogram1D * seldxHisto;
AIDA::IHistogram1D * seldyHisto;
AIDA::IHistogram1D * selkxHisto;
AIDA::IHistogram1D * selkyHisto;

AIDA::IHistogram1D * seldx1Histo;
AIDA::IHistogram1D * seldy1Histo;
AIDA::IHistogram1D * seldx3Histo;
AIDA::IHistogram1D * seldy3Histo;
AIDA::IHistogram1D * seldx4Histo;
AIDA::IHistogram1D * seldy4Histo;
AIDA::IHistogram1D * seldx5Histo;
AIDA::IHistogram1D * seldy5Histo;
AIDA::IHistogram1D * seldx6Histo;
AIDA::IHistogram1D * seldy6Histo;

AIDA::IHistogram1D * gblndfHisto;
AIDA::IHistogram1D * gblchi2aHisto;
AIDA::IHistogram1D * gblchi2bHisto;
AIDA::IHistogram1D * gblprbHisto;

AIDA::IHistogram1D * badxHisto;
AIDA::IHistogram1D * badyHisto;
AIDA::IHistogram1D * badaxHisto;
AIDA::IHistogram1D * badayHisto;
AIDA::IHistogram1D * baddxHisto;
AIDA::IHistogram1D * baddyHisto;
AIDA::IHistogram1D * badkxHisto;
AIDA::IHistogram1D * badkyHisto;

AIDA::IHistogram1D * baddx1Histo;
AIDA::IHistogram1D * baddy1Histo;
AIDA::IHistogram1D * baddx3Histo;
AIDA::IHistogram1D * baddy3Histo;
AIDA::IHistogram1D * baddx4Histo;
AIDA::IHistogram1D * baddy4Histo;
AIDA::IHistogram1D * baddx5Histo;
AIDA::IHistogram1D * baddy5Histo;
AIDA::IHistogram1D * baddx6Histo;
AIDA::IHistogram1D * baddy6Histo;

AIDA::IHistogram1D * goodx1Histo;
AIDA::IHistogram1D * goody1Histo;
AIDA::IHistogram1D * goodx6Histo;
AIDA::IHistogram1D * goody6Histo;

AIDA::IHistogram1D * gblax0Histo;
AIDA::IHistogram1D * gbldx0Histo;
AIDA::IHistogram1D * gblrx0Histo;
AIDA::IHistogram1D * gblpx0Histo;
AIDA::IHistogram1D * gblqx0Histo;

AIDA::IHistogram1D * gblax1Histo;
AIDA::IHistogram1D * gbldx1Histo;
AIDA::IHistogram1D * gblrx1Histo;
AIDA::IHistogram1D * gblpx1Histo;
AIDA::IHistogram1D * gblqx1Histo;

AIDA::IHistogram1D * gblax2Histo;
AIDA::IHistogram1D * gbldx2Histo;
AIDA::IHistogram1D * gblrx2Histo;
AIDA::IHistogram1D * gblpx2Histo;
AIDA::IHistogram1D * gblqx2Histo;

AIDA::IHistogram1D * gblax3Histo;
AIDA::IHistogram1D * gbldx3Histo;
AIDA::IHistogram1D * gblrx3Histo;
AIDA::IHistogram1D * gblpx3Histo;
AIDA::IHistogram1D * gblqx3Histo;

AIDA::IHistogram1D * gblax4Histo;
AIDA::IHistogram1D * gbldx4Histo;
AIDA::IHistogram1D * gblrx4Histo;
AIDA::IHistogram1D * gblpx4Histo;
AIDA::IHistogram1D * gblqx4Histo;

AIDA::IHistogram1D * gblax5Histo;
AIDA::IHistogram1D * gbldx5Histo;
AIDA::IHistogram1D * gblrx5Histo;
AIDA::IHistogram1D * gblpx5Histo;
AIDA::IHistogram1D * gblqx5Histo;

AIDA::IHistogram1D * gblax6Histo;
AIDA::IHistogram1D * gbldx6Histo;
AIDA::IHistogram1D * gbldy6Histo;
AIDA::IHistogram1D * gblrx6Histo;
AIDA::IHistogram1D * gblry6Histo;
AIDA::IHistogram1D * gblpx6Histo;
AIDA::IHistogram1D * gblqx6Histo;

AIDA::IHistogram1D * gblkx1Histo;
AIDA::IHistogram1D * gblkx2Histo;
AIDA::IHistogram1D * gblkx3Histo;
AIDA::IHistogram1D * gblkx4Histo;
AIDA::IHistogram1D * gblkx5Histo;
AIDA::IHistogram1D * gblkx6Histo;

AIDA::IHistogram1D * sixdxcHisto;
AIDA::IHistogram1D * sixdycHisto;
AIDA::IHistogram1D * sixkxcHisto;
AIDA::IHistogram1D * sixkycHisto;
AIDA::IHistogram1D * sixkxsHisto;
AIDA::IHistogram1D * sixkysHisto;

AIDA::IHistogram1D * sixzx3Histo;
AIDA::IHistogram1D * sixzy3Histo;
AIDA::IHistogram1D * sixzx2Histo;
AIDA::IHistogram1D * sixzy2Histo;
AIDA::IHistogram1D * sixzx1Histo;
AIDA::IHistogram1D * sixzy1Histo;

AIDA::IHistogram1D * sixkxzyHisto;
AIDA::IHistogram1D * sixkyzxHisto;
AIDA::IHistogram1D * sixkxzxHisto;
AIDA::IHistogram1D * sixkyzyHisto;

std::string EUTelTestFitter::_linChi2HistoName   = "linChi2";
std::string EUTelTestFitter::_logChi2HistoName   = "logChi2";
std::string EUTelTestFitter::_firstChi2HistoName = "firstChi2";
std::string EUTelTestFitter::_bestChi2HistoName  = "bestChi2";
std::string EUTelTestFitter::_fullChi2HistoName  = "fullChi2";
std::string EUTelTestFitter::_nTrackHistoName    = "nTrack";
std::string EUTelTestFitter::_nAllHitHistoName   = "nAllHit";
std::string EUTelTestFitter::_nAccHitHistoName   = "nAccHit";
std::string EUTelTestFitter::_nHitHistoName      = "nHit";
std::string EUTelTestFitter::_nBestHistoName     = "nBest";
std::string EUTelTestFitter::_hitAmbiguityHistoName  = "nAmbig";
#endif
int nclst = 0;
int nmtch = 0;
int naldut = 0;
int nmille = 0;
bool ldb = 0;

MilleBinary * mille; // for producing MillePede-II binary file


EUTelTestFitter::EUTelTestFitter() : Processor("EUTelTestFitter") {

  // modify processor description
  _description = "Analytical track fitting processor for EUDET telescope";


  // register steering parameters:
  //       name, description, class-variable, default value

  // input collection first:

  registerInputCollection( LCIO::TRACKERHIT,
                           "InputCollectionName" ,
                           "Name of the input TrackerHit collection"  ,
                           _inputColName ,
                           std::string("meshit") );

  // output collection:

  registerOutputCollection(LCIO::TRACK,"OutputTrackCollectionName",
                           "Collection name for fitted tracks",
                           _outputTrackColName, string( "testfittracks"));

  registerOutputCollection(LCIO::TRACKERHIT,"CorrectedHitCollectionName",
                           "Collection name for corrected particle positions",
                           _correctedHitColName, string( "corrfithits"));

  registerOutputCollection(LCIO::TRACKERHIT,"OutputHitCollectionName",
                           "Collection name for fitted particle hits( positions)",
                           _outputHitColName, string( "testfithits"));

  // alignment collections( might be important ):

  EVENT::StringVec	_alignmentCollectionExample;
  _alignmentCollectionExample.push_back("alignment");
  
  registerProcessorParameter( "alignmentCollectionNames",
                            "List of alignment collections which are needed to get track position on a Sensor surface ",
                            _alignmentCollectionNames, _alignmentCollectionExample );

  // compulsory parameters:

  registerProcessorParameter( "InputHitsInTrack",
                              "Flag for storing input( measured) hits in track",
                              _InputHitsInTrack, static_cast < bool >( false));

  registerProcessorParameter( "OutputHitsInTrack",
                              "Flag for storing output( fitted) hits in track",
                              _OutputHitsInTrack, static_cast < bool >( true));

  registerProcessorParameter( "DebugEventCount",
                              "Print out every DebugEnevtCount event",
                              _debugCount, static_cast < int >( 10));

  registerProcessorParameter( "AllowMissingHits",
                              "Allowed number of missing hits in the track",
                              _allowMissingHits, static_cast < int >( 0));

  registerProcessorParameter( "AllowSkipHits",
                              "Allowed number of hits removed from the track",
                              _allowSkipHits, static_cast < int >( 0));

  registerProcessorParameter( "MaxPlaneHits",
                              "Maximum number of considered hits per plane",
                              _maxPlaneHits, static_cast < int >( 100));

  registerProcessorParameter( "MissingHitPenalty",
                              "Chi2 penalty for missing hit in the track",
                              _missingHitPenalty, static_cast < double >( 0.));

  registerProcessorParameter( "SkipHitPenalty",
                              "Chi2 penalty for removing hit from the track",
                              _skipHitPenalty, static_cast < double >( 100.));

  registerProcessorParameter( "Chi2Max",
                              "Maximum Chi2 for accepted track fit",
                              _chi2Max, static_cast < double >( 100.));

  registerProcessorParameter( "UseNominalResolution",
                              "Flag for using nominal resolution instead of position errors",
                              _useNominalResolution, static_cast < bool >( true));

  registerProcessorParameter( "UseDUT",
                              "Flag for including DUT measurement in the fit",
                              _useDUT, static_cast < bool >( false));

  registerProcessorParameter( "Ebeam",
                              "Beam energy [GeV]",
                              _eBeam, static_cast < double >( 4.0));

  registerProcessorParameter("HistoInfoFileName",
                             "Name of the histogram information file",
                             _histoInfoFileName, string( "histoinfo.xml" ) );

  // optional parameters

  // Parameters added to allow correlation band info 02 August 2010 libov@mail.desy.de

  registerOptionalParameter( "nSkip",
                              "Skip n CMS events at begin of file",
                              _nSkip, static_cast < int >(-1) );

  registerOptionalParameter("UseSlope","Use expected track direction to constraint number of considered hit combinations( track preselection).", _UseSlope, true );

  registerOptionalParameter("SlopeXLimit","Limit on track slope change when passing sensor layer( in X direction)", _SlopeXLimit, static_cast <float>( 0.001));

  registerOptionalParameter("SlopeYLimit","Limit on track slope change when passing sensor layer( in Y direction)", _SlopeYLimit, static_cast <float>( 0.001));

  registerOptionalParameter("SlopeDistanceMax","Maximum hit distance from the expected position, used for hit preselection", _SlopeDistanceMax, static_cast <float>( 1.));

  // -------------------------------------------------------------------------------------------------

  std::vector<int > initLayerIDs;
  std::vector<float > initLayerShift;

  registerOptionalParameter( "SkipLayerIDs",
                             "Ids of layers which should NOT be included in the fit",
                             _SkipLayerIDs, initLayerIDs);

  registerOptionalParameter( "PassiveLayerIDs",
                             "Ids of layers which should be treated as passive in the fit",
                             _PassiveLayerIDs, initLayerIDs);

  registerOptionalParameter( "AlignLayerIDs",
                             "Ids of layers for which alignment corrections are given",
                             _AlignLayerIDs, initLayerIDs);

  registerOptionalParameter( "AlignLayerShiftX",
                             "Alignment corrections in X for these layers",
                             _AlignLayerShiftX, initLayerShift);

  registerOptionalParameter( "AlignLayerShiftY",
                             "Alignment corrections in Y for these layers",
                             _AlignLayerShiftY, initLayerShift);

  registerOptionalParameter( "AlignLayerRotZ",
                             "Rotation around Z for layer alignment",
                             _AlignLayerRotZ, initLayerShift);

  registerOptionalParameter( "WindowLayerIDs",
                             "Ids of layers for which position window cut are applied",
                             _WindowLayerIDs, initLayerIDs);

  registerOptionalParameter( "WindowMinX",
                             "Lower window edge in X",
                             _WindowMinX, initLayerShift);

  registerOptionalParameter( "WindowMaxX",
                             "Upper window edge in X",
                             _WindowMaxX, initLayerShift);

  registerOptionalParameter( "WindowMinY",
                             "Lower window edge in Y",
                             _WindowMinY, initLayerShift);

  registerOptionalParameter( "WindowMaxY",
                             "Upper window edge in Y",
                             _WindowMaxY, initLayerShift);


  registerOptionalParameter( "MaskLayerIDs",
                             "Ids of layers for which position masks are applied",
                             _MaskLayerIDs, initLayerIDs);

  registerOptionalParameter( "MaskMinX",
                             "Lower mask edge in X",
                             _MaskMinX, initLayerShift);

  registerOptionalParameter( "MaskMaxX",
                             "Upper mask edge in X",
                             _MaskMaxX, initLayerShift);

  registerOptionalParameter( "MaskMinY",
                             "Lower mask edge in Y",
                             _MaskMinY, initLayerShift);

  registerOptionalParameter( "MaskMaxY",
                             "Upper mask edge in Y",
                             _MaskMaxY, initLayerShift);


  registerOptionalParameter( "UseBeamConstraint",
                             "Flag for using beam direction constraint in the fit",
                             _useBeamConstraint, static_cast < bool >( false));

  registerOptionalParameter( "BeamSpread",
                             "Assumed angular spread of the beam [rad]( for beam constraint)",
                             _beamSpread, static_cast < double >( 0.1) );

  registerOptionalParameter( "BeamSlopeX",
                             "Beam direction tilt in X-Z plane [rad]( for beam constraint)",
                             _beamSlopeX, static_cast < double >( 0.));

  registerOptionalParameter( "BeamSlopeY",
                             "Beam direction tilt in Y-Z plane [rad]( for beam constraint)",
                             _beamSlopeY, static_cast < double >( 0.));

  registerOptionalParameter( "SearchMultipleTracks",
                             "Flag for searching multiple tracks in events with multiple hits",
                             _searchMultipleTracks, static_cast < bool >( true));

  registerOptionalParameter( "AllowAmbiguousHits",
                             "Allow same hit to be used in more than one track",
                             _allowAmbiguousHits, static_cast < bool >( false));

  registerOptionalParameter( "MaximumAmbiguousHits",
                             "Maximum number of hits to be shared by more than one track",
                             _maximumAmbiguousHits, static_cast < int >( 2));

  // initialize all the counters
  _noOfEventWOInputHit = 0;
  _noOfEventWOTrack = 0;
  _noOfTracks       = 0;
}

//------------------------------------------------------------------------------
void EUTelTestFitter::init() {

  // usually a good idea to
  printParameters();

  _nRun = 0;
  _nEvt = 0;

  _isFirstEvent = true;

  // check if Marlin was built with GEAR support or not
#ifndef USE_GEAR

  streamlog_out( ERROR2 ) << "Marlin was not built with GEAR support.\n"
    "You need to install GEAR and recompile Marlin with -DUSE_GEAR before continue." << endl;

  // I'm thinking if this is the case of throwing an exception or
  // not. This is a really error and not something that can
  // exceptionally happens. Still not sure what to do
  exit(-1);

#else

  // check if the GEAR manager pointer is not null!
  if( Global::GEAR == 0x0 ) {
    streamlog_out( ERROR2 ) <<  "The GearMgr is not available, for an unknown reason."  << endl;
    exit(-1);
  }

  streamlog_out(MESSAGE0) << "assumed beam energy " << _eBeam << " GeV" <<  endl;

  // Read geometry information from GEAR

  streamlog_out( MESSAGE0 ) << "Reading telescope geometry description from GEAR " << endl;

  _siPlanesParameters  = const_cast<gear::SiPlanesParameters* >( &(Global::GEAR->getSiPlanesParameters()));
  _siPlanesLayerLayout = const_cast<gear::SiPlanesLayerLayout*>(  &(_siPlanesParameters->getSiPlanesLayerLayout() ));

#endif

  // Test output

  if( _SkipLayerIDs.size() ) {
    streamlog_out( MESSAGE0 )  <<  _SkipLayerIDs.size() << " layers should be skipped" << endl;
  }

  if( _PassiveLayerIDs.size() ) {
    streamlog_out( MESSAGE0 ) <<  _PassiveLayerIDs.size() << " layers should be considered passive" << endl;
  }

  // Take all layers defined in GEAR geometry
  _nTelPlanes = _siPlanesLayerLayout->getNLayers();

  // Check for DUT
  if( _siPlanesParameters->getSiPlanesType()==_siPlanesParameters->TelescopeWithDUT ) {
    _iDUT = _nTelPlanes;
    _nTelPlanes++;
  }
  else {
    _iDUT = -1;
  }

  // _nTelPlanes is the total number of sensitive layers in the setup
  // summing both the telescopes and the DUT 
  // Read position in Z( for sorting), skip layers if requested

  _planeSort = new int[_nTelPlanes];
  _planePosition = new double[_nTelPlanes]; // z pos

  int nSkip = 0;

  for( int ipl = 0; ipl <  _siPlanesLayerLayout->getNLayers(); ipl++ ) {
    _planePosition[ipl-nSkip] = _siPlanesLayerLayout->getLayerPositionZ(ipl);
    _planeSort[ipl-nSkip]=ipl; 

    // Check if not on "skip list"

    int _plID = _siPlanesLayerLayout->getID(ipl);

    for( int spl = 0; spl<( int)_SkipLayerIDs.size(); spl++ ) {
      if( _SkipLayerIDs.at(spl) == _plID) {
        streamlog_out( MESSAGE0 ) <<  "Skipping layer ID " << _plID
                                   << " at Z = " << _siPlanesLayerLayout->getLayerPositionZ(ipl) << endl;
        nSkip++;
        break;
      }
    }
  }

  _nTelPlanes -= nSkip;

  if( _iDUT > 0 )  {
    _planePosition[_iDUT-nSkip] = _siPlanesLayerLayout->getDUTPositionZ();
    _planeSort[_iDUT-nSkip] = _iDUT;
  }

  // Binary sorting
  // Note: all this code can in principle be replaced by a simple sort 

  bool sorted;
  do {
    sorted = false;
    for( int iz = 0; iz < _nTelPlanes-1; iz++ ) {
      if( _planePosition[iz] > _planePosition[iz+1] ) {
        double _posZ = _planePosition[iz];
        _planePosition[iz] = _planePosition[iz+1];
        _planePosition[iz+1] = _posZ;

        int _idZ = _planeSort[iz];
        _planeSort[iz] = _planeSort[iz+1];
        _planeSort[iz+1] = _idZ;

        sorted = true;
      }
    }
  } while(sorted);

  // Book local geometry arrays
  //
  // N.B. the plane ID array is containing the sensorID sorted
  // according to the position along Z. We will be using this array
  // for histogram booking and filling.
  _planeID         = new int[_nTelPlanes];
  _planeShiftX     = new double[_nTelPlanes];
  _planeShiftY     = new double[_nTelPlanes];
  _planeRotZ       = new double[_nTelPlanes];
  _planeThickness  = new double[_nTelPlanes];
  _planeX0         = new double[_nTelPlanes];
  _planeResolution = new double[_nTelPlanes];
  _isActive        = new bool[_nTelPlanes];
  _planeWindowIDs  = new std::vector<int>[_nTelPlanes];
  _planeMaskIDs    = new std::vector<int>[_nTelPlanes];

  _nActivePlanes = 0;

  // Fill arrays with parameters of layer, sorted in Z

  for( int iz = 0; iz < _nTelPlanes; iz++ ) {

    int ipl = _planeSort[iz];

    int iActive;
    double resolution;

    // All dimensions are assumed to be in mm !!!

    if(ipl != _iDUT )    {
      _planeID[iz]=_siPlanesLayerLayout->getID(ipl);
      _planeThickness[iz]=_siPlanesLayerLayout->getLayerThickness(ipl);
      _planeX0[iz]=_siPlanesLayerLayout->getLayerRadLength(ipl);
      resolution = _siPlanesLayerLayout->getSensitiveResolution(ipl);
    }
    else {
      _planeID[iz]=_siPlanesLayerLayout->getDUTID();
      _planeThickness[iz]=_siPlanesLayerLayout->getDUTThickness();
      _planeX0[iz]=_siPlanesLayerLayout->getDUTRadLength();
      resolution = _siPlanesLayerLayout->getDUTSensitiveResolution();
    }

    iActive =( resolution > 0);

    // Check passive layer list
    for( int ppl = 0; ppl<( int)_PassiveLayerIDs.size(); ppl++) {
      if( _PassiveLayerIDs.at(ppl) == _planeID[iz])    {
        streamlog_out( MESSAGE0 ) <<  "Force passive layer ID " << _planeID[iz]
                                   << " at Z = " << _planePosition[iz] << endl;
        iActive = false;
        break;
      }
    }

    if( iActive &&( ipl != _iDUT || _useDUT ) )   {
      _isActive[iz] = true;
      _planeResolution[iz]=resolution;
      _nActivePlanes++;
    }
    else  {
      _isActive[iz] = false;
      _planeResolution[iz] = 0.;
    }

    // Alignment corrections should already be added in HitMaker
    // But look in input options for possible adjustments

    _planeShiftX[iz] = 0.;
    _planeShiftY[iz] = 0.;
    _planeRotZ[iz] = 0.;

    for( int apl = 0; apl<( int)_AlignLayerIDs.size(); apl++) {
      if( _AlignLayerIDs.at(apl) == _planeID[iz]) {
        _planeShiftX[iz]=_AlignLayerShiftX.at(apl);
        _planeShiftY[iz]=_AlignLayerShiftY.at(apl);
        // Rotation can be skipped: check size
        if(apl <( int)_AlignLayerRotZ.size()) {
          _planeRotZ[iz]=_AlignLayerRotZ.at(apl);
        }
        break;
      }
    }

    // Check, if there are additional cuts defined for this plane
    for( int wpl = 0; wpl<( int)_WindowLayerIDs.size(); wpl++) {
      if( _WindowLayerIDs.at(wpl) == _planeID[iz]) {
        _planeWindowIDs[iz].push_back(wpl);
      }
    }
    for( int mpl = 0; mpl<( int)_MaskLayerIDs.size(); mpl++) {
      if( _MaskLayerIDs.at(mpl) == _planeID[iz]) {
        _planeMaskIDs[iz].push_back(mpl);
      }
    }
  } // End of plane loop

  // Get new DUT position( after sorting)
  for( int iz = 0;iz< _nTelPlanes; iz++) {
    if(_planeSort[iz]==_iDUT)  {
      _iDUT=iz;
      break;
    }
  }

  // Print out geometry information
  streamlog_out( MESSAGE2 ) <<  "Telescope configuration with " << _nTelPlanes << " planes" << endl;

  for( int ipl = 0; ipl < _nTelPlanes; ipl++) {
    stringstream ss;
    if( ipl == _iDUT) {
      ss << "D.U.T.  plane";
    }
    else {
      if(_isActive[ipl]) {
        ss << "Active  plane";
      }
      else {
        ss << "Passive plane";
      }
      ss << "  ID = " << _planeID[ipl]
         << "  at Z [mm] = " << _planePosition[ipl]
         << " dZ [um] = " << _planeThickness[ipl]*1000.;

      if( _isActive[ipl] )  {
        ss << "  Res [um] = " << _planeResolution[ipl]*1000.;

        if( _planeShiftX[ipl] != 0. || _planeShiftY[ipl] != 0. ) {
          ss << "\n  alignment corrections:"
             <<  " dX [mm] = " << _planeShiftX[ipl]
             << " dY [mm] = " << _planeShiftY[ipl];
        }
        if( _planeRotZ[ipl] != 0. ) {
          ss << " RotZ [rad] = " << _planeRotZ[ipl];
        }

        for( int wpl = 0; wpl < ( int)_planeWindowIDs[ipl].size(); wpl++ )  {
          int iwin= _planeWindowIDs[ipl].at(wpl);
          ss << "\n accepted window: X = " <<  _WindowMinX.at(iwin)
             << " to " <<  _WindowMaxX.at(iwin)
             << " Y = " <<  _WindowMinY.at(iwin)
             << " to " <<  _WindowMaxY.at(iwin);
        }

        for( int mpl = 0; mpl < ( int)_planeMaskIDs[ipl].size(); mpl++ )  {
          int imsk= _planeMaskIDs[ipl].at(mpl);
          ss << "\n imposed mask: X = " <<  _MaskMinX.at(imsk)
             << " to " <<  _MaskMaxX.at(imsk)
             << " Y = " <<  _MaskMinY.at(imsk)
             << " to " <<  _MaskMaxY.at(imsk);
        }
      }

      streamlog_out( MESSAGE2 ) <<  ss.str() << endl;
    }
  }
  streamlog_out( MESSAGE2 ) << "Total of " << _nActivePlanes << " active sensor planes " << endl;

  // Allocate arrays for track fitting

  _planeHits   = new int[_nTelPlanes];
  _planeChoice = new int[_nTelPlanes];
  _planeMod    = new type_fitcount[_nTelPlanes];

  _planeX  = new double[_nTelPlanes];
  _planeEx = new double[_nTelPlanes];
  _planeY  = new double[_nTelPlanes];
  _planeEy = new double[_nTelPlanes];

  _planeScatAngle  = new double[_nTelPlanes];

  _planeDist = new double[_nTelPlanes];
  _planeScat = new double[_nTelPlanes];

  _fitX  = new double[_nTelPlanes];
  _fitEx = new double[_nTelPlanes];
  _fitY  = new double[_nTelPlanes];
  _fitEy = new double[_nTelPlanes];

  int arrayDim = _nTelPlanes * _nTelPlanes;

  _fitArray = new double[arrayDim];
  _nominalFitArrayX = new double[arrayDim];
  _nominalErrorX = new double[_nTelPlanes];

  _nominalFitArrayY = new double[arrayDim];
  _nominalErrorY = new double[_nTelPlanes];

  // Fill nominal fit matrices and
  // calculate expected precision of track fitting

  // Planes are ordered in position along the beam line !

  double totalScatAngle = 0.;

  for( int ipl = 0; ipl < _nTelPlanes; ipl++ ) {
    if( ipl > 0) {
      _planeDist[ipl-1] = 1. / ( _planePosition[ipl] - _planePosition[ipl-1] );
    }

    _planeScatAngle[ipl]= 0.0136/_eBeam * sqrt(_planeThickness[ipl]/_planeX0[ipl])
      *( 1.+0.038*std::log(_planeThickness[ipl]/_planeX0[ipl]));

    if( ipl == 0 && _useBeamConstraint ) {
      _planeScat[ipl]= 1./(_planeScatAngle[ipl]*_planeScatAngle[ipl]+ _beamSpread*_beamSpread);
    }
    else {
      _planeScat[ipl]= 1./(_planeScatAngle[ipl] * _planeScatAngle[ipl]);
    }

    totalScatAngle+= _planeScatAngle[ipl] * _planeScatAngle[ipl];

    _fitX[ipl] =_fitY[ipl] = 0.;
    _nominalErrorX[ipl]= _planeResolution[ipl];
    _nominalErrorY[ipl]= _planeResolution[ipl];
  }

  totalScatAngle = sqrt(totalScatAngle);

  // Fit with nominal parameters for X direction

  int status = DoAnalFit(_fitX,_nominalErrorX,_beamSlopeX);

  if(status) {
    streamlog_out( ERROR2 ) << "\n Fit in X with nominal geometry failed !?!" << endl;
  }
  // Store fit matrix

  for( int imx = 0; imx<arrayDim; imx++) {
    _nominalFitArrayX[imx] = _fitArray[imx];
  }

  stringstream ss;
  ss << "Expected position resolutions in X [um]: ";
  for( int ipl = 0; ipl < _nTelPlanes; ipl++) {
    ss << _nominalErrorX[ipl]*1000. << "  ";
  }

  ss << endl << "Expected scattering angle [mrad]: ";
  for( int ipl = 0; ipl < _nTelPlanes; ipl++) {
    ss << _planeScatAngle[ipl]*1000. << "  ";
  }

  ss << endl << "Expected total scattering angle [mrad]: " << totalScatAngle*1000.;

  streamlog_out( MESSAGE2 ) << ss.str() << endl;


  // Fit with nominal parameters for Y direction

  status = DoAnalFit(_fitY,_nominalErrorY,_beamSlopeY);

  if(status) {
    streamlog_out( ERROR2 ) << "\n Fit in Y with nominal geometry failed !?!" << endl;
  }
  // Store fit matrix

  for( int imx = 0; imx<arrayDim; imx++) {
    _nominalFitArrayY[imx] = _fitArray[imx];
  }

// Check if slope-based preselection parameter values are not too small

  if( _UseSlope && 
        _SlopeXLimit < 5.*totalScatAngle 
     )
    streamlog_out( ERROR2 ) << "SlopeXLimit cut probably too tight! Check parameters!" << endl; 
  

  if( _UseSlope && 
        _SlopeYLimit < 5.*totalScatAngle 
     )
    streamlog_out( ERROR2 ) << "SlopeXLimit cut probably too tight! Check parameters!" << endl; 
  
  // Take into account beam spread

  totalScatAngle = sqrt(totalScatAngle*totalScatAngle + _beamSpread*_beamSpread);

  if( _UseSlope && 
  _SlopeDistanceMax < 5.*totalScatAngle *( _planePosition[_nTelPlanes-1] - _planePosition[0])
     )
    streamlog_out( ERROR2 ) << "SlopeDistanceMax cut probably too tight! Check parameters!" << endl; 

  // Book histograms:

#if defined(USE_AIDA) || defined(MARLIN_USE_AIDA)
  bookHistos();
#endif

  mille = new MilleBinary::MilleBinary( "mille.bin" );

}//init

//------------------------------------------------------------------------------
void EUTelTestFitter::processRunHeader( LCRunHeader* runHeader) {

  auto_ptr<EUTelRunHeaderImpl> eutelHeader( new EUTelRunHeaderImpl( runHeader ) );
  eutelHeader->addProcessor( type() );

  _nRun++;

  // Decode and print out Run Header information - just a check

  streamlog_out( MESSAGE2 )  << "Processing run header " << _nRun
                             << ", run nr " << runHeader->getRunNumber() << endl;

  const std::string detectorName = runHeader->getDetectorName();
  const std::string detectorDescription = runHeader->getDescription();
  const std::vector<std::string> * subDets = runHeader->getActiveSubdetectors();

  streamlog_out( MESSAGE0 ) << detectorName << " : " << detectorDescription << endl;

  int nDet = subDets->size();

  streamlog_out( MESSAGE0 ) << nDet << " subdetectors defined :" << endl;
  stringstream ss;
  for( int idet = 0;idet<nDet;idet++) {
    streamlog_out( MESSAGE0 )  << idet+1 << " : " << subDets->at(idet) << endl;
  }

} // processRunHeader


//------------------------------------------------------------------------------
TMatrixD Jac5( double ds ) {
  /* for GBL:
     Jacobian for straight line track
     track = q/p, x', y', x, y
             0,   1,  2,  3, 4
  */
  TMatrixD jac(5, 5);
  jac.UnitMatrix();
  jac[3][1] = ds; // x = x0 + xp * ds
  jac[4][2] = ds; // y = y0 + yp * ds
  return jac;
}

//------------------------------------------------------------------------------
void EUTelTestFitter::processEvent( LCEvent * event ) {

  EUTelEventImpl * euEvent = static_cast<EUTelEventImpl*>( event );

  if( euEvent->getEventType() == kEORE ) {
    streamlog_out( DEBUG ) <<  "EORE found: nothing else to do." << endl;
    return;
  }

  bool debug = _nEvt < _debugCount;

  static unsigned long prevtime = event->getTimeStamp(); //init

  static unsigned long time0 = 0;
  if( _isFirstEvent ) time0 = event->getTimeStamp();

  // phase from TLU time between events:

  bool lprnt = 0;

  // TLU time stamp: count 384 MHz clock (384 = 48*8)
  // CMS timestamp: count 40 MHz clock
  // sysdt: TLU/CMS = 1.00007
  // => CMS = TLU/1.00007
  // or TLU/1.00007 - CMS = 0
  // sysddt: <TLU/1.00007344 - CMS> = 0

  //double gTLU = 0.3840282; // [GHz] TB 21
  //double gTLU = 0.3846769; // [GHz] TB 22
  double gTLU = 0.3846740; // [GHz] TB 22

  double fTLU = gTLU * 1E9; // [Hz]

  if( _nEvt < 10 || _nEvt % 100 == 0 ) {

    //lprnt = 1;

    streamlog_out( MESSAGE2 ) << "Processing event "
                              << setw( 7) << setiosflags(ios::right)
			      << event->getEventNumber()
			      << " in run " << setw( 4) << setiosflags(ios::right)
			      << event->getRunNumber()
			      << ", time " << setw(13) << setiosflags(ios::right)
			      << event->getTimeStamp()
			      << " = " << fixed << setprecision(1)
			      << ( event->getTimeStamp() - time0 )/fTLU << " s"
			      << ", dt = " << int( (event->getTimeStamp() - prevtime)/gTLU/1E3 ) << " us"
			      << ", rate "
			      << (_nEvt+1) / ( event->getTimeStamp() - time0 + 0.1 ) * fTLU
			      << " Hz"
                              << endl;
  }//lpr 1'127'591'287'209

  t100Histo->fill( ( event->getTimeStamp() - time0 ) / fTLU ); //event time
  t300Histo->fill( ( event->getTimeStamp() - time0 ) / fTLU ); //event time
  t1000Histo->fill( ( event->getTimeStamp() - time0 ) / fTLU ); //event time
  t3600Histo->fill( ( event->getTimeStamp() - time0 ) / fTLU ); //event time

  if( _nEvt > 0 ) {

    dtHisto->fill( ( event->getTimeStamp() - prevtime ) / 384E0 ); //us
    dtmsHisto->fill( ( event->getTimeStamp() - prevtime ) / 384E3 ); //ms
    logdtHisto->fill( std::log( ( event->getTimeStamp() - prevtime ) / 384E3 ) / std::log(10.0) );

    if( ( event->getTimeStamp() - prevtime ) / fTLU > 10 ) {
      cout << "long gap event "
	   << setw( 6) << setiosflags(ios::right)
	   << event->getEventNumber()
	   << " in run " << setw( 6) << setiosflags(ios::right)
	   << event->getRunNumber()
	   << ", time " << setw(15) << setiosflags(ios::right)
	   << event->getTimeStamp()
	   << ", " << setw(10) << setprecision(9)
	   << ( event->getTimeStamp() - time0 ) /fTLU << " s"
	   << ", dt = " << int( (event->getTimeStamp() - prevtime)/384E0 ) << " us"
	   << endl;
    }//gap

  }//dt plots

  //----------------------------------------------------------------------------
  // CMS pixel:

  // Open the CMS pixel data stream:
  if( _isFirstEvent ) {
    _dutReader.open( event->getRunNumber(), 1 );
    _refReader.open( event->getRunNumber(), 2 );
  }//first event

  vector <cluster> ClustDUT;
  vector <cluster> ClustREF;

  int NpixDUT = 0;
  int NpixREF = 0;

  unsigned long dutTimeStamp = 0;
  unsigned long refTimeStamp = 0;
  if( _dutReader.isOpen() ) ClustDUT = _dutReader.getCMSclust( dutTimeStamp, NpixDUT );
  if( _refReader.isOpen() ) ClustREF = _refReader.getCMSclust( refTimeStamp, NpixREF );

  if( _isFirstEvent ) {

    int nskip = 0;
    if( event->getRunNumber() == 1946 ) nskip = 70;

    if( event->getRunNumber() == 2153 ) nskip = 19;
    if( event->getRunNumber() == 2156 ) nskip = 52;
    if( event->getRunNumber() == 2163 ) nskip =  7;

    if( event->getRunNumber() == 2257 ) nskip = 70;
    if( event->getRunNumber() == 2258 ) nskip = 91;
    if( event->getRunNumber() == 2261 ) nskip =102;
    if( event->getRunNumber() == 2265 ) nskip =  0;
    if( event->getRunNumber() == 2266 ) nskip =  0;
    if( event->getRunNumber() == 2267 ) nskip =  0;
    if( event->getRunNumber() == 2269 ) nskip =107;

    if( event->getRunNumber() == 2293 ) nskip =  1;
    if( event->getRunNumber() == 2295 ) nskip = 22;

    if( event->getRunNumber() == 2314 ) nskip = 25;

    if( event->getRunNumber() == 2322 ) nskip = 39;
    if( event->getRunNumber() == 2325 ) nskip =  0;
    if( event->getRunNumber() == 2326 ) nskip =  1;
    if( event->getRunNumber() == 2328 ) nskip =  2;
    if( event->getRunNumber() == 2330 ) nskip =  3;
    if( event->getRunNumber() == 2331 ) nskip =  3;

    if( event->getRunNumber() ==   77 ) nskip =  0;
    if( event->getRunNumber() ==   81 ) nskip =294;

    if( event->getRunNumber() ==   92 ) nskip =  0;

    if( event->getRunNumber() ==  125 ) nskip =  0;

    if( event->getRunNumber() ==  133 ) nskip =  0;
    if( event->getRunNumber() ==  134 ) nskip =  0;

    if( event->getRunNumber() ==  143 ) nskip =  0;
    if( event->getRunNumber() ==  146 ) nskip =  0;
    if( event->getRunNumber() ==  147 ) nskip =  0;

    if( event->getRunNumber() ==  158 ) nskip =  0;

    if( event->getRunNumber() ==  164 ) nskip =  0;
    if( event->getRunNumber() ==  166 ) nskip =  0;
    if( event->getRunNumber() ==  169 ) nskip =  0;
    if( event->getRunNumber() ==  171 ) nskip =  0;
    if( event->getRunNumber() ==  173 ) nskip =  0;
    if( event->getRunNumber() ==  175 ) nskip =  0;
    if( event->getRunNumber() ==  177 ) nskip =  0;
    if( event->getRunNumber() ==  179 ) nskip =  0;

    if( _nSkip > -1 ) nskip = _nSkip; // from input option

    cout << "skip first " << nskip << " CMS events" << endl;

    if( _dutReader.isOpen() ){
      for( int iskip = 0; iskip < nskip; ++iskip ){
	ClustDUT = _dutReader.getCMSclust( dutTimeStamp, NpixDUT ); // skip first CMS events
      }
    }

    if( _refReader.isOpen() ){
      for( int iskip = 0; iskip < nskip; ++iskip ){
	ClustREF = _refReader.getCMSclust( refTimeStamp, NpixREF ); // skip first CMS events
      }
    }

  }//first

  if( lprnt ) {
    
    cout << setw(6) << setiosflags(ios::right) << event->getEventNumber();
    cout << ". DUT pixels " << NpixDUT;
    if( NpixDUT > 0 ) {
      cout << ", clusters at";
      for( vector<cluster>::iterator c = ClustDUT.begin(); c != ClustDUT.end(); c++ ){
	cout << "  (" << c->col << ", " << c->row << ")";
      }
    }
    cout << endl;

  }//lprnt

  static unsigned long dutTimeStamp0 = dutTimeStamp; // init
  static unsigned long prevdutTimeStamp = dutTimeStamp; // init

  double dtns = ( event->getTimeStamp() - prevtime ) /gTLU - ( dutTimeStamp - prevdutTimeStamp) / 0.040; // +-40 ns

  if( _dutReader.isOpen() ) {

    if( dutTimeStamp > dutTimeStamp0 ) {

      cmsdtHisto->fill( ( dutTimeStamp - prevdutTimeStamp ) / 40E0 );
      sysrtHisto->fill( (event->getTimeStamp()-time0)/fTLU / (dutTimeStamp-dutTimeStamp0)*40E6 );  // mean = 1.00007
      sysrdtHisto->fill( (event->getTimeStamp()-prevtime)/fTLU / (dutTimeStamp-prevdutTimeStamp)*40E6 );

      // time between events, TLU - DUT:

      sysdtHisto->fill( dtns );  // +-40 ns

      sysdtvsdt->fill( (event->getTimeStamp()-time0)/gTLU/1E3, dtns );

    }// > time0

  }//DUT open

  if( NpixDUT > 0 ) {

    for( vector<cluster>::iterator c = ClustDUT.begin(); c != ClustDUT.end(); c++ ){

      nclst++;
      cmscolHisto->fill( c->col );
      cmsrowHisto->fill( c->row );
      cmsnpxHisto->fill( c->size );
      cmsadcHisto->fill( c->charge );

    }//DUT clusters

    if( _nEvt > 0 ) {
      logdtcmsHisto->fill( std::log( ( event->getTimeStamp() - prevtime ) / gTLU/1E6 ) / std::log(10.0) );
    }

  }//pix

  if( NpixREF > 0 ) {

    for( vector<cluster>::iterator c = ClustREF.begin(); c != ClustREF.end(); c++ ){

      nclst++;
      refcolHisto->fill( c->col );
      refrowHisto->fill( c->row );
      refnpxHisto->fill( c->size );
      refadcHisto->fill( c->charge );

    }//REF clusters

  }//pixREF

  prevdutTimeStamp = dutTimeStamp;

  _nEvt++;

  prevtime = event->getTimeStamp(); // remember

  //----------------------------------------------------------------------------

  if( _isFirstEvent ) {
 
    // apply all GEAR/alignment offsets to get corrected X,Y,Z position of the
    // sensor center

    for( int iplane = 0; iplane < _siPlanesLayerLayout->getNLayers(); iplane++ ) {

      printf(":\n");

      map< unsigned int , double > _planeCenter;
      map< unsigned int , double > _planeNormal;
 
      int sensorID = _siPlanesLayerLayout->getID(iplane);

      // start filling the map with Gear values:

      _planeCenter[ 0 ] =  _siPlanesLayerLayout->getLayerPositionX(iplane); // X
      _planeCenter[ 1 ] =  _siPlanesLayerLayout->getLayerPositionY(iplane); // Y
      _planeCenter[ 2 ] =  _siPlanesLayerLayout->getLayerPositionZ(iplane); // Z

      _planeNormal[ 0 ] =  0.; // X
      _planeNormal[ 1 ] =  0.; // Y
      _planeNormal[ 2 ] =  1.; // Z

      TVector3  _normalTVec( _planeNormal[0], _planeNormal[1], _planeNormal[2]); 

      // do initial rotation( from GEAR)
      try{
 
	double gRotation[3] = { 0., 0., 0.}; // not rotated
                         
	gRotation[0] = _siPlanesLayerLayout->getLayerRotationXY(iplane); // Euler alpha
	gRotation[1] = _siPlanesLayerLayout->getLayerRotationZX(iplane); // Euler alpha
	gRotation[2] = _siPlanesLayerLayout->getLayerRotationZY(iplane); // Euler alpha
                          
	// input angles are in DEGREEs !!!
	// translate into radians
	gRotation[0] =  gRotation[0]*3.1415926/180.; // 
	gRotation[1] =  gRotation[1]*3.1415926/180.; //
	gRotation[2] =  gRotation[2]*3.1415926/180.; //

	TRotation	r;
	r.RotateX( gRotation[2] );                          
	r.RotateY( gRotation[1] );                          
	r.RotateZ( gRotation[0] );                          

	_normalTVec.Transform( r );
                                  
	_planeNormal[0] = _normalTVec[0];
	_planeNormal[1] = _normalTVec[1];
	_planeNormal[2] = _normalTVec[2];
      }
      catch(...) {
	printf(" no sensor rotation is given in the GEAR steering file, assume NONE \n" );
      }
 
      printf("sensorID: %5d    %8.3f %8.3f %8.3f \n", sensorID, _planeCenter[ 0 ], _planeCenter[ 1 ], _planeCenter[ 2 ] );
             
      if( _alignmentCollectionNames.size() > 0 ) {

	for( unsigned int i = 0; i < _alignmentCollectionNames.size(); i++ ) {

	  // if offsets are given via the alignment collections apply:

	  try{
	    LCCollectionVec * alignmentCollectionVec = dynamic_cast < LCCollectionVec * >( event->getCollection( _alignmentCollectionNames[i] ) );
	    // next, find the alignment constant corresponding to the DUT

	    EUTelAlignmentConstant * c = NULL;
	    for( size_t iPos = 0; iPos < alignmentCollectionVec->size(); ++iPos ) {
	      c = static_cast< EUTelAlignmentConstant * >( alignmentCollectionVec->getElementAt( iPos ) );
	      if( c -> getSensorID() == sensorID  ) {                
		break;
	      }
	    }

	    _planeCenter[ 0 ] += c->getXOffset();
	    _planeCenter[ 1 ] += c->getYOffset();
	    _planeCenter[ 2 ] += c->getZOffset();

	    double gRotation[3] = { 0., 0., 0.}; // not rotated

	    try{
	      gRotation[0] = c->getGamma();
	      gRotation[1] = c->getBeta();
	      gRotation[2] = c->getAlpha();

	      TRotation	r;
	      r.RotateX( gRotation[2] );
	      r.RotateY( gRotation[1] );
	      r.RotateZ( gRotation[0] );

	      _normalTVec.Transform( r );

	    }
	    catch(...) {
	      printf(" no sensor rotation is given in the GEAR steering file, assume NONE \n" );
	    }

	    _planeNormal[0] = _normalTVec[0];
	    _planeNormal[1] = _normalTVec[1];
	    _planeNormal[2] = _normalTVec[2];
                      
	    printf("%-20s %5d( %8.3f %8.3f %8.3f) %8.3f %8.3f %8.3f  normal( %8.3f %8.3f %8.3f) %8.3f %8.3f %8.3f \n", 
		   _alignmentCollectionNames[i].c_str(), sensorID, 
		   c->getXOffset(), c->getYOffset(), c->getZOffset(),
		   _planeCenter[0], _planeCenter[1], _planeCenter[2],
		   gRotation[0], gRotation[1], gRotation[2], 
		   _planeNormal[0], _planeNormal[1], _planeNormal[2] 
		   );
	  }
	  catch(...) {
	    printf("collection %s not found \n", _alignmentCollectionNames[i].c_str() );
	  }

	  printf(":\n");
	}

      }
 
      _siPlaneCenter[ sensorID]=  _planeCenter;
      _siPlaneNormal[ sensorID]=  _planeNormal;
      //         _siPlaneCenter.insert( make_pair( sensorID, _planeCenter ));
      //          _siPlaneNormal.insert( make_pair( sensorID, _planeNormal ));
    }

    if( isFirstEvent() ) _isFirstEvent = false;

  }//isFirstEvent

  //----------------------------------------------------------------------------
  // check input collection (aligned hits):

  LCCollection* col;
  try {
    col = event->getCollection( _inputColName );
  }
  catch( lcio::DataNotAvailableException& e) {
    streamlog_out( DEBUG ) << "Not able to get collection "
			   << _inputColName
			   << "\nfrom event " << event->getEventNumber()
			   << " in run " << event->getRunNumber()  << endl;

#if defined(USE_AIDA) || defined(MARLIN_USE_AIDA)
   ( dynamic_cast<AIDA::IHistogram1D*>( _aidaHistoMap[_nTrackHistoName]))->fill(0);
#endif
    ++_noOfEventWOInputHit;
    return;
  }

  //----------------------------------------------------------------------------
  // Copy hits to local table
  // Assign hits to sensor planes

  int nHit = col->getNumberOfElements();  //  from event collection

  if( debug ) {
    streamlog_out( TESTFITTERMESSAGE )  << "Total of " << nHit
					<< " tracker hits in input collection " << endl;
  }

#if defined(USE_AIDA) || defined(MARLIN_USE_AIDA)
  ( dynamic_cast<AIDA::IHistogram1D*>( _aidaHistoMap[_nAllHitHistoName]))->fill(nHit);
#endif

  if( nHit + _allowMissingHits < _nActivePlanes ) {

    streamlog_out( TESTFITTERMESSAGE ) << "Not enough hits to perform the fit, exiting... " << endl;

#if defined(USE_AIDA) || defined(MARLIN_USE_AIDA)
    ( dynamic_cast<AIDA::IHistogram1D*>( _aidaHistoMap[_nTrackHistoName]))->fill(0);
#endif

    _noOfEventWOTrack++;
    return;

  }

  //----------------------------------------------------------------------------

  double * hitX  = new double[nHit];
  double * hitEx = new double[nHit];
  double * hitY  = new double[nHit];
  double * hitEy = new double[nHit];
  double * hitZ  = new double[nHit];
  int    * hitPlane = new int[nHit];
  int    * hitFits = new int[nHit];

  IntVec * planeHitID = new IntVec[_nTelPlanes];

  // Loop over telescope hits:

  int nGoodHit = 0;

  for( int ihit = 0; ihit < nHit; ihit++ ) {

    TrackerHit * meshit = dynamic_cast<TrackerHit*>( col->getElementAt(ihit) );

    // Telescope hit position:

    const double * pos = meshit->getPosition();

    hitZ[ihit] = pos[2];

    hitFits[ihit] = -1;

    // We have to find Plane ID of the hit by looking at the Z position

    double distMin = 99; // [mm]
    hitPlane[ihit] = -1;

    for( int ipl = 0; ipl < _nTelPlanes; ipl++ ) {

      double dist =  hitZ[ihit] - _planePosition[ipl];
      if( dist < 0 ) dist = -dist; // always positive defined !

      if( dist < distMin ) {
        hitPlane[ihit] = ipl;
        distMin        = dist;
      }
    }

    // Ignore hits not matched to any plane:

    if( hitPlane[ihit] < 0 ) {
      if( _isActive[hitPlane[ihit]] ) {
	streamlog_out( WARNING0 )  << "Reconstructed hit outside sensor planes z [mm] = "  << hitZ[ihit] << endl;
      }
      continue;
    }

    // Ignore hit, if plane not declared as active (i.e. not used in the fit):

    if( ! _isActive[hitPlane[ihit]] ) {
      hitPlane[ihit] = -1;
      continue;
    }

    // Ignore hit also, if maximum number of hits already matched to this plane:

    if( _maxPlaneHits>0 && _maxPlaneHits - planeHitID[hitPlane[ihit]].size() <= 0 ) {
      hitPlane[ihit] = -1;
      continue;
    }

    // Hit will be used: correct X and Y position for extra user plane alignment:

    hitX[ihit]
      = pos[0] * cos(_planeRotZ[hitPlane[ihit]] )
      + pos[1] * sin(_planeRotZ[hitPlane[ihit]] )
      - _planeShiftX[hitPlane[ihit]];

    hitY[ihit]
      = pos[1] * cos(_planeRotZ[hitPlane[ihit]] )
      - pos[0] * sin(_planeRotZ[hitPlane[ihit]] )
      - _planeShiftY[hitPlane[ihit]];

    // DP: trust millepede:

    hitX[ihit] = pos[0];
    hitY[ihit] = pos[1];

    // Check Window and Mask cuts, if defined
    bool hitcut = false;

    for( int wpl = 0; wpl < ( int)_planeWindowIDs[hitPlane[ihit]].size(); wpl++ ) {
      int iwin = _planeWindowIDs[hitPlane[ihit]].at(wpl);

      if( hitX[ihit] <  _WindowMinX.at(iwin) ) hitcut = true;
      if( hitX[ihit] >  _WindowMaxX.at(iwin) ) hitcut = true;
      if( hitY[ihit] <  _WindowMinY.at(iwin) ) hitcut = true;
      if( hitY[ihit] >  _WindowMaxY.at(iwin) ) hitcut = true;
    }

    if( hitcut ) {
      hitPlane[ihit] = -1;
      continue;
    }

    for( int mpl = 0; mpl<( int)_planeMaskIDs[hitPlane[ihit]].size(); mpl++) {
      int imsk= _planeMaskIDs[hitPlane[ihit]].at(mpl);

      if(
	 hitX[ihit] >  _MaskMinX.at(imsk) &&  hitX[ihit] <  _MaskMaxX.at(imsk) &&
	 hitY[ihit] >  _MaskMinY.at(imsk) &&  hitY[ihit] <  _MaskMaxY.at(imsk)
	 ) {
        hitcut = true;
      }
    }

    if( hitcut ) {
      hitPlane[ihit] = -1;
      continue;
    }

    // Add hit to hit list for given plane - to be used in track selection:

    planeHitID[hitPlane[ihit]].push_back(ihit);
    nGoodHit++;

    // Position uncertainty. Use nominal resolution if not properly defined:

    const EVENT::FloatVec cov = meshit->getCovMatrix();

    if( cov.at(0) > 0. ) {
      hitEx[ihit]=sqrt(cov.at(0));
    }
    else {
      hitEx[ihit]=_planeResolution[hitPlane[ihit]];
    }
    if( cov.at(2) > 0. ) {
      hitEy[ihit]=sqrt(cov.at(2));
    }
    else {
      hitEy[ihit]=_planeResolution[hitPlane[ihit]];
    }

    hitFits[ihit] = 0;

    if( debug ) {
      streamlog_out( TESTFITTERMESSAGE ) << "Hit " << ihit
					 << "   X = " << hitX[ihit] << " +/- " << hitEx[ihit]
					 << "   Y = " << hitY[ihit] << " +/- " << hitEy[ihit]
					 << "   Z = " << hitZ[ihit] << "( plane " << hitPlane[ihit] << ")" << endl;
    }

  }//loop over hits

  if( debug ) {
    cout << "event " << event->getEventNumber();
    cout << ", hits " << nHit;
    cout << ", good " << nGoodHit;
    cout << endl;
  }

  if( lprnt ) {
    cout << "                                           ";
  }

#if defined(USE_AIDA) || defined(MARLIN_USE_AIDA)
  ( dynamic_cast<AIDA::IHistogram1D*>( _aidaHistoMap[_nAccHitHistoName]))->fill(nGoodHit);
#endif

  hits0Histo->fill( planeHitID[0].size() );
  hits1Histo->fill( planeHitID[1].size() );
  hits2Histo->fill( planeHitID[2].size() );
  hits3Histo->fill( planeHitID[3].size() );
  hits4Histo->fill( planeHitID[4].size() );
  hits5Histo->fill( planeHitID[5].size() );

  //----------------------------------------------------------------------------
  // loop over telescope hit pairs:

  for( int ihit = 0; ihit < nHit; ihit++ ) {

    int ipl = hitPlane[ihit];

    for( int jhit = 0; jhit < nHit; jhit++ ) {

      int jpl = hitPlane[jhit];

      double dx = hitX[jhit] - hitX[ihit];
      double dy = hitY[jhit] - hitY[ihit];

#if defined(USE_AIDA) || defined(MARLIN_USE_AIDA)

      if( ipl == 0 ) {

	if( jpl == 1 ) {
	  dx01Histo->fill( dx );
	  dy01Histo->fill( dy );
	  if( abs(dy) < 1 ) du01Histo->fill( dx );
	}
	if( jpl == 2 ) {
	  dx02Histo->fill( dx );
	}
	if( jpl == 3 ) {
	  dx03Histo->fill( dx );
	}
	if( jpl == 4 ) {
	  dx04Histo->fill( dx );
	}
	if( jpl == 5 ) {
	  dx05Histo->fill( dx );
	}

      }//ipl==0

      if( ipl == 1 ) {

	if( jpl == 2 ) {
	  dx12Histo->fill( dx );
	  dy12Histo->fill( dy );
	  if( abs(dy) < 1 ) du12Histo->fill( dx );
	}
      }

      if( ipl == 2 ) {

	if( jpl == 3 ) {
	  dx23Histo->fill( dx );
	  dy23Histo->fill( dy );
	  if( abs(dy) < 1 ) du23Histo->fill( dx );
	}
      }

      if( ipl == 3 ) {

	if( jpl == 4 ) {
	  dx34Histo->fill( dx );
	  dy34Histo->fill( dy );
	  if( abs(dy) < 1 ) du34Histo->fill( dx );
	}
      }

      if( ipl == 4 ) {

	if( jpl == 5 ) {
	  dx45Histo->fill( dx );
	  dy45Histo->fill( dy );
	  if( abs(dy) < 1 ) du45Histo->fill( dx );
	}
      }

#endif

    }//loop over telescope hits

  }// loop over telescope hits

  //----------------------------------------------------------------------------
  // DUT alignment: relative to _planePosition[2]
  // REF alignment: relative to _planePosition[5]

  bool lDUT = 1; // DUT present
  double tilt = 0;
  double DUTz = 20 + _planePosition[2]; // default
  double DUTalignx =  3.33; //align 2153 Wed 14.12.
  double DUTaligny = -0.90; //align 2153 Wed 14.12.
  double DUTrot = 0.0;
  double wt = atan(1.0) / 45.0; // pi/180 deg
  double DUTX0 = 1; // init

  // Fri 10.2.2012:
  double REFz = 430 + _planePosition[5]; // timing reference plane
  double REFalignx =  2.863; // [mm] run 133 Fri 10.2.2012
  double REFaligny = -2.353; // [mm] run 133 Fri 10.2.2012
  double REFrot = 0.0;

  // Feb 2012:

  if( event->getRunNumber() ==   77 || event->getRunNumber() ==   78 || event->getRunNumber() ==   81 ) { // chip 11, 0 deg

    tilt =  0;
    DUTalignx =  3.687; //from cmssxa
    DUTaligny = -1.633; //from cmsdya
    DUTrot =-0.010; //from cmsdyvsx
    DUTz = 32 + _planePosition[2]; // from cmsdyvsty, was 28
  }

  // chip 10, 0 deg:

  if( event->getRunNumber() ==   92 || event->getRunNumber() ==   93 ) {
    tilt =  0;
    DUTalignx =  3.425; //from cmssxa
    DUTaligny = -1.660; //from cmsdya
    DUTrot =-0.018; //from cmsdyvsx
    DUTz = 36 + _planePosition[2]; // from cmsdyvsty, was 28
  }

  // chip 11, 20 deg:

  if( event->getRunNumber() ==  122 || event->getRunNumber() ==  125 ) {
    tilt = 20;
    DUTalignx =  3.086; //from cmssxa
    DUTaligny =  0.236; //from cmsdya
    DUTrot =-0.010; //from cmsdyvsx
    DUTz = 32 + _planePosition[2]; // from cmsdyvsty
  }

  // Fri 10.2.2012: chip 11, 20 deg

  if( event->getRunNumber() ==  133 || event->getRunNumber() ==  134 || event->getRunNumber() ==  135 ) {
    tilt = 20.2; // from cmsdyvsy
    DUTalignx =  3.064; //from cmssxa
    DUTaligny =  0.133; //from cmsdya
    DUTrot =-0.003; //from cmsdyvsx, was -0.000
    DUTz = 23 + _planePosition[2]; // from cmsdyvsty, was 32
  }

  // Mon 13.2.2012: chip 11, 20 deg

  if( event->getRunNumber() >= 137 ) { // Mon 13.2.2012
    DUTalignx =  3.372; //from cmssxa
    DUTaligny =  0.118; //from cmsdya, cmsdyfcntq
    DUTrot =-0.003; //from cmsdyvsx
    DUTz = 25 + _planePosition[2]; // from cmsdyvsty, was 23
    tilt = 20.2; // from cmsdyvsy

    REFalignx =  2.643; // [mm] run 143 Mon 13.2.2012, refsxa
    REFaligny =  0.827; // [mm] run 143 Mon 13.2.2012, refdya
    REFz = 197 + _planePosition[5];
  }

  if( event->getRunNumber() >= 145 ) { // Mon 13.2.2012
    DUTalignx =  3.350; //from cmssxa
    DUTaligny =  0.130; //from cmsdya, cmsdyfcntq
    DUTrot =-0.003; //from cmsdyvsx

    REFalignx =  3.487; // [mm]  refsxa
    REFaligny = -0.649; // [mm]  refdya
  }

  if( event->getRunNumber() >= 147 ) { // Mon 13.2.2012
    REFalignx =  3.418; // [mm]  refsxa
    REFaligny = -0.823; // [mm]  refdya
  }

  // Mon 20.2.2012: eff

  if( event->getRunNumber() == 158 ) {
    DUTalignx =  3.344; //from cmssxa
    DUTaligny =  0.230; //from cmsdya, cmsdyfcntq
    DUTrot = -0.002; //from cmsdyvsx
    tilt = 20.2; // from cmsdyvsy, was 20.0
    DUTz = 22 + _planePosition[2]; // from cmsdyvsty

    REFalignx =  3.444; // [mm], refsxa
    REFaligny = -0.705; // [mm], refdya
    REFrot = 0.011; //from refdyvsx
  }

  // eff runs from Tue 21.2. with alignment from 177:

  if( event->getRunNumber() >= 164 ) {
    DUTalignx =  3.301; //from cmssxa
    DUTaligny =  0.212; //from cmsdya, cmsdyfcntq
    DUTaligny =  0.196; //from cmsdya, cmsdyfcntq
    DUTrot = -0.003; //from cmsdyvsx
    tilt = 19.8; // from cmsdyvsy
    DUTz = 25 + _planePosition[2]; // from cmsdyvsty, was 24
    //DUTz += 2.2689; // pede

    //DUTalignx += 0.150; // test
    //DUTaligny += 0.050; // test
    //DUTz += 4.0; // test
    //DUTrot += 0.005; // test
    //tilt += 2.0; // test

    REFalignx =  3.204; // [mm], refsxa
    REFaligny = -0.601; // [mm], refdya
  }

  // scattering studies without CMS pixel:

  if( event->getRunNumber() >= 192 ) { // empty telescope
    DUTz = 0.5 * ( _planePosition[2] + _planePosition[3] ); // mid air
    lDUT = 0; // no DUT
  }

  if( event->getRunNumber() >= 202 ) { // 1 mm Al plate
    DUTz = 71 + _planePosition[2];
    lDUT = 1;
    DUTX0 = 0.7 / 88.9; // from GBL/ex150
    tilt = 0;
  }

  if( event->getRunNumber() >= 217 ) { // 1 mm Cu plate
    DUTz = 69 + _planePosition[2];
    lDUT = 1;
    DUTX0 = 0.75 / 14.35; // from GBL/ex150
    DUTX0 = 0.60 / 14.35; // for flat prob
    tilt = 0;
  }

  if( event->getRunNumber() >= 230 ) { // 1 mm Cu plate
    DUTz = 84 + _planePosition[2];
    lDUT = 1;
    DUTX0 = 0.70 * 1.0 / 14.35; // from GBL/ex38
    tilt = 0;
  }

  if( event->getRunNumber() >= 255 ) { // 0.4 mm Cu plate
    DUTz = 84 + _planePosition[2];
    lDUT = 1;
    DUTX0 = 0.68 * 0.4 / 14.35; // from GBL/ex38
    tilt = 0;
  }

  if( event->getRunNumber() >= 263 ) { // 4 mm Al plate
    DUTz = 85 + _planePosition[2];
    lDUT = 1;
    DUTX0 = 0.65 * 4.0 / 88.9; // from GBL/ex38
    tilt = 0;
  }

  // Dec 2011:

  if( event->getRunNumber() == 1946 ) { //Mon 12.12.
    DUTalignx =  0.641;
    DUTaligny = -0.772;
    lDUT = 1;
  }

  if( event->getRunNumber() == 2257 ) {
    tilt = 10;
  }

  if( event->getRunNumber() == 2258 ) {
    tilt = 10;
  }

  if( event->getRunNumber() == 2261 ) {
    tilt = 10;
    DUTalignx =  3.685;
    DUTaligny = -2.75 ;
  }

  if( event->getRunNumber() == 2265 || event->getRunNumber() == 2266 ) {
    tilt = 20;
    DUTalignx =  3.658; //align 2266 Thu 16.12.
    DUTaligny = -4.445; //align 2266 Thu 16.12.
  }

  if( event->getRunNumber() == 2267 ) {
    tilt = 15;
  }

  if( event->getRunNumber() == 2269 ) {
    tilt = 30;
  }

  if( event->getRunNumber() == 2293 ) {
    tilt = 20;
    DUTalignx =  3.689; //align 2293 Fri 16.12.
    DUTaligny = -4.491; //align 2293 Fri 16.12.
  }

  if( event->getRunNumber() == 2295 ) {
    tilt = 20;
  }

  if( event->getRunNumber() == 2314 ) {
    tilt = 20;
  }

  if( event->getRunNumber() == 2322 ) { // chip 7, 20 deg
    tilt = 19.9; //was 20.0, 19.8
    //tilt = 20.5; //test
    DUTalignx =  3.299; //from cmssxa
    DUTaligny = -0.844; //from cmsdya
    DUTrot = 0.0030; //from cmsdyvsx
    //DUTrot = 0.0020; //test
    DUTz = 22 + _planePosition[2]; // from cmsdyvsty
    //DUTz = 25 + _planePosition[2]; // test
  }

  if( event->getRunNumber() == 2325 ) { // chip 7, 25 deg
    tilt = 24.4; // cmsdyvsy
    DUTalignx =  3.232; //from cmssxa
    DUTaligny = -0.629; //from cmsdya. was -0.631
    DUTrot =-0.0027; //from cmsdyvsx
    DUTz = 22 + _planePosition[2]; // from cmsdyvsty
  }

  if( event->getRunNumber() == 2326 ) { // chip 7, 30 deg
    tilt = 29.5;
    DUTalignx =  3.203; //from cmssxa
    DUTaligny = -1.451; //from cmsdya, was -1.454
    DUTrot =-0.0035; //from cmsdyvsx
    DUTz = 22 + _planePosition[2]; // from cmsdyvsty
  }

  if( event->getRunNumber() == 2328 ) { // chip 7, 15 deg
    tilt = 16.0; // from cmsdyvsy
    DUTalignx =  3.096; //from cmssxa
    DUTaligny =  0.804; //from cmsdya
    DUTrot =-0.0018; //from cmsdyvsx
    DUTz = 20 + _planePosition[2]; // from cmsdyvsty
  }

  if( event->getRunNumber() == 2330 ) { // chip 7, 10 deg
    tilt = 10.5;
    DUTalignx =  3.102; //from cmssxa
    DUTaligny =  1.841; //from cmsdya
    DUTrot =-0.0010; //from cmsdyvsx
    DUTz = 20 + _planePosition[2]; // from cmsdyvsty
  }

  if( event->getRunNumber() == 2331 ) { // chip 7, 0 deg
    tilt =  0;
    DUTalignx =  3.081; //from cmssxa
    DUTaligny =  3.658; //from cmsdya
    DUTrot =-0.0007; //from cmsdyvsx, was -0.0012
    DUTz = 20 + _planePosition[2]; // from cmsdyvsty
  }

  // Mar 2012 TB 22 EUDET

  if( event->getRunNumber() == 2415 ||
      event->getRunNumber() == 2416 ||
      event->getRunNumber() == 2419 ||
      event->getRunNumber() == 2438 ||
      event->getRunNumber() == 2439 ||
      event->getRunNumber() == 2440 ) {
    lDUT = 0; // no DUT, empty telescope
  }

  // Apr 2012, TB 22:

  if( event->getRunNumber() == 2844 ) { // chip 10, 0 deg
    tilt =  0;
    DUTalignx = -5.302; //from cmssxa
    DUTaligny =  1.994; //from cmsdya
    DUTrot = 0.0025; //from cmsdyvsx
    DUTz = 28 + _planePosition[2]; // from cmsdyvsty

    REFalignx = -3.553; // [mm], refsxa
    REFaligny =  2.469; // [mm], refdya
    REFrot = 0.000; //from refdyvsx
    REFz = 105 + _planePosition[5];
  }

  if( event->getRunNumber() == 2878 ) { // chip 10, 20 deg
    DUTalignx = -2.350; //from cmssxa
    DUTaligny = -0.405; //from cmsdya
    DUTrot = -0.0042; //from cmsdyvsx
    tilt = 19.6;
    DUTz = 30 + _planePosition[2]; // from cmsdyvsty

    REFalignx = -0.894; // [mm], refsxa
    REFaligny = -1.870; // [mm], refdya
    REFrot = 0.024; //from refdyvsx
    REFz = 105 + _planePosition[5];
  }

  if( event->getRunNumber() == 2880 ) { // chip 10, 20 deg
    DUTalignx = -2.373; // from cmssxa
    DUTaligny = -0.427; // from cmsdya
    DUTrot = -0.0045; // from cmsdyvsx
    tilt = 19.5; // from cmsdyvsy
    DUTz = 30 + _planePosition[2]; // from cmsdyvsty, was 32

    REFalignx = -0.993; // [mm], refsxa
    REFaligny = -1.930; // [mm], refdya
    REFrot = 0.024; //from refdyvsx
    REFz = 105 + _planePosition[5];
  }

  double pitchy = 0.10 * cos( tilt * wt ); // [mm] in telescope system
  double pitchx = 0.15; //[mm]

  if( DUTX0 > 0.99 ) { // not yet set
    //DUTX0 = 0.08 / cos( tilt * wt ); // CMS Pixel thickness
    DUTX0 = 0.07 / cos( tilt * wt ); // for flat prob in run 158
  }

  //----------------------------------------------------------------------------
  // 3-4-5 downstream driplets:

  double driCut = 0.1; // [mm]

  int ndri = 0;
  int ndricms = 0;
  double xmB[99];
  double sxB[99];
  double ymB[99];
  double syB[99];
  double zmB = 0.5 * ( _planePosition[3] + _planePosition[5] );
  double dzB = _planePosition[5] - _planePosition[3];
  bool lkB[99];
  int hts[6][99]; // 6 planes

  for( int ihit = 0; ihit < nHit; ihit++ ) {

    if( hitPlane[ihit] != 3 ) continue; // want plane 3: noisy

    for( int jhit = 0; jhit < nHit; jhit++ ) {

      if( hitPlane[jhit] != 5 ) continue; //want last plane 5

      double dx35 = hitX[jhit] - hitX[ihit];
      double dy35 = hitY[jhit] - hitY[ihit];

      if( abs(dy35) < 1 ) dx35Histo->fill( dx35 );
      if( abs(dx35) < 1 ) dy35Histo->fill( dy35 );

      if( abs(dy35) > 0.010 * dzB ) continue; // [mm]
      if( abs(dx35) > 0.010 * dzB ) continue; // cuts on track angle:+-10 mrad

      double avx = 0.5 * ( hitX[ihit] + hitX[jhit] );
      double avy = 0.5 * ( hitY[ihit] + hitY[jhit] );
      double avz = 0.5 * ( hitZ[jhit] + hitZ[ihit] ); // mid z

      double tx = ( hitX[jhit] - hitX[ihit] ) / dzB; // angle theta x
      double ty = ( hitY[jhit] - hitY[ihit] ) / dzB; // angle theta y

      //extrapolate to CMS timing plane: [mm]

      double zs = REFz - avz; // z REF from mid of driplet
      double xs = avx + tx * zs; // driplet impact point on REF
      double ys = avy + ty * zs;

      double xmod = fmod( 20 + xs, 2*0.15 ) * 1E3; // [0,300] um, 2 pixel wide
      double ymod = fmod( 20 + ys, 2*0.10 ) * 1E3; // [0,200] um

      // middle plane k for driplets:

      for( int khit = 0; khit < nHit; khit++ ) {

	if( hitPlane[khit] != 4 ) continue; //want middle plane 4

	// driplet residual: equidistant planes

	double dx = hitX[khit] - avx;
	double dy = hitY[khit] - avy;

	if( abs(dy) < driCut ) dridxHisto->fill( dx*1E3 ); //sig = 7 um at 5 GeV
	if( abs(dx) < driCut ) dridyHisto->fill( dy*1E3 );

	// alignment check:

	if( abs(dy) < driCut ) {
	  dridxvsx->fill( avx, dx*1E3 ); // check for rot
	  dridxvsy->fill( avy, dx*1E3 );
	  dridxvstx->fill( tx*1E3, dx*1E3 ); // check for z shift
	  dridxvsty->fill( ty*1E3, dx*1E3 );
	}
	if( abs(dx) < driCut ) {
	  dridyvsx->fill( avx, dy*1E3 );
	  dridyvsy->fill( avy, dy*1E3 );
	  dridyvstx->fill( tx*1E3, dy*1E3 );
	  dridyvsty->fill( ty*1E3, dy*1E3 );
	}

	// driplet found:

	if( abs(dx) < driCut && abs(dy) < driCut ){

	  if( ndri < 99 ) {
	    xmB[ndri] = avx;
	    ymB[ndri] = avy;
	    sxB[ndri] = tx;
	    syB[ndri] = ty;
	    lkB[ndri] = 0; //no link to CMS yet
	    hts[3][ndri] = ihit;
	    hts[5][ndri] = jhit;
	    hts[4][ndri] = khit;
	  }

	  ndri++;

	  drixHisto->fill( hitX[khit] );
	  driyHisto->fill( hitY[khit] );
	  drixyHisto->fill( hitX[khit], hitY[khit] );
	  dritxHisto->fill( tx*1E3 );
	  drityHisto->fill( ty*1E3 );

	  refxHisto->fill( xs );
	  refyHisto->fill( ys );
	  refxyHisto->fill( xs, ys );

	  // correlate with REF:

	  if( NpixREF > 0 ){

	    ndricms++;

	    for( vector<cluster>::iterator c = ClustREF.begin(); c != ClustREF.end(); c++ ){

	      refxxHisto->fill( c->col, xs ); // anti-correlation
	      refyyHisto->fill( c->row, ys ); // correlation

	      double refx = ( c->col - 26 ) * 0.15; // -3.9..3.9 mm
	      double refy = ( c->row - 40 ) * 0.10; // -4..4 mm

	      refsxaHisto->fill( refx + xs );
	      refdyaHisto->fill( refy - ys );

	      double refsx = refx + REFrot*refy - REFalignx + xs; // residual x
	      double refdy = refy - REFrot*refx - REFaligny - ys; // residual y

	      refsxHisto->fill( refsx );
	      refdyHisto->fill( refdy );

	      if( abs( refdy ) < 1 ) refsxcHisto->fill( refsx );
	      if( abs( refsx ) < 1 ) refdycHisto->fill( refdy );

	      if( abs( refdy ) < 1 && abs( refsx ) < 1 ) {

		drixcmsHisto->fill( hitX[ihit] );
		driycmsHisto->fill( hitY[ihit] );

		refxcmsHisto->fill( xs );
		refycmsHisto->fill( ys );
		refxycmsHisto->fill( xs, ys );
		refpixvsxy->fill( xmod, ymod ); // occupancy map

		refqHisto->fill( c->charge );
		refqvsxy->fill( xmod, ymod, c->charge ); // cluster charge profile

		refdyvsx->fill( refx, refdy*1E3 ); // dy vs x: rot?
		refdyvsy->fill( refy, refdy*1E3 ); // dy vs y: tilt?
		refdyvsty->fill( ty*1E3, refdy*1E3 ); // dy vs tety: z shift?

	      }

	      if( abs( refdy ) < 0.3 && abs( refsx ) < 0.3 ) {
		if( ndri < 99 ) lkB[ndri-1] = 1; // link to CMS
	      }

	    }//loop over CMS clusters

	  }//have CMS cluster

	  lkBvst->fill( (event->getTimeStamp()-time0)/fTLU, lkB[ndri-1] );

	}//have telescope driplet

	// correlate 3-4-5 driplet with DUT:

	if( NpixDUT > 0 && abs(dx) < driCut && abs(dy) < driCut ){

	  //extrapolate to DUT timing plane: [mm]

	  double zs = DUTz - avz; // z REF from mid of driplet
	  double xs = avx + tx * zs; // driplet impact point on REF
	  double ys = avy + ty * zs;

	  for( vector<cluster>::iterator c = ClustDUT.begin(); c != ClustDUT.end(); c++ ){

	    double bacx = ( c->col - 26 ) * pitchx; // -3.9..3.9 mm
	    double bacy = ( c->row - 40 ) * pitchy; // -4..4 mm

	    bacsxaHisto->fill( bacx + xs );
	    bacdyaHisto->fill( bacy - ys );

	    double bacsx = bacx + DUTrot*bacy - DUTalignx + xs; // residual x
	    double bacdy = bacy - DUTrot*bacx - DUTaligny - ys; // residual y

	    if( abs( bacdy ) < 1 ) bacsxcHisto->fill( bacsx*1E3 );
	    if( abs( bacsx ) < 1 ) bacdycHisto->fill( bacdy*1E3 );

	    if( c->charge > 18 ){
	      if( abs( bacdy ) < 1 ) bacsxcqHisto->fill( bacsx*1E3 );
	      if( abs( bacsx ) < 1 ) bacdycqHisto->fill( bacdy*1E3 ); // 18 um @ 5 GeV
	    }

	  }//loop over BAC clusters

	}//have telescope driplet

      }//loop over hits

    }//loop over hits

  }// loop over hits

  ndriHisto->fill( ndri );
  if( NpixDUT > 0 ) ndricmsHisto->fill( ndricms );

  //----------------------------------------------------------------------------
  // 0-1-2 triplets:

  double triCut = 0.1; // [mm]

  int ntri = 0;
  double xmA[99];
  double sxA[99];
  double ymA[99];
  double syA[99];
  double zmA = 0.5 * ( _planePosition[0] + _planePosition[2] );
  double dzA = _planePosition[2] - _planePosition[0];
  bool lkA[99];
  double DUTsxA[99];
  double DUTdyA[99];

  for( int ihit = 0; ihit < nHit; ihit++ ) { //hit in 0

    if( hitPlane[ihit] > 0 ) continue; //want plane 0

    for( int jhit = 0; jhit < nHit; jhit++ ) { // hit in 2

      if( hitPlane[jhit] != 2 ) continue; //want plane 2

      double dx02 = hitX[jhit] - hitX[ihit];
      double dy02 = hitY[jhit] - hitY[ihit];

      if( abs(dy02) < 1 ) da02Histo->fill( dx02 );
      if( abs(dx02) < 1 ) db02Histo->fill( dy02 );

      if( abs(dy02) > 0.005 * dzA ) continue; // [mm]
      if( abs(dx02) > 0.005 * dzA ) continue; // cuts on track angle:+-5 mrad

      double avx = 0.5 * ( hitX[ihit] + hitX[jhit] ); // average
      double avy = 0.5 * ( hitY[ihit] + hitY[jhit] ); // = mid point
      double avz = 0.5 * ( hitZ[jhit] + hitZ[ihit] ); // mid z

      double tx = ( hitX[jhit] - hitX[ihit] ) / dzA; // angle theta x
      double ty = ( hitY[jhit] - hitY[ihit] ) / dzA; // angle theta x

      //extrapolate to CMS: [mm]

      double zs = DUTz - avz; // z CMS from mid of triplet
      double xs = avx + tx * zs; // triplet impact point on CMS
      double ys = avy + ty * zs;
      ys = ys / ( 1.0 - ty * tan(tilt*wt) ); // tilt correction

      double xmod = fmod( 20 + xs, 2*pitchx ) * 1E3; // [0,300] um, 2 pixel wide
      double ymod = fmod( 20 + ys, 2*pitchy ) * 1E3; // [0,200] um

      // middle plane k for triplets:

      for( int khit = 0; khit < nHit; khit++ ) { //hit in 1

	if( hitPlane[khit] != 1 ) continue; //want middle plane

	// triplet residual: equidistant planes

	double dx = hitX[khit] - avx;
	double dy = hitY[khit] - avy;

	if( abs(dy) < triCut ) tridxHisto->fill( dx*1E3 );
	if( abs(dx) < triCut ) tridyHisto->fill( dy*1E3 );

	if( abs(dy) < triCut ) tridx1Histo->fill( dx*1E0 );
	if( abs(dx) < triCut ) tridy1Histo->fill( dy*1E0 );

	// alignment check:

	if( abs(dy) < triCut ) {
	  tridxvsx->fill( avx, dx*1E3 ); // check for rot
	  tridxvsy->fill( avy, dx*1E3 );
	  tridxvstx->fill( tx*1E3, dx*1E3 ); // check for z shift
	  tridxvsty->fill( ty*1E3, dx*1E3 );
	}
	if( abs(dx) < triCut ) {
	  tridyvsx->fill( avx, dy*1E3 );
	  tridyvsy->fill( avy, dy*1E3 );
	  tridyvstx->fill( tx*1E3, dy*1E3 );
	  tridyvsty->fill( ty*1E3, dy*1E3 );
	}

	// triplet found:

	bool ltri = 0;

	if( abs(dx) < triCut && abs(dy) < triCut ) {

	  if( ntri < 99 ) {
	    ltri = 1;
	    xmA[ntri] = avx;
	    ymA[ntri] = avy;
	    sxA[ntri] = tx;
	    syA[ntri] = ty;
	    lkA[ntri] = 0; //no link to CMS yet
	    hts[0][ntri] = ihit;
	    hts[2][ntri] = jhit;
	    hts[1][ntri] = khit;
	  }

	  ntri++;

	  trixHisto->fill( hitX[khit] );
	  triyHisto->fill( hitY[khit] );
	  trixyHisto->fill( hitX[khit], hitY[khit] );
	  tritxHisto->fill( tx*1E3 );
	  trityHisto->fill( ty*1E3 );

	  dutxHisto->fill( xs );
	  dutyHisto->fill( ys );
	  dutxyHisto->fill( xs, ys );

	  // scale hit to CMS pixel:
	  // sizeX="21.2"  sizeY="10.6"
	  // npixelX="1152"  npixelY="576" 

	  if( lprnt ) {
	    cout << "  x " << int( ( hitX[jhit] + 10.6 ) / 0.15 )
		 << ", y " << int( ( hitY[jhit] +  5.3 ) / 0.10 );
	  }

	  // extrapolate to planes 3,4,5: resolution study

	  for( int lhit = 0; lhit < nHit; lhit++ ) {

	    if( hitPlane[lhit] <= 2 ) continue; //want 3,4, or 5

	    double dz = hitZ[lhit] - avz;
	    double xs = avx + tx * dz; // triplet impact point on CMS
	    double ys = avy + ty * dz;
	    double dx = hitX[lhit] - xs;
	    double dy = hitY[lhit] - ys;

	    if( hitPlane[lhit] == 3 ) {
	      tridx3Histo->fill( dx*1E0 );
	      tridy3Histo->fill( dy*1E0 ); // 65 um at 4.7 GeV with CMS
	      tridx3bHisto->fill( dx*1E3 ); // finer binning
	      tridy3bHisto->fill( dy*1E3 ); // 
	    }
	    else if( hitPlane[lhit] == 4 ) {
	      tridx4Histo->fill( dx*1E0 );
	      tridy4Histo->fill( dy*1E0 ); //174 um at 4.7 GeV
	      tridx4bHisto->fill( dx*1E3 ); // finer binning
	      tridy4bHisto->fill( dy*1E3 ); // 
	    }
	    else if( hitPlane[lhit] == 5 ) {
	      tridx5Histo->fill( dx*1E0 );
	      tridy5Histo->fill( dy*1E0 ); //273 um at 4.7 GeV
	      tridx5bHisto->fill( dx*1E3 ); // finer binning
	      tridy5bHisto->fill( dy*1E3 ); // 
	    }

	  }//lhit in 3,4,5

	  //----------------------------------------------------------------------
	  // correlate with CMS DUT:
	  if( NpixDUT > 0 ){

	    // CMS pixel clusters:

	    for( vector<cluster>::iterator c = ClustDUT.begin(); c != ClustDUT.end(); c++ ){

	      cmsxxHisto->fill( c->col, xs ); // anti-correlation x-x
	      cmsyyHisto->fill( c->row, ys ); // correlation y-y

	      // pix in clus:

	      int colmin = 99;
	      int colmax = -1;
	      int rowmin = 99;
	      int rowmax = -1;

	      double arow[80];
	      for( int irow = 0; irow < 80; ++irow ) arow[irow] = 0;

	      for( vector<pixel>::iterator px = c->vpix.begin(); px != c->vpix.end(); px++ ){

		arow[px->row] += px->anaVcal; // project cluster onto rows
		if( px->col < colmin ) colmin = px->col;
		if( px->col > colmax ) colmax = px->col;
		if( px->row < rowmin ) rowmin = px->row;
		if( px->row > rowmax ) rowmax = px->row;

	      }//pix

	      bool fiducial = 1;
	      if( rowmin ==  0 ) fiducial = 0;
	      if( rowmax == 79 ) fiducial = 0;
	      if( colmin ==  0 ) fiducial = 0;
	      if( colmax == 51 ) fiducial = 0;
	      if( _dutReader.getChip() == 7 ) {
		if( colmax == 25 ) fiducial = 0; // 26 dead
		if( colmin == 28 ) fiducial = 0; // 27 dead
		if( colmax == 35 ) fiducial = 0; // 36 dead
		if( colmin == 38 ) fiducial = 0; // 37 dead
	      }

	      int ncol = colmax - colmin + 1;
	      int nrow = rowmax - rowmin + 1;

	      cmsncolHisto->fill( ncol );
	      cmsnrowHisto->fill( nrow );

	      double a1 = 0; //1st
	      double a2 = 0; //2nd
	      int i1 = 99;
	      int i2 = 99;

	      for( int irow = 0; irow < 80; ++irow ) {
		if( arow[irow] > a1 ) {
		  a2 = a1;
		  a1 = arow[irow];
		  i2 = i1;
		  i1 = irow;
		}
		else if( arow[irow] > a2 ) {
		  a2 = arow[irow];
		  i2 = irow;
		}
	      }
	      double a12 = a1 + a2;
	      double eta = 0;
	      if( a12 > 1 ) eta = ( a1 - a2 ) / a12;
	      if( i1 > i2 ) eta = -eta;

	      if( nrow == 2 ) cmsetaHisto->fill( eta );

	      double cmsx = ( c->col - 26 ) * pitchx; // -3.9..3.9 mm
	      double cmsy = ( c->row - 40 ) * pitchy; // -4..4 mm

	      cmssxaHisto->fill( cmsx + xs );
	      cmsdyaHisto->fill( cmsy - ys );

	      double cmssx = cmsx + DUTrot*cmsy - DUTalignx + xs; // residual x
	      double cmsdy = cmsy - DUTrot*cmsx - DUTaligny - ys; // residual y

	      cmssxHisto->fill( cmssx*1E3 );
	      cmsdyHisto->fill( cmsdy*1E3 );

	      if( fiducial ) {

		cmssxfHisto->fill( cmssx*1E3 );
		cmsdyfHisto->fill( cmsdy*1E3 );

		if( abs( cmsdy ) < 0.10 ) cmssxfcHisto->fill( cmssx*1E3 );
		if( abs( cmssx ) < 0.15 ) cmsdyfcHisto->fill( cmsdy*1E3 );

		if( nrow == 1 )
		  cmsdy1Histo->fill( cmsdy*1E3 );
		else if( nrow == 2 )
		  cmsdy2Histo->fill( cmsdy*1E3 );
		else
		  cmsdy3Histo->fill( cmsdy*1E3 );

	      }//CMS fiducial

	      // accumulate cuts:

	      if( fiducial && nrow < 3 ) {

		if( abs( cmsdy ) < 0.10 ) cmssxfcnHisto->fill( cmssx*1E3 );

		if( abs( cmssx ) < 0.15 ) {

		  cmsdyfcnHisto->fill( cmsdy*1E3 );

		  if(      c->charge < 18 )
		    cmsdyq0Histo->fill( cmsdy*1E3 );
		  else if( c->charge < 40 )
		    cmsdyq1Histo->fill( cmsdy*1E3 );
		  else
		    cmsdyq2Histo->fill( cmsdy*1E3 );

		  //if( abs( ty-0.00056 ) < 0.002 &&  abs( tx-0.00016 ) < 0.002 )
		  if( abs( ty-0.000 ) < 0.002 &&  abs( tx-0.000 ) < 0.002 )
		    cmsdyfcntHisto->fill( cmsdy*1E3 );

		  if( abs( ty-0.000 ) < 0.002 &&
		      abs( tx-0.000 ) < 0.002 &&
		      c->charge > 18 ){
		    cmsdyfcntqHisto->fill( cmsdy*1E3 ); // 9.2 um @ 6 GeV, 20 deg
		    if( c->charge * cos(tilt*wt) < 40 ) cmsdyfcntq1Histo->fill( cmsdy*1E3 ); // 8.8 um @ 6 GeV, 20 deg, more Gaussian
		  }
		}

	      }//CMS fiducial

	      // match CMS cluster and telescope track:

	      if( abs( cmssx ) < 0.3 && abs( cmsdy ) < 0.3 ){

		nmtch++;

		cmsqHisto->fill( c->charge );
		if( fiducial ) cmsqfHisto->fill( c->charge );
		cmsxHisto->fill( xs );
		cmsyHisto->fill( ys );
		cmsxyHisto->fill( xs, ys );
		cmsdyvsx->fill( cmsx, cmsdy*1E3 ); // dy vs x: rot?
		cmsdyvsy->fill( cmsy, cmsdy*1E3 ); // dy vs y: tilt?
		cmsdyvsty->fill( ty*1E3, cmsdy*1E3 ); // dy vs tety: z shift?
		cmsdyvsxHisto->fill( cmsx, cmsdy*1E3 ); // resolution map [um]

		if( fiducial ) {

		  cmspixvsxy->fill( xmod, ymod ); // occupancy map

		  cmsqvsxy->fill( xmod, ymod, c->charge ); // cluster charge profile
		  cmsqvsdt->fill( dtns, c->charge ); // cluster charge vs phase
		  if( dtns < -20 )
		    cmsqdt0Histo->fill( c->charge );
		  else if( dtns < -10 )
		    cmsqdt1Histo->fill( c->charge );
		  else if( dtns <   0 )
		    cmsqdt2Histo->fill( c->charge );
		  else if( dtns <  10 )
		    cmsqdt3Histo->fill( c->charge );
		  else if( dtns <  20 )
		    cmsqdt4Histo->fill( c->charge );
		  else
		    cmsqdt5Histo->fill( c->charge );

		  cmsrmsxvsxm->fill( xmod, abs(cmssx)*1E3 ); //resolution within pixel
		  cmsrmsyvsym->fill( ymod, abs(cmsdy)*1E3 ); //resolution within pixel
		  cmsncolvsxm->fill( xmod, ncol ); //within pixel
		  cmsnrowvsym->fill( ymod, nrow ); //within pixel
		  cmsetavsym->fill( ymod, eta ); //eta within pixel
		  if( nrow == 2 ) cmsetavsym2->fill( ymod, eta ); //eta within pixel
		  if( nrow == 1 ) cmsym1Histo->fill( ymod ); //where are 1-rows?

		}//fiducial

	      }//CMS - triplet match

	      if( ltri && abs( cmssx ) < 0.5 && abs( cmsdy ) < 0.5 ){ // link to CMS
		lkA[ntri-1] = 1; // ntri was already incremented
		DUTsxA[ntri-1] = cmssx;
		DUTdyA[ntri-1] = cmsdy;
	      }

	    } // loop over CMS clusters

	  }// have some CMS hit

	  lkAvst->fill( (event->getTimeStamp()-time0)/fTLU, lkA[ntri-1] );

	} // have telescope triplet

      }//loop over plane 1 hits k

    }//loop over plane 2 hits j

  }// loop over plane 0 hits i

  if( lprnt ) {
    cout << endl;
  }

  ntriHisto->fill( ntri );

  //----------------------------------------------------------------------------
  // telescope triplet efficiency w.r.t. CMS pixel:

  if( NpixDUT > 0 ){

    // CMS pixel clusters:

    for( vector<cluster>::iterator c = ClustDUT.begin(); c != ClustDUT.end(); c++ ){

      double cmsx = ( c->col - 26 ) * pitchx; // -3.9..3.9 mm
      double cmsy = ( c->row - 40 ) * pitchy; // -4..4 mm

      double cmsxe = cmsx + DUTrot*cmsy - DUTalignx; // aligned w.r.t. telescope
      double cmsye = cmsy - DUTrot*cmsx - DUTaligny; // 

      cmsxeHisto->fill( cmsxe );
      cmsyeHisto->fill( cmsye );

      // triplets:

      int  nm = 0;
      bool im = 0;

      for( int i = 0; i < ntri && i < 99; ++i ) {

	double xA = xmA[i] + sxA[i] * ( DUTz - zmA ); // A at DUT
	double yA = ymA[i] + syA[i] * ( DUTz - zmA );

	double sx = cmsxe + xA;
	double dy = cmsye - yA;

	cmssxeHisto->fill( sx*1E3 );
	cmsdyeHisto->fill( dy*1E3 );

	if( abs(sx) < 0.300 && abs(dy) < 0.200 ) {
	  nm++;
	  im = 1;
	}

      }//triplets

      cmsnmHisto->fill( nm );
      trieffvsxy->fill( cmsxe, cmsye, im ); //efficiency profile

    }//clusters

  }//have CMS data

  //----------------------------------------------------------------------------
  // kinks: triplets A vs driplets B
  // scattering point = DUT:
  // matching and GBL fit

  double sixCut = 0.1; // [mm]

  for( int iA = 0; iA < ntri && iA < 99; ++iA ) { // i = A = upstream

    double xA = xmA[iA] + sxA[iA] * ( DUTz - zmA ); // A at DUTz
    double yA = ymA[iA] + syA[iA] * ( DUTz - zmA );
    yA = yA / ( 1.0 - syA[iA] * tan( tilt*wt ) ); // tilt correction

    double xmod = fmod( 20 + xA, 2*pitchx ) * 1E3; // [0,300] um, 2 pixel wide
    double ymod = fmod( 20 + yA, 2*pitchy ) * 1E3; // [0,200] um

    for( int jB = 0; jB < ndri && jB < 99; ++jB ) { // j = B = downstream

      double kx = sxB[jB] - sxA[iA]; //kink
      double ky = syB[jB] - syA[iA];
      double k2 = kx*kx + ky*ky;
      double xB = xmB[jB] + sxB[jB] * ( DUTz - zmB ); //B at DUTz
      double yB = ymB[jB] + syB[jB] * ( DUTz - zmB );
      yB = yB / ( 1.0 - syB[jB] * tan( tilt*wt ) ); // tilt correction
      double dx = xB - xA;
      double dy = yB - yA;

      double xR = xmB[jB] + sxB[jB] * ( REFz - zmB ); //B at REFz
      double yR = ymB[jB] + syB[jB] * ( REFz - zmB );

      sixkxHisto->fill( kx*1E3 );
      sixkyHisto->fill( ky*1E3 );
      sixdxHisto->fill( dx );
      sixdyHisto->fill( dy );
      if( abs(dy) < 0.4 ) sixdxcHisto->fill( dx*1E3 ); // sig = 17 um at 5 GeV
      if( abs(dx) < 0.4 ) sixdycHisto->fill( dy*1E3 );

      double probchi = 0;

      // match driplet and triplet:

      if( abs(dx) < sixCut && abs(dy) < sixCut ) {

	sixkxcHisto->fill( kx*1E3 );
	sixkycHisto->fill( ky*1E3 );
	sixxyHisto->fill( xA, yA );

	// large kink map:

	if( abs( kx ) > 0.003 || abs( ky ) > 0.003 ) sixxycHisto->fill( xA, yA ); //islands

	kinkvsxy->fill( xA, yA, k2*1E6 ); //<kink^2> [mrad^2]

	// GBL with triplet A as seed:

	GblTrajectory traj( false ); // curvature = false

	// build up trajectory:

	vector<double> sPoint;

	// plane 0:

	double s = 0;

	TMatrixD proL2m(2,2);
	proL2m.UnitMatrix();

	TVectorD meas(2);

	double res = 3.5E-3; // [mm] Anemone telescope intrinsic resolution
	//res = 4.5E-3; // EUDET

	TVectorD measPrec(2); // precision = 1/resolution^2
	measPrec[0] = 1.0 / res / res;
	measPrec[1] = 1.0 / res / res;

	// scatter:
	TVectorD scat(2);
	scat.Zero(); //mean is zero

	double p = _eBeam; // beam momentum
	double X0Si = 65e-3 / 94; // Si + Kapton
	double tetSi = 0.0136 * sqrt(X0Si) / p * ( 1 + 0.038*std::log(X0Si) );

	TVectorD wscatSi(2);
	wscatSi[0] = 1.0 / ( tetSi * tetSi ); //weight
	wscatSi[1] = 1.0 / ( tetSi * tetSi );

	TMatrixD alDer( 2, 3 ); // alignment derivatives
	alDer[0][0] = 1.0; // dx/dx
	alDer[1][0] = 0.0; // dy/dx
	alDer[0][1] = 0.0; // dx/dy
	alDer[1][1] = 1.0; // dy/dy

	unsigned int ilab[8]; // 0-5 = telescope, 6 = DUT, 7 = REF

	// plane 0-5:

	double rx[6];
	double ry[6];
	double zprev = _planePosition[0];

	for( int ipl = 0; ipl < 6; ++ipl ){

	  int lhit;
	  if( ipl < 3 )
	    lhit = hts[ipl][iA];
	  else
	    lhit = hts[ipl][jB];

	  double step = _planePosition[ipl] - zprev;
	  zprev = _planePosition[ipl];

	  TMatrixD jacPointToPoint( 5, 5 );
	  jacPointToPoint = Jac5( step );

	  GblPoint * point = new GblPoint( jacPointToPoint );
	  s += step;

	  double dz = hitZ[lhit] - zmA;
	  double xs = xmA[iA] + sxA[iA] * dz; // Ax at plane
	  double ys = ymA[iA] + syA[iA] * dz; // Ay at plane

	  rx[ipl] = hitX[lhit] - xs;
	  ry[ipl] = hitY[lhit] - ys;

	  meas[0] = rx[ipl];
	  meas[1] = ry[ipl];

	  point->addMeasurement( proL2m, meas, measPrec );

	  point->addScatterer( scat, wscatSi );

	  std::vector<int> globalLabels(3);
	  globalLabels[0] = 10 + ipl; // dx
	  globalLabels[1] = 20 + ipl; // dy
	  globalLabels[2] = 40 + ipl; // rot
	  alDer[0][2] = -ys; // dx/rot
	  alDer[1][2] =  xs; // dy/rot

	  point->addGlobals( globalLabels, alDer ); // for MillePede alignment

	  unsigned int iLabel = traj.addPoint(*point);
	  ilab[ipl] = iLabel;
	  sPoint.push_back( s );

	  delete point;

	  if( lDUT && ipl == 2 ) { // insert DUT

	    step = DUTz - zprev;
	    zprev = DUTz;

	    jacPointToPoint = Jac5( step );
	    point = new GblPoint( jacPointToPoint );
	    s += step;

	    double tetDUT = 0.0136 * sqrt(DUTX0) / p * ( 1 + 0.038*std::log(DUTX0) );

	    TVectorD wscatDUT(2);
	    wscatDUT[0] = 1.0 / ( tetDUT * tetDUT ); //weight
	    wscatDUT[1] = 1.0 / ( tetDUT * tetDUT );

	    point->addScatterer( scat, wscatDUT );

	    // DUT measurement:

	    if( lkA[iA] ) {

	      naldut++;

	      meas[0] = DUTsxA[iA];
	      meas[1] = DUTdyA[iA];

	      double resx = 53E-3; // [mm] CMS col resolution
	      double resy =  8E-3; // [mm] CMS row resolution at 20 deg tilt

	      TVectorD measWeight(2);
	      measWeight[0] = 1.0 / resx / resx; // weight = 1/resolution^2
	      measWeight[1] = 1.0 / resy / resy;

	      point->addMeasurement( proL2m, meas, measWeight );

	      // DUT alignment: dx, dy, dz, drot, dtilt

	      int nAlignmentParametersForDUT = 4;
	      if( abs( tilt ) > 1.1 ) nAlignmentParametersForDUT = 5; // tilt [deg]

	      TMatrixD DUTDer( 2, nAlignmentParametersForDUT ); // alignment derivatives
	      DUTDer[0][0] = 1.0; // dx/dx
	      DUTDer[1][0] = 0.0;

	      DUTDer[0][1] = 0.0;
	      DUTDer[1][1] = 1.0; // dy/dy

	      DUTDer[0][2] = sxA[iA]; // dx/dz
	      DUTDer[1][2] = syA[iA]; // dy/dz

	      DUTDer[0][3] = yA; // dx/drot
	      DUTDer[1][3] =-xA; // dy/drot

	      if( nAlignmentParametersForDUT > 4 ) {
		DUTDer[0][4] = 0.0; // dx/dtilt
		DUTDer[1][4] = yA*sin(tilt*wt); // dy/dtilt
	      }

	      // global labels for Pede:

	      std::vector<int> DUTLabels( nAlignmentParametersForDUT );
	      DUTLabels[0] = 16; // dx
	      DUTLabels[1] = 26; // dy
	      DUTLabels[2] = 36; // dz
	      DUTLabels[3] = 46; // drot
	      if( nAlignmentParametersForDUT > 4 )
		DUTLabels[4] = 56; // dtilt

	      point->addGlobals( DUTLabels, DUTDer ); // for MillePede alignment

	    }//lkA linked hit in DUT

	    iLabel = traj.addPoint(*point);

	    ilab[6] = iLabel;

	    delete point;

	  }//DUT present

	} // loop over planes

	// REF:

	// monitor what we put into GBL:

	selxHisto->fill( xA ); // triplet at DUT
	selyHisto->fill( yA );
	selaxHisto->fill( sxA[iA]*1E3 );
	selayHisto->fill( syA[iA]*1E3 );
	seldxHisto->fill( dx*1E3 ); // triplet-driplet match
	seldyHisto->fill( dy*1E3 );
	selkxHisto->fill( kx*1E3 ); // triplet-driplet kink
	selkyHisto->fill( ky*1E3 );

	seldx1Histo->fill( rx[1]*1E3 ); // triplet interpol
	seldy1Histo->fill( ry[1]*1E3 );
	seldx3Histo->fill( rx[3]*1E3 ); // triplet extrapol
	seldy3Histo->fill( ry[3]*1E3 );
	seldx4Histo->fill( rx[4]*1E3 );
	seldy4Histo->fill( ry[4]*1E3 );
	seldx5Histo->fill( rx[5]*1E3 );
	seldy5Histo->fill( ry[5]*1E3 );
	if( lkA[iA] ) {
	  seldx6Histo->fill( DUTsxA[iA]*1E3 );
	  seldy6Histo->fill( DUTdyA[iA]*1E3 );
	}

	// debug:
	/*
	cout << "traj with " << traj.getNumPoints() << " points:" << endl;
	for( int ipl = 0; ipl < 6; ++ipl ){
	  cout << "  plane " << ipl << ", lab " << ilab[ipl];
	  cout << "  z " << sPoint[ipl];
	  cout << " dx " << rx[ipl]*1E3;
	  cout << " dy " << ry[ipl]*1E3;
	  cout << endl;
	}
	*/
	double Chi2;
	int Ndf;
	double lostWeight;

	traj.fit( Chi2, Ndf, lostWeight );

	//cout << " chi2 " << Chi2 << ", ndf " << Ndf << endl;

	gblndfHisto->fill( Ndf );
	if( Ndf == 8 ) 
	  gblchi2aHisto->fill( Chi2 );
	else
	  gblchi2bHisto->fill( Chi2 );
	probchi = TMath::Prob( Chi2, Ndf );
	gblprbHisto->fill( probchi );

	// bad fits:

	if( probchi < 0.01 ) {

	  badxHisto->fill( xA ); // triplet at DUT
	  badyHisto->fill( yA );
	  badaxHisto->fill( sxA[iA]*1E3 );
	  badayHisto->fill( syA[iA]*1E3 );
	  baddxHisto->fill( dx*1E3 ); // triplet-driplet match
	  baddyHisto->fill( dy*1E3 );
	  badkxHisto->fill( kx*1E3 ); // triplet-driplet kink
	  badkyHisto->fill( ky*1E3 );

	  baddx1Histo->fill( rx[1]*1E3 ); // triplet interpol
	  baddy1Histo->fill( ry[1]*1E3 );
	  baddx3Histo->fill( rx[3]*1E3 ); // triplet extrapol
	  baddy3Histo->fill( ry[3]*1E3 );
	  baddx4Histo->fill( rx[4]*1E3 );
	  baddy4Histo->fill( ry[4]*1E3 );
	  baddx5Histo->fill( rx[5]*1E3 );
	  baddy5Histo->fill( ry[5]*1E3 );
	  if( lkA[iA] ) {
	    baddx6Histo->fill( DUTsxA[iA]*1E3 );
	    baddy6Histo->fill( DUTdyA[iA]*1E3 );
	  }

	}// bad fit

	else {

	  goodx1Histo->fill( rx[1]*1E3 ); // triplet interpol
	  goody1Histo->fill( ry[1]*1E3 );
	  if( lkA[iA] ) {
	    goodx6Histo->fill( DUTsxA[iA]*1E3 );
	    goody6Histo->fill( DUTdyA[iA]*1E3 );
	  }
	} // OK fit

	// look at fit:

	TVectorD aCorrection(5);
	TMatrixDSym aCovariance(5);

	double ax[8];
	double ay[8];
	unsigned int k = 0;

	// at plane 0:

	int ipos = ilab[0];
	traj.getResults( ipos, aCorrection, aCovariance );

	unsigned int ndim = 2;
	TVectorD aResiduals(ndim);
	TVectorD aMeasErrors(ndim);
	TVectorD aResErrors(ndim);
	TVectorD aDownWeights(ndim);

	traj.getMeasResults( (unsigned int) ipos, ndim, aResiduals, aMeasErrors, aResErrors, aDownWeights );

	TVectorD aKinks(ndim);
	TVectorD aKinkErrors(ndim);
	TVectorD kResErrors(ndim);
	TVectorD kDownWeights(ndim);
	traj.getScatResults( (unsigned int) ipos, ndim, aKinks, aKinkErrors, kResErrors, kDownWeights );

	//track = q/p, x', y', x, y
	//        0,   1,  2,  3, 4

	gblax0Histo->fill( aCorrection[1]*1E3 ); // angle x [mrad]
	gbldx0Histo->fill( aCorrection[3]*1E3 ); // shift x [um]
	gblrx0Histo->fill( ( rx[0] - aCorrection[3] ) * 1E3 ); // residual x [um]
	gblpx0Histo->fill( aResiduals[0] / aResErrors[0] ); // pull
	gblqx0Histo->fill( aKinks[0]*1E3 ); // kink
	ax[k] = aCorrection[1]; // angle correction at plane, for kinks
	ay[k] = aCorrection[2]; // angle correction at plane, for kinks
	k++;

	ipos = ilab[1];
	traj.getResults( ipos, aCorrection, aCovariance );
	traj.getMeasResults( (unsigned int) ipos, ndim, aResiduals, aMeasErrors, aResErrors, aDownWeights );
	traj.getScatResults( (unsigned int) ipos, ndim, aKinks, aKinkErrors, kResErrors, kDownWeights );
	gblax1Histo->fill( aCorrection[1]*1E3 ); // angle x [mrad]
	gbldx1Histo->fill( aCorrection[3]*1E3 ); // shift x [um]
	gblrx1Histo->fill( ( rx[1] - aCorrection[3] ) * 1E3 ); // residual x [um]
	gblpx1Histo->fill( aResiduals[0] / aResErrors[0] ); // pull
	gblqx1Histo->fill( aKinks[0]*1E3 ); // kink
	ax[k] = aCorrection[1]; // angle correction at plane, for kinks
	ay[k] = aCorrection[2]; // angle correction at plane, for kinks
	k++;

	ipos = ilab[2];
	traj.getResults( ipos, aCorrection, aCovariance );
	traj.getMeasResults( (unsigned int) ipos, ndim, aResiduals, aMeasErrors, aResErrors, aDownWeights );
	traj.getScatResults( (unsigned int) ipos, ndim, aKinks, aKinkErrors, kResErrors, kDownWeights );
	gblax2Histo->fill( aCorrection[1]*1E3 ); // angle x [mrad]
	gbldx2Histo->fill( aCorrection[3]*1E3 ); // shift x [um]
	gblrx2Histo->fill( ( rx[2] - aCorrection[3] ) * 1E3 ); // residual x [um]
	gblpx2Histo->fill( aResiduals[0] / aResErrors[0] ); // pull
	gblqx2Histo->fill( aKinks[0]*1E3 ); // kink
	ax[k] = aCorrection[1]; // angle correction at plane, for kinks
	ay[k] = aCorrection[2]; // angle correction at plane, for kinks
	k++;

	if( lDUT ) {
	  ipos = ilab[6]; // 6 = DUT
	  traj.getResults( ipos, aCorrection, aCovariance );
	  traj.getScatResults( (unsigned int) ipos, ndim, aKinks, aKinkErrors, kResErrors, kDownWeights );
	  gblax6Histo->fill( aCorrection[1]*1E3 ); // angle x [mrad]
	  gbldx6Histo->fill( aCorrection[3]*1E3 ); // shift x [um]
	  gbldy6Histo->fill( aCorrection[4]*1E3 ); // shift x [um]
	  gblqx6Histo->fill( aKinks[0]*1E3 ); // kink
	  if( lkA[iA] ) {
	    traj.getMeasResults( (unsigned int) ipos, ndim, aResiduals, aMeasErrors, aResErrors, aDownWeights );
	    gblrx6Histo->fill( ( DUTsxA[iA] - aCorrection[3] ) * 1E3 ); // residual x [um]
	    gblry6Histo->fill( ( DUTdyA[iA] - aCorrection[4] ) * 1E3 ); // residual y [um]
	    gblpx6Histo->fill( aResiduals[0] / aResErrors[0] ); // pull
	  }
	  ax[k] = aCorrection[1]; // angle correction at plane, for kinks
	  ay[k] = aCorrection[2]; // angle correction at plane, for kinks
	  k++;
	}

	ipos = ilab[3];
	traj.getResults( ipos, aCorrection, aCovariance );
	traj.getMeasResults( (unsigned int) ipos, ndim, aResiduals, aMeasErrors, aResErrors, aDownWeights );
	traj.getScatResults( (unsigned int) ipos, ndim, aKinks, aKinkErrors, kResErrors, kDownWeights );
	gblax3Histo->fill( aCorrection[1]*1E3 ); // angle x [mrad]
	gbldx3Histo->fill( aCorrection[3]*1E3 ); // shift x [um]
	gblrx3Histo->fill( ( rx[3] - aCorrection[3] ) * 1E3 ); // residual x [um]
	gblpx3Histo->fill( aResiduals[0] / aResErrors[0] ); // pull
	gblqx3Histo->fill( aKinks[0]*1E3 ); // kink
	ax[k] = aCorrection[1]; // angle correction at plane, for kinks
	ay[k] = aCorrection[2]; // angle correction at plane, for kinks
	k++;

	ipos = ilab[4];
	traj.getResults( ipos, aCorrection, aCovariance );
	traj.getMeasResults( (unsigned int) ipos, ndim, aResiduals, aMeasErrors, aResErrors, aDownWeights );
	traj.getScatResults( (unsigned int) ipos, ndim, aKinks, aKinkErrors, kResErrors, kDownWeights );
	gblax4Histo->fill( aCorrection[1]*1E3 ); // angle x [mrad]
	gbldx4Histo->fill( aCorrection[3]*1E3 ); // shift x [um]
	gblrx4Histo->fill( ( rx[4] - aCorrection[3] ) * 1E3 ); // residual x [um]
	gblpx4Histo->fill( aResiduals[0] / aResErrors[0] ); // pull
	gblqx5Histo->fill( aKinks[0]*1E3 ); // kink
	ax[k] = aCorrection[1]; // angle correction at plane, for kinks
	ay[k] = aCorrection[2]; // angle correction at plane, for kinks
	k++;

	ipos = ilab[5];
	traj.getResults( ipos, aCorrection, aCovariance );
	traj.getMeasResults( (unsigned int) ipos, ndim, aResiduals, aMeasErrors, aResErrors, aDownWeights );
	traj.getScatResults( (unsigned int) ipos, ndim, aKinks, aKinkErrors, kResErrors, kDownWeights );
	gblax5Histo->fill( aCorrection[1]*1E3 ); // angle x [mrad]
	gbldx5Histo->fill( aCorrection[3]*1E3 ); // shift x [um]
	gblrx5Histo->fill( ( rx[5] - aCorrection[3] ) * 1E3 ); // residual x [um]
	gblpx5Histo->fill( aResiduals[0] / aResErrors[0] ); // pull
	gblqx5Histo->fill( aKinks[0]*1E3 ); // kink
	ax[k] = aCorrection[1]; // angle correction at plane, for kinks
	ay[k] = aCorrection[2]; // angle correction at plane, for kinks
	k++;

	// kinks: 1,2 = tele, 3 = DUT, 4,5 = tele

	gblkx1Histo->fill( (ax[1] - ax[0])*1E3 ); // kink at 1 [mrad]
	gblkx2Histo->fill( (ax[2] - ax[1])*1E3 ); // kink at 2 [mrad]
	gblkx3Histo->fill( (ax[3] - ax[2])*1E3 ); // kink at 3 [mrad]
	gblkx4Histo->fill( (ax[4] - ax[3])*1E3 ); // kink at 4 [mrad]
	gblkx5Histo->fill( (ax[5] - ax[4])*1E3 ); // kink at 5 [mrad]
	gblkx6Histo->fill( (ax[6] - ax[5])*1E3 ); // kink at 6 [mrad]

	traj.milleOut( *mille );

	nmille++;

      } // driplet-triplet match

      //------------------------------------------------------------------------
      // CMS DUT efficiency:

      if( abs(dx) < 0.1 && abs(dy) < 0.1 && lkB[jB] ) {

	sixxylkHisto->fill( xA, yA ); // CMS DUT efficiency profile

	bool nm = 0;
	if( lkA[iA] ) nm = 1;

	// transform into DUT system:

	double xAm =-xA + DUTalignx; // invert x, shift to mid
	double yAm = yA + DUTaligny; // shift to mid
	double xAr = xAm - DUTrot*yAm; // rotate CMS pixel
	double yAr = yAm + DUTrot*xAm;
	double xAt = xAr + 0.075; // shift by half bin?
	double yAt = yAr / cos(tilt*wt); // tilt in y

	effxyHisto->fill( xAt, yAt );
	effvsxy->fill( xAt, yAt, nm ); // CMS DUT efficiency profile

	double fidxmax =  12;
	double fidxmin = -12;
	double fidymax =   5;
	double fidymin =  -6;

	if(  event->getRunNumber() < 188 ) { // Feb 2012
	  fidymin = -3;
	  fidymax =  3;
	  fidxmin = -3;
	  fidxmax =  4;
	}
	if(  event->getRunNumber() > 2842 ) { // Apr 2012 TB 22
	  fidymin = -2.7;
	  fidymax =  3.5;
	  fidxmin = -3.8;
	  fidxmax =  2.5;
	}
	if(  event->getRunNumber() > 2877 ) { // Apr 2012 TB 22
	  fidymin = -3.8;
	  fidymax =  3.8;
	  fidxmin = -4.0;
	  fidxmax =  3.0;
	}

	if( yAt > fidymin && yAt < fidymax ) {
	  effvsx->fill( xAt, nm ); // CMS DUT efficiency profile
	  if( probchi > 0.01 ) effvsxg->fill( xAt, nm );
	}

	if( xAt > fidxmin && xAt < fidxmax ) {
	  effvsy->fill( yAt, nm ); // CMS DUT efficiency profile
	}

	if( xAt > fidxmin && xAt < fidxmax && yAt > fidymin && yAt < fidymax ) {
	  effvst1->fill( (event->getTimeStamp()-time0)/fTLU, nm );
	  effvst3->fill( (event->getTimeStamp()-time0)/fTLU, nm );
	  effvsdt->fill( dtns, nm );
	  effvsmm->fill( xmod, ymod, nm ); // CMS DUT efficiency profile
	}

      }// triplet-driplet match

      //------------------------------------------------------------------------
      // eff(REF) with DUT as timing plane:

      if( abs(dx) < 0.1 && abs(dy) < 0.1 && lkA[iA] ) {

	bool nm = 0;
	if( lkB[jB] ) nm = 1;

	rffvsxy->fill( xR, yR, nm ); // CMS REF efficiency profile

	if( abs( yR ) < 3 ) {
	  rffvsx->fill( xR, nm ); // CMS DUT efficiency profile
	}

      }// triplet-driplet match

      //------------------------------------------------------------------------

      if( xA > -1 && // CMS pixel sensor for run 2331
	  xA <  7 &&
	  yA > -5 &&
	  yA <  1 &&
	  abs( sxA[iA] ) < 0.002 && // no upstream scattering
	  abs( syA[iA] ) < 0.002
	  ){
	sixkxsHisto->fill( kx*1E3 );
	sixkysHisto->fill( ky*1E3 );
      }

      //------------------------------------------------------------------------
      // intersect point in z:

      if( abs(dy) < 0.1 ) { // no cut on dx
	if( abs( kx ) > 0.003 ) { // 
	  double zx = ( xmA[iA] - sxA[iA]*zmA - xmB[jB] + sxB[jB]*zmB ) / kx; //z intersect
	  sixzx3Histo->fill( zx - _planePosition[2] );
	}
	if( abs( kx ) > 0.002 ) { // 
	  double zx = ( xmA[iA] - sxA[iA]*zmA - xmB[jB] + sxB[jB]*zmB ) / kx; //z intersect
	  sixzx2Histo->fill( zx - _planePosition[2] );
	}
      }

      if( abs(dx) < 0.1 ) { // no cut on dy
	if( abs( ky ) > 0.003 ) { // 
	  double zy = ( ymA[iA] - syA[iA]*zmA - ymB[jB] + syB[jB]*zmB ) / ky;
	  sixzy3Histo->fill( zy - _planePosition[2] );
	}
	if( abs( ky ) > 0.002 ) { // 
	  double zy = ( ymA[iA] - syA[iA]*zmA - ymB[jB] + syB[jB]*zmB ) / ky;
	  sixzy2Histo->fill( zy - _planePosition[2] );
	}
      }

      //------------------------------------------------------------------------
      // z intersect:

      if( abs(dx) < 0.2 && abs(dy) < 0.2 ) { // looser cut allows more z range

	if( abs( kx ) > 0.001 ) {
	  double zx = ( xmA[iA] - sxA[iA]*zmA - xmB[jB] + sxB[jB]*zmB ) / kx; //z intersect
	  sixzx1Histo->fill( zx - _planePosition[2] );
	}

	if( abs( ky ) > 0.001 ) {
	  double zy = ( ymA[iA] - syA[iA]*zmA - ymB[jB] + syB[jB]*zmB ) / ky;
	  sixzy1Histo->fill( zy - _planePosition[2] );
	}

	// measure scattering angle x after cuts in y:
	// cut on ky creates bias in kx

	if( abs( ky ) > 0.001 ) {
	  double zy = ( ymA[iA] - syA[iA]*zmA - ymB[jB] + syB[jB]*zmB ) / ky;
	  if( abs( zy - DUTz ) < 30 ) {
	    sixkxzyHisto->fill( kx*1E3 );
	  }
	}

	if( abs( kx ) > 0.001 ) {
	  double zx = ( xmA[iA] - sxA[iA]*zmA - xmB[jB] + sxB[jB]*zmB ) / kx;
	  if( abs( zx - DUTz ) < 30 ) {
	    sixkyzxHisto->fill( ky*1E3 );
	  }
	}

	// kx at DUT:

	if( abs( kx ) > 0.001 ) {
	  double zx = ( xmA[iA] - sxA[iA]*zmA - xmB[jB] + sxB[jB]*zmB ) / kx;
	  if( abs( zx - DUTz ) < 30 ) {
	    sixkxzxHisto->fill( kx*1E3 ); // plot with gap, fittp0g.C("sixkxzx")
	  }
	}

	// ky at DUT:

	if( abs( ky ) > 0.001 ) {
	  double zy = ( ymA[iA] - syA[iA]*zmA - ymB[jB] + syB[jB]*zmB ) / ky;
	  if( abs( zy - DUTz ) < 30 ) {
	    sixkyzyHisto->fill( ky*1E3 ); // plot with gap
	  }
	}

      }//match

      //------------------------------------------------------------------------

    }//driplets B j

  }//triplets A i

  // Clear all working arrays

  delete [] planeHitID;
  delete [] hitPlane;
  delete [] hitFits;
  delete [] hitZ;
  delete [] hitEy;
  delete [] hitY;
  delete [] hitEx;
  delete [] hitX;

  return;

}


//------------------------------------------------------------------------------
void EUTelTestFitter::check( LCEvent * /* evt */  ) {
  // nothing to check here - could be used to fill checkplots in reconstruction processor
}


//------------------------------------------------------------------------------
void EUTelTestFitter::end(){

  // Print the summary:

  streamlog_out( MESSAGE ) << "Total number of processed events:    "
			   << setw(10) << setiosflags(ios::right)
			   << _nEvt << resetiosflags(ios::right) << endl
                           << "Total number of event w/o input hit: "
			   << setw(10) << setiosflags(ios::right) << _noOfEventWOInputHit 
                           << resetiosflags(ios::right) << endl
                           << "Total number of event w/o track:     "
			   << setw(10) << setiosflags(ios::right) << _noOfEventWOTrack
                           << resetiosflags(ios::right) << endl
                           << "Total number of reconstructed tracks "
			   << setw(10) << setiosflags(ios::right)
			   << _noOfTracks << resetiosflags(ios::right) << endl
			   << "CMS clusters " << nclst << endl
			   << "CMS matches  " << nmtch << endl
                           << endl;

  cout << cmsdyaHisto->title() << endl;
  cout << "Entries " << cmsdyaHisto->entries() << endl;
  cout << "Maximum " << cmsdyaHisto->maxBinHeight() << endl;

  // Clean memory:

  delete [] _planeSort;
  delete [] _planeID;
  delete [] _planeShiftX;
  delete [] _planeShiftY;
  delete [] _planeRotZ;
  delete [] _planePosition;
  delete [] _planeThickness;
  delete [] _planeX0;
  delete [] _planeResolution;
  delete [] _planeDist;
  delete [] _planeScat;
  delete [] _planeScatAngle;
  delete [] _isActive;

  delete [] _planeHits;
  delete [] _planeChoice;
  delete [] _planeMod;

  delete [] _planeX;
  delete [] _planeEx;
  delete [] _planeY;
  delete [] _planeEy;

  delete [] _fitX;
  delete [] _fitEx;
  delete [] _fitY;
  delete [] _fitEy;
  delete [] _fitArray;

  delete [] _nominalFitArrayX;
  delete [] _nominalErrorX;

  delete [] _nominalFitArrayY;
  delete [] _nominalErrorY;

  // Pede:

  // close the output file:
  delete mille;

  streamlog_out( MESSAGE4 ) << endl << "Generating the steering file for the pede program..." << endl;

  streamlog_out( MESSAGE2 ) << "have " << nmille << " align tracks" << endl;
  streamlog_out( MESSAGE2 ) << "have " << naldut << " DUT links" << endl;

  bool ldut = naldut > 99;

  ofstream steerFile;
  steerFile.open( "steerPede.txt" );

  if( steerFile.is_open()) {

    steerFile << "! generated by EUTelTestFitter" << endl;
    steerFile << "Cfiles" << endl;
    steerFile << "mille.bin" << endl;
    steerFile << endl;

    steerFile << "Parameter" << endl;

    // loop over all planes:

    for( int ipl = 0; ipl < _nTelPlanes; ipl++) {

      if( ipl == 2 ) { // fixed
	steerFile << 10 + ipl << "  0.0 -1.0" << endl;
	steerFile << 20 + ipl << "  0.0 -1.0" << endl;
	steerFile << 40 + ipl << "  0.0 -1.0" << endl;
      }
      else{
	steerFile << 10 + ipl << "  0.0  0.0" << endl; // dx
	steerFile << 20 + ipl << "  0.0  0.0" << endl; // dy
	steerFile << 40 + ipl << "  0.0  0.0" << endl; // rot
      }// not fixed

    }// loop over planes

    // DUT is 6:

    if( ldut ) {
      steerFile << 16 << "  0.0  0.0" << endl; // dx
      steerFile << 26 << "  0.0  0.0" << endl; // dy
      steerFile << 36 << "  0.0  0.0" << endl; // dz
      steerFile << 46 << "  0.0  0.0" << endl; // rot
      steerFile << 56 << "  0.0  0.0" << endl; // dtilt
    }

    steerFile << "Constraint 0 ! sum dx = 0" << endl;
    steerFile << "10  1.0" << endl;
    steerFile << "11  1.0" << endl;
    steerFile << "12  1.0" << endl;
    steerFile << "13  1.0" << endl;
    steerFile << "14  1.0" << endl;
    steerFile << "15  1.0" << endl;
    if( ldut )
      steerFile << "16  1.0" << endl;

    steerFile << "Constraint 0 ! sum dy = 0" << endl;
    steerFile << "20  1.0" << endl;
    steerFile << "21  1.0" << endl;
    steerFile << "22  1.0" << endl;
    steerFile << "23  1.0" << endl;
    steerFile << "24  1.0" << endl;
    steerFile << "25  1.0" << endl;
    if( ldut )
      steerFile << "26  1.0" << endl;

    steerFile << "Constraint 0 ! sum dphi = 0" << endl;
    steerFile << "40  1.0" << endl;
    steerFile << "41  1.0" << endl;
    steerFile << "42  1.0" << endl;
    steerFile << "43  1.0" << endl;
    steerFile << "44  1.0" << endl;
    steerFile << "45  1.0" << endl;
    if( ldut )
      steerFile << "46  1.0" << endl;

    steerFile << endl;
    steerFile << "! chiscut 5.0 2.5" << endl;
    steerFile << "! outlierdownweighting 4" << endl;
    steerFile << endl;
    steerFile << "method inversion 10  0.1" << endl;
    steerFile << endl;
    steerFile << "histprint" << endl;
    steerFile << endl;
    steerFile << "end" << endl;

    steerFile.close();

    streamlog_out( MESSAGE2 ) << "Pede steer file written." << endl;

    // before starting pede, let's check if it is in the path
    bool isPedeInPath = true;

    // create a new process
    redi::ipstream which("which pede");

    // wait for the process to finish
    which.close();

    // get the status
    // if it 255 then the program wasn't found in the path
    isPedeInPath = !( which.rdbuf()->status() == 255 );

    if( !isPedeInPath ) {
      streamlog_out( ERROR ) << "Pede cannot be executed because not found in the path" << endl;
    }
    else {

      std::string command = "pede  steerPede.txt";

      streamlog_out( MESSAGE2 ) << endl;
      streamlog_out( MESSAGE2 ) << "Starting pede..." << endl;
      streamlog_out( MESSAGE2 ) << command.c_str() << endl;

      redi::ipstream pede( command.c_str() );
      string output;
      while ( getline( pede, output ) ) {
	streamlog_out( MESSAGE2 ) << output << endl;
      }

      // wait for the pede execution to finish
      pede.close();

      // check the exit value of pede
      if( pede.rdbuf()->status() == 0 )
	streamlog_out( MESSAGE2 ) << "Pede successfully finished" << endl;

      // reading back the millepede.res file:

      string millepedeResFileName = "millepede.res";

      streamlog_out( MESSAGE2 ) << "Reading back the " << millepedeResFileName << endl;

      // open the millepede ASCII output file
      ifstream millepede( millepedeResFileName.c_str() );

      if( millepede.bad() || !millepede.is_open() ) {
	streamlog_out( ERROR4 )
	  << "Error opening the " << millepedeResFileName << endl
	  << "The alignment file cannot be saved" << endl;
      }
      else {
	vector<double > tokens;
	stringstream tokenizer;
	string line;
	double buffer;

	// get the first line and throw it away since it is a comment!

	getline( millepede, line );
	std::cout << "line: " <<  line  << std::endl;

	unsigned int numplane = 6; // align pars per plane in Pede
	if( ldut ) numplane = 7;
	unsigned int numpars = 3*numplane; // dx dy rot
	if( ldut ) numpars += 2; // DUT has dz and tilt

	map< unsigned int, double > alpar; // map = associative array

	for( unsigned int ipar = 0; ipar < numpars; ++ipar ) {

	  if( millepede.eof() ) break;

	  getline( millepede, line );

	  if( line.empty() ) continue;

	  tokens.clear();
	  tokenizer.clear();
	  tokenizer.str( line );

	  while( tokenizer >> buffer ) {
	    tokens.push_back( buffer );
	  }

	  int lpar = (int) ( tokens[0] + 0.5 ); // par label

	  bool isFixed = ( tokens.size() == 3 );

	  if( isFixed )
	    cout << "Parameter " << lpar
		 << " is at " << tokens[1]
		 << " (fixed)"  << endl;
	  else
	    cout << "Parameter " << lpar
		 << " is at " << tokens[1]
		 << " +/- " << tokens[4] << endl;

	  alpar[lpar] = tokens[1];

	}//loop param

	if( ldut ){
	  cout << "DUT alignment corrections:" << endl;
	  cout << "dx " << alpar[16] << " mm" << endl;
	  cout << "dy " << alpar[26] << " mm" << endl;
	  cout << "dz " << alpar[36] << " mm" << endl;
	  cout << "df " << alpar[46] << " rad" << endl;
	  cout << "dt " << alpar[56]*180/3.141592654 << " deg" << endl;
	}
      }//millepede OK

      millepede.close();

    }//PedeInPath

  }// pede steer file open

  else {
    streamlog_out( ERROR2 ) << "Could not open steering file." << endl;
  }

} // end end


//------------------------------------------------------------------------------
void EUTelTestFitter::bookHistos()
{

#if defined(USE_AIDA) || defined(MARLIN_USE_AIDA)

  streamlog_out( MESSAGE2 ) <<  "Booking histograms " << endl;
  streamlog_out( MESSAGE2 ) << "Histogram information searched in " << _histoInfoFileName << endl;

  auto_ptr<EUTelHistogramManager> histoMgr( new EUTelHistogramManager( _histoInfoFileName ));
  EUTelHistogramInfo    * histoInfo;
  bool                    isHistoManagerAvailable;

  try {
    isHistoManagerAvailable = histoMgr->init();
  }
  catch( ios::failure& e) {
    streamlog_out( ERROR ) << "I/O problem with " << _histoInfoFileName << "\n"
                           << "Continuing without histogram manager"    << endl;
    isHistoManagerAvailable = false;
  }
  catch( ParseException& e ) {
    streamlog_out( ERROR ) << e.what() << "\n"
                           << "Continuing without histogram manager" << endl;
    isHistoManagerAvailable = false;
  }

  // event time

  t100Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "t100", 100, 0, 100 );
  t100Histo->setTitle( "telescope event time;telescope event time [s];events/s" );

  t300Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "t300", 300, 0, 300 );
  t300Histo->setTitle( "telescope event time;telescope event time [s];events/s" );

  t1000Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "t1000", 100, 0, 1000 );
  t1000Histo->setTitle( "telescope event time;telescope event time [s];events/10s" );

  t3600Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "t3600", 360, 0, 3600 );
  t3600Histo->setTitle( "telescope event time;telescope event time [s];events/10s" );

  dtHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "dt", 500, 0, 5000 );
  dtHisto->setTitle( "telescope time between events;telescope time between events [#mus];events" );

  dtmsHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "dtms", 500, 0, 100 );
  dtmsHisto->setTitle( "telescope time between events;telescope time between events [ms];events" );

  logdtHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "logdt", 500, -1, 4 );
  logdtHisto->setTitle( "telescope time between events;telescope time between events log_{10}(#Deltat [ms]);events" );

  logdtcmsHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "logdtcms", 500, -1, 4 );
  logdtcmsHisto->setTitle( "telescope time between events with DUT;telescope time between events log_{10}(#Deltat [ms]);events" );

  // DUT clusters:

  cmscolHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "cmscol", 52, -0.5, 51.5 );
  cmscolHisto->setTitle( "DUT column;DUT cluster col;DUT clusters" );

  cmsrowHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "cmsrow", 80, -0.5, 79.5 );
  cmsrowHisto->setTitle( "DUT row;DUT cluster row;DUT clusters" );

  cmsnpxHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "cmsnpx", 11, -0.5, 10.5 );
  cmsnpxHisto->setTitle( "DUT cluster size;DUT pixel per cluster;DUT clusters" );

  cmsadcHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "cmsadc", 100, 0, 100 );
  cmsadcHisto->setTitle( "DUT cluster charge;DUT cluster charge [ke];DUT clusters" );


  refcolHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "refcol", 52, -0.5, 51.5 );
  refcolHisto->setTitle( "REF column;REF cluster col;REF clusters" );

  refrowHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "refrow", 80, -0.5, 79.5 );
  refrowHisto->setTitle( "REF row;REF cluster row;REF clusters" );

  refnpxHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "refnpx", 11, -0.5, 10.5 );
  refnpxHisto->setTitle( "REF cluster size;REF pixel per cluster;REF clusters" );

  refadcHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "refadc", 100, 0, 100 );
  refadcHisto->setTitle( "REF cluster charge;REF cluster charge [ke];REF clusters" );


  cmsncolHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "cmsncol", 11, -0.5, 10.5 );
  cmsncolHisto->setTitle( "DUT cluster col size;columns per cluster;DUT clusters" );

  cmsnrowHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "cmsnrow", 11, -0.5, 10.5 );
  cmsnrowHisto->setTitle( "DUT cluster row size;rows per cluster;DUT clusters" );

  cmsetaHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "cmseta", 100, -1, 1 );
  cmsetaHisto->setTitle( "DUT row eta;DUT cluster row eta;DUT clusters" );

  // telescope hits per plane:

  hits0Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "hits0", 101, -0.5, 100.5 );
  hits0Histo->setTitle( "hits in plane 0;hits in plane 0;events" );

  hits1Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "hits1", 101, -0.5, 100.5 );
  hits1Histo->setTitle( "hits in plane 1;hits in plane 1;events" );

  hits2Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "hits2", 101, -0.5, 100.5 );
  hits2Histo->setTitle( "hits in plane 2;hits in plane 2;events" );

  hits3Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "hits3", 101, -0.5, 100.5 );
  hits3Histo->setTitle( "hits in plane 3;hits in plane 3;events" );

  hits4Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "hits4", 101, -0.5, 100.5 );
  hits4Histo->setTitle( "hits in plane 4;hits in plane 4;events" );

  hits5Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "hits5", 101, -0.5, 100.5 );
  hits5Histo->setTitle( "hits in plane 5;hits in plane 5;events" );

  // telescope dx:

  dx01Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "dx01", 100, -1, 1 );
  dx01Histo->setTitle( "x1-x0;x_{1}-x_{0} [mm];hit pairs" );

  dy01Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "dy01", 100, -1, 1 );
  dy01Histo->setTitle( "y1-y0;y_{1}-y_{0} [mm];hit pairs" );

  du01Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "du01", 100, -1, 1 );
  du01Histo->setTitle( "x1-x0, |dy| < 1;x_{1}-x_{0} [mm];hit pairs" );

  dx02Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "dx02", 100, -1, 1 );
  dx02Histo->setTitle( "x2-x0;x_{2}-x_{0} [mm];hit pairs" );

  dx03Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "dx03", 100, -1, 1 );
  dx03Histo->setTitle( "x3-x0;x_{3}-x_{0} [mm];hit pairs" );

  dx04Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "dx04", 100, -1, 1 );
  dx04Histo->setTitle( "x4-x0;x_{4}-x_{0} [mm];hit pairs" );

  dx05Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "dx05", 100, -1, 1 );
  dx05Histo->setTitle( "x5-x0;x_{5}-x_{0} [mm];hit pairs" );

  dx12Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "dx12", 100, -1, 1 );
  dx12Histo->setTitle( "x2-x1;x_{2}-x_{1} [mm];hit pairs" );

  dy12Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "dy12", 100, -1, 1 );
  dy12Histo->setTitle( "y2-y1;y_{2}-y_{1} [mm];hit pairs" );

  du12Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "du12", 100, -1, 1 );
  du12Histo->setTitle( "x2-x1, |dy| < 1;x_{2}-x_{1} [mm];hit pairs" );

  dx23Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "dx23", 100, -1, 1 );
  dx23Histo->setTitle( "x3-x2;x_{3}-x_{2} [mm];hit pairs" );

  dy23Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "dy23", 100, -1, 1 );
  dy23Histo->setTitle( "y3-y2;y_{3}-y_{2} [mm];hit pairs" );

  du23Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "du23", 100, -1, 1 );
  du23Histo->setTitle( "x3-x2, |dy| < 1;x_{3}-x_{2} [mm];hit pairs" );

  dx34Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "dx34", 100, -1, 1 );
  dx34Histo->setTitle( "x4-x3;x_{4}-x_{3} [mm];hit pairs" );

  dy34Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "dy34", 100, -1, 1 );
  dy34Histo->setTitle( "y4-y3;y_{4}-y_{3} [mm];hit pairs" );

  du34Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "du34", 100, -1, 1 );
  du34Histo->setTitle( "x4-x3, |dy| < 1;x_{4}-x_{3} [mm];hit pairs" );

  dx45Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "dx45", 100, -1, 1 );
  dx45Histo->setTitle( "x5-x4;x_{5}-x_{4} [mm];hit pairs" );

  dy45Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "dy45", 100, -1, 1 );
  dy45Histo->setTitle( "y5-y4;y_{5}-y_{4} [mm];hit pairs" );

  du45Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "du45", 100, -1, 1 );
  du45Histo->setTitle( "x5-x4, |dy| < 1;x_{5}-x_{4} [mm];hit pairs" );

  // triplets:

  double dz = _planePosition[2] - _planePosition[0];
  double dx = 0.005 * dz;

  da02Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "da02", 100, -dx, dx );
  da02Histo->setTitle( "x2-x0, |dy| < 2 mm;x_{2}-x_{0} [mm];hit pairs" );

  db02Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "db02", 100, -dx, dx );
  db02Histo->setTitle( "y2-y0, |dx| < 2 mm;y_{2}-y_{0} [mm];hit pairs" );

  tridxHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "tridx", 100, -100, 100 );
  tridxHisto->setTitle( "triplet dx;x_{1}-x_{m} [#mum];telescope triplets" );

  tridyHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "tridy", 100, -100, 100 );
  tridyHisto->setTitle( "triplet dy;y_{1}-y_{m} [#mum];telescope triplets" );

  tridx1Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "tridx1", 100, -1, 1 );
  tridx1Histo->setTitle( "triplet dx;x_{1}-x_{t} [mm];telescope triplets" );

  tridy1Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "tridy1", 100, -1, 1 );
  tridy1Histo->setTitle( "triplet dy;y_{1}-y_{t} [mm];telescope triplets" );

  tridxvsx = AIDAProcessor::histogramFactory(this)->
    createProfile1D( "tridxvsx", 100, -10, 10, -100, 100 );
  tridxvsx->setTitle( "triplet x resid vs x;x [mm];triplet <#Deltax> [#mum]" );

  tridxvsy = AIDAProcessor::histogramFactory(this)->
    createProfile1D( "tridxvsy", 50, -5, 5, -100, 100 );
  tridxvsy->setTitle( "triplet x resid vs y;y [mm];triplet <#Deltax> [#mum]" );

  tridxvstx = AIDAProcessor::histogramFactory(this)->
    createProfile1D( "tridxvstx", 80, -2, 2, -100, 100 );
  tridxvstx->setTitle( "triplet x resid vs tx;t_{x} [mrad];triplet <#Deltax> [#mum]" );

  tridxvsty = AIDAProcessor::histogramFactory(this)->
    createProfile1D( "tridxvsty", 80, -2, 2, -100, 100 );
  tridxvsty->setTitle( "triplet x resid vs ty;t_{y} [mrad];triplet <#Deltax> [#mum]" );

  tridyvsx = AIDAProcessor::histogramFactory(this)->
    createProfile1D( "tridyvsx", 100, -10, 10, -100, 100 );
  tridyvsx->setTitle( "triplet y resid vs x;x [mm];triplet <#Deltay> [#mum]" );

  tridyvsy = AIDAProcessor::histogramFactory(this)->
    createProfile1D( "tridyvsy", 50, -5, 5, -100, 100 );
  tridyvsy->setTitle( "triplet y resid vs y;y [mm];triplet <#Deltay> [#mum]" );

  tridyvstx = AIDAProcessor::histogramFactory(this)->
    createProfile1D( "tridyvstx", 80, -2, 2, -100, 100 );
  tridyvstx->setTitle( "triplet y resid vs tx;t_{x} [mrad];triplet <#Deltay> [#mum]" );

  tridyvsty = AIDAProcessor::histogramFactory(this)->
    createProfile1D( "tridyvsty", 80, -2, 2, -100, 100 );
  tridyvsty->setTitle( "triplet y resid vs ty;t_{y} [mrad];triplet <#Deltay> [#mum]" );

  tridx3Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "tridx3", 100, -1, 1 );
  tridx3Histo->setTitle( "triplet dx;x_{3}-x_{t} [mm];telescope triplets" );

  tridy3Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "tridy3", 100, -1, 1 );
  tridy3Histo->setTitle( "triplet dy;y_{3}-y_{t} [mm];telescope triplets" );

  tridx3bHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "tridx3b", 100, -250, 250 );
  tridx3bHisto->setTitle( "triplet dx;x_{3}-x_{t} [um];telescope triplets" );

  tridy3bHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "tridy3b", 100, -250, 250 );
  tridy3bHisto->setTitle( "triplet dy;y_{3}-y_{t} [um];telescope triplets" );

  tridx4Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "tridx4", 100, -2, 2 );
  tridx4Histo->setTitle( "triplet dx;x_{4}-x_{t} [mm];telescope triplets" );

  tridy4Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "tridy4", 100, -2, 2 );
  tridy4Histo->setTitle( "triplet dy;y_{4}-y_{t} [mm];telescope triplets" );

  tridx4bHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "tridx4b", 100, -400, 400 );
  tridx4bHisto->setTitle( "triplet dx;x_{4}-x_{t} [um];telescope triplets" );

  tridy4bHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "tridy4b", 100, -400, 400 );
  tridy4bHisto->setTitle( "triplet dy;y_{4}-y_{t} [um];telescope triplets" );

  tridx5Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "tridx5", 100, -3, 3 );
  tridx5Histo->setTitle( "triplet dx;x_{5}-x_{t} [mm];telescope triplets" );

  tridy5Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "tridy5", 100, -3, 3 );
  tridy5Histo->setTitle( "triplet dy;y_{5}-y_{t} [mm];telescope triplets" );

  tridx5bHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "tridx5b", 100, -1000, 1000 );
  tridx5bHisto->setTitle( "triplet dx;x_{5}-x_{t} [um];telescope triplets" );

  tridy5bHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "tridy5b", 100, -1000, 1000 );
  tridy5bHisto->setTitle( "triplet dy;y_{5}-y_{t} [um];telescope triplets" );

  trixHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "trix", 240, -12, 12 );
  trixHisto->setTitle( "triplet x1;x1 [mm];telescope triplets" );

  triyHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "triy", 80, -8, 8 );
  triyHisto->setTitle( "triplet y1;y1 [mm];telescope triplets" );

  trixyHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram2D( "trixy", 240, -12, 12, 160, -8, 8 );
  trixyHisto->setTitle( "triplet y1 vs x1;x1 [mm];y1 [mm];telescope triplets" );

  tritxHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "tritx", 100, -10, 10 );
  tritxHisto->setTitle( "triplet slope x;#theta_{x} [mrad];telescope triplets" );

  trityHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "trity", 100, -10, 10 );
  trityHisto->setTitle( "triplet slope y;#theta_{y} [mrad];telescope triplets" );

  dutxHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "dutx", 120, -12, 12 );
  dutxHisto->setTitle( "triplet x at DUT;x_{DUT} [mm];telescope triplets" );

  dutyHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "duty", 80, -8, 8 );
  dutyHisto->setTitle( "triplet y at DUT;y_{DUT} [mm];telescope triplets" );

  dutxyHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram2D( "dutxy", 120, -12, 12, 80, -8, 8 );
  dutxyHisto->setTitle( "triplet x-y at DUT;x_{DUT} [mm];y_{DUT} [mm];telescope triplets" );

  // DUT pixel vs triplets:

  cmssxaHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "cmssxa", 400, -10, 10 );
  cmssxaHisto->setTitle( "Pixel + Telescope x;cluster + triplet #Sigmax [mm];clusters" );

  cmsdyaHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "cmsdya", 200, -5, 5 );
  cmsdyaHisto->setTitle( "Pixel - telescope y;cluster - triplet #Deltay [mm];clusters" );

  cmssxHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "cmssx", 200, -500, 500 );
  cmssxHisto->setTitle( "Pixel + Telescope x;cluster + triplet #Sigmax [#mum];clusters" );

  cmsdyHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "cmsdy", 500, -500, 500 );
  cmsdyHisto->setTitle( "Pixel - telescope y;cluster - triplet #Deltay [#mum];clusters" );

  cmssxfHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "cmssxf", 200, -500, 500 );
  cmssxfHisto->setTitle( "fiducial Pixel - telescope x;fiducial cluster - triplet #Deltax [#mum];fiducial clusters" );

  cmsdyfHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "cmsdyf", 500, -500, 500 );
  cmsdyfHisto->setTitle( "fiducial Pixel - telescope y;fiducial cluster - triplet #Deltay [#mum];fiducial clusters" );

  cmssxfcHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "cmssxfc", 200, -500, 500 );
  cmssxfcHisto->setTitle( "fiducial Pixel - telescope x;fiducial cluster - triplet #Deltax [#mum];fiducial clusters" );

  cmsdyfcHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "cmsdyfc", 500, -500, 500 );
  cmsdyfcHisto->setTitle( "fiducial Pixel - telescope y;fiducial cluster - triplet #Deltay [#mum];fiducial clusters" );

  cmsdy1Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "cmsdy1", 500, -500, 500 );
  cmsdy1Histo->setTitle( "Pixel - telescope y;1-row cluster - triplet #Deltay [#mum];clusters" );

  cmsdy2Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "cmsdy2", 500, -500, 500 );
  cmsdy2Histo->setTitle( "Pixel - telescope y;2-row cluster - triplet #Deltay [#mum];clusters" );

  cmsdy3Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "cmsdy3", 500, -500, 500 );
  cmsdy3Histo->setTitle( "Pixel - telescope y;3-row cluster - triplet #Deltay [#mum];clusters" );

  cmssxfcnHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "cmssxfcn", 200, -500, 500 );
  cmssxfcnHisto->setTitle( "fiducial Pixel - telescope x;fiducial cluster - triplet #Deltax [#mum];fiducial clusters" );

  cmsdyfcnHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "cmsdyfcn", 500, -500, 500 );
  cmsdyfcnHisto->setTitle( "fiducial Pixel - telescope y;fiducial cluster - triplet #Deltay [#mum];fiducial clusters" );

  cmsdyfcntHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "cmsdyfcnt", 500, -500, 500 );
  cmsdyfcntHisto->setTitle( "fiducial Pixel - telescope y;fiducial cluster - triplet #Deltay [#mum];fiducial clusters" );

  cmsdyfcntqHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "cmsdyfcntq", 500, -500, 500 );
  cmsdyfcntqHisto->setTitle( "fiducial Pixel - telescope y;fiducial cluster - triplet #Deltay [#mum];fiducial clusters" );

  cmsdyfcntq1Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "cmsdyfcntq1", 500, -500, 500 );
  cmsdyfcntq1Histo->setTitle( "fiducial Pixel - telescope y;fiducial cluster - triplet #Deltay [#mum];fiducial clusters" );

  cmsdyq0Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "cmsdyq0", 500, -500, 500 );
  cmsdyq0Histo->setTitle( "Pixel - telescope y, Q < 18;cluster - triplet #Deltay [#mum];low dE/dx clusters" );

  cmsdyq1Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "cmsdyq1", 500, -500, 500 );
  cmsdyq1Histo->setTitle( "Pixel - telescope y, 18 < Q < 40;cluster - triplet #Deltay [#mum];low dE/dx clusters" );

  cmsdyq2Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "cmsdyq2", 500, -500, 500 );
  cmsdyq2Histo->setTitle( "Pixel - telescope y, Q > 40;cluster - triplet #Deltay [#mum];high dE/dx clusters" );

  cmsdyvsx = AIDAProcessor::histogramFactory(this)->
    createProfile1D( "cmsdyvsx", 52, -26*0.15, 26*0.15, -500, 500 );
  cmsdyvsx->setTitle( "DUT y resid vs x;x [mm];<DUT cluster - telescope triplet #Deltay> [#mum]" );

  cmsdyvsy = AIDAProcessor::histogramFactory(this)->
    createProfile1D( "cmsdyvsy", 80, -4, 4, -500, 500 );
  cmsdyvsy->setTitle( "DUT y resid vs y;y [mm];<DUT cluster - telescope triplet #Deltay> [#mum]" );

  cmsdyvsty = AIDAProcessor::histogramFactory(this)->
    createProfile1D( "cmsdyvsty", 80, -2, 2, -500, 500 );
  cmsdyvsty->setTitle( "DUT y resid vs tet y;#theta_{y} [mrad];<DUT cluster - telescope triplet #Deltay> [#mum]" );

  cmspixvsxy = AIDAProcessor::histogramFactory(this)->
    createHistogram2D( "cmspixvsxy", 60, 0, 300, 40, 0, 200 );
  cmspixvsxy->setTitle( "DUT pixel occupancy;x mod 300 #mum;y mod 200 #mum;clusters" );

  cmsqvsxy = AIDAProcessor::histogramFactory(this)->
    createProfile2D( "cmsqvsxy", 60, 0, 300, 40, 0, 200, 0, 250 );
  cmsqvsxy->setTitle( "DUT cluster charge map;x mod 300 #mum;y mod 200 #mum;<cluster charge> [ke]" );

  cmsqvsdt = AIDAProcessor::histogramFactory(this)->
    createProfile1D( "cmsqvsdt", 40, -40, 40, 0, 250 );
  cmsqvsdt->setTitle( "DUT cluster charge vs phase;TLU-DUT #Deltat [ns];<cluster charge> [ke]" );

  cmsrmsxvsxm = AIDAProcessor::histogramFactory(this)->
    createProfile1D( "cmsrmsxvsxm", 60, 0, 300, 0, 250 );
  cmsrmsxvsxm->setTitle( "DUT x resolution;telescope x mod 300 [#mum];RMS(#Deltax) [#mum]" );

  cmsrmsyvsym = AIDAProcessor::histogramFactory(this)->
    createProfile1D( "cmsrmsyvsym", 40, 0, 200, 0, 250 );
  cmsrmsyvsym->setTitle( "DUT y resolution;telescope y mod 200 [#mum];RMS(#Deltay) [#mum]" );

  cmsetavsym = AIDAProcessor::histogramFactory(this)->
    createProfile1D( "cmsetavsym", 40, 0, 200, -1.5, 1.5 );
  cmsetavsym->setTitle( "DUT eta vs y;telescope y mod 200 [#mum];<eta>" );

  cmsetavsym2 = AIDAProcessor::histogramFactory(this)->
    createProfile1D( "cmsetavsym2", 40, 0, 200, -1.5, 1.5 );
  cmsetavsym2->setTitle( "DUT eta 2-row vs y;telescope y mod 200 [#mum];<eta> 2-row" );

  cmsncolvsxm = AIDAProcessor::histogramFactory(this)->
    createProfile1D( "cmsncolvsxm", 60, 0, 300, 0, 5 );
  cmsncolvsxm->setTitle( "DUT cols vs x;telescope x mod 300 [#mum];<cols/cluster>" );

  cmsnrowvsym = AIDAProcessor::histogramFactory(this)->
    createProfile1D( "cmsnrowvsym", 40, 0, 200, 0, 5 );
  cmsnrowvsym->setTitle( "DUT rows vs y;telescope y mod 200 [#mum];<rows/cluster>" );

  cmsym1Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "cmsym1", 40, 0, 200 );
  cmsym1Histo->setTitle( "y at DUT 1-row clus;telescope y mod 200 [#mum];1-row clusters" );

  cmsdyvsxHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram2D( "cmsdyvsxh", 52, -26*0.15, 26*0.15, 50, -250, 250 );
  cmsdyvsxHisto->setTitle( "DUT resolution;x [mm];#Deltay [#mum];clusters" );


  cmsxxHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram2D( "cmsxx", 52, -0.5, 51.5, 110, -11, 11 );
  cmsxxHisto->setTitle( "x correlation;DUT cluster col;telescope triplet x [mm];clusters" );

  cmsyyHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram2D( "cmsyy", 80, -0.5, 79.5, 55, -5.5, 5.5 );
  cmsyyHisto->setTitle( "y correlation;DUT cluster row;telescope triplet y [mm];clusters" );

  cmsqHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "cmsq", 100, 0, 100 );
  cmsqHisto->setTitle( "DUT cluster charge linked;DUT cluster charge [ke];DUT linked clusters" );

  cmsqfHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "cmsqf", 100, 0, 100 );
  cmsqfHisto->setTitle( "DUT cluster charge linked fiducial;DUT cluster charge [ke];DUT linked fiducial clusters" );

  cmsqdt0Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "cmsqdt0", 100, 0, 100 );
  cmsqdt0Histo->setTitle( "DUT cluster charge #Deltat < -20 ns;DUT cluster charge [ke];DUT linked clusters" );

  cmsqdt1Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "cmsqdt1", 100, 0, 100 );
  cmsqdt1Histo->setTitle( "DUT cluster charge #Deltat < -10 ns;DUT cluster charge [ke];DUT linked clusters" );

  cmsqdt2Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "cmsqdt2", 100, 0, 100 );
  cmsqdt2Histo->setTitle( "DUT cluster charge #Deltat <   0 ns;DUT cluster charge [ke];DUT linked clusters" );

  cmsqdt3Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "cmsqdt3", 100, 0, 100 );
  cmsqdt3Histo->setTitle( "DUT cluster charge #Deltat <  10 ns;DUT cluster charge [ke];DUT linked clusters" );

  cmsqdt4Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "cmsqdt4", 100, 0, 100 );
  cmsqdt4Histo->setTitle( "DUT cluster charge #Deltat <  20 ns;DUT cluster charge [ke];DUT linked clusters" );

  cmsqdt5Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "cmsqdt5", 100, 0, 100 );
  cmsqdt5Histo->setTitle( "DUT cluster charge #Deltat >  20 ns;DUT cluster charge [ke];DUT linked clusters" );

  cmsxHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "cmsx", 120, -12, 12 );
  cmsxHisto->setTitle( "DUT - triplet link;triplet x [mm];clusters" );

  cmsyHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "cmsy", 80, -8, 8 );
  cmsyHisto->setTitle( "DUT - triplet link;triplet y [mm];clusters" );

  cmsxyHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram2D( "cmsxy", 120, -12, 12, 80, -8, 8 );
  cmsxyHisto->setTitle( "DUT - triplet link;triplet x [mm];triplet y [mm];clusters" );


  effxyHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram2D( "effxy", 60, -4.5, 4.5, 90, -4.5, 4.5 );
  effxyHisto->setTitle( "DUT - triplet link;triplet xp [mm];triplet yp [mm];clusters" );

  effvsxy = AIDAProcessor::histogramFactory(this)->
    createProfile2D( "effvsxy", 60, -4.5, 4.5, 90, -4.5, 4.5, -1, 2 );
  effvsxy->setTitle( "DUT efficiency;telescope xp [mm];telescope yp [mm];CMS DUT efficiency" );

  effvsx = AIDAProcessor::histogramFactory(this)->
    createProfile1D( "effvsx", 60, -4.5, 4.5, -1, 2 );
  effvsx->setTitle( "DUT efficiency;telescope x [mm];CMS DUT efficiency" );

  effvsxg = AIDAProcessor::histogramFactory(this)->
    createProfile1D( "effvsxg", 60, -4.5, 4.5, -1, 2 );
  effvsxg->setTitle( "DUT efficiency;telescope x [mm];CMS DUT efficiency" );

  effvsy = AIDAProcessor::histogramFactory(this)->
    createProfile1D( "effvsy", 90, -4.5, 4.5, -1, 2 );
  effvsy->setTitle( "DUT efficiency;telescope y [mm];<CMS DUT efficiency>" );

  effvst1 = AIDAProcessor::histogramFactory(this)->
    createProfile1D( "effvst1", 100, 0, 1000, -1, 2 );
  effvst1->setTitle( "DUT efficiency;time [s];<CMS DUT efficiency>" );

  effvst3 = AIDAProcessor::histogramFactory(this)->
    createProfile1D( "effvst3", 360, 0, 3600, -1, 2 );
  effvst3->setTitle( "DUT efficiency;time [s];<CMS DUT efficiency>" );

  effvsdt = AIDAProcessor::histogramFactory(this)->
    createProfile1D( "effvsdt", 40, -40, 40, -1, 2 );
  effvsdt->setTitle( "DUT efficiency;TLU-DUT #Deltat [ns];<CMS DUT efficiency>" );

  effvsmm = AIDAProcessor::histogramFactory(this)->
    createProfile2D( "effvsmm", 60, 0, 300, 40, 0, 200, -1, 2 );
  effvsmm->setTitle( "DUT efficiency;telescope x mod 300 [um];telescope y mod 200 [um];efficiency" );


  rffvsxy = AIDAProcessor::histogramFactory(this)->
    createProfile2D( "rffvsxy", 120, -12, 12, 80, -8, 8, -1, 2 );
  rffvsxy->setTitle( "REF efficiency;telescope x [mm];telescope y [mm];DUT REF efficiency" );

  rffvsx = AIDAProcessor::histogramFactory(this)->
    createProfile1D( "rffvsx", 100, -2, 8, -1, 2 );
  rffvsx->setTitle( "REF efficiency;telescope x [mm];DUT REF efficiency" );


  lkAvst = AIDAProcessor::histogramFactory(this)->
    createProfile1D( "lkAvst", 360, 0, 3600, -1, 2 );
  lkAvst->setTitle( "DUT link fraction;time [s];<CMS DUT linking>" );

  ntriHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "ntri", 31, -0.5, 30.5 );
  ntriHisto->setTitle( "telescope triplets;0-1-2 triplets;events" );

  // triplet efficiency w.r.t. DUT:

  cmsxeHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "cmsxe", 240, -12, 12 );
  cmsxeHisto->setTitle( "DUT hits;aligned x [mm];clusters" );

  cmsyeHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "cmsye", 240, -6, 6 );
  cmsyeHisto->setTitle( "DUT hits;aligned y [mm];clusters" );

  cmssxeHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "cmssxe", 200, -500, 500 );
  cmssxeHisto->setTitle( "Pixel + Triplet x;cluster + triplet #Sigmax [#mum];clusters" );

  cmsdyeHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "cmsdye", 500, -500, 500 );
  cmsdyeHisto->setTitle( "Pixel - triplet y;cluster - triplet #Deltay [#mum];clusters" );

  cmsnmHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "cmsnm", 11, -0.5, 10.5 );
  cmsnmHisto->setTitle( "triplets linked to clusters;linked triplets;clusters" );

  trieffvsxy = AIDAProcessor::histogramFactory(this)->
    createProfile2D( "trieffvsxy", 120, -12, 12, 80, -8, 8, -1, 2 );
  trieffvsxy->setTitle( "triplet efficiency;DUT cluster x [mm];DUT cluster y [mm];efficiency" );

  // driplets:

  dz = _planePosition[5] - _planePosition[3];
  dx = 0.01 * dz;

  dx35Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "dx35", 100, -dx, dx );
  dx35Histo->setTitle( "x5-x3, |dy| < 2 mm;x_{5}-x_{3} [mm];hit pairs" );

  dy35Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "dy35", 100, -dx, dx );
  dy35Histo->setTitle( "y5-y3, |dx| < 2 mm;y_{5}-y_{3} [mm];hit pairs" );

  dridxHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "dridx", 100, -100, 100 );
  dridxHisto->setTitle( "driplet dx;x_{4}-x_{m} [#mum];driplets" );

  dridyHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "dridy", 100, -100, 100 );
  dridyHisto->setTitle( "driplet dy;y_{4}-y_{m} [#mum];driplets" );

  dridxvsx = AIDAProcessor::histogramFactory(this)->
    createProfile1D( "dridxvsx", 100, -10, 10, -100, 100 );
  dridxvsx->setTitle( "driplet x resid vs x;x [mm];driplet <#Deltax> [#mum]" );

  dridxvsy = AIDAProcessor::histogramFactory(this)->
    createProfile1D( "dridxvsy", 50, -5, 5, -100, 100 );
  dridxvsy->setTitle( "driplet x resid vs y;y [mm];driplet <#Deltax> [#mum]" );

  dridxvstx = AIDAProcessor::histogramFactory(this)->
    createProfile1D( "dridxvstx", 80, -2, 2, -100, 100 );
  dridxvstx->setTitle( "driplet x resid vs tx;t_{x} [mrad];driplet <#Deltax> [#mum]" );

  dridxvsty = AIDAProcessor::histogramFactory(this)->
    createProfile1D( "dridxvsty", 80, -2, 2, -100, 100 );
  dridxvsty->setTitle( "driplet x resid vs ty;t_{y} [mrad];driplet <#Deltax> [#mum]" );

  dridyvsx = AIDAProcessor::histogramFactory(this)->
    createProfile1D( "dridyvsx", 100, -10, 10, -100, 100 );
  dridyvsx->setTitle( "driplet y resid vs x;x [mm];driplet <#Deltay> [#mum]" );

  dridyvsy = AIDAProcessor::histogramFactory(this)->
    createProfile1D( "dridyvsy", 50, -5, 5, -100, 100 );
  dridyvsy->setTitle( "driplet y resid vs y;y [mm];driplet <#Deltay> [#mum]" );

  dridyvstx = AIDAProcessor::histogramFactory(this)->
    createProfile1D( "dridyvstx", 80, -2, 2, -100, 100 );
  dridyvstx->setTitle( "driplet y resid vs tx;t_{x} [mrad];driplet <#Deltay> [#mum]" );

  dridyvsty = AIDAProcessor::histogramFactory(this)->
    createProfile1D( "dridyvsty", 80, -2, 2, -100, 100 );
  dridyvsty->setTitle( "driplet y resid vs ty;t_{y} [mrad];driplet <#Deltay> [#mum]" );

  drixHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "drix", 120, -12, 12 );
  drixHisto->setTitle( "driplet x4;x4 [mm];driplets" );

  driyHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "driy", 80, -8, 8 );
  driyHisto->setTitle( "driplet y4;y4 [mm];driplets" );

  drixyHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram2D( "drixy", 240, -12, 12, 160, -8, 8 );
  drixyHisto->setTitle( "driplet y4 vs x4;x4 [mm];y4 [mm];driplets" );

  dritxHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "dritx", 100, -10, 10 );
  dritxHisto->setTitle( "driplet slope x;#theta_{x} [mrad];driplets" );

  drityHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "drity", 100, -10, 10 );
  drityHisto->setTitle( "driplet slope y;#theta_{y} [mrad];driplets" );


  refxHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "refx", 120, -12, 12 );
  refxHisto->setTitle( "driplet x at REF;x_{REF} [mm];driplets" );

  refyHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "refy", 80, -8, 8 );
  refyHisto->setTitle( "driplet y at REF;y_{REF} [mm];driplets" );

  refxyHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram2D( "refxy", 120, -12, 12, 80, -8, 8 );
  refxyHisto->setTitle( "driplet x-y at REF;x_{REF} [mm];y_{REF} [mm];driplets" );


  drixcmsHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "drixcms", 240, -12, 12 );
  drixcmsHisto->setTitle( "driplet x5 with REF link;x5 [mm];linked driplets" );

  driycmsHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "driycms", 120, -6, 6 );
  driycmsHisto->setTitle( "driplet y5 with REF link;y5 [mm];linked driplets" );

  refxcmsHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "refxcms", 120, -12, 12 );
  refxcmsHisto->setTitle( "driplet x at REF with hit;x_{REF} [mm];linked driplets" );

  refycmsHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "refycms", 120, -12, 12 );
  refycmsHisto->setTitle( "driplet y at REF with hit;y_{REF} [mm];linked driplets" );

  refxycmsHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram2D( "refxycms", 120, -12, 12, 80, -8, 8 );
  refxycmsHisto->setTitle( "driplet x-y at REF with hit;x_{REF} [mm];y_{REF} [mm];linked driplets" );

  refpixvsxy = AIDAProcessor::histogramFactory(this)->
    createHistogram2D( "refpixvsxy", 60, 0, 300, 40, 0, 200 );
  refpixvsxy->setTitle( "REF pixel occupancy;x mod 300 #mum;y mod 200 #mum;clusters" );

  refqHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "refq", 100, 0, 100 );
  refqHisto->setTitle( "REF cluster charge linked;REF cluster charge [ke];REF linked clusters" );

  refqvsxy = AIDAProcessor::histogramFactory(this)->
    createProfile2D( "refqvsxy", 60, 0, 300, 40, 0, 200, 0, 250 );
  refqvsxy->setTitle( "REF cluster charge map;x mod 300 #mum;y mod 200 #mum;<cluster charge> [ke]" );


  refxxHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram2D( "refxx", 52, -0.5, 51.5, 110, -11, 11 );
  refxxHisto->setTitle( "x correlation;REF pixel cluster col;driplet x [mm];REF clusters" );

  refyyHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram2D( "refyy", 80, -0.5, 79.5, 55, -5.5, 5.5 );
  refyyHisto->setTitle( "y correlation;REF pixel cluster row;driplet y [mm];REF clusters" );

  refsxaHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "refsxa", 100, -10, 10 );
  refsxaHisto->setTitle( "REF + driplet x;cluster + driplet #Sigmax [mm];REF clusters" );

  refdyaHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "refdya", 100, -10, 10 );
  refdyaHisto->setTitle( "REF - driplet y;cluster - driplet #Deltay [mm];REF clusters" );

  refsxHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "refsx", 100, -5, 5 );
  refsxHisto->setTitle( "REF + driplet x;cluster + driplet #Sigmax [mm];REF clusters" );

  refdyHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "refdy", 100, -5, 5 );
  refdyHisto->setTitle( "REF - driplet y;cluster - driplet #Deltay [mm];REF clusters" );

  refsxcHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "refsxc", 100, -1, 1 );
  refsxcHisto->setTitle( "REF + driplet x;cluster + driplet #Sigmax [mm];REF clusters" );

  refdycHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "refdyc", 100, -1, 1 );
  refdycHisto->setTitle( "REF - driplet y;cluster - driplet #Deltay [mm];REF clusters" );

  refdyvsx = AIDAProcessor::histogramFactory(this)->
    createProfile1D( "refdyvsx", 52, -26*0.15, 26*0.15, -500, 500 );
  refdyvsx->setTitle( "REF y resid vs x;x [mm];<REF cluster - driplet #Deltay> [#mum]" );

  refdyvsy = AIDAProcessor::histogramFactory(this)->
    createProfile1D( "refdyvsy", 80, -4, 4, -500, 500 );
  refdyvsy->setTitle( "REF y resid vs y;y [mm];<REF cluster - driplet #Deltay> [#mum]" );

  refdyvsty = AIDAProcessor::histogramFactory(this)->
    createProfile1D( "refdyvsty", 80, -2, 2, -500, 500 );
  refdyvsty->setTitle( "REF y resid vs tet y;#theta_{y} [mrad];<REF cluster - driplet #Deltay> [#mum]" );


  bacsxaHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "bacsxa", 100, -10, 10 );
  bacsxaHisto->setTitle( "DUT + driplet x;DUT cluster + driplet #Sigmax [mm];DUT clusters" );

  bacdyaHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "bacdya", 100, -10, 10 );
  bacdyaHisto->setTitle( "DUT - driplet y;DUT cluster - driplet #Deltay [mm];DUT clusters" );

  bacsxcHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "bacsxc", 200, -500, 500 );
  bacsxcHisto->setTitle( "DUT + driplet x;DUT cluster + driplet #Sigmax [mm];DUT clusters" );

  bacdycHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "bacdyc", 200, -500, 500 );
  bacdycHisto->setTitle( "DUT - driplet y;DUT cluster - driplet #Deltay [mm];DUT clusters" );

  bacsxcqHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "bacsxcq", 200, -500, 500 );
  bacsxcqHisto->setTitle( "DUT + driplet x;DUT cluster + driplet #Sigmax [mm];DUT clusters" );

  bacdycqHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "bacdycq", 100, -500, 500 );
  bacdycqHisto->setTitle( "DUT - driplet y;DUT cluster - driplet #Deltay [mm];DUT clusters" );


  ndriHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "ndri", 31, -0.5, 30.5 );
  ndriHisto->setTitle( "telescope driplets;3-4-5 driplets;events" );

  ndricmsHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "ndricms", 31, -0.5, 30.5 );
  ndricmsHisto->setTitle( "telescope driplets and DUT hit in event;3-4-5 driplets and DUT;events" );

  lkBvst = AIDAProcessor::histogramFactory(this)->
    createProfile1D( "lkBvst", 360, 0, 3600, -1, 2 );
  lkBvst->setTitle( "REF link fraction;time [s];<CMS REF linking>" );

  //driplets-triplets

  sixkxHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "sixkx", 100, -10, 10 );
  sixkxHisto->setTitle( "kink x;kink x [mrad];triplet-driplet pairs" );

  sixkyHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "sixky", 100, -10, 10 );
  sixkyHisto->setTitle( "kink y;kink y [mrad];triplet-driplet pairs" );

  sixdxHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "sixdx", 100, -1, 1 );
  sixdxHisto->setTitle( "six match x;match x [mm];triplet-driplet pairs" );

  sixdyHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "sixdy", 100, -1, 1 );
  sixdyHisto->setTitle( "six match y;match y [mm];triplet-driplet pairs" );

  sixdxcHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "sixdxc", 100, -250, 250 );
  sixdxcHisto->setTitle( "six match x;track #Deltax[#mum];triplet-driplet pairs" );

  sixdycHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "sixdyc", 100, -250, 250 );
  sixdycHisto->setTitle( "six match y;track #Deltay[#mum];triplet-driplet pairs" );

  sixkxcHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "sixkxc", 100, -10, 10 );
  sixkxcHisto->setTitle( "kink x, x-y matched;kink x [mrad];triplet-driplet pairs" );

  sixkycHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "sixkyc", 100, -10, 10 );
  sixkycHisto->setTitle( "kink y, x-y matched;kink y [mrad];triplet-driplet pairs" );

  sixkxsHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "sixkxs", 100, -10, 10 );
  sixkxsHisto->setTitle( "kink x at DUT sensor;kink x [mrad];triplet-driplet pairs" );

  sixkysHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "sixkys", 100, -10, 10 );
  sixkysHisto->setTitle( "kink y at DUT sensor;kink y [mrad];triplet-driplet pairs" );

  sixxyHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram2D( "sixxy", 120, -12, 12, 80, -8, 8 );
  sixxyHisto->setTitle( "match x-y;x(zDUT) [mm];y(zDUT) [mm]" );

  sixxycHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram2D( "sixxyc", 120, -12, 12, 80, -8, 8 );
  sixxycHisto->setTitle( "match x-y large kink;x(zDUT) [mm];y(zDUT) [mm]" );

  sixxylkHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram2D( "sixxylk", 120, -12, 12, 80, -8, 8 );
  sixxylkHisto->setTitle( "match x-y REF link;x(zDUT) [mm];y(zDUT) [mm]" );

  kinkvsxy = AIDAProcessor::histogramFactory(this)->
    createProfile2D( "kinkvsxy", 120, -12, 12, 80, -8, 8, 0, 100 );
  kinkvsxy->setTitle( "kink;x(zDUT) [mm];y(zDUT) [mm];<kink^{2}> [mrad^{2}]" );

  // GBL:

  selxHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "selx", 120, -12, 12 );
  selxHisto->setTitle( "x at DUT, sel GBL;x [mm];tracks" );

  selyHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "sely", 80, -8, 8 );
  selyHisto->setTitle( "y at DUT, sel GBL;y [mm];tracks" );

  selaxHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "selax", 100, -5, 5 );
  selaxHisto->setTitle( "track angle x, sel GBL;x angle [mrad];tracks" );

  selayHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "selay", 100, -5, 5 );
  selayHisto->setTitle( "track angle y, sel GBL;y angle [mrad];tracks" );

  seldxHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "seldx", 100, -150, 150 );
  seldxHisto->setTitle( "track match x, sel GBL;#Deltax [#mum];tracks" );

  seldyHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "seldy", 100, -150, 150 );
  seldyHisto->setTitle( "track match y, sel GBL;#Deltay [#mum];tracks" );

  selkxHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "selkx", 100, -10, 10 );
  selkxHisto->setTitle( "kink x, sel GBL;kink x [mrad];tracks" );

  selkyHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "selky", 100, -10, 10 );
  selkyHisto->setTitle( "kink y, sel GBL;kink y [mrad];tracks" );

  seldx1Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "seldx1", 100, -100, 100 );
  seldx1Histo->setTitle( "triplet resid x at 1, sel GBL;#Deltax [#mum];tracks" );

  seldy1Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "seldy1", 100, -100, 100 );
  seldy1Histo->setTitle( "triplet resid y at 1, sel GBL;#Deltay [#mum];tracks" );

  seldx3Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "seldx3", 100, -1000, 1000 );
  seldx3Histo->setTitle( "triplet resid x at 3, sel GBL;#Deltax [#mum];tracks" );

  seldy3Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "seldy3", 100, -1000, 1000 );
  seldy3Histo->setTitle( "triplet resid y at 3, sel GBL;#Deltay [#mum];tracks" );

  seldx4Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "seldx4", 100, -1500, 1500 );
  seldx4Histo->setTitle( "triplet resid x at 4, sel GBL;#Deltax [#mum];tracks" );

  seldy4Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "seldy4", 100, -1500, 1500 );
  seldy4Histo->setTitle( "triplet resid y at 4, sel GBL;#Deltay [#mum];tracks" );

  seldx5Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "seldx5", 100, -3000, 3000 );
  seldx5Histo->setTitle( "triplet resid x at 5, sel GBL;#Deltax [#mum];tracks" );

  seldy5Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "seldy5", 100, -3000, 3000 );
  seldy5Histo->setTitle( "triplet resid y at 5, sel GBL;#Deltay [#mum];tracks" );

  seldx6Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "seldx6", 100, -250, 250 );
  seldx6Histo->setTitle( "triplet resid x at DUT, sel GBL;#Deltax [#mum];tracks" );

  seldy6Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "seldy6", 100, -250, 250 );
  seldy6Histo->setTitle( "triplet resid y at DUT, sel GBL;#Deltay [#mum];tracks" );

  gblndfHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "gblndf", 16, -0.5, 15.5 );
  gblndfHisto->setTitle( "GBL fit NDF;GBL NDF;tracks" );

  gblchi2aHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "gblchi2a", 100, 0, 100 );
  gblchi2aHisto->setTitle( "GBL fit chi2, DoF 8;GBL chi2;tracks" );

  gblchi2bHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "gblchi2b", 100, 0, 100 );
  gblchi2bHisto->setTitle( "GBL fit chi2;GBL chi2;tracks" );

  gblprbHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "gblprb", 100, 0, 1 );
  gblprbHisto->setTitle( "GBL fit probability;GBL fit probability;tracks" );

  // bad fits:

  badxHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "badx", 120, -12, 12 );
  badxHisto->setTitle( "x at DUT, bad GBL;x [mm];tracks" );

  badyHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "bady", 80, -8, 8 );
  badyHisto->setTitle( "y at DUT, bad GBL;y [mm];tracks" );

  badaxHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "badax", 100, -5, 5 );
  badaxHisto->setTitle( "track angle x, bad GBL;x angle [mrad];tracks" );

  badayHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "baday", 100, -5, 5 );
  badayHisto->setTitle( "track angle y, bad GBL;y angle [mrad];tracks" );

  baddxHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "baddx", 100, -150, 150 );
  baddxHisto->setTitle( "track match x, bad GBL;#Deltax [#mum];tracks" );

  baddyHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "baddy", 100, -150, 150 );
  baddyHisto->setTitle( "track match y, bad GBL;#Deltay [#mum];tracks" );

  badkxHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "badkx", 100, -10, 10 );
  badkxHisto->setTitle( "kink x, bad GBL;kink x [mrad];tracks" );

  badkyHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "badky", 100, -10, 10 );
  badkyHisto->setTitle( "kink y, bad GBL;kink y [mrad];tracks" );

  baddx1Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "baddx1", 100, -100, 100 );
  baddx1Histo->setTitle( "triplet resid x at 1, bad GBL;#Deltax [#mum];tracks" );

  baddy1Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "baddy1", 100, -100, 100 );
  baddy1Histo->setTitle( "triplet resid y at 1, bad GBL;#Deltay [#mum];tracks" );

  baddx3Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "baddx3", 100, -1000, 1000 );
  baddx3Histo->setTitle( "triplet resid x at 3, bad GBL;#Deltax [#mum];tracks" );

  baddy3Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "baddy3", 100, -1000, 1000 );
  baddy3Histo->setTitle( "triplet resid y at 3, bad GBL;#Deltay [#mum];tracks" );

  baddx4Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "baddx4", 100, -1500, 1500 );
  baddx4Histo->setTitle( "triplet resid x at 4, bad GBL;#Deltax [#mum];tracks" );

  baddy4Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "baddy4", 100, -1500, 1500 );
  baddy4Histo->setTitle( "triplet resid y at 4, bad GBL;#Deltay [#mum];tracks" );

  baddx5Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "baddx5", 100, -3000, 3000 );
  baddx5Histo->setTitle( "triplet resid x at 5, bad GBL;#Deltax [#mum];tracks" );

  baddy5Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "baddy5", 100, -3000, 3000 );
  baddy5Histo->setTitle( "triplet resid y at 5, bad GBL;#Deltay [#mum];tracks" );

  baddx6Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "baddx6", 100, -250, 250 );
  baddx6Histo->setTitle( "triplet resid x at DUT, bad GBL;#Deltax [#mum];tracks" );

  baddy6Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "baddy6", 100, -250, 250 );
  baddy6Histo->setTitle( "triplet resid y at DUT, bad GBL;#Deltay [#mum];tracks" );

  // good fits:

  goodx1Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "goodx1", 100, -100, 100 );
  goodx1Histo->setTitle( "triplet resid x at 1, good GBL;#Deltax [#mum];tracks" );

  goody1Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "goody1", 100, -100, 100 );
  goody1Histo->setTitle( "triplet resid y at 1, good GBL;#Deltay [#mum];tracks" );

  goodx6Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "goodx6", 100, -250, 250 );
  goodx6Histo->setTitle( "triplet resid x at 6, good GBL;#Deltax [#mum];tracks" );

  goody6Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "goody6", 100, -250, 250 );
  goody6Histo->setTitle( "triplet resid y at 6, good GBL;#Deltay [#mum];tracks" );

  // look at fit:

  gblax0Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "gblax0", 100, -1, 1 );
  gblax0Histo->setTitle( "GBL angle at plane 0;x angle at plane 0 [mrad];tracks" );

  gbldx0Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "gbldx0", 100, -10, 10 );
  gbldx0Histo->setTitle( "GBL shift at plane 0;x shift at plane 0 [#mum];tracks" );

  gblrx0Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "gblrx0", 100, -25, 25 );
  gblrx0Histo->setTitle( "GBL resid at plane 0;x resid at plane 0 [#mum];tracks" );

  gblpx0Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "gblpx0", 100, -10, 10 );
  gblpx0Histo->setTitle( "GBL pull at plane 0;x pull at plane 0 [#sigma];tracks" );

  gblqx0Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "gblqx0", 100, -1, 1 );
  gblqx0Histo->setTitle( "GBL kink at plane 0;x kink at plane 0 [mrad];tracks" );


  gblax1Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "gblax1", 100, -1, 1 );
  gblax1Histo->setTitle( "GBL angle at plane 1;x angle at plane 1 [mrad];tracks" );

  gbldx1Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "gbldx1", 100, -100, 100 );
  gbldx1Histo->setTitle( "GBL shift at plane 1;x shift at plane 1 [#mum];tracks" );

  gblrx1Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "gblrx1", 100, -25, 25 );
  gblrx1Histo->setTitle( "GBL resid at plane 1;x resid at plane 1 [#mum];tracks" );

  gblpx1Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "gblpx1", 100, -10, 10 );
  gblpx1Histo->setTitle( "GBL pull at plane 1;x pull at plane 1 [#sigma];tracks" );

  gblqx1Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "gblqx1", 100, -1, 1 );
  gblqx1Histo->setTitle( "GBL kink at plane 1;x kink at plane 1 [mrad];tracks" );


  gblax2Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "gblax2", 100, -1, 1 );
  gblax2Histo->setTitle( "GBL angle at plane 2;x angle at plane 2 [mrad];tracks" );

  gbldx2Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "gbldx2", 100, -20, 20 );
  gbldx2Histo->setTitle( "GBL shift at plane 2;x shift at plane 2 [#mum];tracks" );

  gblrx2Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "gblrx2", 100, -25, 25 );
  gblrx2Histo->setTitle( "GBL resid at plane 2;x resid at plane 2 [#mum];tracks" );

  gblpx2Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "gblpx2", 100, -10, 10 );
  gblpx2Histo->setTitle( "GBL pull at plane 2;x pull at plane 2 [#sigma];tracks" );

  gblqx2Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "gblqx2", 100, -1, 1 );
  gblqx2Histo->setTitle( "GBL kink at plane 2;x kink at plane 2 [mrad];tracks" );


  gblax3Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "gblax3", 100, -1, 1 );
  gblax3Histo->setTitle( "GBL angle at plane 3;x angle at plane 3 [mrad];tracks" );

  gbldx3Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "gbldx3", 100, -250, 250 );
  gbldx3Histo->setTitle( "GBL shift at plane 3;x shift at plane 3 [#mum];tracks" );

  gblrx3Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "gblrx3", 100, -25, 25 );
  gblrx3Histo->setTitle( "GBL resid at plane 3;x resid at plane 3 [#mum];tracks" );

  gblpx3Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "gblpx3", 100, -10, 10 );
  gblpx3Histo->setTitle( "GBL pull at plane 3;x pull at plane 3 [#sigma];tracks" );

  gblqx3Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "gblqx3", 100, -1, 1 );
  gblqx3Histo->setTitle( "GBL kink at plane 3;x kink at plane 3 [mrad];tracks" );


  gblax4Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "gblax4", 100, -1, 1 );
  gblax4Histo->setTitle( "GBL angle at plane 4;x angle at plane 4 [mrad];tracks" );

  gbldx4Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "gbldx4", 100, -500, 500 );
  gbldx4Histo->setTitle( "GBL shift at plane 4;x shift at plane 4 [#mum];tracks" );

  gblrx4Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "gblrx4", 100, -25, 25 );
  gblrx4Histo->setTitle( "GBL resid at plane 4;x resid at plane 4 [#mum];tracks" );

  gblpx4Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "gblpx4", 100, -10, 10 );
  gblpx4Histo->setTitle( "GBL pull at plane 4;x pull at plane 4 [#sigma];tracks" );

  gblqx4Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "gblqx4", 100, -1, 1 );
  gblqx4Histo->setTitle( "GBL kink at plane 4;x kink at plane 4 [mrad];tracks" );


  gblax5Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "gblax5", 100, -1, 1 );
  gblax5Histo->setTitle( "GBL angle at plane 5;x angle at plane 5 [mrad];tracks" );

  gbldx5Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "gbldx5", 100, -1000, 1000 );
  gbldx5Histo->setTitle( "GBL shift at plane 5;x shift at plane 5 [#mum];tracks" );

  gblrx5Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "gblrx5", 100, -25, 25 );
  gblrx5Histo->setTitle( "GBL resid at plane 5;x resid at plane 5 [#mum];tracks" );

  gblpx5Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "gblpx5", 100, -10, 10 );
  gblpx5Histo->setTitle( "GBL pull at plane 5;x pull at plane 5 [#sigma];tracks" );

  gblqx5Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "gblqx5", 100, -1, 1 );
  gblqx5Histo->setTitle( "GBL kink at plane 5;x kink at plane 5 [mrad];tracks" );


  gblax6Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "gblax6", 100, -1, 1 );
  gblax6Histo->setTitle( "GBL angle at DUT;x angle at DUT [mrad];tracks" );

  gbldx6Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "gbldx6", 100, -1000, 1000 );
  gbldx6Histo->setTitle( "GBL shift at DUT;x shift at DUT [#mum];tracks" );

  gbldy6Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "gbldy6", 100, -1000, 1000 );
  gbldy6Histo->setTitle( "GBL shift at DUT;y shift at DUT [#mum];tracks" );

  gblrx6Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "gblrx6", 100, -250, 250 );
  gblrx6Histo->setTitle( "GBL resid at DUT;x resid at DUT [#mum];tracks" );

  gblry6Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "gblry6", 100, -100, 100 );
  gblry6Histo->setTitle( "GBL resid at DUT;y resid at DUT [#mum];tracks" );

  gblpx6Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "gblpx6", 100, -10, 10 );
  gblpx6Histo->setTitle( "GBL pull at DUT;x pull at DUT [#sigma];tracks" );

  gblqx6Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "gblqx6", 100, -10, 10 );
  gblqx6Histo->setTitle( "GBL kink at DUT;x kink at DUT [mrad];tracks" );


  gblkx1Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "gblkx1", 100, -1, 1 );
  gblkx1Histo->setTitle( "GBL kink angle at plane 1;plane 1 kink [mrad];tracks" );

  gblkx2Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "gblkx2", 100, -1, 1 );
  gblkx2Histo->setTitle( "GBL kink angle at plane 2;plane 2 kink [mrad];tracks" );

  gblkx3Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "gblkx3", 100, -10, 10 );
  gblkx3Histo->setTitle( "GBL kink angle at plane 3;plane 3 kink [mrad];tracks" );

  gblkx4Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "gblkx4", 100, -1, 1 );
  gblkx4Histo->setTitle( "GBL kink angle at plane 4;plane 4 kink [mrad];tracks" );

  gblkx5Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "gblkx5", 100, -1, 1 );
  gblkx5Histo->setTitle( "GBL kink angle at plane 5;plane 5 kink [mrad];tracks" );

  gblkx6Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "gblkx6", 100, -1, 1 );
  gblkx6Histo->setTitle( "GBL kink angle at plane 6;plane 6 kink [mrad];tracks" );

  // z intersect:

  sixzx3Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "sixzx3", 100, -50, 250 );
  sixzx3Histo->setTitle( "intersect z-x, kink > 3 mrad;intersect z(x) - z_{2} [mm];tracks" );

  sixzy3Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "sixzy3", 100, -50, 250 );
  sixzy3Histo->setTitle( "intersect z-y, kink > 3 mrad;intersect z(y) - z_{2} [mm];tracks" );

  sixzx2Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "sixzx2", 100, -50, 250 );
  sixzx2Histo->setTitle( "intersect z-x, kink > 2 mrad;intersect z(x) - z_{2} [mm];tracks" );

  sixzy2Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "sixzy2", 100, -50, 250 );
  sixzy2Histo->setTitle( "intersect z-y, kink > 2 mrad;intersect z(y) - z_{2} [mm];tracks" );

  sixzx1Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "sixzx1", 100, -50, 250 );
  sixzx1Histo->setTitle( "intersect z-x, kink > 1 mrad;intersect z(x) - z_{2} [mm];tracks" );

  sixzy1Histo = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "sixzy1", 100, -50, 250 );
  sixzy1Histo->setTitle( "intersect z-y, kink > 1 mrad;intersect z(y) - z_{2} [mm];tracks" );

  sixkxzyHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "sixkxzy", 200, -20, 20 );
  sixkxzyHisto->setTitle( "kink x at DUT;kink x [mrad];tracks" );

  sixkyzxHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "sixkyzx", 200, -20, 20 );
  sixkyzxHisto->setTitle( "kink y at DUT;kink y [mrad];tracks" );

  sixkxzxHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "sixkxzx", 200, -20, 20 );
  sixkxzxHisto->setTitle( "kink x at DUT;kink x [mrad];tracks" );

  sixkyzyHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "sixkyzy", 200, -20, 20 );
  sixkyzyHisto->setTitle( "kink y at DUT;kink y [mrad];tracks" );

  // system times:

  cmsdtHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "cmsdt", 500, 0, 5000 );
  cmsdtHisto->setTitle( "DUT time between events;DUT time between events [us];events" );

  sysrtHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "sysrt", 200, 0.9998, 1.0002 );
  sysrtHisto->setTitle( "TLU time / DUT time;event time ratio;events" );

  sysrdtHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "sysrdt", 200, 0.999, 1.001 );
  sysrdtHisto->setTitle( "TLU time / DUT time;time between events ratio;events" );

  sysdtHisto = AIDAProcessor::histogramFactory(this)->
    createHistogram1D( "sysdt", 100, -100, 100 );
  sysdtHisto->setTitle( "TLU - DUT time;TLU - DUT #Deltat [ns];events" );

  sysdtvsdt = AIDAProcessor::histogramFactory(this)->
    createProfile1D( "sysdtvsdt", 100, 0, 10000, -1E9, 1E9 );
  sysdtvsdt->setTitle( "TLU - DUT time;time between events [us];<#Deltat> [ns]" );

  // Chi2 distribution for all accepted tracks( linear scale)

  int    chi2NBin  = 1000;
  double chi2Min   = 0.;
  double chi2Max   = 200.;
  string chi2Title = "Chi2 distribution";

  if( isHistoManagerAvailable ) {
    histoInfo = histoMgr->getHistogramInfo(_linChi2HistoName);
    if( histoInfo ) {
      streamlog_out( DEBUG )  <<( * histoInfo ) << endl;
      chi2NBin = histoInfo->_xBin;
      chi2Min  = histoInfo->_xMin;
      chi2Max  = histoInfo->_xMax;
      if( histoInfo->_title != "" ) chi2Title = histoInfo->_title;
    }
  }

  AIDA::IHistogram1D * linChi2Histo = AIDAProcessor::histogramFactory(this)->createHistogram1D( _linChi2HistoName.c_str(),chi2NBin,chi2Min,chi2Max);
  linChi2Histo->setTitle(chi2Title.c_str());
  _aidaHistoMap.insert(make_pair(_linChi2HistoName, linChi2Histo));

  // log(Chi2) distribution for all accepted tracks

  chi2NBin  = 100;
  chi2Min   = -2.;
  chi2Max   = 8.;
  chi2Title = "log(Chi2) distribution";

  if( isHistoManagerAvailable ) {
    histoInfo = histoMgr->getHistogramInfo(_logChi2HistoName);
    if( histoInfo ) {
      streamlog_out( DEBUG )   <<( * histoInfo ) << endl;
      chi2NBin = histoInfo->_xBin;
      chi2Min  = histoInfo->_xMin;
      chi2Max  = histoInfo->_xMax;
      if( histoInfo->_title != "" ) chi2Title = histoInfo->_title;
    }
  }

  AIDA::IHistogram1D * logChi2Histo = AIDAProcessor::histogramFactory(this)->createHistogram1D( _logChi2HistoName.c_str(),chi2NBin,chi2Min,chi2Max);
  logChi2Histo->setTitle(chi2Title.c_str());
  _aidaHistoMap.insert(make_pair(_logChi2HistoName, logChi2Histo));

  // Additional Chi2 histogram for first track candidate( without chi2
  // cut)

  string firstchi2Title = chi2Title + ", first candidate( before cut)";
  AIDA::IHistogram1D * firstChi2Histo = AIDAProcessor::histogramFactory(this)->createHistogram1D( _firstChi2HistoName.c_str(),chi2NBin,chi2Min,chi2Max);
  firstChi2Histo->setTitle(firstchi2Title.c_str());
  _aidaHistoMap.insert(make_pair(_firstChi2HistoName, firstChi2Histo));

  // Chi2 histogram for best tracks in an event - use same binning
  if(_searchMultipleTracks) {
    string bestchi2Title = chi2Title + ", best fit";

    AIDA::IHistogram1D * bestChi2Histo = AIDAProcessor::histogramFactory(this)->createHistogram1D( _bestChi2HistoName.c_str(),chi2NBin,chi2Min,chi2Max);

    bestChi2Histo->setTitle(bestchi2Title.c_str());

    _aidaHistoMap.insert(make_pair(_bestChi2HistoName, bestChi2Histo));
  }

  // Another Chi2 histogram for all full tracks in an event - use same binning
  if(_allowMissingHits) {
    string fullchi2Title = chi2Title + ", full tracks";
    AIDA::IHistogram1D * fullChi2Histo = AIDAProcessor::histogramFactory(this)->createHistogram1D( _fullChi2HistoName.c_str(),chi2NBin,chi2Min,chi2Max);
    fullChi2Histo->setTitle(fullchi2Title.c_str());
    _aidaHistoMap.insert(make_pair(_fullChi2HistoName, fullChi2Histo));
  }

  // Distribution of number of tracks

  int    trkNBin  = 21;
  double trkMin   = -0.5;
  double trkMax   = 20.5;
  string trkTitle = "Number of reconstructed tracks in an event";

  if( isHistoManagerAvailable ) {
    histoInfo = histoMgr->getHistogramInfo(_nTrackHistoName);
    if( histoInfo ) {
      streamlog_out( DEBUG ) <<( * histoInfo ) << endl;
      trkNBin = histoInfo->_xBin;
      trkMin  = histoInfo->_xMin;
      trkMax  = histoInfo->_xMax;
      if( histoInfo->_title != "" ) trkTitle = histoInfo->_title;
    }
  }

  AIDA::IHistogram1D * nTrackHisto = AIDAProcessor::histogramFactory(this)->createHistogram1D( _nTrackHistoName.c_str(),trkNBin,trkMin,trkMax);
  nTrackHisto->setTitle(trkTitle.c_str());
  _aidaHistoMap.insert(make_pair(_nTrackHistoName, nTrackHisto));

  // Number of hits in input collection
  int    hitNBin  = 200;
  double hitMin   = 0.;
  double hitMax   = 200.;
  string hitTitle = "Number of hits in input collection";

  if( isHistoManagerAvailable ) {
    histoInfo = histoMgr->getHistogramInfo(_nAllHitHistoName);
    if( histoInfo ) {
      streamlog_out( DEBUG ) <<( * histoInfo ) << endl;
      hitNBin = histoInfo->_xBin;
      hitMin  = histoInfo->_xMin;
      hitMax  = histoInfo->_xMax;
      if( histoInfo->_title != "" ) hitTitle = histoInfo->_title;
    }
  }

  AIDA::IHistogram1D * nAllHitHisto = AIDAProcessor::histogramFactory(this)->createHistogram1D( _nAllHitHistoName.c_str(),hitNBin,hitMin,hitMax);
  nAllHitHisto->setTitle(hitTitle.c_str());
  _aidaHistoMap.insert(make_pair(_nAllHitHistoName, nAllHitHisto));

  // Number of accepted hits
  hitNBin  = 200;
  hitMin   = 0.;
  hitMax   = 200.;
  hitTitle = "Number of hits accepted from input collection";

  if( isHistoManagerAvailable ) {
    histoInfo = histoMgr->getHistogramInfo(_nAccHitHistoName);
    if( histoInfo ) {
      streamlog_out( DEBUG ) <<( * histoInfo ) << endl;
      hitNBin = histoInfo->_xBin;
      hitMin  = histoInfo->_xMin;
      hitMax  = histoInfo->_xMax;
      if( histoInfo->_title != "" ) hitTitle = histoInfo->_title;
    }
  }

  AIDA::IHistogram1D * nAccHitHisto = AIDAProcessor::histogramFactory(this)->createHistogram1D( _nAccHitHistoName.c_str(),hitNBin,hitMin,hitMax);
  nAccHitHisto->setTitle(hitTitle.c_str());
  _aidaHistoMap.insert(make_pair(_nAccHitHistoName, nAccHitHisto));

  // Number of hits per track
  hitNBin  = 11;
  hitMin   = -0.5;
  hitMax   = 10.5;
  hitTitle = "Number of hits used to reconstruct the track";

  if( isHistoManagerAvailable ) {
    histoInfo = histoMgr->getHistogramInfo(_nHitHistoName);
      if( histoInfo ) {
	streamlog_out( DEBUG ) <<( * histoInfo ) << endl;
	hitNBin = histoInfo->_xBin;
	hitMin  = histoInfo->_xMin;
	hitMax  = histoInfo->_xMax;
	if( histoInfo->_title != "" ) hitTitle = histoInfo->_title;
      }
  }

  AIDA::IHistogram1D * nHitHisto = AIDAProcessor::histogramFactory(this)->createHistogram1D( _nHitHistoName.c_str(),hitNBin,hitMin,hitMax);
  nHitHisto->setTitle(hitTitle.c_str());
  _aidaHistoMap.insert(make_pair(_nHitHistoName, nHitHisto));

  // Additional hit number histogram for best tracks in an event - use same binning
  if(_searchMultipleTracks) {
    string bestTitle = hitTitle + ", best fit";
    AIDA::IHistogram1D * BestHisto = AIDAProcessor::histogramFactory(this)->createHistogram1D( _nBestHistoName.c_str(),hitNBin,hitMin,hitMax);
    BestHisto->setTitle(bestTitle.c_str());
    _aidaHistoMap.insert(make_pair(_nBestHistoName, BestHisto));
  }

  // Additional histogram for number of tracks fitted to given hit - use same binning
  if(_allowAmbiguousHits ) {
    string ambigTitle = "Number of tracks containing given hit( ambiguity)";
    AIDA::IHistogram1D * AmbigHisto = AIDAProcessor::histogramFactory(this)->createHistogram1D( _hitAmbiguityHistoName.c_str(),hitNBin,hitMin,hitMax);
    AmbigHisto->setTitle(ambigTitle.c_str());
    _aidaHistoMap.insert(make_pair(_hitAmbiguityHistoName, AmbigHisto));
  }

  // List all booked histogram - check of histogram map filling
  streamlog_out( DEBUG ) <<  _aidaHistoMap.size() << " histograms booked" << endl;

  map<string, AIDA::IBaseHistogram *>::iterator mapIter;
  for(mapIter = _aidaHistoMap.begin(); mapIter != _aidaHistoMap.end(); mapIter++ )
    streamlog_out( DEBUG ) <<  mapIter->first << " : " << ( mapIter->second)->title()  << endl;

  streamlog_out( DEBUG ) << "Histogram booking completed \n\n" << endl;

#else

  streamlog_out( MESSAGE4 ) << "No histogram produced because Marlin doesn't use AIDA" << endl;

#endif

  return;
}
//
// ===============================================================================
//
//  Private function members
//

double EUTelTestFitter::MatrixFit()
{
  for( int ipl = 0; ipl < _nTelPlanes; ipl++ ) {
    _fitX[ipl] = _planeX[ipl];
    _fitEx[ipl] = _planeEx[ipl];
    _fitY[ipl] = _planeY[ipl];
    _fitEy[ipl] = _planeEy[ipl];
  }

  int status = DoAnalFit( _fitX, _fitEx, _beamSlopeX );

  if( status ) return -1.;

  status = DoAnalFit( _fitY, _fitEy, _beamSlopeY );

  if( status ) return -1.;

  double chi2 = GetFitChi2();

  return chi2;
}


double EUTelTestFitter::SingleFit()
{
  for( int ipl = 0; ipl < _nTelPlanes; ipl++ ) {
    _fitX[ipl] = _planeX[ipl];
    _fitEx[ipl] = _planeEx[ipl];
  }

  int status = DoAnalFit( _fitX, _fitEx, _beamSlopeX );

  if( status ) return -1.;

  // Use same matrix to solve equation in Y

  for( int ipl = 0; ipl < _nTelPlanes; ipl++ ) {
    _fitEy[ipl] = _fitEx[ipl];

    _fitY[ipl] =  0.;
    for( int jpl = 0; jpl < _nTelPlanes; jpl++ )
      if( _planeEy[jpl]> 0. )
	_fitY[ipl] += _fitArray[ipl+jpl*_nTelPlanes]*_planeY[jpl]/_planeEy[jpl]/_planeEy[jpl];
  }

  double chi2 = GetFitChi2();

  return chi2;
}


double EUTelTestFitter::NominalFit()
{
  for( int ipl = 0; ipl < _nTelPlanes; ipl++ ) {
    _fitEx[ipl] = _nominalErrorX[ipl];
    _fitEy[ipl] = _nominalErrorY[ipl];

    _fitX[ipl] = 0.;
    _fitY[ipl] = 0.;

    for( int jpl = 0; jpl<_nTelPlanes;jpl++) {
      if( _planeEx[jpl] > 0. )
	_fitX[ipl] += _nominalFitArrayX[ipl+jpl*_nTelPlanes]*_planeX[jpl]/_planeEx[jpl]/_planeEx[jpl];
      if( _planeEy[jpl] > 0. )
	_fitY[ipl] += _nominalFitArrayY[ipl+jpl*_nTelPlanes]*_planeY[jpl]/_planeEy[jpl]/_planeEy[jpl];
    }

    // Correction for beam slope

    if( _useBeamConstraint && _beamSlopeX != 0. ) {
      _fitX[ipl] -= _nominalFitArrayX[ipl]*_beamSlopeX*_planeDist[0]*_planeScat[0];
      _fitX[ipl] += _nominalFitArrayX[ipl+_nTelPlanes]*_beamSlopeX*_planeDist[0]*_planeScat[0];
    }

    if( _useBeamConstraint && _beamSlopeY != 0. ) {
      _fitY[ipl] -= _nominalFitArrayY[ipl]*_beamSlopeY*_planeDist[0]*_planeScat[0];
      _fitY[ipl] += _nominalFitArrayY[ipl+_nTelPlanes]*_beamSlopeY*_planeDist[0]*_planeScat[0];
    }

  }

  double chi2 = GetFitChi2();

  return chi2;
}


int EUTelTestFitter::DoAnalFit( double * pos, double *err, double slope )
{
  for( int ipl = 0; ipl < _nTelPlanes; ipl++ ) {
    if( _isActive[ipl] && err[ipl] > 0 )
      err[ipl] = 1./err[ipl]/err[ipl];
    else
      err[ipl] = 0.;

    pos[ipl] *= err[ipl];
  }

  // To take into account beam tilt

  if( _useBeamConstraint && slope != 0. ) {
    pos[0] -= slope*_planeDist[0]*_planeScat[0];
    pos[1] += slope*_planeDist[0]*_planeScat[0];
  }

  for( int ipl = 0; ipl < _nTelPlanes; ipl++ ) {
    for( int jpl = 0; jpl < _nTelPlanes; jpl++ ) {
      int imx = ipl + jpl * _nTelPlanes;

      _fitArray[imx] = 0.;

      if( jpl == ipl-2 )
	_fitArray[imx] += _planeDist[ipl-2]*_planeDist[ipl-1]*_planeScat[ipl-1];

      if( jpl == ipl+2 )
	_fitArray[imx] += _planeDist[ipl]*_planeDist[ipl+1]*_planeScat[ipl+1];

      if( jpl == ipl-1 ) {
	if( ipl>0 &&  ipl < _nTelPlanes-1 )
	  _fitArray[imx] -= _planeDist[ipl-1]*(_planeDist[ipl]+_planeDist[ipl-1])*_planeScat[ipl];
	if( ipl > 1 )
	  _fitArray[imx] -= _planeDist[ipl-1]*(_planeDist[ipl-1]+_planeDist[ipl-2])*_planeScat[ipl-1];
      }

      if( jpl == ipl + 1 ) {
	if( ipl > 0 && ipl < _nTelPlanes-1 )
	  _fitArray[imx] -= _planeDist[ipl]*(_planeDist[ipl]+_planeDist[ipl-1])*_planeScat[ipl];
	if( ipl < _nTelPlanes-2 )
	  _fitArray[imx] -= _planeDist[ipl]*(_planeDist[ipl+1]+_planeDist[ipl])*_planeScat[ipl+1];
      }

      if( jpl == ipl ) {
	_fitArray[imx] += err[ipl];

	if(ipl > 0 && ipl < _nTelPlanes-1 )
	  _fitArray[imx] += _planeScat[ipl]*(_planeDist[ipl]+_planeDist[ipl-1])*(_planeDist[ipl]+_planeDist[ipl-1]);

	if( ipl > 1 )
	  _fitArray[imx] += _planeScat[ipl-1]*_planeDist[ipl-1]*_planeDist[ipl-1];

	if( ipl < _nTelPlanes-2 )
	  _fitArray[imx] += _planeScat[ipl+1]*_planeDist[ipl]*_planeDist[ipl];
      }

      // For beam constraint

      if( ipl == jpl && ipl < 2 && _useBeamConstraint )
	_fitArray[imx] += _planeScat[0]*_planeDist[0]*_planeDist[0];

      if( ipl+jpl == 1 && _useBeamConstraint )
	_fitArray[imx] -= _planeScat[0]*_planeDist[0]*_planeDist[0];        
    }
  }

  int status = GaussjSolve( _fitArray,pos, _nTelPlanes );

  if( status ) {
    cerr << "Singular matrix in track fitting algorithm ! " << endl;
    for( int ipl = 0; ipl < _nTelPlanes; ipl++ )
      err[ipl] = 0.;
  }
  else
    for( int ipl = 0; ipl < _nTelPlanes; ipl++ )
      err[ipl] = sqrt( _fitArray[ipl + ipl*_nTelPlanes] );

  return status;
}


double EUTelTestFitter::GetFitChi2()
{
  double chi2 = 0.;

  // Measurements

  for( int ipl = 0; ipl<_nTelPlanes;ipl++)
    if( _isActive[ipl] ) {
      if( _planeEx[ipl] > 0. )
	chi2 += (_fitX[ipl]-_planeX[ipl])*(_fitX[ipl]-_planeX[ipl])/_planeEx[ipl]/_planeEx[ipl];

      if( _planeEy[ipl] > 0. )
	chi2 += (_fitY[ipl]-_planeY[ipl])*(_fitY[ipl]-_planeY[ipl])/_planeEy[ipl]/_planeEy[ipl];

    }

  // Scattering angles
  // Use approximate formulas, corresponding to the approximation
  // used in fitting algorithm

  for( int ipl = 1; ipl<_nTelPlanes-1; ipl++ ) {
    double th1,th2,dth;

    th2 = (_fitX[ipl+1]-_fitX[ipl])*_planeDist[ipl];
    th1 = (_fitX[ipl]-_fitX[ipl-1])*_planeDist[ipl-1];
    //    dth=atan(th2)-atan(th1);
    dth = th2 - th1;
    chi2 += _planeScat[ipl] * dth * dth;

    th2 = (_fitY[ipl+1]-_fitY[ipl])*_planeDist[ipl];
    th1 = (_fitY[ipl]-_fitY[ipl-1])*_planeDist[ipl-1];
    //    dth=atan(th2)-atan(th1);
    dth = th2 - th1;
    chi2 += _planeScat[ipl] * dth * dth;
  }

  // Beam constraint
  // Beam slope taken into account.

  if( _useBeamConstraint ) {
    double dth;

    // Use small angle approximation: atan(x) = x
    // Should be:
    // dth=atan((_fitX[1]-_fitX[0])*_planeDist[0]);
      
    dth = (_fitX[1]-_fitX[0])*_planeDist[0] - _beamSlopeX;
    chi2 += _planeScat[0] * dth * dth;

    // dth=atan((_fitY[1]-_fitY[0])*_planeDist[0]);

    dth=(_fitY[1]-_fitY[0])*_planeDist[0] - _beamSlopeY;
    chi2 += _planeScat[0] * dth * dth;
  }

  return chi2;
}


int EUTelTestFitter::GaussjSolve( double *alfa, double *beta, int n )
{
  int *ipiv;
  int *indxr;
  int *indxc;
  int i,j,k;
  int irow = 0;
  int icol = 0;
  double abs,big,help,pivinv;

  ipiv = new int[n];
  indxr = new int[n];
  indxc = new int[n];

  for( i = 0;i<n;i++)ipiv[i] = 0;

  for( i = 0;i<n;i++) {
    big = 0.;
    for(j = 0;j<n;j++) {
      if(ipiv[j]==1)continue;
      for(k = 0;k<n;k++) {
	if(ipiv[k]!= 0)continue;
	abs=fabs(alfa[n*j+k]);
	if(abs>big) {
	  big=abs;
	  irow=j;
	  icol=k;
	}
      }
    }
    ipiv[icol]++;

    if(ipiv[icol]>1)
      return 1;

    if(irow!=icol) {
      help=beta[irow];
      beta[irow]=beta[icol];
      beta[icol]=help;
      for(j = 0;j<n;j++) {
	help=alfa[n*irow+j];
	alfa[n*irow+j]=alfa[n*icol+j];
	alfa[n*icol+j]=help;
      }
    }
    indxr[i]=irow;
    indxc[i]=icol;

    if(alfa[n*icol+icol] == 0.)
      return 1;

    help=alfa[n*icol+icol];
    pivinv=1./help;
    alfa[n*icol+icol]=1.;
    for(j = 0;j<n;j++) alfa[n*icol+j]*=pivinv;

    beta[icol]*=pivinv;

    for(j = 0;j<n;j++) {
      if(j==icol)continue;
      help=alfa[n*j+icol];
      alfa[n*j+icol] = 0.;
      for(k = 0;k<n;k++)
	alfa[n*j+k]-=alfa[n*icol+k]*help;
      beta[j]-=beta[icol]*help;
    }
  }

  for( i = n-1; i >= 0; i-- ) {
    if( indxr[i] == indxc[i] )continue;
    for( j = 0; j < n; j++ ) {
      help = alfa[ n*j + indxr[i] ];
      alfa[n*j+indxr[i]] = alfa[n*j+indxc[i]];
      alfa[n*j+indxc[i]] = help;
    }
  }

  delete [] ipiv;
  delete [] indxr;
  delete [] indxc;

  return 0;
}


void EUTelTestFitter::getFastTrackImpactPoint(double & x, double & y, double & z, Track * tr, LCEvent * ev) {

  // given maps:
  //            _siPlaneCenter.insert( make_pair( sensorID, _planeCenter ));
  //            _siPlaneNormal.insert( make_pair( sensorID, _planeNormal ));

  // this is a simplified version of getTrackImpactPoint
  // 
  int sensorID = guessSensorID( x, y, z );

  //
  // assume no slope for a track
  //
  double slopeX = 0.;
  double slopeY = 0.;

  double offsetX = x;
  double offsetY = y;   

  // this is a normal vector to the plane
  TVectorD	NormalVector(3);
  NormalVector(0) = _siPlaneNormal[sensorID][0];
  NormalVector(1) = _siPlaneNormal[sensorID][1];
  NormalVector(2) = _siPlaneNormal[sensorID][2];

  // this is a vector to the point in the plane
  TVectorD    r0Vector(3);
  r0Vector(0) = _siPlaneCenter[sensorID][0];
  r0Vector(1) = _siPlaneCenter[sensorID][1];
  r0Vector(2) = _siPlaneCenter[sensorID][2];

  // now have to solve the equation
  TVectorD	trackImpact(3);
  TMatrixD	equationMatrix(3,3);
  TVectorD	b(3);

  equationMatrix(0, 0) = NormalVector( 0);
  equationMatrix(0, 1) = NormalVector( 1);
  equationMatrix(0, 2) = NormalVector( 2);
  equationMatrix(1, 0) = 1;
  equationMatrix(1, 1) = 0;
  equationMatrix(1, 2) =( -1)*slopeX;
  equationMatrix(2, 0) = 0;
  equationMatrix(2, 1) = 1;
  equationMatrix(2, 2) =( -1)*slopeY;

  b(0) = r0Vector(0) * NormalVector( 0) + r0Vector(1) * NormalVector( 1) + r0Vector(2) * NormalVector( 2);
  b(1) = offsetX;
  b(2) = offsetY;

  trackImpact = equationMatrix.Invert() * b;

/*
  cout << trackImpact(0) << " ** " << x << endl;
  cout << trackImpact(1) << " ** " << y << endl;
  cout << trackImpact(2) << " ** " << z << endl;
*/

  x = trackImpact(0);
  y = trackImpact(1);
  z = trackImpact(2);

}


void EUTelTestFitter::getTrackImpactPoint(double & x, double & y, double & z, Track * tr, LCEvent * ev) {

    // what should be here:
    // get the center of sensor( X,Y,Z) and normal vector( a,b,c)
    // based on Gear information and all available alignment collections 
    // at the moment Millepede "alignment" and offset "preAlignment" collections
    //
    // solve linear algebra
    //  intersection of sensor Plane and Track
    //
    //  input: 
    //        Track   *tr
    //        LCEvent *ev
    //  output:
    //        double x,y,z
    //
    //
    
    int iplane = guessSensorID( tr );

    // angles only!
    TVector3 _GEAR_(0.,0.,0.);
    TVector3 _PreAlignment_(0.,0.,0.);
    TVector3 _Alignment_(0.,0.,0.);

    _GEAR_.SetX(   _siPlanesLayerLayout->getLayerRotationXY(iplane) ); // alfa
    _GEAR_.SetY(   _siPlanesLayerLayout->getLayerRotationZX(iplane) ); // beta
    _GEAR_.SetZ(   _siPlanesLayerLayout->getLayerRotationZY(iplane) ); // gamma

    // code taken from EUTelAPIXHistograms
    // small fixes to make it compile,
    // definitely not the way to write code
    // 
    double _zDUT           = 0.;
    double _zDUTneighbour  = 0.;
    int    _indexDUT       = 0;
    int    _manualDUTid    = 0;

    std::vector<EVENT::TrackerHit*>  trackhits = tr->getTrackerHits();
    int nHit =   trackhits.size();

    // first, get the track impact points at the DUT and the layer closest to it

    // coordinates of the track impact point at the closest to DUT layer
    double x1 = 0, y1 = 0, z1 = 0;
    // coordinates of the track impact point at the DUT layer
    double x2 = 0, y2 = 0, z2 = 0;

    // for sanity check
    bool foundHitDUT = false;
    bool foundHitNeighbour = false;

    double dist;
    for( int ihit = 0; ihit <  nHit; ihit++ ) {
      TrackerHit * meshit = trackhits.at(ihit);

      // Look at fitted hits only!
      if( meshit->getType() < 32 ) continue;
      // Hit position
      const double * pos = meshit->getPosition();

      // look for a hit at DUT
      dist =  pos[2] - _zDUT;
      if( dist * dist < 1 ) {
	if( foundHitDUT ) {
	  cout << "hit at DUT layer already found! Terminating" << endl;
	  abort();
	}
	x2 = pos[0];
	y2 = pos[1];
	z2 = pos[2];
	foundHitDUT = true;
      }
      // look for a hit at DUT's neighbour
      dist =  pos[2] - _zDUTneighbour;
      if( dist * dist < 1 ) {
	if( foundHitNeighbour ) {
	  cout << "hit at DUT neigbour layer already found! Terminating" << endl;
	  abort();
	}
	x1 = pos[0];
	y1 = pos[1];
	z1 = pos[2];
      }
    }

    if( !( foundHitDUT && foundHitNeighbour) ) {
      cout << "Was not possible to find hits at dut and next-to-dut layers. Terminating." << endl;
      abort();
    }

    /*cout << "track impact points:" << endl;
      cout << x1 << " , " << z1 << endl;
      cout << x2 << " , " << z2 << endl;*/

    // now proceed to calculation of the intersection point
    // of this track and the rotated DUT layer

    // for the formulas see logbook 20/01/2011
    // track parametrization: X = offsetX + slopeX * Z, Y = offsetY + slopeY * Z,

    double slopeX =( x2 - x1) /( z2-z1);
    double offsetX = x1 - z1 * slopeX;

    double slopeY =( y2 - y1) /( z2-z1);
    double offsetY = y1 - z1 * slopeY;

    (dynamic_cast<AIDA::IHistogram1D*>( _aidaHistoMap["TrackSlopeX"]))->fill(slopeX);
    (dynamic_cast<AIDA::IHistogram1D*>( _aidaHistoMap["TrackOffsetX"]))->fill(offsetX);
    (dynamic_cast<AIDA::IHistogram1D*>( _aidaHistoMap["TrackSlopeY"]))->fill(slopeY);
    (dynamic_cast<AIDA::IHistogram1D*>( _aidaHistoMap["TrackOffsetY"]))->fill(offsetY);

    // now determine the total shift in X due to alignment
    //double	deltaX = 0;

    // this is a normal vector to the plane
    TVectorD	NormalVector(3);
    NormalVector(0) = 0.;
    NormalVector(1) = 0.;
    NormalVector(2) = 1.;

    double z_sensor = _siPlanesLayerLayout->getSensitivePositionZ(_indexDUT)+ 0.5 * _siPlanesLayerLayout->getSensitiveThickness( _indexDUT );

    // this is a vector to the point in the plane
    TVectorD r0Vector(3);
    r0Vector(0) = 0;
    r0Vector(1) = 0;
    r0Vector(2)=z_sensor;

    // this is a vector, pointing to the( 0, 0, z_sensor )
    // it's needed to translate the coordinate system to
    // and to perform the rotation around that point, not the( 0, 0, 0)
    // as it's done in the alignment/applyAlignment steps.
    TVectorD auxVector(3);
    auxVector(0) = 0;
    auxVector(1) = 0;
    auxVector(2) = z_sensor;

    for( unsigned i = 0; i < _alignmentCollectionNames.size(); i++ ) {

      LCCollectionVec * alignmentCollectionVec = dynamic_cast < LCCollectionVec * >( ev->getCollection(_alignmentCollectionNames[i]));
      // next, find the alignment constant corresponding to the DUT
      EUTelAlignmentConstant * c=NULL;
      for( size_t iPos = 0; iPos < alignmentCollectionVec->size(); ++iPos ) {

	c = static_cast< EUTelAlignmentConstant * >( alignmentCollectionVec->getElementAt( iPos ) );
	if( c -> getSensorID() == _manualDUTid ) break;	// this means we found the alignment constant corresponding
	// to the DUT; the pointer to it is now stored in c and can
	// be further used
      }
      if( c == NULL ) {
	cout << "Was not possible to found alignment constant, terminating" << endl;
	abort();
      }

      //------------------------------------------------------------------------------------
      TMatrixD	RotationMatrix(3,3);
      RotationMatrix(0,0) = 1;
      RotationMatrix(0,1) = c -> getGamma();
      RotationMatrix(0,2) = c -> getBeta();
      RotationMatrix(1,0) =( -1) * c -> getGamma();
      RotationMatrix(1,1) = 1;
      RotationMatrix(1,2) = c -> getAlpha();
      RotationMatrix(2,0) =( -1) * c -> getBeta();
      RotationMatrix(2,1) =( -1) * c -> getAlpha();
      RotationMatrix(2,2) = 1;

      // transform the normal vector( only the rotation)
      NormalVector = RotationMatrix * NormalVector;

      // transform the vector to the plane point( rotation+translations)
      r0Vector = RotationMatrix * ( r0Vector - auxVector) + auxVector;
      r0Vector(0) -= c -> getXOffset();
      r0Vector(1) -= c -> getYOffset();
      r0Vector(2) -= c -> getZOffset();

      //------------------------------------------------------------------------------------

      //deltaX -= c->getXOffset();
    }

    //r0Vector.Print();
    //NormalVector.Print();

    // now have to solve the equation
    TVectorD	trackImpact(3);
    TMatrixD	equationMatrix(3,3);
    TVectorD	b(3);

    equationMatrix(0, 0) = NormalVector( 0);
    equationMatrix(0, 1) = NormalVector( 1);
    equationMatrix(0, 2) = NormalVector( 2);
    equationMatrix(1, 0) = 1;
    equationMatrix(1, 1) = 0;
    equationMatrix(1, 2) =( -1)*slopeX;
    equationMatrix(2, 0) = 0;
    equationMatrix(2, 1) = 1;
    equationMatrix(2, 2) =( -1)*slopeY;

    b(0) = r0Vector(0) * NormalVector( 0) + r0Vector(1) * NormalVector( 1) + r0Vector(2) * NormalVector( 2);
    b(1) = offsetX;
    b(2) = offsetY;

    trackImpact = equationMatrix.Invert() * b;

    /*
    // very very naive approach
    // finally calculate track impact point
    z =( z2 /tan((-1)*_beta) + deltaX - offsetX) /( slopeX + 1./tan((-1)*_beta));
    x = slopeX * z + offsetX;
    // in this simple approximation, no correction for y
    y = y2;*/

    /*cout << trackImpact(0) << " ** " << x << endl;
      cout << trackImpact(1) << " ** " << y << endl;
      cout << trackImpact(2) << " ** " << z << endl;*/

    x = trackImpact(0);
    y = trackImpact(1);
    z = trackImpact(2);

}


int EUTelTestFitter::guessSensorID( Track * track ) {

  int sensorID = -1;
  /*
  double minDistance =  numeric_limits< double >::max();
  double * hitPosition = const_cast<double * >( hit->getPosition());

  for( int iPlane = 0; iPlane < _siPlanesLayerLayout->getNLayers(); ++iPlane ) 
  {
      printf("iPlane %5d   hitPos:  %8.3f  siZpos: %8.3f \n", iPlane, hitPosition[2] , _siPlaneZPosition[ iPlane ] );
      double distance = std::abs( hitPosition[2] - _siPlaneZPosition[ iPlane ] );
      if( distance < minDistance ) 
      {
          minDistance = distance;
          sensorID = _siPlanesLayerLayout->getID( iPlane );
      }
  }
  if( _siPlanesParameters->getSiPlanesType() == _siPlanesParameters->TelescopeWithDUT ) 
  {
      double distance = std::abs( hitPosition[2] - _siPlanesLayerLayout->getDUTPositionZ() );
      if( distance < minDistance )
      {
        minDistance = distance;
        sensorID = _siPlanesLayerLayout->getDUTID();
      }
  }
  if( minDistance > 10  // mm  ) 
  {
    // advice the user that the guessing wasn't successful 
    streamlog_out( WARNING3 ) << "A hit was found " << minDistance << " mm far from the nearest plane\n"
      "Please check the consistency of the data with the GEAR file " << endl;
    throw SkipEventException(this);
  }
*/
  return sensorID;
}


int EUTelTestFitter::guessSensorID( double & x, double & y, double & z) 
{

  int sensorID = -1;
  double minDistance =  numeric_limits< double >::max();

  for( int iPlane = 0; iPlane < _siPlanesLayerLayout->getNLayers(); ++iPlane ) {
    int _id = _siPlanesLayerLayout->getID(iPlane);
    // printf("iPlane %5d   hitPos:  %8.3f  siZpos: %8.3f \n", iPlane, z , _siPlaneCenter[ _id  ][2] );
    double distance = std::abs( z - _siPlaneCenter[ _id ][2] );
    if( distance < minDistance ) {
      minDistance = distance;
      sensorID = _id;
    }
  }

  if( minDistance > 10 /* mm */ ) {
    // inform the user that the guessing wasn't successful 
    streamlog_out( WARNING3 ) << "A hit was found " << minDistance << " mm far from the nearest plane\n"
      "Please check the consistency of the data with the GEAR file " << endl;
    throw SkipEventException(this);
  }

  return sensorID;
}

#include "CMSpixelFiles.cc"
#include "CMSpixelClust.cc"

//#include "../GeneralBrokenLines/trunk/cpp/GblTrajectory.cpp"
//#include "../GeneralBrokenLines/trunk/cpp/GblPoint.cpp"
//#include "../GeneralBrokenLines/trunk/cpp/GblData.cpp"
//#include "../GeneralBrokenLines/trunk/cpp/MilleBinary.cpp"
//#include "../GeneralBrokenLines/trunk/cpp/VMatrix.cpp"
//#include "../GeneralBrokenLines/trunk/cpp/BorderedBandMatrix.cpp"

#endif