#ifndef __CINT__
  #include <iostream.h>
  #include "AliTRDtracker.h"
  #include "AliTRDcluster.h" 
  #include "AliTRDv1.h"
  #include "AliTRDgeometry.h"    
  #include "AliTRDparameter.h"    
  #include "alles.h"  
  #include "AliTRDmcTrack.h"
  #include "AliTRDtrack.h"
  
  #include "TFile.h"
  #include "TParticle.h"
  #include "TStopwatch.h"
#endif


void AliTRDseedsAnalysis() {

  const Int_t nPrimaries = (Int_t) 3*86030/4;
  //  const Int_t nPrimaries = 500;

  Bool_t page[6]; for(Int_t i=0; i<6; i++) page[i]=kTRUE;

  const Int_t nPtSteps = 30;
  const Float_t maxPt = 3.;
  const Float_t minPt = 0.;
  
  const Int_t nEtaSteps = 40;
  const Float_t maxEta = 1.;
  const Float_t minEta = -1.;
  
  // page[0]
  TH1F *hNcl=new TH1F("hNcl","Seeds No of Clusters",255,-9.5,500.5); 
  hNcl->SetFillColor(4);
  hNcl->SetXTitle("No of clusters"); 
  hNcl->SetYTitle("counts"); 
  TH1F *hNtb=new TH1F("hNtb","Seeds No of TB with clusters",160,-9.5,150.5); 
  hNtb->SetFillColor(4);
  hNtb->SetXTitle("No of timebins with clusters"); 
  hNtb->SetYTitle("counts"); 
  TH1F *hNep=new TH1F("hNep","Seeds end point",160, -9.5, 150.5); 
  hNep->SetFillColor(4);
  hNep->SetXTitle("outermost timebin with cluster"); 
  hNep->SetYTitle("counts"); 
  TH1F *hNmg=new TH1F("hNmg","Seeds max gap",160, -9.5, 150.5); 
  hNmg->SetFillColor(4);
  hNmg->SetXTitle("max gap (tb)"); 
  hNmg->SetYTitle("counts"); 

  // page[1]
  Float_t fMin = -5.5, fMax = 134.5;
  Int_t   iChan = 140; 
  TH2F *h2ep = new TH2F("h2ep","MC vs RT (end point)",iChan,fMin,fMax,iChan,fMin,fMax); 
  h2ep->SetMarkerColor(4);
  h2ep->SetMarkerSize(2);
  TH2F *h2ntb = new TH2F("h2ntb","MC vs RT (TB with Clusters)",iChan,fMin,fMax,iChan,fMin,fMax); 
  h2ntb->SetMarkerColor(4);
  h2ntb->SetMarkerSize(2);
  TH2F *h2mg = new TH2F("h2mg","MC vs RT (max gap)",iChan/2,fMin,fMax/2,iChan/2,fMin,fMax/2); 
  h2mg->SetMarkerColor(4);
  h2mg->SetMarkerSize(2);

  // page[2]
  TH1F *hPt_all=new TH1F("hPt_all","Seeds Pt",nPtSteps,minPt,maxPt);
  hPt_all->SetLineColor(4);
  hPt_all->SetXTitle("Pt (GeV/c)"); 
  hPt_all->SetYTitle("counts"); 

  TH1F *hPt_short=new TH1F("hPt_short","Short seeds Pt",nPtSteps,minPt,maxPt);
  hPt_short->SetLineColor(8);
  TH1F *hPt_long=new TH1F("hPt_long","Long seeds Pt",nPtSteps,minPt,maxPt);
  hPt_long->SetLineColor(2);
  
  TH1F *hrtPt_short=new TH1F("hrtPt_short","RT from short seeds",nPtSteps,minPt,maxPt);
  hrtPt_short->SetLineColor(8);
  TH1F *hrtPt_long=new TH1F("hrtPt_long","RT from long seeds",nPtSteps,minPt,maxPt);
  hrtPt_long->SetLineColor(2);

  TH1F *hEta_all=new TH1F("hEta_all","Seeds Eta",nEtaSteps,minEta,maxEta);
  hEta_all->SetLineColor(4);
  hEta_all->SetXTitle("Eta"); 
  hEta_all->SetYTitle("counts"); 
  TH1F *hEta_short=new TH1F("hEta_short","Short seeds Eta",nEtaSteps,minEta,maxEta);
  hEta_short->SetLineColor(8);
  TH1F *hEta_long=new TH1F("hEta_long","Long seeds Eta",nEtaSteps,minEta,maxEta);
  hEta_long->SetLineColor(2);
  
  TH1F *hrtEta_short=new TH1F("hrtEta_short","RT from short seeds",nEtaSteps,minEta,maxEta);
  hrtEta_short->SetLineColor(8);
  TH1F *hrtEta_long=new TH1F("hrtEta_long","RT from long seeds",nEtaSteps,minEta,maxEta);
  hrtPt_long->SetLineColor(2);

  // page[3]
  TH1F *hEff_long = new TH1F("hEff_long","Efficiency vs Pt",nPtSteps,minPt,maxPt); 
  hEff_long->SetLineColor(2); hEff_long->SetLineWidth(2);  
  hEff_long->SetXTitle("Pt, GeV/c"); 
  hEff_long->SetYTitle("Efficiency"); 
  
  TH1F *hEff_short = new TH1F("hEff_short","Efficiency short",nPtSteps,minPt,maxPt); 
  hEff_short->SetLineColor(8); hEff_short->SetLineWidth(2);  
  hEff_short->SetXTitle("Pt, GeV/c"); 
  hEff_short->SetYTitle("Efficiency"); 
  
  TH1F *hEff_long_eta = new TH1F("hEff_long_eta","Efficiency vs Eta",nEtaSteps,minEta,maxEta); 
  hEff_long_eta->SetLineColor(2); hEff_long_eta->SetLineWidth(2);  
  hEff_long_eta->SetXTitle("Pt, GeV/c"); 
  hEff_long_eta->SetYTitle("Efficiency"); 
  
  TH1F *hEff_short_eta = new TH1F("hEff_short_eta","Efficiency short",nEtaSteps,minEta,maxEta); 
  hEff_short_eta->SetLineColor(8); hEff_short_eta->SetLineWidth(2);  
  hEff_short_eta->SetXTitle("Pt, GeV/c"); 
  hEff_short_eta->SetYTitle("Efficiency"); 
  
  // page[4]
  TH1F *hY=new TH1F("hY","Y resolution",50,-0.5,0.5);hY->SetLineColor(4);
  hY->SetLineWidth(2);
  hY->SetXTitle("(cm)");
  TH1F *hZ=new TH1F("hZ","Z resolution",80,-4,4);hZ->SetLineColor(4);
  hZ->SetLineWidth(2);
  hZ->SetXTitle("(cm)");
  TH1F *hPhi=new TH1F("hPhi","PHI resolution",100,-20.,20.); 
  hPhi->SetFillColor(4);
  hPhi->SetXTitle("(mrad)");
  TH1F *hLambda=new TH1F("hLambda","Lambda resolution",100,-100,100);
  hLambda->SetFillColor(17);
  hLambda->SetXTitle("(mrad)");
  TH1F *hPt=new TH1F("hPt","Relative Pt resolution",30,-10.,10.);  
  
  // page[5]
  TH1F *hY_n=new TH1F("hY_n","Normalized Y resolution",50,-0.5,0.5); hY_n->SetLineColor(4);
  hY_n->SetLineWidth(2);
  hY_n->SetXTitle("  ");
  TH1F *hZ_n=new TH1F("hZ_n","Normalized Z resolution",50,-0.5,0.5); hZ_n->SetLineColor(2);
  hZ_n->SetLineWidth(2);
  hZ_n->SetXTitle("   ");
  TH1F *hLambda_n=new TH1F("hLambda_n","Normalized Lambda resolution",50,-0.5,0.5);
  hLambda_n->SetFillColor(17);
  TH1F *hPt_n=new TH1F("hPt_n","Normalized Pt resolution",50,-0.5,0.5); 
  

  Int_t nEvent = 0;

  // Load Seeds saved as MC tracks 
  TObjArray mctarray(2000);
  TFile *mctf=TFile::Open("AliTRDtrackableSeeds.root");
  if (!mctf->IsOpen()) {cerr<<"Can't open AliTRDtrackableSeeds.root !\n"; return;}
  TTree *mctracktree=(TTree*)mctf->Get("MCtracks");
  TBranch *mctbranch=mctracktree->GetBranch("MCtracks");
  Int_t const nMCtracks = (Int_t) mctracktree->GetEntries();
  cerr<<"Found "<<nMCtracks<<" entries in the MC tracks tree"<<endl;
  for (Int_t i=0; i<nMCtracks; i++) {
    AliTRDmcTrack *ioMCtrack=new AliTRDmcTrack();
    mctbranch->SetAddress(&ioMCtrack);
    mctracktree->GetEvent(i);
    mctarray.AddLast(ioMCtrack);
  }
  mctf->Close();                 

  TFile *tf=TFile::Open("AliTRDtracks.root");

  if (!tf->IsOpen()) {cerr<<"Can't open AliTRDtracks.root !\n"; return;}
  TObjArray rtarray(2000);

  char   tname[100];
  sprintf(tname,"TRDb_%d",nEvent);     
  TTree *tracktree=(TTree*)tf->Get(tname);

  TBranch *tbranch=tracktree->GetBranch("tracks");

  Int_t nRecTracks = (Int_t) tracktree->GetEntries();
  cerr<<"Found "<<nRecTracks<<" entries in the track tree"<<endl;

  for (Int_t i=0; i<nRecTracks; i++) {
    AliTRDtrack *iotrack=new AliTRDtrack();
    tbranch->SetAddress(&iotrack);
    tracktree->GetEvent(i);
    rtarray.AddLast(iotrack);
  }
  tf->Close();                 

  //  return;

  AliTRDmcTrack *mct;
  AliTRDtrack *rt;
  Int_t rtIndex[nMCtracks];
  for(Int_t j = 0; j < nMCtracks; j++) {
    mct = (AliTRDmcTrack*) mctarray.UncheckedAt(j);
    rtIndex[j] = 0;
    for (Int_t i=0; i<nRecTracks; i++) {
      rt = (AliTRDtrack*) rtarray.UncheckedAt(i);
      Int_t label = rt->GetSeedLabel();
      if(mct->GetTrackIndex() == label) rtIndex[j] = i;
    }
  } 

  // Load clusters

  TFile *geofile =TFile::Open("AliTRDclusters.root");   
  AliTRDtracker *Tracker = new AliTRDtracker(geofile);
  Tracker->SetEventNumber(nEvent);

  AliTRDgeometry *fGeom   = (AliTRDgeometry*) geofile->Get("TRDgeometry"); 
  AliTRDparameter *fPar   = (AliTRDparameter*) geofile->Get("TRDparameter"); 
  
  Char_t *alifile = "AliTRDclusters.root";
  TObjArray carray(2000);
  TObjArray *ClustersArray = &carray;
  Tracker->ReadClusters(ClustersArray,alifile);   


  // Connect the AliRoot file containing Geometry, Kine, Hits, and Digits
  alifile = "galice.root";
  TFile *gafl = (TFile*) gROOT->GetListOfFiles()->FindObject(alifile);
  if (!gafl) {
    cout << "Open the ALIROOT-file " << alifile << endl;
    gafl = new TFile(alifile);
  }
  else {
    cout << alifile << " is already open" << endl;
  }
  gAlice = (AliRun*) gafl->Get("gAlice");
  if (gAlice)
    cout << "AliRun object found on file" << endl;
  else
    gAlice = new AliRun("gAlice","Alice test program");

  AliTRDv1       *TRD = (AliTRDv1*) gAlice->GetDetector("TRD");    
  const Int_t nparticles = gAlice->GetEvent(nEvent);
  if (nparticles <= 0) return;

  TParticle *p;
  Bool_t electrons[300000];
  for(Int_t i = 0; i < 300000; i++) electrons[i] = kFALSE;

  Bool_t mark_electrons = kFALSE;
  if(mark_electrons) {
    printf("mark electrons\n");

    for(Int_t i = nPrimaries; i < nparticles; i++) {
      p = gAlice->Particle(i);
      if(p->GetMass() > 0.01) continue;
      if(p->GetMass() < 0.00001) continue;
      electrons[i] = kTRUE;
    }
  }

  AliTRDcluster *cl = 0;
  Int_t nw, label, index, ti[3];
  Int_t mct_tbwc, rt_tbwc, mct_max_gap, rt_max_gap, mct_sector, rt_sector; 
  Int_t mct_max_tb, rt_max_tb, mct_min_tb, rt_min_tb;

  Int_t det, plane, ltb, gtb, gap, max_gap, sector;
  Double_t Pt, Px, Py, Pz;
  Double_t mcPt, mcPx, mcPy, mcPz;
  Double_t x,y,z, mcX;
  Int_t rtClusters, rtCorrect;

  Int_t nToFind_long, nFound_long, nToFind_short, nFound_short;

  printf("\n");

  Double_t dxAmp = (Double_t) fGeom->CamHght();   // Amplification region
  Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region
  Double_t dx = (Double_t) fPar->GetTimeBinSize();   

  Int_t tbAmp = fPar->GetTimeBefore();
  Int_t maxAmp = 0;
  Int_t tbDrift = fPar->GetTimeMax();
  Int_t maxDrift = (Int_t) ((dxDrift+0.000001)/dx);

  tbDrift = TMath::Min(tbDrift,maxDrift);
  tbAmp = TMath::Min(tbAmp,maxAmp);

  const Int_t nPlanes = fGeom->Nplan();
  const Int_t tbpp = Tracker->GetTimeBinsPerPlane();
  const Int_t nTB = tbpp * nPlanes;
  Int_t mask[nTB][2];

  nToFind_long = 0; nFound_long = 0;
  nToFind_short = 0; nFound_short = 0;

  for(Int_t i=0; i < nMCtracks; i++) {
    mct = (AliTRDmcTrack *) mctarray.UncheckedAt(i);
    label = mct->GetTrackIndex();

    // Waveform of the MC track
    for(nw = 0; nw < nTB; nw++) { mask[nw][0] = -1; mask[nw][1] = -1; }
    Int_t mct_ncl = mct->GetNumberOfClusters();

    for(ltb = 0; ltb < kMAX_TB; ltb++) {
      for(plane = 0; plane < nPlanes; plane++) {
	for(Int_t n = 0; n < 2; n++) {
	  index = mct->GetClusterIndex(ltb,plane,n);
	  if(index < 0) continue; 
	  cl = (AliTRDcluster *) ClustersArray->UncheckedAt(index);

	  for(nw = 0; nw < 3; nw++) ti[nw] = cl->GetLabel(nw);

	  if((ti[0] != label) && (ti[1] != label) &&  (ti[2] != label)) {
	    printf("mc wrong track label: %d, %d, %d != %d \n",
		   ti[0], ti[1], ti[2], label);
	  } 
	  det=cl->GetDetector();
	  if(fGeom->GetPlane(det) != plane) 
	    printf("cluster plane = %d != %d expected plane\n", 
		   fGeom->GetPlane(det), plane);
	  if(cl->GetLocalTimeBin() != ltb) 
	    printf("cluster ltb = %d != %d expected ltb\n", 
		   cl->GetLocalTimeBin(), ltb);
	  gtb = Tracker->GetGlobalTimeBin(0,plane,ltb);
	  mask[gtb][n] = index;
	}
      }
    }
    
    for(plane = 0; plane < nPlanes; plane++) {
      for(ltb = tbDrift-1; ltb >= -tbAmp; ltb--) {
	gtb = Tracker->GetGlobalTimeBin(0,plane,ltb);
	if(mask[gtb][0] > -1) break;
      }
      if(ltb >= -tbAmp) break;
    }  
    if((plane == nPlanes) && (ltb == -tbAmp-1)) {
      // printf("warning: for track %d min tb is not found and set to %d!\n",
      //	 label, nTB-1);
      mct_min_tb = nTB-1;
      // for(Int_t tb = 0; tb<nTB; tb++) printf("gtb %d; cl index %d\n",tb,mask[tb]);
    }
    else {
      mct_min_tb = Tracker->GetGlobalTimeBin(0,plane,ltb);
    }

    for(plane = nPlanes-1 ; plane>=0; plane--) {
      for(ltb = -tbAmp; ltb < tbDrift; ltb++) {
	gtb = Tracker->GetGlobalTimeBin(0,plane,ltb);
	if(mask[gtb][0] > -1) break;
      }
      if(ltb < tbDrift) break;
    }  
    if((plane == -1) && (ltb == tbDrift)) {
      //      printf("warning: for track %d max tb is not found and set to 0!\n",label);
      //      for(Int_t tb = 0; tb<nTB; tb++) printf("gtb %d; cl index %d\n",tb,mask[tb]);
      mct_max_tb = 0;
      //      mcX = Tracker->GetX(0,0,tbDrift-1);
    }
    else {
      mct_max_tb = Tracker-> GetGlobalTimeBin(0,plane,ltb);
      //      mcX = Tracker->GetX(0,plane,ltb);
      cl = (AliTRDcluster *) ClustersArray->UncheckedAt(mask[gtb][0]);
      det=cl->GetDetector();
      sector = fGeom->GetSector(det);  
      mct_sector = sector;
    }

    mct_tbwc = 0;
    mct_max_gap = 0;
    gap = -1;
    
    for(nw = mct_min_tb; nw < mct_max_tb+1; nw++) {
      gap++;
      if(mask[nw][0] > -1) {
	mct_tbwc++;
	if(gap > mct_max_gap) mct_max_gap = gap;
	gap = 0;
      }
    }

    //  Waveform of the reconstructed track
    if(rtIndex[i] >= 0) rt = (AliTRDtrack *) rtarray.UncheckedAt(rtIndex[i]);

    for(nw = 0; nw < nTB; nw++) { mask[nw][0] = -1; mask[nw][1] = -1; }
    Int_t rt_ncl = rt->GetNumberOfClusters();

    for(Int_t n = 0; n < rt_ncl; n++) {
      index = rt->GetClusterIndex(n);
      cl = (AliTRDcluster *) ClustersArray->UncheckedAt(index);
      
      for(nw = 0; nw < 3; nw++) ti[nw] = cl->GetLabel(nw);
      
      if((ti[0] != label) && (ti[1] != label) &&  (ti[2] != label)) {
	// printf("rt wrong track label: %d, %d, %d != %d \n", ti[0], ti[1], ti[2], label);
	continue;
      } 
      
      det=cl->GetDetector();
      plane = fGeom->GetPlane(det); 
      ltb = cl->GetLocalTimeBin(); 
      gtb = Tracker->GetGlobalTimeBin(0,plane,ltb);
      mask[gtb][0] = index;
    }

    for(plane = 0; plane < nPlanes; plane++) {
      for(ltb = tbDrift-1; ltb >= -tbAmp; ltb--) {
	gtb = Tracker->GetGlobalTimeBin(0,plane,ltb);
	if(mask[gtb][0] > -1) break;
      }
      if(ltb >= -tbAmp) break;
    }  
    if((plane == nPlanes) && (ltb == -tbAmp-1)) {
      // printf("warning: for track %d min tb is not found and set to %d!\n",
      //	     label, nTB-1);
      rt_min_tb = nTB-1;
      // for(Int_t tb = 0; tb<nTB; tb++) printf("gtb %d; cl index %d\n",tb,mask[tb]);
    }
    else {
      rt_min_tb = Tracker->GetGlobalTimeBin(0,plane,ltb);
    }

    for(plane = nPlanes-1 ; plane>=0; plane--) {
      for(ltb = -tbAmp; ltb < tbDrift; ltb++) {
	gtb = Tracker->GetGlobalTimeBin(0,plane,ltb);
	if(mask[gtb][0] > -1) break;
      }
      if(ltb < tbDrift) break;
    }  
    if((plane == -1) && (ltb == tbDrift)) {
      // printf("warning: for track %d max tb is not found and set to 0!\n",label);
      //      for(Int_t tb = 0; tb<nTB; tb++) printf("gtb %d; cl index %d\n",tb,mask[tb]);
      rt_max_tb = 0;
      //      mcX = Tracker->GetX(0,0,tbDrift-1);
    }
    else {
      rt_max_tb = Tracker-> GetGlobalTimeBin(0,plane,ltb);
      //      mcX = Tracker->GetX(0,plane,ltb);
      cl = (AliTRDcluster *) ClustersArray->UncheckedAt(mask[gtb][0]);
      det=cl->GetDetector();
      sector = fGeom->GetSector(det);  
      rt_sector = sector;
    }

    rt_tbwc = 0;
    rt_max_gap = 0;
    gap = -1;
    
    for(nw = rt_min_tb; nw < rt_max_tb+1; nw++) {
      gap++;
      if(mask[nw][0] > -1) {
	rt_tbwc++;
	if(gap > rt_max_gap) rt_max_gap = gap;
	gap = 0;
      }
    }

    // Fill the histoes

    if(page[0]) {
      hNcl->Fill((Float_t) mct_ncl);
      hNtb->Fill((Float_t) mct_tbwc);
      hNep->Fill((Float_t) mct_max_tb);
      hNmg->Fill((Float_t) mct_max_gap);
    }
    if(page[1]) {
      h2ep->Fill((Float_t) rt_max_tb, (Float_t) mct_max_tb);
      h2ntb->Fill((Float_t) rt_tbwc, (Float_t) mct_tbwc);
      h2mg->Fill((Float_t) rt_max_gap, (Float_t) mct_max_gap);
    }
    if(page[2]) {
      p = gAlice->Particle(label);
      hPt_all->Fill(p->Pt());
      hEta_all->Fill(p->Eta());
      if(mct_max_tb > 60) {
	nToFind_long++;
	hPt_long->Fill(p->Pt());
	hEta_long->Fill(p->Eta());
	if(((mct_max_tb - rt_max_tb) < 10) && 
	   (((Float_t) rt_tbwc) / ((Float_t) mct_tbwc) > 0.7)) {
	  nFound_long++;
	  hrtPt_long->Fill(p->Pt());
	  hrtEta_long->Fill(p->Eta());
	}	  
      }
      if((mct_max_tb < 60) && (mct_max_tb > 10)) {
	nToFind_short++;
	hPt_short->Fill(p->Pt());
	hEta_short->Fill(p->Eta());
	if(((mct_max_tb - rt_max_tb) < 10) && 
	   (((Float_t) rt_tbwc) / ((Float_t) mct_tbwc) > 0.7)) {
	  nFound_short++;
	  hrtPt_short->Fill(p->Pt());
	  hrtEta_short->Fill(p->Eta());
	}	  
      }
    }
    if(page[4] && page[5]) {
      if((mct_tbwc > 50) && (rt_tbwc > 50)) {
	if(rt->GetSeedLabel() != label) {
	  printf("mct vs rt index mismatch: %d != %d \n",
		 label, rt->GetSeedLabel());
	  return;
	}
	Pt = TMath::Abs(rt->GetPt()); 
	Double_t cc = TMath::Abs(rt->GetSigmaC2()); 
	mct->GetPxPyPzXYZ(mcPx,mcPy,mcPz,x,y,z,-1);
	rt->GetPxPyPz(Px,Py,Pz);      

	printf("\n\ntrack %d \n", label);
	printf("rt Px, Py, Pz: %f, %f, %f \n", Px, Py, Pz); 
	printf("mc Px, Py, Pz: %f, %f, %f \n", mcPx, mcPy, mcPz); 
    
	mcPt = TMath::Sqrt(mcPx * mcPx + mcPy * mcPy);
	if(mcPt < 0.0001) mcPt = 0.0001;
    
	Float_t lamg=TMath::ATan2(mcPz,mcPt);
	Float_t lam =TMath::ATan2(Pz,Pt);
	if(TMath::Abs(mcPt) < 0.0001) printf("attempt to divide by mcPt = %f\n",mcPt);
	else hPt_n->Fill((0.3*0.4/100*(1/Pt - 1/mcPt))/TMath::Sqrt(cc));
    
	Float_t phig=TMath::ATan2(mcPy,mcPx);
	Float_t phi =TMath::ATan2(Py,Px);    
	//	if(!(rt->PropagateTo(x))) continue;
	//    if((mcPt > Pt_min) && (mcPt < Pt_max)) {
	hLambda->Fill((lam - lamg)*1000.);
	hLambda_n->Fill((lam - lamg)/TMath::Sqrt(rt->GetSigmaTgl2()));
	if(TMath::Abs(mcPt) < 0.0001) printf("attempt to divide by mcPt = %f\n",mcPt);
	else hPt->Fill((1/Pt - 1/mcPt)/(1/mcPt)*100.); 
	hPhi->Fill((phi - phig)*1000.);
	hY->Fill(rt->GetY() - y);
	hZ->Fill(rt->GetZ() - z);
	hY_n->Fill((rt->GetY() - y)/TMath::Sqrt(rt->GetSigmaY2()));
	hZ_n->Fill((rt->GetZ() - z)/TMath::Sqrt(rt->GetSigmaZ2()));
      }
    }  
  }

  // plot pictures

  if(page[0]) {
    TCanvas *c0=new TCanvas("c0","",0,0,700,850);
    gStyle->SetOptStat(111110);
    c0->Divide(2,2);
    c0->cd(1); gPad->SetLogy(); hNcl->Draw();
    c0->cd(2); hNtb->Draw();
    c0->cd(3); hNep->Draw();
    c0->cd(4); hNmg->Draw();
    c0->Print("c0.ps","ps"); 
  }
  if(page[1]) {
    TCanvas *c1=new TCanvas("c1","",0,0,700,850);
    gStyle->SetOptStat(0);
    c1->Divide(2,2);
    c1->cd(1); h2ep->Draw();
    c1->cd(2); h2ntb->Draw();
    c1->cd(3); h2mg->Draw();
    //    c1->cd(4); hNmg->Draw();
    c1->Print("c1.ps","ps"); 
  }
  if(page[2]) {
    TCanvas *c2=new TCanvas("c2","",0,0,700,850);
    gStyle->SetOptStat(0);
    c2->Divide(2,2);
    c2->cd(1); hPt_all->Draw(); hPt_long->Draw("same"); hPt_short->Draw("same"); 
    c2->cd(2); hEta_all->Draw(); hEta_long->Draw("same"); hEta_short->Draw("same"); 
    //    c2->cd(3); h2mg->Draw();
    //    c2->cd(4); hNmg->Draw();
    c2->Print("c2.ps","ps"); 
  }
  if(page[3]) {
    TCanvas *c3=new TCanvas("c3","",0,0,700,850);
    gStyle->SetOptStat(0);
    c3->Divide(1,2);
    c3->cd(1);
    hrtPt_long->Sumw2(); hPt_long->Sumw2(); 
    hEff_long->Divide(hrtPt_long,hPt_long,1,1.,"b");
    hEff_long->SetMaximum(1.4);
    hEff_long->SetYTitle("Matching Efficiency");
    hEff_long->SetXTitle("Pt (GeV/c)");
    hEff_long->Draw();          
    hrtPt_short->Sumw2(); hPt_short->Sumw2(); 
    hEff_short->Divide(hrtPt_short,hPt_short,1,1.,"b");
    hEff_short->SetMaximum(1.4);
    hEff_short->SetYTitle("Matching Efficiency");
    hEff_short->SetXTitle("Pt (GeV/c)");
    hEff_short->Draw("same");          

    c3->cd(2);
    hrtEta_long->Sumw2(); hEta_long->Sumw2(); 
    hEff_long_eta->Divide(hrtEta_long,hEta_long,1,1.,"b");
    hEff_long_eta->SetMaximum(1.4);
    hEff_long_eta->SetYTitle("Matching Efficiency");
    hEff_long_eta->SetXTitle("Eta");
    hEff_long_eta->Draw();          
    hrtEta_short->Sumw2(); hEta_short->Sumw2(); 
    hEff_short_eta->Divide(hrtEta_short,hEta_short,1,1.,"b");
    hEff_short_eta->SetMaximum(1.4);
    hEff_short_eta->SetYTitle("Matching Efficiency");
    hEff_short_eta->SetXTitle("Eta");
    hEff_short_eta->Draw("same");          
    c3->Print("c3.ps","ps"); 
  }
  if(page[4]) {
    TCanvas *c4=new TCanvas("c4","",0,0,700,850);
    gStyle->SetOptStat(111110);
    c4->Divide(2,3);
    c4->cd(1); hY->Draw();
    c4->cd(2); hZ->Draw();
    c4->cd(3); hPhi->Draw();
    c4->cd(4); hLambda->Draw();
    c4->cd(5); hPt->Draw();
    c4->Print("c4.ps","ps"); 
  }
  if(page[5]) {
    TCanvas *c5=new TCanvas("c5","",0,0,700,850);
    gStyle->SetOptStat(111110);
    c5->Divide(2,3);
    c5->cd(1); hY_n->Draw();
    c5->cd(2); hZ_n->Draw();
    c5->cd(3); hLambda_n->Draw();
    c5->cd(4); hPt_n->Draw();
    c5->Print("c5.ps","ps"); 
  }
  printf("Efficiency (long) = %d/%d = %f\n",nFound_long,nToFind_long,
	 ((Float_t) nFound_long / (Float_t) nToFind_long));
  printf("Efficiency (shrt) = %d/%d = %f\n",nFound_short,nToFind_short,
	 ((Float_t) nFound_short / (Float_t) nToFind_short));
}