New comparison macros (M.Ivanov)
authorhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 21 Jul 2004 14:28:41 +0000 (14:28 +0000)
committerhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 21 Jul 2004 14:28:41 +0000 (14:28 +0000)
STEER/AliESDComparisonMI.C [new file with mode: 0644]
STEER/AliESDComparisonMI.h [new file with mode: 0644]
STEER/AliGenInfo.C [new file with mode: 0644]
STEER/AliGenInfo.h [new file with mode: 0644]
STEER/AliITSclusterComparison.C [new file with mode: 0644]
STEER/AliITSclusterComparison.h [new file with mode: 0644]
STEER/AliTPCComparisonMI.C [new file with mode: 0644]
STEER/AliTPCComparisonMI.h [new file with mode: 0644]

diff --git a/STEER/AliESDComparisonMI.C b/STEER/AliESDComparisonMI.C
new file mode 100644 (file)
index 0000000..bc333e9
--- /dev/null
@@ -0,0 +1,722 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+
+///////////////////////////////////////////////////////////////////////////////
+//                                                                           //
+//  Time Projection Chamber                                                  //
+//  Comparison macro for ESD                                                 //
+//  responsible: 
+//  marian.ivanov@cern.ch                                                    //
+//
+//
+
+/*
+marian.ivanov@cern.ch
+Usage:
+
+
+.L $ALICE_ROOT/STEER/AliGenInfo.C+
+//be sure you created genTracks file before
+.L $ALICE_ROOT/STEER/AliESDComparisonMI.C+
+
+//
+ESDCmpTr *t2 = new ESDCmpTr("genTracks.root","cmpESDTracks.root","galice.root",-1,1,0);
+t2->Exec();
+
+//
+//some cuts definition
+TCut cprim("cprim","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<0.00005&&abs(MC.fVDist[2])<0.0005")
+//TCut cprim("cprim","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<0.5&&abs(MC.fVDist[2])<0.5")
+TCut citsin("citsin","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<3.9");
+
+TCut csec("csec","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)>0.5");
+
+
+TCut crec("crec","fReconstructed==1");
+TCut cteta1("cteta1","abs(MC.fTrackRef.Theta()/3.1415-0.5)<0.25");
+TCut cpos1("cpos1","abs(MC.fParticle.fVz/sqrt(MC.fParticle.fVx*MC.fParticle.fVx+MC.fParticle.fVy*MC.fParticle.fVy))<1");
+TCut csens("csens","abs(sqrt(fVDist[0]**2+fVDist[1]**2)-170)<50");
+TCut cmuon("cmuon","abs(MC.fParticle.fPdgCode==-13)");
+TCut cchi2("cchi2","fESDTrack.fITSchi2MIP[0]<7.&&fESDTrack.fITSchi2MIP[1]<5.&&fESDTrack.fITSchi2MIP[2]<7.&&fESDTrack.fITSchi2MIP[3]<7.5&&fESDTrack.fITSchi2MIP[4]<6.")
+
+AliESDComparisonDraw comp;  
+comp.SetIO(); 
+
+//
+//example
+comp.fTree->SetAlias("radius","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)");
+
+TH1F his("his","his",100,0,20);
+TH1F hpools("hpools","hpools",100,-7,7);
+TH1F hfake("hfake","hfake",1000,0,150);
+TProfile profp0("profp0","profp0",20,-0.4,0.9)
+
+comp.DrawLogXY("fTPCinP0[3]","fTPCDelta[4]/fTPCinP1[3]","fReconstructed==1&&abs(fPdg)==211"+cprim,"1",4,0.2,1.5,-0.06,0.06)
+comp.fRes->Draw();
+comp.fMean->Draw();  
+comp.Eff("fTPCinP0[3]","fRowsWithDigits>120&&abs(fPdg)==211"+cteta1+cpos1+cprim,"fTPCOn",20,0.2,1.5)
+comp.fRes->Draw();
+
+comp.fTree->Draw("fESDTrack.fITSsignal/fESDTrack.fTPCsignal","fITSOn&&fTPCOn&&fESDTrack.fITSFakeRatio==0") 
+
+TH1F his("his","his",100,0,20);
+TH1F hpools("hpools","hpools",100,-7,7);
+
+TH2F * hdedx0 = new TH2F("dEdx0","dEdx0",100, 0,2,200,0,250); hdedx0->SetMarkerColor(2);
+TH2F * hdedx1 = new TH2F("dEdx1","dEdx1",100, 0,2,200,0,250); hdedx1->SetMarkerColor(3);
+TH2F * hdedx2 = new TH2F("dEdx2","dEdx2",100, 0,2,200,0,250); hdedx2->SetMarkerColor(4);
+TH2F * hdedx3 = new TH2F("dEdx3","dEdx3",100, 0,2,200,0,250); hdedx3->SetMarkerColor(5);
+comp.fTree->Draw("fESDTrack.fITSsignal:MC.fParticle.P()>>dEdx0","fITSOn&&MC.fParticle.P()<2&&abs(fPdg)==211&&fITStrack.fN==6"+cprim) 
+comp.fTree->Draw("fESDTrack.fITSsignal:MC.fParticle.P()>>dEdx1","fITSOn&&MC.fParticle.P()<2&&abs(fPdg)==2212&&fITStrack.fN==6"+cprim) 
+comp.fTree->Draw("fESDTrack.fITSsignal:MC.fParticle.P()>>dEdx2","fITSOn&&MC.fParticle.P()<2&&abs(fPdg)==321&&fITStrack.fN==6"+cprim) 
+comp.fTree->Draw("fESDTrack.fITSsignal:MC.fParticle.P()>>dEdx3","fITSOn&&MC.fParticle.P()<2&&abs(fPdg)==11&&fITStrack.fN==6"+cprim) 
+hdedx0->Draw(); hdedx1->Draw("same"); hdedx2->Draw("same"); hdedx3->Draw("same");
+
+comp.DrawXY("fITSinP0[3]","fITSPools[4]","fReconstructed==1&&fPdg==-211&&fITSOn"+cprim,"1",4,0.2,1.0,-8,8)
+
+*/
+
+
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include <stdio.h>
+#include <string.h>
+//ROOT includes
+#include "Rtypes.h"
+#include "TFile.h"
+#include "TTree.h"
+#include "TChain.h"
+#include "TCut.h"
+#include "TString.h"
+#include "TBenchmark.h"
+#include "TStopwatch.h"
+#include "TParticle.h"
+#include "TSystem.h"
+#include "TTimer.h"
+#include "TVector3.h"
+#include "TPad.h"
+#include "TCanvas.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TF1.h"
+#include "TText.h"
+#include "Getline.h"
+#include "TStyle.h"
+
+//ALIROOT includes
+#include "AliRun.h"
+#include "AliStack.h"
+#include "AliESDtrack.h"
+#include "AliSimDigits.h"
+#include "AliTPCParam.h"
+#include "AliTPC.h"
+#include "AliTPCLoader.h"
+#include "AliDetector.h"
+#include "AliTrackReference.h"
+#include "AliRun.h"
+#include "AliTPCParamSR.h"
+#include "AliTracker.h"
+#include "AliComplexCluster.h"
+#include "AliMagF.h"
+#include "AliESD.h"
+#include "AliESDtrack.h"
+#include "AliITStrackV2.h"
+#include "AliTRDtrack.h"
+#endif
+#include "AliGenInfo.h"
+#include "AliESDComparisonMI.h"
+
+
+
+//
+//
+void AliESDRecInfo::Update(AliMCInfo* info,AliTPCParam * par, Bool_t reconstructed)
+{
+  //
+  //
+  //calculates derived variables
+  //  
+  //
+  AliTrackReference * ref = &(info->fTrackRef);
+  fTPCinR0[0] = info->fTrackRef.X();   
+  fTPCinR0[1] = info->fTrackRef.Y();   
+  fTPCinR0[2] = info->fTrackRef.Z();
+  fTPCinR0[3] = TMath::Sqrt(fTPCinR0[0]*fTPCinR0[0]+fTPCinR0[1]*fTPCinR0[1]);
+  fTPCinR0[4] = TMath::ATan2(fTPCinR0[1],fTPCinR0[0]);
+  //
+  fTPCinP0[0] = ref->Px();
+  fTPCinP0[1] = ref->Py();
+  fTPCinP0[2] = ref->Pz();
+  fTPCinP0[3] = ref->Pt();
+  fTPCinP0[4] = ref->P();
+  fDeltaP     = (ref->P()-info->fParticle.P())/info->fParticle.P();
+  //
+  //
+  if (fTPCinP0[3]>0.0000001){
+    //
+    fTPCAngle0[0] = TMath::ATan2(fTPCinP0[1],fTPCinP0[0]);
+    fTPCAngle0[1] = TMath::ATan(fTPCinP0[2]/fTPCinP0[3]);
+  }
+  //
+  //
+  fITSinP0[0]=info->fParticle.Px();
+  fITSinP0[1]=info->fParticle.Py();
+  fITSinP0[2]=info->fParticle.Pz();
+  fITSinP0[3]=info->fParticle.Pt();    
+  //
+  fITSinR0[0]=info->fParticle.Vx();
+  fITSinR0[1]=info->fParticle.Vy();
+  fITSinR0[2]=info->fParticle.Vz();
+  fITSinR0[3] = TMath::Sqrt(fITSinR0[0]*fITSinR0[0]+fITSinR0[1]*fITSinR0[1]);
+  fITSinR0[4] = TMath::ATan2(fITSinR0[1],fITSinR0[0]);
+  //
+  //
+  if (fITSinP0[3]>0.0000001){
+    fITSAngle0[0] = TMath::ATan2(fITSinP0[1],fITSinP0[0]);
+    fITSAngle0[1] = TMath::ATan(fITSinP0[2]/fITSinP0[3]);
+  }
+  //
+  if (reconstructed==kFALSE){
+    fReconstructed = kFALSE;
+    fTPCOn = kFALSE;
+    fITSOn = kFALSE;
+    fTRDOn = kFALSE;
+    return;
+  }
+  fReconstructed = kTRUE;
+  fTPCOn = ((fESDTrack.GetStatus()&AliESDtrack::kTPCrefit)>0) ? kTRUE : kFALSE;
+  fITSOn = ((fESDTrack.GetStatus()&AliESDtrack::kITSrefit)>0) ? kTRUE : kFALSE;
+  fTRDOn = ((fESDTrack.GetStatus()&AliESDtrack::kTRDrefit)>0) ? kTRUE : kFALSE;
+  
+  if (fTPCOn){
+    //TPC
+    fESDTrack.GetInnerXYZ(fTPCinR1);
+    fTPCinR1[3] = TMath::Sqrt(fTPCinR1[0]*fTPCinR1[0]+fTPCinR1[1]*fTPCinR1[1]);
+    fTPCinR1[4] = TMath::ATan2(fTPCinR1[1],fTPCinR1[0]);       
+    fESDTrack.GetInnerPxPyPz(fTPCinP1);
+    fTPCinP1[3] = TMath::Sqrt(fTPCinP1[0]*fTPCinP1[0]+fTPCinP1[1]*fTPCinP1[1]);
+    fTPCinP1[4] = TMath::Sqrt(fTPCinP1[3]*fTPCinP1[3]+fTPCinP1[2]*fTPCinP1[2]);
+    //
+    //
+    if (fTPCinP1[3]>0.0000001){
+      fTPCAngle1[0] = TMath::ATan2(fTPCinP1[1],fTPCinP1[0]);
+      fTPCAngle1[1] = TMath::ATan(fTPCinP1[2]/fTPCinP1[3]);  
+    }    
+    Double_t cov[15], param[5],x;
+    fESDTrack.GetInnerExternalCovariance(cov);
+    fESDTrack.GetInnerExternalParameters(x,param);    
+    //
+    fTPCDelta[0] = (fTPCinR0[4]-fTPCinR1[4])*fTPCinR1[3];  //delta rfi
+    fTPCPools[0] = fTPCDelta[0]/TMath::Sqrt(cov[0]);
+    fTPCDelta[1] = (fTPCinR0[2]-fTPCinR1[2]);              //delta z
+    fTPCPools[1] = fTPCDelta[1]/TMath::Sqrt(cov[2]);
+    fTPCDelta[2] = (fTPCAngle0[0]-fTPCAngle1[0]);
+    fTPCPools[2] = fTPCDelta[2]/TMath::Sqrt(cov[5]);
+    fTPCDelta[3] = (TMath::Tan(fTPCAngle0[1])-TMath::Tan(fTPCAngle1[1]));
+    fTPCPools[3] = fTPCDelta[3]/TMath::Sqrt(cov[9]);
+    fTPCDelta[4] = (fTPCinP0[3]-fTPCinP1[3]);
+    Double_t sign = (param[4]>0)? 1.:-1; 
+    fTPCPools[4] = sign*(1./fTPCinP0[3]-1./fTPCinP1[3])/TMath::Sqrt(cov[14]);    
+  }
+  if (fITSOn){
+    // ITS 
+    Double_t param[5],x;
+    //fESDTrack.GetExternalParameters(x,param);   
+    fESDTrack.GetConstrainedExternalParameters(x,param);   
+    Double_t cov[15];
+    fESDTrack.GetExternalCovariance(cov);
+    fESDTrack.GetConstrainedExternalCovariance(cov);
+    if (TMath::Abs(param[4])<0.0000000001) return;
+
+    fESDTrack.GetXYZ(fITSinR1);
+    fESDTrack.GetPxPyPz(fITSinP1);
+    fITSinP1[3] = TMath::Sqrt(fITSinP1[0]*fITSinP1[0]+fITSinP1[1]*fITSinP1[1]);
+    //
+    fITSinR1[3] = TMath::Sqrt(fITSinR1[0]*fITSinR1[0]+fITSinR1[1]*fITSinR1[1]);
+    fITSinR1[4] = TMath::ATan2(fITSinR1[1],fITSinR1[0]);
+    //
+    //
+    if (fITSinP1[3]>0.0000001){
+      fITSAngle1[0] = TMath::ATan2(fITSinP1[1],fITSinP1[0]);
+      fITSAngle1[1] = TMath::ATan(fITSinP1[2]/fITSinP1[3]);  
+    }
+    //
+    //
+    fITSDelta[0] = (fITSinR0[4]-fITSinR1[4])*fITSinR1[3];  //delta rfi
+    fITSPools[0] = fITSDelta[0]/TMath::Sqrt(cov[0]);
+    fITSDelta[1] = (fITSinR0[2]-fITSinR1[2]);              //delta z
+    fITSPools[1] = fITSDelta[1]/TMath::Sqrt(cov[2]);
+    fITSDelta[2] = (fITSAngle0[0]-fITSAngle1[0]);
+    fITSPools[2] = fITSDelta[2]/TMath::Sqrt(cov[5]);
+    fITSDelta[3] = (TMath::Tan(fITSAngle0[1])-TMath::Tan(fITSAngle1[1]));
+    fITSPools[3] = fITSDelta[3]/TMath::Sqrt(cov[9]);
+    fITSDelta[4] = (fITSinP0[3]-fITSinP1[3]);    
+    Double_t sign = (param[4]>0) ? 1:-1; 
+    fITSPools[4] = sign*(1./fITSinP0[3]-1./fITSinP1[3])/TMath::Sqrt(cov[14]);    
+  }
+  
+}
+
+
+void  AliESDRecV0Info::Update(Float_t vertex[3])
+{
+  Float_t v[3] = {fXr[0]-vertex[0],fXr[1]-vertex[1],fXr[2]-vertex[2]};
+  Float_t p[3] = {fPdr[0]+fPm[0], fPdr[1]+fPm[1],fPdr[2]+fPm[2]};
+
+  Float_t vnorm2 = v[0]*v[0]+v[1]*v[1];
+  Float_t vnorm3 = TMath::Sqrt(v[2]*v[2]+vnorm2);
+  vnorm2 = TMath::Sqrt(vnorm2);
+  Float_t pnorm2 = p[0]*p[0]+p[1]*p[1];
+  Float_t pnorm3 = TMath::Sqrt(p[2]*p[2]+pnorm2);
+  pnorm2 = TMath::Sqrt(pnorm2);
+  
+  fPointAngleFi = (v[0]*p[0]+v[1]*p[1])/(vnorm2*pnorm2);
+  fPointAngleTh = (v[2]*p[2]+vnorm2*pnorm2)/(vnorm3*pnorm3);  
+  fPointAngle   = (v[0]*p[0]+v[1]*p[1]+v[2]*p[2])/(vnorm3*pnorm3);
+
+}
+
+
+////////////////////////////////////////////////////////////////////////
+ESDCmpTr::ESDCmpTr()
+{
+  Reset();
+}
+
+////////////////////////////////////////////////////////////////////////
+ESDCmpTr::ESDCmpTr(const char* fnGenTracks,
+                  const char* fnCmp,
+                  const char* fnGalice, Int_t direction,
+                  Int_t nEvents, Int_t firstEvent)
+{
+  Reset();
+  //  fFnGenTracks = fnGenTracks;
+  //  fFnCmp = fnCmp;
+  sprintf(fFnGenTracks,"%s",fnGenTracks);
+  sprintf(fFnCmp,"%s",fnCmp);
+
+  fFirstEventNr = firstEvent;
+  fEventNr = firstEvent;
+  fNEvents = nEvents;
+  fDirection = direction;
+  //
+  fLoader = AliRunLoader::Open(fnGalice);
+  if (gAlice){
+    //delete gAlice->GetRunLoader();
+    delete gAlice;
+    gAlice = 0x0;
+  }
+  if (fLoader->LoadgAlice()){
+    cerr<<"Error occured while l"<<endl;
+  }
+  Int_t nall = fLoader->GetNumberOfEvents();
+  if (nEvents==0) {
+    nEvents =nall;
+    fNEvents=nall;
+    fFirstEventNr=0;
+  }    
+
+  if (nall<=0){
+    cerr<<"no events available"<<endl;
+    fEventNr = 0;
+    return;
+  }
+  if (firstEvent+nEvents>nall) {
+    fEventNr = nall-firstEvent;
+    cerr<<"restricted number of events availaible"<<endl;
+  }
+  AliMagF * magf = gAlice->Field();
+  AliTracker::SetFieldMap(magf);
+
+}
+
+
+////////////////////////////////////////////////////////////////////////
+ESDCmpTr::~ESDCmpTr()
+{
+  if (fLoader) {
+    delete fLoader;
+  }
+}
+
+//////////////////////////////////////////////////////////////
+Int_t ESDCmpTr::SetIO()
+{
+  //
+  // 
+  CreateTreeCmp();
+  if (!fTreeCmp) return 1;
+  fParamTPC = GetTPCParam();
+  //
+  if (!ConnectGenTree()) {
+    cerr<<"Cannot connect tree with generated tracks"<<endl;
+    return 1;
+  }
+  return 0;
+}
+
+//////////////////////////////////////////////////////////////
+
+Int_t ESDCmpTr::SetIO(Int_t eventNr)
+{
+  //
+  // 
+  // SET INPUT
+  //
+  TFile f("AliESDs.root");
+  //
+  TTree* tree = (TTree*) f.Get("esdTree");
+  if (!tree) { 
+    Char_t ename[100]; 
+    sprintf(ename,"%d",eventNr);
+    fEvent = (AliESD*)f.Get(ename);
+    if (!fEvent){
+      sprintf(ename,"ESD%d",eventNr);
+      fEvent = (AliESD*)f.Get(ename);
+    }
+  }
+  else{
+    tree->SetBranchAddress("ESD", &fEvent);
+    tree->GetEntry(eventNr);
+  }
+
+
+  /*
+  Char_t ename[100]; 
+  sprintf(ename,"%d",eventNr);
+  fEvent = (AliESD*)f.Get(ename);
+  if (!fEvent){
+    sprintf(ename,"ESD%d",eventNr);
+    fEvent = (AliESD*)f.Get(ename);
+  }
+  
+  TTree* tree = (TTree*) f.Get("esdTree");
+  if (!tree) {
+    Error("CheckESD", "no ESD tree found");
+    return kFALSE;
+  }
+  tree->SetBranchAddress("ESD", &fEvent);
+  tree->GetEntry(eventNr);
+  */
+
+  if (!fEvent) return 1;
+
+  return 0;
+}
+
+
+
+////////////////////////////////////////////////////////////////////////
+void ESDCmpTr::Reset()
+{
+  fEventNr = 0;
+  fNEvents = 0;
+  fTreeCmp = 0;
+  //  fFnCmp = "cmpTracks.root";
+  fFileGenTracks = 0;
+  fDebug = 0;
+  //
+  fParamTPC = 0;
+  fEvent =0;
+}
+
+////////////////////////////////////////////////////////////////////////
+Int_t ESDCmpTr::Exec(Int_t nEvents, Int_t firstEventNr)
+{
+  fNEvents = nEvents;
+  fFirstEventNr = firstEventNr;
+  return Exec();
+}
+
+////////////////////////////////////////////////////////////////////////
+Int_t ESDCmpTr::Exec()
+{
+  TStopwatch timer;
+  timer.Start();
+
+  if (SetIO()==1) 
+    return 1;
+   
+  fNextTreeGenEntryToRead = 0;
+  cerr<<"fFirstEventNr, fNEvents: "<<fFirstEventNr<<" "<<fNEvents<<endl;
+  for (Int_t eventNr = fFirstEventNr; eventNr < fFirstEventNr+fNEvents;
+       eventNr++) {
+    fEventNr = eventNr;
+    SetIO(fEventNr);
+    fNParticles = gAlice->GetEvent(fEventNr);    
+
+    fIndexRecTracks = new Int_t[fNParticles*4];  //write at maximum 4 tracks corresponding to particle
+    fFakeRecTracks = new Int_t[fNParticles];
+    fMultiRecTracks = new Int_t[fNParticles];
+    for (Int_t i = 0; i<fNParticles; i++) {
+      fIndexRecTracks[4*i] = -1;
+      fIndexRecTracks[4*i+1] = -1;
+      fIndexRecTracks[4*i+2] = -1;
+      fIndexRecTracks[4*i+3] = -1;
+
+      fFakeRecTracks[i] = 0;
+      fMultiRecTracks[i] = 0;
+    }
+  
+    cout<<"Start to process event "<<fEventNr<<endl;
+    cout<<"\tfNParticles = "<<fNParticles<<endl;
+    if (fDebug>2) cout<<"\tStart loop over TreeT"<<endl;
+    if (TreeTLoop()>0) return 1;
+
+    if (fDebug>2) cout<<"\tStart loop over tree genTracks"<<endl;
+    if (TreeGenLoop(eventNr)>0) return 1;
+    if (fDebug>2) cout<<"\tEnd loop over tree genTracks"<<endl;
+
+    delete [] fIndexRecTracks;
+    delete [] fFakeRecTracks;
+    delete [] fMultiRecTracks;
+  }
+
+  CloseOutputFile();
+
+  cerr<<"Exec finished"<<endl;
+  timer.Stop();
+  timer.Print();
+  return 0;
+
+}
+////////////////////////////////////////////////////////////////////////
+Bool_t ESDCmpTr::ConnectGenTree()
+{
+//
+// connect all branches from the genTracksTree
+// use the same variables as for the new cmp tree, it may work
+//
+  fFileGenTracks = TFile::Open(fFnGenTracks,"READ");
+  if (!fFileGenTracks) {
+    cerr<<"Error in ConnectGenTree: cannot open file "<<fFnGenTracks<<endl;
+    return kFALSE;
+  }
+  fTreeGenTracks = (TTree*)fFileGenTracks->Get("genTracksTree");
+  if (!fTreeGenTracks) {
+    cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file "
+       <<fFnGenTracks<<endl;
+    return kFALSE;
+  }
+  //
+  fMCInfo = new AliMCInfo;
+  fTreeGenTracks->SetBranchAddress("MC",&fMCInfo);
+  //
+  if (fDebug > 1) {
+    cout<<"Number of gen. tracks with TR: "<<fTreeGenTracks->GetEntries()<<endl;
+  }
+  return kTRUE;
+}
+
+
+////////////////////////////////////////////////////////////////////////
+void ESDCmpTr::CreateTreeCmp() 
+{
+  fFileCmp = TFile::Open(fFnCmp,"RECREATE");
+  if (!fFileCmp) {
+    cerr<<"Error in CreateTreeCmp: cannot open file "<<fFnCmp<<endl;
+    return;
+  }
+
+
+  fTreeCmp    = new TTree("ESDcmpTracks","ESDcmpTracks");
+  //
+  fMCInfo = new AliMCInfo;
+  fRecInfo = new AliESDRecInfo;
+  //
+  AliESDtrack * esdTrack = new AliESDtrack;
+  AliITStrackV2 * itsTrack = new AliITStrackV2;  
+  //
+  fTreeCmp->Branch("MC","AliMCInfo",&fMCInfo);
+  fTreeCmp->Branch("RC","AliESDRecInfo",&fRecInfo);
+  fTreeCmp->Branch("fESDTrack","AliESDtrack",&esdTrack);
+  //  fTreeCmp->Branch("ITS","AliITStrackV2",&itsTrack);
+  delete esdTrack;
+  //
+  fTreeCmp->AutoSave(); 
+
+}
+////////////////////////////////////////////////////////////////////////
+void ESDCmpTr::CloseOutputFile()  
+{
+  if (!fFileCmp) {
+    cerr<<"File "<<fFnCmp<<" not found as an open file."<<endl;
+    return;
+  }
+  fFileCmp->cd();
+  fTreeCmp->Write();    
+  delete fTreeCmp;
+  
+  fFileCmp->Close();
+  delete fFileCmp;
+  return;
+}
+////////////////////////////////////////////////////////////////////////
+
+TVector3 ESDCmpTr::TR2Local(AliTrackReference *trackRef,
+                           AliTPCParam *paramTPC) {
+
+  Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
+  Int_t index[4];
+  paramTPC->Transform0to1(x,index);
+  paramTPC->Transform1to2(x,index);
+  return TVector3(x);
+}
+////////////////////////////////////////////////////////////////////////
+
+Int_t ESDCmpTr::TreeTLoop()
+{
+  //
+  // loop over all ESD reconstructed tracks and store info in memory
+  //
+  TStopwatch  timer;
+  timer.Start();
+  //  
+  Int_t nEntries = (Int_t)fEvent->GetNumberOfTracks();  
+  //
+  fTracks      = new TObjArray(nEntries);
+  //
+  //load tracks to the memory
+  for (Int_t i=0; i<nEntries;i++){
+    AliESDtrack * track =fEvent->GetTrack(i);
+    fTracks->AddAt(track,i);
+  }
+  //
+  //
+  AliESDtrack * track=0;
+  for (Int_t iEntry=0; iEntry<nEntries;iEntry++){
+    track = (AliESDtrack*)fTracks->UncheckedAt(iEntry);
+    //
+    Int_t label = track->GetLabel();
+    Int_t absLabel = abs(label);
+    if (absLabel < fNParticles) {
+      //      fIndexRecTracks[absLabel] =  iEntry;
+      if (label < 0) fFakeRecTracks[absLabel]++;      
+      if (fMultiRecTracks[absLabel]>0){
+       if (fMultiRecTracks[absLabel]<4)
+         fIndexRecTracks[absLabel*4+fMultiRecTracks[absLabel]] =  iEntry;      
+      }
+      else      
+       fIndexRecTracks[absLabel*4] =  iEntry;
+      fMultiRecTracks[absLabel]++;
+    }
+  }  
+  printf("Time spended in TreeTLoop\n");
+  timer.Print();
+  
+  if (fDebug > 2) cerr<<"end of TreeTLoop"<<endl;  
+  return 0;
+}
+
+////////////////////////////////////////////////////////////////////////
+Int_t ESDCmpTr::TreeGenLoop(Int_t eventNr)
+{
+//
+// loop over all entries for a given event, find corresponding 
+// rec. track and store in the fTreeCmp
+//
+  TStopwatch timer;
+  timer.Start();
+  Int_t entry = fNextTreeGenEntryToRead;
+  Double_t nParticlesTR = fTreeGenTracks->GetEntriesFast();
+  cerr<<"fNParticles, nParticlesTR, fNextTreeGenEntryToRead: "<<fNParticles<<" "
+      <<nParticlesTR<<" "<<fNextTreeGenEntryToRead<<endl;
+  TBranch * branch = fTreeCmp->GetBranch("RC");
+  branch->SetAddress(&fRecInfo); // set all pointers
+
+  while (entry < nParticlesTR) {
+    fTreeGenTracks->GetEntry(entry);
+    entry++;
+    if (eventNr < fMCInfo->fEventNr) continue;
+    if (eventNr > fMCInfo->fEventNr) continue;;
+    //
+    fNextTreeGenEntryToRead = entry-1;
+    if (fDebug > 2 && fMCInfo->fLabel < 10) {
+      cerr<<"Fill track with a label "<<fMCInfo->fLabel<<endl;
+    }
+    //
+    fRecInfo->Reset();
+    AliESDtrack * track=0;
+    fRecInfo->fReconstructed =0;
+    TVector3 local = TR2Local(&(fMCInfo->fTrackRef),fParamTPC);
+    local.GetXYZ(fRecInfo->fTRLocalCoord);     
+    //
+    if (fIndexRecTracks[fMCInfo->fLabel*4] >= 0) {
+      track= (AliESDtrack*)fTracks->UncheckedAt(fIndexRecTracks[fMCInfo->fLabel*4]);
+      //
+      //
+      // find nearest track if multifound
+      if (fIndexRecTracks[fMCInfo->fLabel*4]+1){
+       Double_t xyz[3];
+       track->GetInnerXYZ(xyz);
+       Float_t dz = TMath::Abs(local.Z()-xyz[2]);
+       for (Int_t i=1;i<4;i++){
+         if (fIndexRecTracks[fMCInfo->fLabel*4+i]>=0){
+           AliESDtrack * track2 = (AliESDtrack*)fTracks->UncheckedAt(fIndexRecTracks[fMCInfo->fLabel*4+i]);
+           track2->GetInnerXYZ(xyz);
+           Float_t dz2 = TMath::Abs(local.Z()-xyz[2]);
+           if  (TMath::Abs(dz2)<dz)
+             track = track2;              
+         }
+       }
+      }        
+      //
+      fRecInfo->fESDTrack =*track; 
+      if (track->GetITStrack())
+       fRecInfo->fITStrack = *((AliITStrackV2*)track->GetITStrack());
+      if (track->GetTRDtrack()){
+       fRecInfo->fTRDtrack = *((AliTRDtrack*)track->GetTRDtrack());
+      }
+      else{
+       fRecInfo->fTRDtrack.SetdEdx(-1);
+      }
+      fRecInfo->fReconstructed = 1;
+      fRecInfo->fFake     = fFakeRecTracks[fMCInfo->fLabel];
+      fRecInfo->fMultiple = fMultiRecTracks[fMCInfo->fLabel];
+      //
+      fRecInfo->Update(fMCInfo,fParamTPC,kTRUE);          
+    }
+    else{
+      fRecInfo->Update(fMCInfo,fParamTPC,kFALSE);
+    }
+
+    fTreeCmp->Fill();
+  }
+  fTreeCmp->AutoSave();
+  fTracks->Delete();
+  printf("Time spended in TreeGenLoop\n");
+  timer.Print();
+  if (fDebug > 2) cerr<<"end of TreeGenLoop"<<endl;
+
+  return 0;
+}
+////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////
+
+void AliESDComparisonDraw::SetIO(const char *fname)
+{
+  //
+   TFile* file = TFile::Open(fname);
+   //
+   fTree = (TTree*) file->Get("ESDcmpTracks");
+   if (!fTree) {
+    printf("no track comparison tree found\n");
+    file->Close();
+    delete file;
+  }
+}
+
+
diff --git a/STEER/AliESDComparisonMI.h b/STEER/AliESDComparisonMI.h
new file mode 100644 (file)
index 0000000..10899a2
--- /dev/null
@@ -0,0 +1,165 @@
+
+class AliESDComparisonDraw: public AliComparisonDraw{
+public:
+  AliESDComparisonDraw(){fTree = 0;}
+  void SetIO(const char *fname = "cmpESDTracks.root");  
+  ClassDef(AliESDComparisonDraw,1)
+};
+ClassImp(AliESDComparisonDraw)
+
+
+/////////////////////////////////////////////////////////////////////////
+class AliESDRecInfo: public TObject {
+
+public:
+  AliESDRecInfo(){}
+  ~AliESDRecInfo(){}
+  //
+  void Update(AliMCInfo* info,AliTPCParam * par, Bool_t reconstructed);
+  Double_t fTPCinR0[5];   //generated position of the track at inner tpc - radius [3] and fi [4]
+  Double_t fTPCinR1[5];   //reconstructed postion of the track           - radius [3] and fi [
+  Double_t fTPCinP0[5];   //generated position of the track at inner tpc
+  Double_t fTPCinP1[5];   //reconstructed postion of the track
+  Double_t fTPCAngle0[2]; // generated angle 
+  Double_t fTPCAngle1[2]; //refconstructed angle 
+  Double_t fTPCDelta[5];  // deltas
+  Double_t fTPCPools[5];  // pools
+  Double_t fITSinR0[5];   //generated position of the track at inner tpc
+  Double_t fITSinR1[5];   //reconstructed postion of the track
+  Double_t fITSinP0[5];   //generated position of the track at inner tpc
+  Double_t fITSinP1[5];   //reconstructed postion of the track
+  Double_t fITSAngle0[2]; // generated angle 
+  Double_t fITSAngle1[2]; //refconstructed angle
+  Double_t fITSDelta[5];  // deltas
+  Double_t fITSPools[5];  // pools
+  AliESDtrack fESDTrack;          // tpc track
+  AliITStrackV2 fITStrack;        //its track
+  AliTRDtrack fTRDtrack;        //its track
+  Float_t fTRLocalCoord[3];       //local coordinates of the track ref.
+  Int_t   fReconstructed;         //flag if track was reconstructed
+  Int_t fFake;             // fake track
+  Int_t fMultiple;         // number of reconstructions
+  Bool_t fTPCOn;           // TPC refitted inward
+  Bool_t fITSOn;           // ITS refitted inward
+  Bool_t fTRDOn;           // ITS refitted inward
+  Float_t fDeltaP;          //delta of momenta
+  void Reset();
+  ClassDef(AliESDRecInfo,1)  // container for 
+};
+ClassImp(AliESDRecInfo)
+
+void AliESDRecInfo::Reset()
+{
+  fMultiple =0; 
+  fFake     =0;
+  fReconstructed=0;
+}
+
+
+/////////////////////////////////////////////////////////
+/////////////////////////////////////////////////////////
+/////////////////////////////////////////////////////////
+
+
+class AliESDRecV0Info: public TObject {
+public:
+  void Update(Float_t vertex[3]);
+  AliESDRecInfo  fT1;      //track1
+  AliESDRecInfo  fT2;      //track2  
+  Double_t       fDist1;    //info about closest distance according closest MC - linear DCA
+  Double_t       fDist2;    //info about closest distance parabolic DCA
+  Double_t       fInvMass;  //reconstructed invariant mass -
+  //
+  Double_t       fPdr[3];    //momentum at vertex daughter  - according approx at DCA
+  Double_t       fXr[3];     //rec. position according helix
+  //
+  Double_t       fPm[3];    //momentum at the vertex mother
+  Double_t       fAngle[3]; //three angles
+  Double_t       fRr;       // rec position of the vertex 
+  Int_t          fLab[2];   //MC label of the partecle
+  Float_t        fPointAngleFi; //point angle fi
+  Float_t        fPointAngleTh; //point angle theta
+  Float_t        fPointAngle;   //point angle full
+  ClassDef(AliESDRecV0Info,1)   // container for  
+};
+
+ClassImp(AliESDRecV0Info)
+
+
+
+////////////////////////////////////////////////////////////////////////
+// 
+// Start of implementation of the class ESDCmpTr
+//
+////////////////////////////////////////////////////////////////////////
+
+class ESDCmpTr {
+
+public:
+  ESDCmpTr();
+  ESDCmpTr(const char* fnGenTracks,
+          const char* fnCmpRes      ="cmpTracks.root", 
+          const char* fnGalice      ="galice.root", Int_t direction=0,
+          Int_t nEvents=1, Int_t firstEvent=0);
+  virtual ~ESDCmpTr();
+  void Reset();
+  Int_t Exec();
+  Int_t Exec(Int_t nEvents, Int_t firstEventNr);
+  Int_t SetIO();
+  Int_t SetIO(Int_t eventNr );
+  void CreateTreeCmp();
+  void CloseOutputFile();
+  Bool_t ConnectGenTree();
+  Int_t TreeGenLoop(Int_t eventNr);
+  Int_t TreeTLoop();
+  void SetFirstEventNr(Int_t i) {fFirstEventNr = i;}
+  void SetNEvents(Int_t i) {fNEvents = i;}
+  void SetDebug(Int_t level) {fDebug = level;}
+
+// tmp method, should go to TrackReferenceESD
+  static TVector3 TR2Local(AliTrackReference *trackRef,
+                   AliTPCParam *paramTPC);
+
+private:
+
+  Int_t fEventNr;                 //! current event number
+  Int_t fNEvents;                 //! number of events to process
+  Int_t fFirstEventNr;            //! first event to process
+  //
+  char  fFnCmp[1000];                   //! output file name with cmp tracks
+  TFile *fFileCmp;                //! output file with cmp tracks
+  TTree *fTreeCmp;                //! output tree with cmp tracks
+  //
+  char  fFnGenTracks[1000];             //! input file name with gen tracks
+  TFile *fFileGenTracks;
+  TTree *fTreeGenTracks;
+  //
+  //
+  Int_t  fDirection;
+  //
+  AliRunLoader * fLoader;         //! pointer to the run loader
+  //TTree *fTreeRecTracks;          //! tree with reconstructed tracks
+  //
+  Int_t *fIndexRecTracks;         //! index of particle label in the TreeT_ESD
+  Int_t *fFakeRecTracks;          //! number of fake tracks
+  Int_t *fMultiRecTracks;         //! number of multiple reconstructions
+  //
+  TObjArray * fTracks;            //!container with tracks 
+  AliESD *fEvent;                 //!event
+
+  //
+  AliTPCParam* fParamTPC;         //! AliTPCParam
+  Int_t fNParticles;              //! number of particles in the input tree genTracks
+  Int_t fDebug;                   //! debug flag  
+  Int_t fNextTreeGenEntryToRead;    //! last entry already read from genTracks tree
+  //
+  AliMCInfo*  fMCInfo;           //! MC information writen per particle
+  AliESDRecInfo*  fRecInfo;          //! Rec. information writen per particle
+  //
+
+  ClassDef(ESDCmpTr,1)    // class which creates and fills tree with ESDGenTrack objects
+};
+ClassImp(ESDCmpTr)
+
+
+
diff --git a/STEER/AliGenInfo.C b/STEER/AliGenInfo.C
new file mode 100644 (file)
index 0000000..ab21613
--- /dev/null
@@ -0,0 +1,1044 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+
+///////////////////////////////////////////////////////////////////////////
+/*
+
+Origin: marian.ivanov@cern.ch
+
+Macro to generate comples MC information - used for Comparison later on
+How to use it?
+
+.L $ALICE_ROOT/STEER/AliGenInfo.C+
+AliGenInfoMaker *t = new AliGenInfoMaker("galice.root","genTracks.root",1,0)
+t->Exec();
+
+*/
+
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include <stdio.h>
+#include <string.h>
+//ROOT includes
+#include "Rtypes.h"
+#include "TFile.h"
+#include "TTree.h"
+#include "TChain.h"
+#include "TCut.h"
+#include "TString.h"
+#include "TBenchmark.h"
+#include "TStopwatch.h"
+#include "TParticle.h"
+#include "TSystem.h"
+#include "TTimer.h"
+#include "TVector3.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TCanvas.h"
+#include "TPad.h"
+#include "TF1.h"
+
+//ALIROOT includes
+#include "AliRun.h"
+#include "AliStack.h"
+#include "AliSimDigits.h"
+#include "AliTPCParam.h"
+#include "AliTPC.h"
+#include "AliTPCLoader.h"
+#include "AliDetector.h"
+#include "AliTrackReference.h"
+#include "AliTPCParamSR.h"
+#include "AliTracker.h"
+#include "AliMagF.h"
+#endif
+#include "AliGenInfo.h" 
+//
+//
+
+AliTPCParam * GetTPCParam(){
+  AliTPCParamSR * par = new AliTPCParamSR;
+  par->Update();
+  return par;
+}
+
+
+//_____________________________________________________________________________
+Float_t TPCBetheBloch(Float_t bg)
+{
+  //
+  // Bethe-Bloch energy loss formula
+  //
+  const Double_t kp1=0.76176e-1;
+  const Double_t kp2=10.632;
+  const Double_t kp3=0.13279e-4;
+  const Double_t kp4=1.8631;
+  const Double_t kp5=1.9479;
+
+  Double_t dbg = (Double_t) bg;
+
+  Double_t beta = dbg/TMath::Sqrt(1.+dbg*dbg);
+
+  Double_t aa = TMath::Power(beta,kp4);
+  Double_t bb = TMath::Power(1./dbg,kp5);
+
+  bb=TMath::Log(kp3+bb);
+  
+  return ((Float_t)((kp2-aa-bb)*kp1/aa));
+}
+
+
+
+
+
+////////////////////////////////////////////////////////////////////////
+AliMCInfo::AliMCInfo()
+{
+  fTPCReferences  = new TClonesArray("AliTrackReference",10);
+  fITSReferences  = new TClonesArray("AliTrackReference",10);
+  fTRDReferences  = new TClonesArray("AliTrackReference",10);
+  fTRdecay.SetTrack(-1);
+}
+
+AliMCInfo::~AliMCInfo()
+{
+  if (fTPCReferences) {
+    delete fTPCReferences;
+  }
+  if (fITSReferences){
+    delete fITSReferences;
+  }
+  if (fTRDReferences){    
+    delete fTRDReferences;  
+  }
+}
+
+
+
+void AliMCInfo::Update()
+{
+  //
+  //
+  fMCtracks =1;
+  if (!fTPCReferences) {
+    fNTPCRef =0;
+    return;
+  }
+  Float_t direction=1;
+  //Float_t rlast=0;
+  fNTPCRef = fTPCReferences->GetEntriesFast();
+  fNITSRef = fITSReferences->GetEntriesFast();
+
+  for (Int_t iref =0;iref<fTPCReferences->GetEntriesFast();iref++){
+    AliTrackReference * ref = (AliTrackReference *) fTPCReferences->At(iref);
+    //Float_t r = (ref->X()*ref->X()+ref->Y()*ref->Y());
+    Float_t newdirection = ref->X()*ref->Px()+ref->Y()*ref->Py(); //inside or outside
+    if (iref==0) direction = newdirection;
+    if ( newdirection*direction<0){
+      //changed direction
+      direction = newdirection;
+      fMCtracks+=1;
+    }
+    //rlast=r;                     
+  }
+  //
+  // decay info
+  fTPCdecay=kFALSE;
+  if (fTRdecay.GetTrack()>0){
+    fDecayCoord[0] = fTRdecay.X();
+    fDecayCoord[1] = fTRdecay.Y();
+    fDecayCoord[2] = fTRdecay.Z();
+    if ( (fTRdecay.R()<250)&&(fTRdecay.R()>85) && (TMath::Abs(fTRdecay.Z())<250) ){
+      fTPCdecay=kTRUE;     
+    }
+    else{
+      fDecayCoord[0] = 0;
+      fDecayCoord[1] = 0;
+      fDecayCoord[2] = 0;
+    }
+  }
+  //
+  //
+  //digits information update
+  fRowsWithDigits    = fTPCRow.RowsOn();    
+  fRowsWithDigitsInn = fTPCRow.RowsOn(63); // 63 = number of inner rows
+  fRowsTrackLength   = fTPCRow.Last() - fTPCRow.First();
+  //
+  //
+  // calculate primary ionization per cm  
+  if (fParticle.GetPDG()){
+    fMass = fParticle.GetMass();  
+    fCharge = fParticle.GetPDG()->Charge();
+    if (fTPCReferences->GetEntriesFast()>0){
+      fTrackRef = *((AliTrackReference*)fTPCReferences->At(0));
+    }
+    if (fMass>0){
+      Float_t p = TMath::Sqrt(fTrackRef.Px()*fTrackRef.Px()+
+                             fTrackRef.Py()*fTrackRef.Py()+
+                             fTrackRef.Pz()*fTrackRef.Pz());    
+      if (p>0.001){
+       Float_t betagama = p /fMass;
+       fPrim = TPCBetheBloch(betagama);
+      }else fPrim=0;
+    }
+  }else{
+    fMass =0;
+    fPrim =0;
+  }  
+}
+
+/////////////////////////////////////////////////////////////////////////////////
+void AliGenV0Info::Update()
+{
+  fMCPd[0] = fMCd.fParticle.Px();
+  fMCPd[1] = fMCd.fParticle.Py();
+  fMCPd[2] = fMCd.fParticle.Pz();
+  fMCPd[3] = fMCd.fParticle.P();
+  fMCX[0]  = fMCd.fParticle.Vx();
+  fMCX[1]  = fMCd.fParticle.Vy();
+  fMCX[2]  = fMCd.fParticle.Vz();
+  fMCR       = TMath::Sqrt( fMCX[0]*fMCX[0]+fMCX[1]*fMCX[1]);
+  fPdg[0]    = fMCd.fParticle.GetPdgCode();
+  fPdg[1]    = fMCm.fParticle.GetPdgCode();
+  //
+  fLab[0]    = fMCd.fParticle.GetUniqueID();
+  fLab[1]    = fMCm.fParticle.GetUniqueID();
+
+}
+
+
+
+
+  
+////////////////////////////////////////////////////////////////////////
+
+////////////////////////////////////////////////////////////////////////
+//
+// End of implementation of the class AliMCInfo
+//
+////////////////////////////////////////////////////////////////////////
+
+
+
+////////////////////////////////////////////////////////////////////////
+digitRow::digitRow()
+{
+  Reset();
+}
+////////////////////////////////////////////////////////////////////////
+digitRow & digitRow::operator=(const digitRow &digOld)
+{
+  for (Int_t i = 0; i<kgRowBytes; i++) fDig[i] = digOld.fDig[i];
+  return (*this);
+}
+////////////////////////////////////////////////////////////////////////
+void digitRow::SetRow(Int_t row) 
+{
+  if (row >= 8*kgRowBytes) {
+    cerr<<"digitRow::SetRow: index "<<row<<" out of bounds."<<endl;
+    return;
+  }
+  Int_t iC = row/8;
+  Int_t iB = row%8;
+  SETBIT(fDig[iC],iB);
+}
+
+////////////////////////////////////////////////////////////////////////
+Bool_t digitRow::TestRow(Int_t row)
+{
+//
+// return kTRUE if row is on
+//
+  Int_t iC = row/8;
+  Int_t iB = row%8;
+  return TESTBIT(fDig[iC],iB);
+}
+////////////////////////////////////////////////////////////////////////
+Int_t digitRow::RowsOn(Int_t upto)
+{
+//
+// returns number of rows with a digit  
+// count only rows less equal row number upto
+//
+  Int_t total = 0;
+  for (Int_t i = 0; i<kgRowBytes; i++) {
+    for (Int_t j = 0; j < 8; j++) {
+      if (i*8+j > upto) return total;
+      if (TESTBIT(fDig[i],j))  total++;
+    }
+  }
+  return total;
+}
+////////////////////////////////////////////////////////////////////////
+void digitRow::Reset()
+{
+//
+// resets all rows to zero
+//
+  for (Int_t i = 0; i<kgRowBytes; i++) {
+    fDig[i] <<= 8;
+  }
+}
+////////////////////////////////////////////////////////////////////////
+Int_t digitRow::Last()
+{
+//
+// returns the last row number with a digit
+// returns -1 if now digits 
+//
+  for (Int_t i = kgRowBytes-1; i>=0; i--) {
+    for (Int_t j = 7; j >= 0; j--) {
+      if TESTBIT(fDig[i],j) return i*8+j;
+    }
+  }
+  return -1;
+}
+////////////////////////////////////////////////////////////////////////
+Int_t digitRow::First()
+{
+//
+// returns the first row number with a digit
+// returns -1 if now digits 
+//
+  for (Int_t i = 0; i<kgRowBytes; i++) {
+    for (Int_t j = 0; j < 8; j++) {
+      if (TESTBIT(fDig[i],j)) return i*8+j;
+    }
+  }
+  return -1;
+}
+
+////////////////////////////////////////////////////////////////////////
+//
+// end of implementation of a class digitRow
+//
+////////////////////////////////////////////////////////////////////////
+  
+////////////////////////////////////////////////////////////////////////
+AliGenInfoMaker::AliGenInfoMaker()
+{
+  Reset();
+}
+
+////////////////////////////////////////////////////////////////////////
+AliGenInfoMaker::AliGenInfoMaker(const char * fnGalice, const char* fnRes,
+                                  Int_t nEvents, Int_t firstEvent)
+{
+  Reset();
+  fFirstEventNr = firstEvent;
+  fEventNr = firstEvent;
+  fNEvents = nEvents;
+  //   fFnRes = fnRes;
+  sprintf(fFnRes,"%s",fnRes);
+  //
+  fLoader = AliRunLoader::Open(fnGalice);
+  if (gAlice){
+    delete gAlice->GetRunLoader();
+    delete gAlice;
+    gAlice = 0x0;
+  }
+  if (fLoader->LoadgAlice()){
+    cerr<<"Error occured while l"<<endl;
+  }
+  Int_t nall = fLoader->GetNumberOfEvents();
+  if (nEvents==0) {
+    nEvents =nall;
+    fNEvents=nall;
+    fFirstEventNr=0;
+  }    
+
+  if (nall<=0){
+    cerr<<"no events available"<<endl;
+    fEventNr = 0;
+    return;
+  }
+  if (firstEvent+nEvents>nall) {
+    fEventNr = nall-firstEvent;
+    cerr<<"restricted number of events availaible"<<endl;
+  }
+  AliMagF * magf = gAlice->Field();
+  AliTracker::SetFieldMap(magf);
+}
+
+
+AliMCInfo * AliGenInfoMaker::MakeInfo(UInt_t i)
+{
+  // 
+  if (i<fNParticles) {
+    if (fGenInfo[i]) return  fGenInfo[i];
+    fGenInfo[i] = new AliMCInfo;  
+    fNInfos++;
+    return fGenInfo[i];
+  }
+  else 
+    return 0;  
+}
+
+////////////////////////////////////////////////////////////////////////
+void AliGenInfoMaker::Reset()
+{
+  fEventNr = 0;
+  fNEvents = 0;
+  fTreeGenTracks = 0;
+  fFileGenTracks = 0;
+  fGenInfo = 0;
+  fNInfos  = 0;
+  //
+  //
+  fDebug = 0;
+  fVPrim[0] = -1000.;
+  fVPrim[1] = -1000.;
+  fVPrim[2] = -1000.;
+  fParamTPC = 0;
+}
+////////////////////////////////////////////////////////////////////////
+AliGenInfoMaker::~AliGenInfoMaker()
+{
+  
+  if (fLoader){
+    fLoader->UnloadgAlice();
+    gAlice = 0;
+    delete fLoader;
+  }
+}
+
+Int_t  AliGenInfoMaker::SetIO()
+{
+  //
+  // 
+  CreateTreeGenTracks();
+  if (!fTreeGenTracks) return 1;
+  //  AliTracker::SetFieldFactor(); 
+  fParamTPC = GetTPCParam();
+  //
+  return 0;
+}
+
+////////////////////////////////////////////////////////////////////////
+Int_t AliGenInfoMaker::SetIO(Int_t eventNr)
+{
+  //
+  // 
+  // SET INPUT
+  fLoader->SetEventNumber(eventNr);
+  //
+  fLoader->LoadHeader();
+  fLoader->LoadKinematics();  
+  fStack = fLoader->Stack();
+  //
+  fLoader->LoadTrackRefs();
+  fTreeTR = fLoader->TreeTR();
+  //
+  AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
+  tpcl->LoadDigits();
+  fTreeD = tpcl->TreeD();  
+  return 0;
+}
+
+Int_t AliGenInfoMaker::CloseIOEvent()
+{
+  fLoader->UnloadHeader();
+  fLoader->UnloadKinematics();
+  fLoader->UnloadTrackRefs();
+  AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
+  tpcl->UnloadDigits();
+  return 0;
+}
+
+Int_t AliGenInfoMaker::CloseIO()
+{
+  fLoader->UnloadgAlice();
+  return 0;
+}
+
+
+
+////////////////////////////////////////////////////////////////////////
+Int_t AliGenInfoMaker::Exec(Int_t nEvents, Int_t firstEventNr)
+{
+  fNEvents = nEvents;
+  fFirstEventNr = firstEventNr;
+  return Exec();
+}
+
+////////////////////////////////////////////////////////////////////////
+Int_t AliGenInfoMaker::Exec()  
+{
+  TStopwatch timer;
+  timer.Start();
+  Int_t status =SetIO();
+  if (status>0) return status;
+  //
+
+  for (fEventNr = fFirstEventNr; fEventNr < fFirstEventNr+fNEvents;
+       fEventNr++) {
+    SetIO(fEventNr);
+    fNParticles = fStack->GetNtrack();
+    //
+    fGenInfo = new AliMCInfo*[fNParticles];
+    for (UInt_t i = 0; i<fNParticles; i++) {
+      fGenInfo[i]=0; 
+    }
+    //
+    cout<<"Start to process event "<<fEventNr<<endl;
+    cout<<"\tfNParticles = "<<fNParticles<<endl;
+    if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeTR"<<endl;
+    if (TreeTRLoop()>0) return 1;
+    //
+    if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeD"<<endl;
+    if (TreeDLoop()>0) return 1;
+    //
+    if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeK"<<endl;
+    if (TreeKLoop()>0) return 1;
+    if (fDebug>2) cout<<"\tEnd loop over TreeK"<<endl;
+    for (UInt_t i = 0; i<fNParticles; i++) {
+      if (fGenInfo[i]) delete fGenInfo[i]; 
+    }
+    delete []fGenInfo;
+    CloseIOEvent();
+  }
+  //
+  CloseIO();
+  CloseOutputFile();
+
+  cerr<<"Exec finished"<<endl;
+
+  timer.Stop();
+  timer.Print();
+  return 0;
+}
+////////////////////////////////////////////////////////////////////////
+void AliGenInfoMaker::CreateTreeGenTracks() 
+{
+  fFileGenTracks = TFile::Open(fFnRes,"RECREATE");
+  if (!fFileGenTracks) {
+    cerr<<"Error in CreateTreeGenTracks: cannot open file "<<fFnRes<<endl;
+    return;
+  }
+  fTreeGenTracks = new TTree("genTracksTree","genTracksTree");  
+  AliMCInfo * info = new AliMCInfo;
+  //
+  fTreeGenTracks->Branch("MC","AliMCInfo",&info);
+  delete info; 
+  fTreeGenTracks->AutoSave();
+}
+////////////////////////////////////////////////////////////////////////
+void AliGenInfoMaker::CloseOutputFile() 
+{
+  if (!fFileGenTracks) {
+    cerr<<"File "<<fFnRes<<" not found as an open file."<<endl;
+    return;
+  }
+  fFileGenTracks->cd();
+  fTreeGenTracks->Write();  
+  delete fTreeGenTracks;
+  fFileGenTracks->Close();
+  delete fFileGenTracks;
+  return;
+}
+
+////////////////////////////////////////////////////////////////////////
+Int_t AliGenInfoMaker::TreeKLoop()
+{
+//
+// open the file with treeK
+// loop over all entries there and save information about some tracks
+//
+
+  AliStack * stack = fStack;
+  if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
+  
+  if (fDebug > 0) {
+    cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
+       <<fEventNr<<endl;
+  }  
+  Int_t  ipdg = 0;
+  TParticlePDG *ppdg = 0;
+  // not all generators give primary vertex position. Take the vertex
+  // of the particle 0 as primary vertex.
+  TDatabasePDG  pdg; //get pdg table  
+  //thank you very much root for this
+  TBranch * br = fTreeGenTracks->GetBranch("MC");
+  TParticle *particle = stack->ParticleFromTreeK(0);
+  fVPrim[0] = particle->Vx();
+  fVPrim[1] = particle->Vy();
+  fVPrim[2] = particle->Vz();
+  for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
+    // load only particles with TR
+    AliMCInfo * info = GetInfo(iParticle);
+    if (!info) continue;
+    //////////////////////////////////////////////////////////////////////
+    info->fLabel = iParticle;
+    //
+    info->fParticle = *(stack->Particle(iParticle));
+    info->fVDist[0] = info->fParticle.Vx()-fVPrim[0];
+    info->fVDist[1] = info->fParticle.Vy()-fVPrim[1];
+    info->fVDist[2] = info->fParticle.Vz()-fVPrim[2]; 
+    info->fVDist[3] = TMath::Sqrt(info->fVDist[0]*info->fVDist[0]+
+                                 info->fVDist[1]*info->fVDist[1]+info->fVDist[2]*info->fVDist[2]);
+    //
+    //
+    ipdg = info->fParticle.GetPdgCode();
+    info->fPdg = ipdg;
+    ppdg = pdg.GetParticle(ipdg);         
+    info->fEventNr = fEventNr;
+    info->Update();
+    //////////////////////////////////////////////////////////////////////    
+    br->SetAddress(&info);    
+    fTreeGenTracks->Fill();
+  }
+  fTreeGenTracks->AutoSave();
+  if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl;
+  return 0;
+}
+
+
+
+
+////////////////////////////////////////////////////////////////////////
+Int_t AliGenInfoMaker::TreeDLoop()
+{
+  //
+  // open the file with treeD
+  // loop over all entries there and save information about some tracks
+  //
+  
+  Int_t nInnerSector = fParamTPC->GetNInnerSector();
+  Int_t rowShift = 0;
+  Int_t zero=fParamTPC->GetZeroSup()+6;  
+  //
+  //
+  AliSimDigits digitsAddress, *digits=&digitsAddress;
+  fTreeD->GetBranch("Segment")->SetAddress(&digits);
+  
+  Int_t sectorsByRows=(Int_t)fTreeD->GetEntries();
+  if (fDebug > 1) cout<<"\tsectorsByRows = "<<sectorsByRows<<endl;
+  for (Int_t i=0; i<sectorsByRows; i++) {
+    if (!fTreeD->GetEvent(i)) continue;
+    Int_t sec,row;
+    fParamTPC->AdjustSectorRow(digits->GetID(),sec,row);
+    if (fDebug > 1) cout<<sec<<' '<<row<<"                          \r";
+    // here I expect that upper sectors follow lower sectors
+    if (sec > nInnerSector) rowShift = fParamTPC->GetNRowLow();
+    //
+    digits->ExpandTrackBuffer();
+    digits->First();        
+    do {
+      Int_t iRow=digits->CurrentRow();
+      Int_t iColumn=digits->CurrentColumn();
+      Short_t digitValue = digits->CurrentDigit();
+      if (digitValue >= zero) {
+       Int_t label;
+       for (Int_t j = 0; j<3; j++) {
+         //      label = digits->GetTrackID(iRow,iColumn,j);
+         label = digits->GetTrackIDFast(iRow,iColumn,j)-2;       
+         if (label >= (Int_t)fNParticles) { //don't label from bakground event
+           continue;
+         }
+         if (label >= 0 && label <= (Int_t)fNParticles) {
+           if (fDebug > 6 ) {
+             cout<<"Inner loop: sector, iRow, iColumn, label, value, row "
+                 <<sec<<" "
+                 <<iRow<<" "<<iColumn<<" "<<label<<" "<<digitValue
+                 <<" "<<row<<endl;
+           }   
+           AliMCInfo * info = GetInfo(label);
+           if (info){
+             info->fTPCRow.SetRow(row+rowShift);
+           }
+         }
+       }
+      }
+    } while (digits->Next());
+  }
+  
+  if (fDebug > 2) cerr<<"end of TreeDLoop"<<endl;  
+  return 0;
+}
+
+
+////////////////////////////////////////////////////////////////////////
+Int_t AliGenInfoMaker::TreeTRLoop()
+{
+  //
+  // loop over TrackReferences and store the first one for each track
+  //  
+  TTree * treeTR = fTreeTR;
+  Int_t nPrimaries = (Int_t) treeTR->GetEntries();
+  if (fDebug > 1) cout<<"There are "<<nPrimaries<<" entries in TreeTR"<<endl;
+  //
+  //
+  //track references for TPC
+  TClonesArray* TPCArrayTR = new TClonesArray("AliTrackReference");
+  TClonesArray* ITSArrayTR = new TClonesArray("AliTrackReference");
+  TClonesArray* TRDArrayTR = new TClonesArray("AliTrackReference");
+  TClonesArray* RunArrayTR = new TClonesArray("AliTrackReference");
+  //
+  if (treeTR->GetBranch("TPC"))    treeTR->GetBranch("TPC")->SetAddress(&TPCArrayTR);
+  if (treeTR->GetBranch("ITS"))    treeTR->GetBranch("ITS")->SetAddress(&ITSArrayTR);
+  if (treeTR->GetBranch("TRD"))    treeTR->GetBranch("TRD")->SetAddress(&TRDArrayTR);
+  if (treeTR->GetBranch("AliRun")) treeTR->GetBranch("AliRun")->SetAddress(&RunArrayTR);
+  //
+  //
+  //
+  for (Int_t iPrimPart = 0; iPrimPart<nPrimaries; iPrimPart++) {
+    treeTR->GetEntry(iPrimPart);
+    //
+    // Loop over TPC references
+    //
+    for (Int_t iTrackRef = 0; iTrackRef < TPCArrayTR->GetEntriesFast(); iTrackRef++) {
+      AliTrackReference *trackRef = (AliTrackReference*)TPCArrayTR->At(iTrackRef);            
+      //
+      if (trackRef->Pt()<fgTPCPtCut) continue;
+      Int_t label = trackRef->GetTrack();      
+      AliMCInfo * info = GetInfo(label);
+      if (!info) info = MakeInfo(label);
+      if (!info) continue;
+      TClonesArray & arr = *(info->fTPCReferences);
+      new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);     
+    }
+    //
+    // Loop over ITS references
+    //
+    for (Int_t iTrackRef = 0; iTrackRef < ITSArrayTR->GetEntriesFast(); iTrackRef++) {
+      AliTrackReference *trackRef = (AliTrackReference*)ITSArrayTR->At(iTrackRef);            
+      //
+      if (trackRef->Pt()<fgTPCPtCut) continue;
+      Int_t label = trackRef->GetTrack();      
+      AliMCInfo * info = GetInfo(label);
+      if ( (!info) && trackRef->Pt()<fgITSPtCut) continue;
+      if (!info) info = MakeInfo(label);
+      if (!info) continue;
+      TClonesArray & arr = *(info->fITSReferences);
+      new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);     
+    }
+    //
+    // Loop over TRD references
+    //
+    for (Int_t iTrackRef = 0; iTrackRef < TRDArrayTR->GetEntriesFast(); iTrackRef++) {
+      AliTrackReference *trackRef = (AliTrackReference*)TRDArrayTR->At(iTrackRef);            
+      //
+      if (trackRef->Pt()<fgTPCPtCut) continue;
+      Int_t label = trackRef->GetTrack();      
+      AliMCInfo * info = GetInfo(label);
+      if ( (!info) && trackRef->Pt()<fgTRDPtCut) continue;
+      if (!info) info = MakeInfo(label);
+      if (!info) continue;
+      TClonesArray & arr = *(info->fTRDReferences);
+      new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);     
+    }
+    //
+    // get dacay position
+    for (Int_t iTrackRef = 0; iTrackRef < RunArrayTR->GetEntriesFast(); iTrackRef++) {
+      AliTrackReference *trackRef = (AliTrackReference*)RunArrayTR->At(iTrackRef);      
+      //
+      Int_t label = trackRef->GetTrack();
+      AliMCInfo * info = GetInfo(label);
+      if (!info) continue;
+      if (!trackRef->TestBit(BIT(2))) continue;  //if not decay
+      //      if (TMath::Abs(trackRef.X());
+      info->fTRdecay = *trackRef;      
+    }
+  }
+  //
+  TPCArrayTR->Delete();
+  delete  TPCArrayTR;
+  TRDArrayTR->Delete();
+  delete  TRDArrayTR;
+  ITSArrayTR->Delete();
+  delete  ITSArrayTR;
+  RunArrayTR->Delete();
+  delete  RunArrayTR;
+  //
+  return 0;
+}
+
+////////////////////////////////////////////////////////////////////////
+Float_t AliGenInfoMaker::TR2LocalX(AliTrackReference *trackRef,
+                                   AliTPCParam *paramTPC) {
+
+  Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
+  Int_t index[4];
+  paramTPC->Transform0to1(x,index);
+  paramTPC->Transform1to2(x,index);
+  return x[0];
+}
+////////////////////////////////////////////////////////////////////////
+
+
+
+TH1F * AliComparisonDraw::DrawXY(const char * chx, const char *chy, const char* selection, 
+               const char * quality, Int_t nbins, Float_t minx, Float_t maxx, Float_t miny, Float_t maxy)
+{
+  //
+  Double_t* bins = CreateLogBins(nbins, minx, maxx);
+  Int_t nBinsRes = 30;
+  TH2F* hRes2 = new TH2F("hRes2", "residuals", nbins, minx, maxx, nBinsRes, miny, maxy);
+  char cut[1000];
+  sprintf(cut,"%s&&%s",selection,quality);
+  char expression[1000];
+  sprintf(expression,"%s:%s>>hRes2",chy,chx);
+  fTree->Draw(expression, cut, "groff");
+  TH1F* hMean=0;
+  TH1F* hRes = CreateResHisto(hRes2, &hMean);
+  AliLabelAxes(hRes, chx, chy);
+  //
+  delete hRes2;
+  delete[] bins;
+  fRes  = hRes;
+  fMean = hMean;
+  return hRes;
+}
+
+TH1F * AliComparisonDraw::DrawLogXY(const char * chx, const char *chy, const char* selection, 
+               const char * quality, Int_t nbins, Float_t minx, Float_t maxx, Float_t miny, Float_t maxy)
+{
+  //
+  Double_t* bins = CreateLogBins(nbins, minx, maxx);
+  Int_t nBinsRes = 30;
+  TH2F* hRes2 = new TH2F("hRes2", "residuals", nbins, bins, nBinsRes, miny, maxy);
+  char cut[1000];
+  sprintf(cut,"%s&&%s",selection,quality);
+  char expression[1000];
+  sprintf(expression,"%s:%s>>hRes2",chy,chx);
+  fTree->Draw(expression, cut, "groff");
+  TH1F* hMean=0;  
+  TH1F* hRes = CreateResHisto(hRes2, &hMean);
+  AliLabelAxes(hRes, chx, chy);
+  //
+  delete hRes2;
+  delete[] bins;
+  fRes  = hRes;
+  fMean = hMean;
+  return hRes;
+}
+
+///////////////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////////////
+TH1F * AliComparisonDraw::Eff(const char *variable, const char* selection, const char * quality, 
+                             Int_t nbins, Float_t min, Float_t max)
+{
+  //
+  //
+  TH1F* hGen = new TH1F("hGen", "gen. tracks", nbins, min, max);
+  TH1F* hRec = new TH1F("hRec", "rec. tracks", nbins, min, max);
+  char inputGen[1000];  
+  sprintf(inputGen,"%s>>hGen", variable);
+  fTree->Draw(inputGen, selection, "groff");
+  char selectionRec[256];
+  sprintf(selectionRec, "%s && %s", selection, quality);
+  char inputRec[1000];  
+  sprintf(inputRec,"%s>>hRec", variable);
+  fTree->Draw(inputRec, selectionRec, "groff");
+  //
+  TH1F* hEff = CreateEffHisto(hGen, hRec);
+  AliLabelAxes(hEff, variable, "#epsilon [%]");
+  fRes = hEff;
+  delete hRec;
+  delete hGen;
+  return hEff;
+}
+
+
+
+///////////////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////////////
+TH1F * AliComparisonDraw::EffLog(const char *variable, const char* selection, const char * quality, 
+                             Int_t nbins, Float_t min, Float_t max)
+{
+  //
+  //
+  Double_t* bins = CreateLogBins(nbins, min, max);
+  TH1F* hGen = new TH1F("hGen", "gen. tracks", nbins, bins);
+  TH1F* hRec = new TH1F("hRec", "rec. tracks", nbins, bins);
+  char inputGen[1000];  
+  sprintf(inputGen,"%s>>hGen", variable);
+  fTree->Draw(inputGen, selection, "groff");
+  char selectionRec[256];
+  sprintf(selectionRec, "%s && %s", selection, quality);
+  char inputRec[1000];  
+  sprintf(inputRec,"%s>>hRec", variable);
+  fTree->Draw(inputRec, selectionRec, "groff");
+  //
+  TH1F* hEff = CreateEffHisto(hGen, hRec);
+  AliLabelAxes(hEff, variable, "#epsilon [%]");
+  fRes = hEff;
+  delete hRec;
+  delete hGen;
+  delete[] bins;
+  return hEff;
+}
+
+
+///////////////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////////////
+
+Double_t* AliComparisonDraw::CreateLogBins(Int_t nBins, Double_t xMin, Double_t xMax)
+{
+  Double_t* bins = new Double_t[nBins+1];
+  bins[0] = xMin;
+  Double_t factor = pow(xMax/xMin, 1./nBins);
+  for (Int_t i = 1; i <= nBins; i++)
+    bins[i] = factor * bins[i-1];
+  return bins;
+}
+
+
+
+
+void AliComparisonDraw::AliLabelAxes(TH1* histo, const char* xAxisTitle, const char* yAxisTitle)
+{
+  //
+  histo->GetXaxis()->SetTitle(xAxisTitle);
+  histo->GetYaxis()->SetTitle(yAxisTitle);
+}
+
+
+TH1F* AliComparisonDraw::CreateEffHisto(TH1F* hGen, TH1F* hRec)
+{
+  //
+  Int_t nBins = hGen->GetNbinsX();
+  TH1F* hEff = (TH1F*) hGen->Clone("hEff");
+  hEff->SetTitle("");
+  hEff->SetStats(kFALSE);
+  hEff->SetMinimum(0.);
+  hEff->SetMaximum(110.);
+  //
+  for (Int_t iBin = 0; iBin <= nBins; iBin++) {
+    Double_t nGen = hGen->GetBinContent(iBin);
+    Double_t nRec = hRec->GetBinContent(iBin);
+    if (nGen > 0) {
+      Double_t eff = nRec/nGen;
+      hEff->SetBinContent(iBin, 100. * eff);
+      Double_t error = sqrt((eff*(1.-eff)+0.01) / nGen);      
+      //      if (error == 0) error = sqrt(0.1/nGen);
+      //
+      if (error == 0) error = 0.0001;
+      hEff->SetBinError(iBin, 100. * error);
+    } else {
+      hEff->SetBinContent(iBin, 100. * 0.5);
+      hEff->SetBinError(iBin, 100. * 0.5);
+    }
+  }
+  return hEff;
+}
+
+
+
+TH1F* AliComparisonDraw::CreateResHisto(TH2F* hRes2, TH1F **phMean,  Bool_t drawBinFits, 
+                    Bool_t overflowBinFits)
+{
+  TVirtualPad* currentPad = gPad;
+  TAxis* axis = hRes2->GetXaxis();
+  Int_t nBins = axis->GetNbins();
+  TH1F* hRes, *hMean;
+  if (axis->GetXbins()->GetSize()){
+    hRes = new TH1F("hRes", "", nBins, axis->GetXbins()->GetArray());
+    hMean = new TH1F("hMean", "", nBins, axis->GetXbins()->GetArray());
+  }
+  else{
+    hRes = new TH1F("hRes", "", nBins, axis->GetXmin(), axis->GetXmax());
+    hMean = new TH1F("hMean", "", nBins, axis->GetXmin(), axis->GetXmax());
+
+  }
+  hRes->SetStats(false);
+  hRes->SetOption("E");
+  hRes->SetMinimum(0.);
+  //
+  hMean->SetStats(false);
+  hMean->SetOption("E");
+  // create the fit function
+  //TKFitGaus* fitFunc = new TKFitGaus("resFunc");
+  //   TF1 * fitFunc = new TF1("G","[3]+[0]*exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))",-3,3);
+  TF1 * fitFunc = new TF1("G","[0]*exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))",-3,3);
+  
+  fitFunc->SetLineWidth(2);
+  fitFunc->SetFillStyle(0);
+  // create canvas for fits
+  TCanvas* canBinFits = NULL;
+  Int_t nPads = (overflowBinFits) ? nBins+2 : nBins;
+  Int_t nx = Int_t(sqrt(nPads-1.));// + 1;
+  Int_t ny = (nPads-1) / nx + 1;
+  if (drawBinFits) {
+    canBinFits = (TCanvas*)gROOT->FindObject("canBinFits");
+    if (canBinFits) delete canBinFits;
+    canBinFits = new TCanvas("canBinFits", "fits of bins", 200, 100, 500, 700);
+    canBinFits->Divide(nx, ny);
+  }
+
+  // loop over x bins and fit projection
+  Int_t dBin = ((overflowBinFits) ? 1 : 0);
+  for (Int_t bin = 1-dBin; bin <= nBins+dBin; bin++) {
+    if (drawBinFits) canBinFits->cd(bin + dBin);
+    TH1D* hBin = hRes2->ProjectionY("hBin", bin, bin);
+    //    
+    if (hBin->GetEntries() > 5) {
+      fitFunc->SetParameters(hBin->GetMaximum(),hBin->GetMean(),hBin->GetRMS());
+      hBin->Fit(fitFunc,"s");
+      Double_t sigma = TMath::Abs(fitFunc->GetParameter(2));
+
+      if (sigma > 0.){
+       hRes->SetBinContent(bin, TMath::Abs(fitFunc->GetParameter(2)));
+       hMean->SetBinContent(bin, fitFunc->GetParameter(1));    
+      }
+      else{
+       hRes->SetBinContent(bin, 0.);
+       hMean->SetBinContent(bin,0);
+      }
+      hRes->SetBinError(bin, fitFunc->GetParError(2));
+      hMean->SetBinError(bin, fitFunc->GetParError(1));
+      
+      //
+      //
+
+    } else {
+      hRes->SetBinContent(bin, 0.);
+      hRes->SetBinError(bin, 0.);
+      hMean->SetBinContent(bin, 0.);
+      hMean->SetBinError(bin, 0.);
+    }
+    
+
+    if (drawBinFits) {
+      char name[256];
+      if (bin == 0) {
+       sprintf(name, "%s < %.4g", axis->GetTitle(), axis->GetBinUpEdge(bin));
+      } else if (bin == nBins+1) {
+       sprintf(name, "%.4g < %s", axis->GetBinLowEdge(bin), axis->GetTitle());
+      } else {
+       sprintf(name, "%.4g < %s < %.4g", axis->GetBinLowEdge(bin),
+               axis->GetTitle(), axis->GetBinUpEdge(bin));
+      }
+      canBinFits->cd(bin + dBin);
+      hBin->SetTitle(name);
+      hBin->SetStats(kTRUE);
+      hBin->DrawCopy("E");
+      canBinFits->Update();
+      canBinFits->Modified();
+      canBinFits->Update();
+    }
+    
+    delete hBin;
+  }
+
+  delete fitFunc;
+  currentPad->cd();
+  *phMean = hMean;
+  return hRes;
+}
+
+
+
diff --git a/STEER/AliGenInfo.h b/STEER/AliGenInfo.h
new file mode 100644 (file)
index 0000000..b80ee58
--- /dev/null
@@ -0,0 +1,207 @@
+
+
+
+////////////////////////////////////////////////////////////////////////
+//
+// Start of implementation of the class digitRow
+//
+////////////////////////////////////////////////////////////////////////
+const Int_t kgRowBytes = 32;
+
+class digitRow: public TObject {
+
+public:
+  digitRow();
+  virtual ~digitRow(){;}
+  void SetRow(Int_t row);
+  Bool_t TestRow(Int_t row);
+  digitRow & operator=(const digitRow &digOld);
+  Int_t RowsOn(Int_t upto=8*kgRowBytes);
+  Int_t Last();
+  Int_t First();
+  void Reset();
+
+//private:
+  UChar_t fDig[kgRowBytes];
+  ClassDef(digitRow,1)  // container for digit pattern
+};
+ClassImp(digitRow)
+
+
+////////////////////////////////////////////////////////////////////////
+//
+// Start of implementation of the class AliMCInfo
+//
+////////////////////////////////////////////////////////////////////////
+
+class AliMCInfo: public TObject {
+
+public:
+  AliMCInfo();
+  ~AliMCInfo();
+  void Update();
+
+  AliTrackReference  fTrackRef;      // track reference saved in the output tree
+  AliTrackReference  fTrackRefOut;   // decay track reference saved in the output tree
+  AliTrackReference  fTRdecay;       // track reference at decay point
+  //
+  TParticle fParticle;           // generated particle 
+  Float_t   fMass;               // mass of the particle
+  Float_t fCharge;               //
+  Int_t fLabel;                   // track label
+  Int_t fEventNr;                 // event number
+  Int_t fMCtracks;                // indication of how many times the track is retuturned back
+  Int_t fPdg;                     //pdg code
+  Float_t fDecayCoord[3];         // position of particle decay
+  Double_t fVDist[4];             //distance of the particle vertex from primary vertex
+  Bool_t fTPCdecay;               //indicates decay in TPC
+  Int_t fRowsWithDigitsInn;    // number of rows with digits in the inner sectors
+  Int_t fRowsWithDigits;       // number of rows with digits in the outer sectors
+  Int_t fRowsTrackLength;      // last - first row with digit
+  Float_t fPrim;               // theoretical dedx in tpc according particle momenta and mass
+  digitRow fTPCRow;               // information about digits row pattern
+  Int_t fNTPCRef;                 // tpc references counter
+  Int_t fNITSRef;                 // tpc references counter
+  TClonesArray * fTPCReferences;  //containner with all track references -in the TPC
+  TClonesArray * fITSReferences;  //container with ITS references
+  TClonesArray * fTRDReferences;  //container with TRD references  
+  //
+  ClassDef(AliMCInfo,3)  // container for 
+};
+ClassImp(AliMCInfo)
+
+
+
+class AliGenV0Info: public TObject {
+public:
+  AliMCInfo fMCd;      //info about daughter particle
+  AliMCInfo fMCm;      //info about mother particle
+  void Update();        // put some derived info to special field 
+  Double_t    fMCDist1;    //info about closest distance according closest MC - linear DCA
+  Double_t    fMCDist2;    //info about closest distance parabolic DCA
+  //
+  Double_t     fMCPdr[3];    //momentum at vertex daughter  - according approx at DCA
+  Double_t     fMCPd[4];     //exact momentum from MC info
+  Double_t     fMCX[3];      //exact position of the vertex
+  Double_t     fMCXr[3];     //rec. position according helix
+  //
+  Double_t     fMCPm[3];    //momentum at the vertex mother
+  Double_t     fMCAngle[3]; //three angels
+  Double_t     fMCRr;       // rec position of the vertex 
+  Double_t     fMCR;        //exact r position of the vertex
+  Int_t        fPdg[2];   //pdg code of mother and daugter particles
+  Int_t        fLab[2];   //MC label of the partecle
+  ClassDef(AliGenV0Info,1)  // container for  
+};
+ClassImp(AliGenV0Info)
+
+
+
+
+////////////////////////////////////////////////////////////////////////
+// 
+// Start of implementation of the class AliGenInfoMaker
+//
+////////////////////////////////////////////////////////////////////////
+
+class AliGenInfoMaker {
+
+public:
+  AliGenInfoMaker();
+  AliGenInfoMaker(const char * fnGalice, const char* fnRes    ="genTracks.root",
+                  Int_t nEvents=1, Int_t firstEvent=0);
+  virtual ~AliGenInfoMaker();
+  void Reset();
+  Int_t Exec();
+  Int_t Exec(Int_t nEvents, Int_t firstEventNr);
+  void CreateTreeGenTracks();
+  void CloseOutputFile();
+  Int_t TreeKLoop();
+  Int_t TreeTRLoop();
+  Int_t TreeDLoop();
+  //void  FillInfo(Int_t iParticle);
+  void SetFirstEventNr(Int_t i) {fFirstEventNr = i;}
+  void SetNEvents(Int_t i) {fNEvents = i;}
+  void SetDebug(Int_t level) {fDebug = level;}
+  Int_t SetIO(Int_t eventNr);
+  Int_t CloseIOEvent();
+  Int_t CloseIO();
+  Int_t SetIO();
+  Float_t TR2LocalX(AliTrackReference *trackRef,
+                   AliTPCParam *paramTPC);
+  AliMCInfo * GetInfo(UInt_t i){return (i<fNParticles)? fGenInfo[i]:0;}
+  AliMCInfo * MakeInfo(UInt_t i);
+
+public:
+  Int_t fDebug;                   //! debug flag  
+  Int_t fEventNr;                 //! current event number
+  Int_t fLabel;                   //! track label
+  Int_t fNEvents;                 //! number of events to process
+  Int_t fFirstEventNr;            //! first event to process
+  UInt_t fNParticles;              //! number of particles in TreeK
+  TTree *fTreeGenTracks;          //! output tree with generated tracks
+  char  fFnRes[1000];                   //! output file name with stored tracks
+  TFile *fFileGenTracks;          //! output file with stored fTreeGenTracks
+  //
+  AliRunLoader * fLoader;         //! pointer to the run loader
+  TTree * fTreeD;                 //! current tree with digits
+  TTree * fTreeTR;                //! current tree with TR
+  AliStack *fStack;               //! current stack
+  // 
+  AliMCInfo **   fGenInfo;    //! array with pointers to gen info
+  Int_t fNInfos;                  //! number of tracks with infos
+  //
+  AliTPCParam* fParamTPC;         //! AliTPCParam
+  Double_t fVPrim[3];             //! primary vertex position
+                                  // the fVDist[3] contains size of the 3-vector
+
+private:
+
+// some constants for the original non-pareller tracking (by Y.Belikov)
+  static const Int_t seedRow11 = 158;  // nRowUp - 1
+  static const Int_t seedRow12 = 139;  // nRowUp - 1 - (Int_t) 0.125*nRowUp
+  static const Int_t seedRow21 = 149;  // seedRow11 - shift
+  static const Int_t seedRow22 = 130;  // seedRow12 - shift
+  static const Double_t kRaddeg = 180./3.14159265358979312;
+  // 
+  static const Double_t fgTPCPtCut = 0.1; // do not store particles with generated pT less than this
+  static const Double_t fgITSPtCut = 0.2; // do not store particles with generated pT less than this
+  static const Double_t fgTRDPtCut = 0.2; // do not store particles with generated pT less than this
+  ClassDef(AliGenInfoMaker,1)    // class which creates and fills tree with TPCGenTrack objects
+};
+ClassImp(AliGenInfoMaker)
+
+
+
+class AliComparisonDraw: public TObject{
+public:
+  TH1F * DrawXY(const char * chx, const char *chy, const char* selection, 
+               const char * quality,Int_t nbins, Float_t minx, Float_t maxx, 
+               Float_t miny, Float_t maxy);
+  TH1F * DrawLogXY(const char * chx, const char *chy, const char* selection, 
+                  const char * quality, Int_t nbins,Float_t minx, Float_t maxx, 
+                  Float_t miny, Float_t maxy); 
+  TH1F * Eff(const char *variable, const char* selection, const char * quality, 
+            Int_t nbins,Float_t min, Float_t max); 
+  TH1F * EffLog(const char *variable, const char* selection, const char * quality, 
+            Int_t nbins,Float_t min, Float_t max);
+  //
+  static void   AliLabelAxes(TH1* histo, const char* xAxisTitle, const char* yAxisTitle);
+  static Double_t* CreateLogBins(Int_t nBins, Double_t xMin, Double_t xMax);
+  //
+  static TH1F*  CreateEffHisto(TH1F* hGen, TH1F* hRec);
+  static TH1F*  CreateResHisto(TH2F* hRes2, TH1F **phMean, 
+                               Bool_t drawBinFits = kTRUE,Bool_t overflowBinFits = kFALSE);
+ public: 
+  TTree * fTree;
+  TH1F  * fRes;  //temporary file
+  TH1F  * fMean;  //temporary file
+  ClassDef(AliComparisonDraw,1)
+};
+ClassImp(AliComparisonDraw)
+
+
+
+AliTPCParam * GetTPCParam();
+Float_t TPCBetheBloch(Float_t bg);
diff --git a/STEER/AliITSclusterComparison.C b/STEER/AliITSclusterComparison.C
new file mode 100644 (file)
index 0000000..52da02f
--- /dev/null
@@ -0,0 +1,892 @@
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include <Riostream.h>
+
+#include "AliITStrackerMI.h"
+#include "AliITSclusterV2.h"
+
+#include "AliITShit.h"
+#include "AliITS.h"
+#include "AliITSgeom.h"
+//#include "AliITSgeometryDetail.h"
+//#include "AliITSparameter.h"
+#include "alles.h"
+#include "TFile.h"
+#include "TStopwatch.h"
+#include "Rtypes.h"
+#include "TTree.h"
+
+
+#include "AliRunLoader.h"
+#include "AliStack.h"
+#include "TF1.h"
+#include "AliTrackReference.h"
+#endif    
+#include "AliITSclusterComparison.h"
+
+
+ClassImp(AliITSCI)
+ClassImp(AliITSClusterErrAnal)
+
+
+
+AliITSClusterErrAnal::AliITSClusterErrAnal(Char_t *chloader )
+{
+  //
+  //SET Input loaders
+  if (gAlice){
+    delete gAlice->GetRunLoader();
+    delete gAlice;
+    gAlice = 0x0;
+  }    
+  fRunLoader = AliRunLoader::Open(chloader);
+  if (fRunLoader == 0x0){
+    cerr<<"Can not open session"<<endl;
+    return;
+  }
+  fITSLoader = fRunLoader->GetLoader("ITSLoader");  
+  if (fITSLoader == 0x0){
+    cerr<<"Can not get ITS Loader"<<endl;
+    return ;
+  }   
+  if (fRunLoader->LoadgAlice()){
+    cerr<<"Error occured while l"<<endl;
+    return;
+  }
+  fRunLoader->CdGAFile();  
+  fITS = (AliITS*)gAlice->GetDetector("ITS");
+  if (!fITS) {
+    cerr<<"AliITSclusterComparison.C : Can not find the ITS detector !"<<endl;
+  }
+  AliITSgeom *geom = fITS->GetITSgeom();
+  //An instance of the ITS tracker
+  fTracker = new AliITStrackerMI(geom);
+  //
+  //
+  AliITSCI * clinfo = new AliITSCI();
+  fFileA  = new TFile("itsclusteranal.root","recreate");
+  fFileA->cd();
+  fTreeA  = new TTree("itscl","itscl");
+  fTreeA->Branch("itscl","AliITSCI",&clinfo);
+
+  AliITSclusterV2 * cl = new AliITSclusterV2;
+  fTreeB  = new TTree("Clusters","Clusters");
+  fTreeB->Branch("cl","AliITSclusterV2",&cl);
+
+  fClusters = 0;
+  delete clinfo;
+}
+
+void AliITSClusterErrAnal::SetIO(Int_t event)
+{
+  //
+  //set input output for given event
+  fRunLoader->SetEventNumber(event);
+  fRunLoader->LoadHeader();
+  fRunLoader->LoadKinematics();
+  fRunLoader->LoadTrackRefs();
+  fITSLoader->LoadHits();
+  fITSLoader->LoadRecPoints("read");
+  //
+  fStack = fRunLoader->Stack();
+  fHitTree = fITSLoader->TreeH();
+  fClusterTree = fITSLoader->TreeR();
+  fReferenceTree = fRunLoader->TreeTR();
+  //
+}
+
+void AliITSClusterErrAnal::LoadClusters()
+{
+  //
+  //
+  // Load clusters
+  if (fClusters) {
+    fClusters->Delete();
+    delete fClusters;
+  }
+  Int_t nparticles = fStack->GetNtrack();
+  fClusters = new TObjArray(nparticles);
+  //
+  TObjArray *ClusterArray = new TClonesArray("AliITSclusterV2");   
+  TObjArray carray(2000);
+  TBranch *branch=fClusterTree->GetBranch("Clusters");
+  if (!branch) {
+    Error("ReadClusters","Can't get the branch !");
+    return;
+  }
+  branch->SetAddress(&ClusterArray);
+  TBranch *brcl = fTreeB->GetBranch("cl");
+  AliITSclusterV2* cluster= new AliITSclusterV2;
+  brcl->SetAddress(&cluster);
+  Int_t nentries = (Int_t)fClusterTree->GetEntries();
+  for (Int_t i=0;i<nentries;i++){
+    fClusterTree->GetEvent(i);
+    Int_t nCluster = ClusterArray->GetEntriesFast();
+    SignDeltas(ClusterArray,-3.8);
+    for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) { 
+      cluster = (AliITSclusterV2*)ClusterArray->UncheckedAt(iCluster);
+      carray.AddLast(new AliITSclusterV2(*cluster));
+      fTreeB->Fill();
+    }
+  }
+  //
+  Int_t nClusters = carray.GetEntriesFast();
+  printf("Total number of clusters %d \n", nClusters);
+  fTracker->LoadClusters(fClusterTree);
+  //
+  //
+  //SORT clusters
+  //
+  Int_t all=0;
+  Int_t noisecl =0;
+  for (Int_t i=0;i<nClusters;i++){
+    AliITSclusterV2 *cl = (AliITSclusterV2 *) carray.UncheckedAt(i);  
+    if (cl->GetLabel(0)<0) noisecl++;
+    //
+    for (Int_t itrack=0; itrack<3;itrack++){
+      Int_t lab = cl->GetLabel(itrack);
+      if (lab>=0){
+       TObjArray * array = (TObjArray*)fClusters->At(lab);
+       if (array==0){
+         array = new TObjArray(20);
+         fClusters->AddAt(array,lab);
+       }
+       array->AddLast(cl);
+       all++;
+      }        
+    }    
+  }
+  printf("Total number of track clusters %d \n", all);
+  printf("Total number of noise clusters %d \n", noisecl);
+}
+
+void AliITSClusterErrAnal::SignDeltas( TObjArray *ClusterArray, Float_t vz)
+{
+  //
+  //  
+  Int_t entries = ClusterArray->GetEntriesFast();
+  if (entries<4) return;
+  AliITSclusterV2* cluster = (AliITSclusterV2*)ClusterArray->At(0);
+  Int_t layer = cluster->GetLayer();
+  if (layer>1) return;
+  Int_t index[10000];
+  Int_t ncandidates=0;
+  Float_t r = (layer>0)? 7:4;
+  // 
+  for (Int_t i=0;i<entries;i++){
+    AliITSclusterV2* cl0 = (AliITSclusterV2*)ClusterArray->At(i);
+    Float_t nz = 1+(cl0->GetZ()-vz)/r;
+    if (cl0->GetNy()+cl0->GetNz()<=5+2*layer+nz) continue;
+    index[ncandidates] = i;  //candidate to belong to delta electron track
+    ncandidates++;
+    if (cl0->GetNy()+cl0->GetNz()>9+2*layer+nz) {
+      cl0->SetDeltaProbability(1);
+    }
+  }
+  //
+  //  
+  //
+  for (Int_t i=0;i<ncandidates;i++){
+    AliITSclusterV2* cl0 = (AliITSclusterV2*)ClusterArray->At(index[i]);
+    if (cl0->GetDeltaProbability()>0.8) continue;
+    // 
+    Int_t ncl = 0;
+    Float_t y[100],z[100],sumy,sumz,sumy2, sumyz, sumw;
+    sumy=sumz=sumy2=sumyz=sumw=0.0;
+    for (Int_t j=0;j<ncandidates;j++){
+      if (i==j) continue;
+      AliITSclusterV2* cl1 = (AliITSclusterV2*)ClusterArray->At(index[j]);
+      //
+      Float_t dz = cl0->GetZ()-cl1->GetZ();
+      Float_t dy = cl0->GetY()-cl1->GetY();
+      if (TMath::Sqrt(dz*dz+dy*dy)<0.2){
+       Float_t weight = cl1->GetNy()+cl1->GetNz()-2;
+       y[ncl] = cl1->GetY();
+       z[ncl] = cl1->GetZ();
+       sumy+= y[ncl]*weight;
+       sumz+= z[ncl]*weight;
+       sumy2+=y[ncl]*y[ncl]*weight;
+       sumyz+=y[ncl]*z[ncl]*weight;
+       sumw+=weight;
+       ncl++;
+      }
+    }
+    if (ncl<4) continue;
+    Float_t det = sumw*sumy2  - sumy*sumy;
+    Float_t delta=1000;
+    if (TMath::Abs(det)>0.01){
+      Float_t z0  = (sumy2*sumz - sumy*sumyz)/det;
+      Float_t k   = (sumyz*sumw - sumy*sumz)/det;
+      delta = TMath::Abs(cl0->GetZ()-(z0+k*cl0->GetY()));
+    }
+    else{
+      Float_t z0  = sumyz/sumy;
+      delta = TMath::Abs(cl0->GetZ()-z0);
+    }
+    if ( delta<0.05) {
+      cl0->SetDeltaProbability(1-20.*delta);
+    }   
+  }
+}
+
+void AliITSClusterErrAnal::LoadParticles()
+{
+  //
+  //
+  //  fReferences = new TObjArray;
+}
+
+void AliITSClusterErrAnal::SortReferences()
+{
+  //
+  //
+  //
+  printf("Sorting references\n");
+  fReferences = new TObjArray;
+  Int_t ntracks = fStack->GetNtrack();
+  fReferences->Expand(ntracks);
+  Int_t nentries = (Int_t)fReferenceTree->GetEntries();
+  TClonesArray * arr = new TClonesArray("AliTrackReference");
+  TBranch * br = fReferenceTree->GetBranch("ITS");
+  br->SetAddress(&arr);
+  //
+  TClonesArray *labarr=0;
+  Int_t nreferences=0;
+  Int_t nreftracks=0;
+  for (Int_t iprim=0;iprim<nentries;iprim++){
+    if (br->GetEntry(iprim)){
+      for (Int_t iref=0;iref<arr->GetEntriesFast();iref++){
+       AliTrackReference *ref =(AliTrackReference*)arr->At(iref);
+       if (!ref) continue;
+       Int_t lab = ref->GetTrack();
+       //if ( (lab<0) || (lab>ntracks)) continue;
+       //
+       if (fReferences->At(lab)==0) {
+         labarr = new TClonesArray("AliTrackReference"); 
+         fReferences->AddAt(labarr,lab);
+         nreftracks++;
+       }
+       TClonesArray &larr = *labarr;
+       new(larr[larr.GetEntriesFast()]) AliTrackReference(*ref);
+       nreferences++;
+      }
+    }
+  }
+  printf("Total number of references = \t%d\n", nreferences);
+  printf("Total number of tracks with references = \t%d\n", nreftracks);
+  printf("End - Sorting references\n");
+  
+}
+
+
+void AliITSClusterErrAnal::GetNTeor(Int_t layer, Float_t theta, Float_t phi, Float_t &ny, Float_t &nz)
+{
+  //
+  //get "mean shape"
+  if (layer==0){
+    ny = 1.+TMath::Abs(phi)*3.2;
+    nz = 1.+TMath::Abs(theta)*0.34;
+    return;
+  }
+  if (layer==1){
+    ny = 1.+TMath::Abs(phi)*3.2;
+    nz = 1.+TMath::Abs(theta)*0.28;
+    return;
+  }
+  
+  if (layer>3){
+    ny = 2.02+TMath::Abs(phi)*1.95;
+    nz = 2.02+TMath::Abs(phi)*2.35;
+    return;
+  }
+  ny  = 6.6-2.7*abs(phi);
+  nz  = 2.8-3.11*TMath::Abs(phi)+0.45*TMath::Abs(theta);
+}
+
+
+Int_t AliITSClusterErrAnal::GetError(Int_t layer, const AliITSclusterV2*cl, Float_t theta, Float_t phi,Float_t expQ, Float_t &erry, Float_t &errz)
+{
+  //calculate cluster position error
+  //
+  Float_t nz,ny;
+  GetNTeor(layer, theta,phi,ny,nz);  
+  erry   = TMath::Sqrt(cl->GetSigmaY2()); 
+  errz   = TMath::Sqrt(cl->GetSigmaZ2()); 
+  //
+  // PIXELS
+  if (layer<2){
+    
+    if (TMath::Abs(ny-cl->GetNy())>0.6)  {
+      if (ny<cl->GetNy()){
+       erry*=0.4+TMath::Abs(ny-cl->GetNy());
+       errz*=0.4+TMath::Abs(ny-cl->GetNy());
+      }else{
+       erry*=0.7+0.5*TMath::Abs(ny-cl->GetNy());
+       errz*=0.7+0.5*TMath::Abs(ny-cl->GetNy());
+      }
+    }
+    if (TMath::Abs(nz-cl->GetNz())>1.)  {
+      erry*=TMath::Abs(nz-cl->GetNz());
+      errz*=TMath::Abs(nz-cl->GetNz());              
+    }
+    erry*=0.85;
+    errz*=0.85;
+    erry= TMath::Min(erry,float(0.005));
+    errz= TMath::Min(errz,float(0.03));
+    return 10;
+  }
+  //STRIPS
+  if (layer>3){ 
+    if (cl->GetNy()==100||cl->GetNz()==100){
+      erry = 0.004;
+      errz = 0.2;
+      return 100;
+    }
+    if (cl->GetNy()+cl->GetNz()>12){
+      erry = 0.06;
+      errz = 0.57;
+      return 100;
+    }
+    Float_t normq = cl->GetQ()/(TMath::Sqrt(1+theta*theta+phi*phi));
+    Float_t chargematch = TMath::Max(double(normq/expQ),2.);
+    //
+    if (cl->GetType()==1 || cl->GetType()==10 ){                                                              
+      if (chargematch<1.0 || (cl->GetNy()+cl->GetNz()<nz+ny+0.5)){
+       errz = 0.043;
+       erry = 0.00094;
+       return 101;
+      }
+      if (cl->GetNy()+cl->GetNz()<nz+ny+1.2){
+       errz = 0.06;
+       erry =0.0013;
+       return 102;
+      }
+      erry = 0.0027;
+      errz = TMath::Min(0.028*(chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.15);
+      return 103;
+    }
+    if (cl->GetType()==2 || cl->GetType()==11 ){ 
+      erry = TMath::Min(0.0010*(1+chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.05);
+      errz = TMath::Min(0.025*(1+chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.5);
+      return 104;
+    }
+    
+    if (cl->GetType()>100 ){                                                                  
+      if ((chargematch+cl->GetNy()+cl->GetNz()-nz-ny<1.5)){
+       errz = 0.05;
+       erry = 0.00096;
+       return 105;
+      }
+      if (cl->GetNy()+cl->GetNz()-nz-ny<1){
+       errz = 0.10;
+       erry = 0.0025;
+       return 106;
+      }
+
+      errz = TMath::Min(0.05*(chargematch+cl->GetNy()+cl->GetNz()-nz-ny),0.4);
+      erry = TMath::Min(0.003*(chargematch+cl->GetNy()+cl->GetNz()-nz-ny),0.05);
+      return 107;
+    }    
+    Float_t diff = cl->GetNy()+cl->GetNz()-ny-nz;
+    if (diff<1) diff=1;
+    if (diff>4) diff=4;
+        
+    if (cl->GetType()==5||cl->GetType()==6||cl->GetType()==7||cl->GetType()==8){
+      errz = 0.14*diff;
+      erry = 0.003*diff;
+      return 108;
+    }  
+    erry = 0.04*diff;
+    errz = 0.06*diff;
+    return 109;
+  }
+
+  //DRIFTS
+  Float_t normq = cl->GetQ()/(TMath::Sqrt(1+theta*theta+phi*phi));
+  Float_t chargematch = normq/expQ;
+  Float_t factorz=1;
+  Int_t   cnz = cl->GetNz()%10;
+  //charge match
+  if (cl->GetType()==1){
+    if (chargematch<1.25){
+      erry =  0.0028*(1.+6./cl->GetQ());  // gold clusters
+    }
+    else{
+      erry = 0.003*chargematch;
+      if (cl->GetNz()==3) erry*=1.5;
+    }
+    if (chargematch<1.0){
+      errz =  0.0011*(1.+6./cl->GetQ());
+    }
+    else{
+      errz = 0.002*(1+2*(chargematch-1.));
+    }
+    if (cnz>nz+0.6) {
+      erry*=(cnz-nz+0.5);
+      errz*=1.4*(cnz-nz+0.5);
+    }
+  }
+  if (cl->GetType()>1){
+    if (chargematch<1){
+      erry =  0.00385*(1.+6./cl->GetQ());  // gold clusters
+      errz =  0.0016*(1.+6./cl->GetQ());
+    }
+    else{
+      errz = 0.0014*(1+3*(chargematch-1.));
+      erry = 0.003*(1+3*(chargematch-1.));
+    } 
+    if (cnz>nz+0.6) {
+      erry*=(cnz-nz+0.5);
+      errz*=1.4*(cnz-nz+0.5);
+    }
+  }
+
+  if (TMath::Abs(cl->GetY())>2.5){
+    factorz*=1+2*(TMath::Abs(cl->GetY())-2.5);
+  }
+  if (TMath::Abs(cl->GetY())<1){
+    factorz*=1.+0.5*TMath::Abs(TMath::Abs(cl->GetY())-1.);
+  }
+  factorz= TMath::Min(factorz,float(4.));  
+  errz*=factorz;
+
+  erry= TMath::Min(erry,float(0.05));
+  errz= TMath::Min(errz,float(0.05));  
+  return 200;
+}
+
+
+AliITSclusterV2 * AliITSClusterErrAnal::FindNearestCluster(AliITSCI * clinfo,  Float_t dmax)
+{
+  //
+  //
+  AliTrackReference * ref = &(clinfo->fRef);
+  TObjArray * clusters = (TObjArray*)(fClusters->At(ref->GetTrack()));
+  if (!clusters) return 0;
+  Int_t nclusters = clusters->GetEntriesFast();  
+  Float_t maxd = dmax;
+  AliITSclusterV2 * cluster=0;
+  Float_t radius = ref->R();
+  clinfo->fNClusters=0;
+  for (Int_t icluster=0;icluster<nclusters;icluster++){
+    AliITSclusterV2 * cl = (AliITSclusterV2*)clusters->At(icluster);
+    if (!cl) continue;
+    Float_t dz = cl->GetZ()-ref->Z();
+    //
+    if (TMath::Abs(dz)>dmax) continue;        
+    Int_t idet=cl->GetDetectorIndex();        
+    Double_t r=-1;
+    Float_t  dy =0;
+    //    for (Int_t i=cl->fLayer;i<=cl->fLayer-1;i++){
+    
+    AliITStrackerMI::AliITSdetector & det = fTracker->GetDetector(cl->fLayer,idet);
+      if (TMath::Abs(radius-det.GetR())<1.5) {
+       Bool_t isOK=kFALSE;
+       for (Int_t itrack=0;itrack<3;itrack++){
+         if (cl->fTracks[itrack]== clinfo->fRef.GetTrack()) {
+           clinfo->fNClusters++;
+           isOK= kTRUE;
+         }
+       }
+       if (!isOK) continue;
+       r = det.GetR();
+       Double_t phicl= det.GetPhi();
+       Double_t ca=TMath::Cos(phicl), sa=TMath::Sin(phicl);
+       Double_t x =  double(ref->X())*ca +double(ref->Y())*sa;
+       Double_t y = -double(ref->X())*sa +double(ref->Y())*ca;
+       if (TMath::Abs(x-r-0.012)>0.1) continue;
+       //
+       Double_t px = ref->Px(), py = ref->Py();
+       Double_t lpx =  double(px)*ca +double(py)*sa;
+       Double_t lpy =  -double(px)*sa +double(py)*ca;
+       Double_t lphi =  TMath::ATan2(lpy,lpx);
+
+       Double_t theta  = ref->Pz()/ref->Pt();
+       Double_t dx = x-r;
+       dy = y       -cl->GetY()-dx*lphi;
+       dz = ref->Z()-cl->GetZ()-dx*theta;
+       //dx-=0.012;
+       Double_t d = TMath::Sqrt(dy*dy+dz*dz);
+       if (d<maxd){
+         cluster=cl;
+         maxd=d;
+         clinfo->Update(fTracker);
+         clinfo->fCl = *cluster;
+         clinfo->fLx    = x;
+         clinfo->fLy    = y;
+         clinfo->fLayer = cl->fLayer;
+         clinfo->fPhiCl = phicl;
+         clinfo->fRCl   = det.GetR();
+         clinfo->fPhiL  = lphi;
+         clinfo->fDz    = dz;
+         clinfo->fDy    = dy;
+         cl->SetSigmaZ2(TMath::Abs(cl->GetSigmaZ2()));
+         cl->SetSigmaY2(TMath::Abs(cl->GetSigmaY2()));
+         clinfo->fPoolZ  = dz/TMath::Sqrt(cl->GetSigmaZ2());
+         clinfo->fPoolY  = dy/TMath::Sqrt(cl->GetSigmaY2());
+         clinfo->fPoolY2  = dy/TMath::Sqrt(cl->GetSigmaY2());    
+         clinfo->fErrY   = TMath::Sqrt(cl->GetSigmaY2()); 
+         clinfo->fErrZ   = TMath::Sqrt(cl->GetSigmaZ2()); 
+         Float_t erry,errz,ny,nz;
+         clinfo->fErrType = GetError(cl->fLayer,cl,theta,lphi,clinfo->fExpectedQ,erry,errz);
+         GetNTeor(cl->fLayer, theta,lphi,ny,nz);  
+         clinfo->fNormQ = cl->GetQ()/TMath::Sqrt(1+theta*theta+lpy*lpy/(lpx*lpx));
+         clinfo->fTNy  = ny;
+         clinfo->fTNz  = nz;
+         clinfo->fErrY = erry;
+         clinfo->fErrZ = errz;  
+         clinfo->fPoolY2 = dy/clinfo->fErrY;
+         clinfo->fPoolZ2 = dz/clinfo->fErrZ;
+         clinfo->fLabPos =-1;
+         Double_t xyz[3];
+         det.GetGlobalXYZ(cl,xyz);
+         clinfo->fGx = xyz[0];
+         clinfo->fGy = xyz[1];
+         for (Int_t i=0;i<3;i++) if (cl->GetLabel(i)==clinfo->fRef.GetTrack()) clinfo->fLabPos = i;
+       }
+      }
+    }
+  //  }
+  return cluster;
+}
+
+
+
+Int_t AliITSClusterErrAnal::Analyze(Int_t trackmax) {    
+  //
+  //
+  // dummy cluster to be fill if not cluster info
+  //
+  AliITSCI * clinfo = new AliITSCI;
+  AliITSclusterV2 dummy;
+  TBranch *  branch = fTreeA->GetBranch("itscl");
+  branch->SetAddress(&clinfo);
+  Int_t nall =0;
+  Int_t withcl =0;
+  Int_t ntracks = fReferences->GetEntriesFast();
+  ntracks = TMath::Min(trackmax,ntracks);
+  for (Int_t itrack=0;itrack<ntracks;itrack++){
+    TClonesArray * references = (TClonesArray*) fReferences->At(itrack);
+    if (!references) continue;
+    Int_t nreferences = references->GetEntriesFast();
+    if (nreferences<4) continue;
+    TObjArray * clarray = (TObjArray*)fClusters->At(itrack);
+    for (Int_t iref=0;iref<nreferences;iref++){
+      AliTrackReference* trackref = (AliTrackReference*)references->At(iref);     
+      if (trackref&& trackref->P()>0.01){
+       TParticle * p = fStack->Particle(trackref->GetTrack());
+       if (p) clinfo->fP = *p;
+       clinfo->fRef = *trackref;
+       clinfo->fStatus = 0;
+       clinfo->Update(fTracker);
+       nall++;
+       AliITSclusterV2 * cluster = FindNearestCluster(clinfo,3.0);
+       if (cluster){
+         //clinfo->fCl = *cluster;
+         clinfo->fStatus=1;
+         withcl++;
+       }else{
+         clinfo->fCl= dummy;
+       }
+       clinfo->Update(fTracker);
+       fTreeA->Fill();
+      }
+    }
+  }
+  
+  fFileA->cd();
+  fTreeA->Write();
+  fTreeB->Write();
+  fFileA->Close();
+  return 0;
+}
+
+
+
+void AliITSClusterErrAnal::MakeSeeds(Double_t zv)
+{
+  //
+  //
+  // 
+
+  AliITStrackerMI::AliITSlayer & layer3 =   fTracker->GetLayer(3); 
+  AliITStrackerMI::AliITSlayer & layer2 =   fTracker->GetLayer(2); 
+  AliITStrackerMI::AliITSlayer & layer1 =   fTracker->GetLayer(1); 
+  AliITStrackerMI::AliITSlayer & layer0 =   fTracker->GetLayer(0); 
+  Int_t ncl3 =  layer3.GetNumberOfClusters();
+  Int_t ncl2 =  layer2.GetNumberOfClusters();
+  Int_t ncl1 =  layer1.GetNumberOfClusters();
+  Int_t ncl0 =  layer0.GetNumberOfClusters();
+  TH1F *hc = new TH1F("hc","hc",100,-0.01,0.01);
+  TH1F *hrc = new TH1F("hrc","hrc",100,-1,1);
+  TH1F *hl = new TH1F("hl","hl",100,-0.05,0.05);
+  TH1F *hlr = new TH1F("hlr","hlr",1000,-30.,30.);
+  TH2F *hcl = new TH2F("hcl","hcl",1000,-0.01,0.01,1000,-0.03,0.03);
+  TH2F *hcc = new TH2F("hcc","hcc",1000,-0.01,0.01,1000,-0.03,0.03);
+  TH1F *hr = new TH1F("hr","hr",100,0,2000.);
+  //
+  TH1F *hpoolc = new TH1F("hpoolc","hpoolc",100,-5,5);
+  TH1F *hpooll = new TH1F("hpooll","hpooll",100,-5,5);
+  //
+  Double_t ratio23 = layer2.GetR()/layer3.GetR();
+  Double_t ratio12 = layer1.GetR()/layer2.GetR();
+  Double_t referenceR =  (layer0.GetR()+layer1.GetR()+layer2.GetR()+layer3.GetR())*0.25;
+  //
+  Double_t *clusterxyzr3 = new Double_t[ncl3*6];
+  Double_t *clusterxyzr2 = new Double_t[ncl2*6];
+  Double_t *clusterxyzr1 = new Double_t[ncl1*6];
+  Double_t *clusterxyzr0 = new Double_t[ncl0*6];
+  //
+  // Get global coordinates 3
+  for (Int_t icl3=0;icl3<ncl3;icl3++){
+    AliITSclusterV2* cl3 = layer3.GetCluster(icl3);
+    Double_t *xyz3 = &clusterxyzr3[6*icl3];
+    AliITStrackerMI::AliITSdetector & det = layer3.GetDetector(cl3->GetDetectorIndex());
+    det.GetGlobalXYZ(cl3,xyz3);
+    xyz3[3] = TMath::Sqrt(xyz3[0]*xyz3[0]+xyz3[1]*xyz3[1]);
+    Double_t fi3 = TMath::ATan2(xyz3[1],xyz3[0]);
+    xyz3[4] = TMath::Cos(fi3);
+    xyz3[5] = TMath::Sin(fi3);
+  }  
+  //
+  // Get global coordinates 2
+  for (Int_t icl2=0;icl2<ncl2;icl2++){
+    AliITSclusterV2* cl2 = layer2.GetCluster(icl2);
+    Double_t *xyz2 = &clusterxyzr2[6*icl2];
+    AliITStrackerMI::AliITSdetector & det = layer2.GetDetector(cl2->GetDetectorIndex());
+    det.GetGlobalXYZ(cl2,xyz2);
+    xyz2[3] = TMath::Sqrt(xyz2[0]*xyz2[0]+xyz2[1]*xyz2[1]);
+    Double_t fi2 = TMath::ATan2(xyz2[1],xyz2[0]);
+    xyz2[4] = TMath::Cos(fi2);
+    xyz2[5] = TMath::Sin(fi2);
+  }  
+  //
+  // Get global coordinates 1
+  for (Int_t icl1=0;icl1<ncl1;icl1++){
+    AliITSclusterV2* cl1 = layer1.GetCluster(icl1);
+    Double_t *xyz1 = &clusterxyzr1[6*icl1];
+    AliITStrackerMI::AliITSdetector & det = layer1.GetDetector(cl1->GetDetectorIndex());
+    det.GetGlobalXYZ(cl1,xyz1);
+    xyz1[3] = TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
+    Double_t fi1 = TMath::ATan2(xyz1[1],xyz1[0]);
+    xyz1[4] = TMath::Cos(fi1);
+    xyz1[5] = TMath::Sin(fi1);
+  }  
+  //
+  // Get global coordinates 0
+  for (Int_t icl0=0;icl0<ncl0;icl0++){
+    AliITSclusterV2* cl0 = layer0.GetCluster(icl0);
+    Double_t *xyz0 = &clusterxyzr0[6*icl0];
+    AliITStrackerMI::AliITSdetector & det = layer0.GetDetector(cl0->GetDetectorIndex());
+    det.GetGlobalXYZ(cl0,xyz0);
+    xyz0[3] = TMath::Sqrt(xyz0[0]*xyz0[0]+xyz0[1]*xyz0[1]);
+    Double_t fi0 = TMath::ATan2(xyz0[1],xyz0[0]);
+    xyz0[4] = TMath::Cos(fi0);
+    xyz0[5] = TMath::Sin(fi0);
+
+  }  
+  //
+  //
+  Int_t ntracks=0;
+  Int_t ntracks0=0;
+  Int_t good0=0;
+  Int_t accepted0=0;
+  //
+  Int_t nfoundtracks=0;
+  Int_t good1 =0;;
+  Int_t accepted1=0;
+  Int_t accepted2=0;
+  //
+  Int_t *trackindex = new Int_t[ncl3*4*300];
+  const AliITSclusterV2 **trackclusters = new const AliITSclusterV2*[ncl3*4*300];
+  //
+  for (Int_t icl3=0;icl3<ncl3;icl3++){
+    AliITSclusterV2* cl3 = layer3.GetCluster(icl3);
+    Double_t *xyz3 = &clusterxyzr3[icl3*6];
+    Double_t &r3   = clusterxyzr3[icl3*6+3];
+    ntracks=0;  //number of possible tracks corresponding to given cluster
+    //
+    //
+    // find pairs in layer 2
+    Float_t zc   = (xyz3[2]-zv)*ratio23+zv;
+    Float_t lphi  = TMath::ATan2(xyz3[1],xyz3[0]);
+    layer2.SelectClusters(zc-1,zc+1, layer2.GetR()*lphi-1, layer2.GetR()*lphi+1.);
+    const AliITSclusterV2 * cl2=0;
+    Int_t npairs2=0;
+    Int_t ci2[100];
+    while ( (cl2=layer2.GetNextCluster(ci2[npairs2]))!=0){
+      Int_t index2 = ci2[npairs2];
+      Double_t *xyz2 = &clusterxyzr2[index2*6];
+      Double_t &r2   = clusterxyzr2[index2*6+3];
+      //      
+      if (TMath::Abs((xyz3[2]-zv)*r2/r3-(xyz2[2]-zv))>0.3) continue;
+      npairs2++;
+    }    
+    //
+    //
+    // find pairs in layer 2
+    for (Int_t jcl2 =0;jcl2<npairs2; jcl2++){
+      Int_t index2 = ci2[jcl2];
+      Double_t *xyz2 = &clusterxyzr2[index2*6];
+      Float_t zc   = (xyz2[2]-zv)*ratio12+zv;
+      Float_t lphi  = TMath::ATan2(xyz2[1],xyz2[0]);
+      cl2 = layer2.GetCluster(index2);
+      //
+      //      
+      layer1.SelectClusters(zc-1,zc+1, layer1.GetR()*lphi-1, layer1.GetR()*lphi+1.);
+      const AliITSclusterV2 * cl1=0;
+      Int_t ci1;
+      //
+      Double_t lr2 = (xyz3[0]-xyz2[0])*(xyz3[0]-xyz2[0])+(xyz3[1]-xyz2[1])*(xyz3[1]-xyz2[1]);
+      Double_t lr  = TMath::Sqrt(lr2); 
+      Double_t sinphi = (xyz3[1]-xyz2[1])/lr;
+      Double_t cosphi = (xyz3[0]-xyz2[0])/lr;
+      Double_t t1 = (xyz3[2]-xyz2[2])/lr;
+      //
+      while ( (cl1=layer1.GetNextCluster(ci1))!=0){
+       if (cl1->GetQ()==0) continue;
+       Int_t index1 = ci1;
+       Double_t *xyz1 = &clusterxyzr1[index1*6];
+       AliITSclusterV2* cl1 = layer1.GetCluster(index1);
+  
+       if (cl3->GetLabel(0)==cl2->GetLabel(0) && cl3->GetLabel(0)==cl1->GetLabel(0)) good0++;
+       Double_t loc1[2];
+       //
+       loc1[0] =  (xyz1[0]-xyz2[0])*cosphi+(xyz1[1]-xyz2[1])*sinphi;
+       loc1[1] = -(xyz1[0]-xyz2[0])*sinphi+(xyz1[1]-xyz2[1])*cosphi;
+       //
+       Double_t c1 = 2.*loc1[1]/(loc1[0]*(loc1[0]-lr));   //curvature approximated using 2 derivation
+       Double_t c=c1;
+       Double_t t2 = (xyz2[2]-xyz1[2])/TMath::Sqrt(loc1[0]*loc1[0]+loc1[1]*loc1[1]); 
+       if (TMath::Abs(c)>0.008) continue;
+       if (TMath::Abs(t2-t1)>0.02) continue;
+       //      
+       trackclusters[ntracks*4]    = 0;
+       trackclusters[ntracks*4+1]  = cl1;
+       trackclusters[ntracks*4+2]  = cl2;
+       trackclusters[ntracks*4+3]  = cl3;
+       trackindex[ntracks*4+0] = 0;
+       trackindex[ntracks*4+1] = index1;
+       trackindex[ntracks*4+2] = index2;
+       trackindex[ntracks*4+3] = icl3;
+       if (ntracks>=300*ncl3){   
+         break;
+       }
+       if (cl3->GetLabel(0)==cl2->GetLabel(0) && cl3->GetLabel(0)==cl1->GetLabel(0)) {
+         //printf("%f\n",c);
+         
+         //hc->Fill(c);
+         //hl->Fill(t2-t1);
+         //hr->Fill(1/c);
+         //hcc->Fill(c,TMath::Abs(c1));
+         //hcl->Fill(c1,t2-t1);
+         
+         accepted0++;
+       }
+       ntracks++;
+      }
+    }   // end - find pairs in layer 2
+    //    
+    ntracks0+=ntracks;
+    //hcc->Draw();
+    Int_t   ntracks2=0;
+    Int_t   cindex2[10000*4];
+    Float_t pools[10000];
+    
+    for (Int_t itrack=0;itrack<ntracks;itrack++){
+      const AliITSclusterV2 *cl0, *cl1,*cl2,*cl3;
+      cl3 = trackclusters[itrack*4+3];
+      cl2 = trackclusters[itrack*4+2];
+      cl1 = trackclusters[itrack*4+1];
+      Int_t index1 = trackindex[itrack*4+1];
+      Int_t index2 = trackindex[itrack*4+2];
+      Int_t index3 = trackindex[itrack*4+3];
+      //
+      Double_t *xyz1 = &clusterxyzr1[index1*6],*xyz2 = &clusterxyzr2[index2*6],*xyz3=&clusterxyzr3[index3*6];
+      Double_t &r1 = clusterxyzr1[index1*6+3],&r2 = clusterxyzr2[index2*6+3], &r3 = clusterxyzr3[index3*6+3];
+      
+      Double_t r0 = layer0.GetR();
+      Float_t zc   = xyz1[2] + (xyz2[2]-xyz1[2])*(r0-r1)/(r2-r1);
+      Float_t lphi  = TMath::ATan2(xyz1[1],xyz1[0]);
+      
+      layer0.SelectClusters(zc-1,zc+1, layer0.GetR()*lphi-0.3, layer0.GetR()*lphi+0.3);
+      Int_t ci0;
+      //
+      //
+      Double_t lr2 = (xyz2[0]-xyz1[0])*(xyz2[0]-xyz1[0])+(xyz2[1]-xyz1[1])*(xyz2[1]-xyz1[1]);
+      Double_t lr  = TMath::Sqrt(lr2); 
+      Double_t sinphi = (xyz2[1]-xyz1[1])/lr;
+      Double_t cosphi = (xyz2[0]-xyz1[0])/lr;
+      //
+      Double_t loc3[2]; // cluster3 in local coordindate 
+      loc3[0] =  (xyz3[0]-xyz1[0])*cosphi+(xyz3[1]-xyz1[1])*sinphi;
+      loc3[1] = -(xyz3[0]-xyz1[0])*sinphi+(xyz3[1]-xyz1[1])*cosphi;
+      
+      Double_t tan1 = (xyz3[2]-xyz2[2])/TMath::Sqrt((xyz3[0]-xyz2[0])*(xyz3[0]-xyz2[0])+(xyz3[1]-xyz2[1])*(xyz3[1]-xyz2[1]));     
+      Double_t c3 = 2.*loc3[1]/(loc3[0]*(loc3[0]-lr));   //curvature approximated using 2 derivation
+      //
+      while ( (cl0=layer0.GetNextCluster(ci0))!=0){ 
+       if (cl0->GetQ()==0) continue;
+       Double_t *xyz0 = &clusterxyzr0[ci0*6];  
+       if (cl3->GetLabel(0)==cl0->GetLabel(0) && cl3->GetLabel(0)==cl2->GetLabel(0) 
+           && cl3->GetLabel(0)==cl1->GetLabel(0)&&cl3->GetLabel(0)==cl0->GetLabel(0)){
+         good1++;
+       }       
+       Double_t loc0[2]; // cluster3 in local coordindate 
+       loc0[0] =  (xyz0[0]-xyz1[0])*cosphi+(xyz0[1]-xyz1[1])*sinphi;
+       loc0[1] = -(xyz0[0]-xyz1[0])*sinphi+(xyz0[1]-xyz1[1])*cosphi;
+       //
+       Double_t c0 = 2.*loc0[1]/(loc0[0]*(loc0[0]-lr));   //curvature approximated using 2 derivation
+       Double_t tan2 = (xyz1[2]-xyz0[2])/TMath::Sqrt((xyz1[0]-xyz0[0])*(xyz1[0]-xyz0[0])+(xyz1[1]-xyz0[1])*(xyz1[1]-xyz0[1]));
+       //if (TMath::Abs(tan1-tan2)>0.05) continue;
+       //if (TMath::Abs(c3-c0)>0.006) continue;
+       Double_t poolphi = (c3-c0)/0.0011; 
+       poolphi*=poolphi;
+       Double_t poolr = ((c3-c0)/(TMath::Abs(c3+c0)+0.001))/0.14;
+       Double_t poolth = (tan1-tan2)/0.0094; 
+       poolth*=poolth;
+       Double_t poolthr = ((tan1-tan2)/(TMath::Abs(c3+c0)+0.005))/0.83;
+       if (poolr*poolr+poolthr*poolthr>40) continue;
+       //if ( TMath::Abs((c3-c0)/(c3+c0))>0.6) continue; 
+       cindex2[ntracks2*4+0]   = ci0;
+       cindex2[ntracks2*4+1]   = index1;
+       cindex2[ntracks2*4+2]   = index2;
+       cindex2[ntracks2*4+3]   = index3;
+       pools[ntracks2] = poolr*poolr+poolthr*poolthr;
+       //
+       ntracks2++;
+       if (cl3->GetLabel(0)==cl0->GetLabel(0) && cl3->GetLabel(0)==cl2->GetLabel(0) 
+           && cl3->GetLabel(0)==cl1->GetLabel(0)){
+         hl->Fill(tan1-tan2);
+         hrc->Fill(poolr);
+         hc->Fill((c3-c0));
+         hlr->Fill(poolthr);
+         hpoolc->Fill(poolr);
+         hpooll->Fill(poolthr);
+         accepted1++;
+       }       
+      }
+    }
+    Int_t sindex[10000];
+    TMath::Sort(ntracks2,pools,sindex,kFALSE);
+    Float_t minpool = pools[sindex[0]];
+    Int_t max = TMath::Min(ntracks2,20);
+    for (Int_t itrack=0;itrack<ntracks2;itrack++){
+      if (itrack>=max) continue;
+      //if (pools[sindex[itrack]]>minpool+5) continue;
+      const AliITSclusterV2 *cl0, *cl1,*cl2,*cl3;
+      cl3 = layer3.GetCluster(cindex2[sindex[itrack]*4+3]);
+      cl2 = layer2.GetCluster(cindex2[sindex[itrack]*4+2]);
+      cl1 = layer1.GetCluster(cindex2[sindex[itrack]*4+1]);
+      cl0 = layer0.GetCluster(cindex2[sindex[itrack]*4+0]);
+      if (cl3->GetLabel(0)==cl0->GetLabel(0) && cl3->GetLabel(0)==cl2->GetLabel(0) 
+         && cl3->GetLabel(0)==cl1->GetLabel(0)){
+       accepted2++;
+      }
+      nfoundtracks++;
+    }
+  }
+  hc->Draw();
+  printf("First pass \t%d\t%d\t%d\n",ntracks0,good0,accepted0);  
+  printf("Found tracks\t%d\t%d\t%d\t%d\n",nfoundtracks,good1, accepted1,accepted2);
+  
+}
diff --git a/STEER/AliITSclusterComparison.h b/STEER/AliITSclusterComparison.h
new file mode 100644 (file)
index 0000000..c6e7af9
--- /dev/null
@@ -0,0 +1,195 @@
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include "AliTrackReference.h"
+#endif
+
+
+class AliITSCI: public TObject {
+  public :
+  AliITSCI(){;}
+  virtual ~AliITSCI(){;}
+  Float_t GetExpectedQ();
+  //
+  AliITSclusterV2 fCl;
+  AliTrackReference fRef;
+  TParticle fP;
+  Char_t fStatus;
+  //
+  Float_t fBeta2;
+  Float_t fVertex;
+  Int_t fNClusters;
+  Int_t   fLabPos;
+  Float_t fExpectedQ;
+  Float_t fNormQ;
+  Int_t   fPDG;
+  Int_t   fNTracks;
+  Float_t fPt;
+  Float_t fPtotal;
+  Float_t fR;
+  Float_t fCharge;
+  Bool_t  fIsPrim;
+  Double_t fDy;
+  Double_t fDz;
+  Double_t fPoolY;
+  Double_t fPoolY2;
+  Double_t fPoolZ;
+  Double_t fPoolZ2;
+  Double_t fErrY;
+  Double_t fErrZ;
+  Int_t    fErrType;
+  Double_t fPhiCl;
+  Double_t fRCl;
+  Double_t fLx;
+  Double_t fLy;
+  Double_t fGx;
+  Double_t fGy;
+  Float_t  fTNy;
+  Float_t  fTNz;
+  Float_t  fTheta;
+  Double_t fPhi;
+  Double_t fPhiL;
+  Int_t fLayer;
+  void Update(AliITStrackerMI*tracker);
+  ClassDef(AliITSCI,1)
+};
+
+class AliITSClusterErrAnal: public TObject{
+public: 
+  AliITSClusterErrAnal(Char_t *chloader  = "galice.root");
+  void SetIO(Int_t event);  
+  Int_t Analyze(Int_t trackmax);
+  void LoadClusters();
+  void LoadParticles();
+  void SignDeltas( TObjArray *ClusterArray, Float_t vz);
+  void SortReferences();
+  AliITSclusterV2 * FindNearestCluster(AliITSCI * clinfo, Float_t dmax=0.5);
+  void GetNTeor(Int_t layer, Float_t theta, Float_t phi, Float_t &ny, Float_t &nz);
+  Int_t GetError(Int_t layer, const AliITSclusterV2*cl, Float_t theta, Float_t phi, Float_t expQ, Float_t &erry, Float_t &errz);
+  void MakeSeeds(Double_t zv);
+public:
+  AliRunLoader * fRunLoader;
+  AliLoader * fITSLoader;
+  TTree * fHitTree;
+  TTree * fClusterTree;
+  TTree * fReferenceTree;
+  AliITS * fITS;
+  //
+  TTree * fTreeA;
+  TTree * fTreeB;
+  TFile * fFileA;
+  AliITStrackerMI *fTracker;
+  AliStack *fStack;
+  TObjArray *fClusters;  //array of clusters
+  TObjArray fExactPoints;
+  TObjArray *fReferences;
+  TObjArray *fParticles;
+  ClassDef(AliITSClusterErrAnal,1)
+};
+
+
+void AliITSCI::Update(AliITStrackerMI*tracker)
+{
+  //
+  //
+  fPt     = fRef.Pt();
+  fR      = fRef.R();
+  fCharge = fP.GetPDG()->Charge();
+  //fIsPrim = (fP.GetMother(0)>0)? kFALSE :kTRUE;
+  Double_t vertex = TMath::Sqrt( fP.Vx()*fP.Vx()
+                               +fP.Vy()*fP.Vy());
+  fVertex = vertex; 
+  fIsPrim = (vertex<1) ? kTRUE:kFALSE;
+  fTheta  = fRef.Pz()/fRef.Pt();
+  fPhi    =  TMath::ATan2(fRef.Py(),fRef.Px());  
+  fPDG    = fP.GetPdgCode();
+  fExpectedQ = GetExpectedQ();
+}
+
+
+Float_t AliITSCI::GetExpectedQ()
+{
+  //
+  //
+  //
+  //fPtotal = fP.P();
+  fPtotal = fRef.P();
+
+  Double_t p2= fPtotal*fPtotal;
+  Float_t mass = fP.GetPDG()->Mass();
+  Double_t beta2=p2/(p2 + mass*mass);
+  fBeta2 = beta2;
+  Float_t d=0.003;
+  //  Double_t dE=0.153e-3/beta2*(log(5940*beta2/(1-beta2)) - beta2)*d;
+  //dE+=0.33e-3/(beta2)*d;
+  //if (beta2/(1-beta2)>3.5*3.5)
+  //  dE=0.153e-3/beta2*(log(3.5*5940)+0.5*log(beta2/(1-beta2)) - beta2)*d;
+  Double_t dE = 0.22*0.153e-3*(39.2-55.6*beta2+28.7*beta2*beta2+27.41/beta2)*d;
+  
+  return dE*10000000.; //normalization  to ADC
+  //return dE; //normalization  to ADC
+}
+
+/*
+
+.L $ALICE_ROOT/COMPARISON/AliITSclusterComparison.C+ 
+.L $ALICE_ROOT/COMPARISON/AliGenInfo.C+ 
+AliITSClusterErrAnal anal;
+anal.SetIO(0); 
+anal.LoadParticles();
+anal.SortReferences();
+anal.LoadClusters();
+
+
+anal.Analyze(30000000);
+
+
+anal.MakeSeeds(-3.63); >seed1.txt
+
+
+TFile f("itsclusteranal.root");
+TTree * tree= (TTree*)f.Get("itscl")
+TTree * tree2= (TTree*)f.Get("Clusters")
+
+AliComparisonDraw comp;
+comp.fTree = tree;
+
+TH1F * hpools = new TH1F("pools","pools",100,-6,6); hpools->SetLineColor(4); hpools->SetFillColor(0);
+TH1F * hpools0 = new TH1F("pools0","pools0",100,-6,6);  hpools0->SetLineColor(2); hpools0->SetFillColor(0);
+
+TH1F * hsigma = new TH1F("hsigma","hsigma",1000,0,0.3)
+
+TCut cprim("cprim","fCl.fTracks[0]<110000&&fCl.fTracks[1]<110000&&fCl.fTracks[2]<110000");
+
+TCut cprim("cprim","fCl.fTracks[0]<24000&&fCl.fTracks[1]<24000&&fCl.fTracks[2]<24000");
+TCut cpt("cpt","sqrt(fRef.fPx**2+fRef.fPy**2)>0.05");
+TCut cphi("cphi","abs(fPhiL)<0.9");
+TCut ctheta("ctheta","abs(fTheta)<1.1");
+TCut cangle = cphi+ctheta;
+TCut cpools5("cpools5","abs(fPoolY)<5&&abs(fPoolZ)<5");
+TCut cpools52("cpools52","abs(fPoolY2)<5&&abs(fPoolZ2)<5");
+TCut cpools62("cpools52","sqrt(fPoolY2**2+fPoolZ2**2)<6");
+TCut cgold("cgold","abs(fCl.fNy+fCl.fNz-fTNy-fTNz)<1.5");
+TCut cback("cback","fRef.fX*fRef.fPx+fRef.fPy*fRef.fY<0");
+
+comp.fTree->Draw("fPoolY2>>pools","fStatus==1&&abs(fLayer-0.5)<0.6&&fIsPrim"+cpt+cangle+cpools62,"")
+
+comp.DrawXY("fErrType","fPoolZ2","fStatus==1&&abs(fLayer-4.5)<0.6","1",8,99.9,108,-5,5)
+
+comp.DrawXY("fCl.fType","fPoolY2","fStatus==1&&abs(fLayer-4.5)<0.6&&fIsPrim","1",4,102.9,106.1,-5,5)
+comp.DrawXY("abs(fCl.fChargeRatio)","fPoolY2","fStatus==1&&abs(fLayer-4.5)<0.6&&fIsPrim","1",4,0.0,1,-5,5)
+
+TProfile prof1("prof1","prof1",30,0.1,1); prof1.SetLineColor(2);
+TProfile prof2("prof2","prof2",30,0.1,1); prof2.SetLineColor(3);
+
+comp.fTree->Draw("fNormQ/fExpectedQ:fBeta2>>prof1","fStatus==1&&fCl.fLayer>1&&fNormQ/fExpectedQ<1.5&&fCl.fNz<4&&abs(fPDG)==211&&fIsPrim")
+comp.fTree->Draw("fNormQ/fExpectedQ:fBeta2>>prof2","fStatus==1&&fCl.fLayer>1&&fNormQ/fExpectedQ<1.5&&fCl.fNz<4&&abs(fPDG)!=211&&fIsPrim")
+
+prof1->Draw();prof2->Draw("same");
+
+
+TH1F * hdeltam   = new TH1F("hdeltam" ,"hdeltam",100,-3,10);  hdeltam->SetLineColor(4); hdeltam->SetFillColor(0);
+TH1F * hdelta0  =  new TH1F("hdelta0" ,"hdelta0",100,-3,10);  hdelta0->SetLineColor(2); hdelta0->SetFillColor(0);
+
+
+
+*/
diff --git a/STEER/AliTPCComparisonMI.C b/STEER/AliTPCComparisonMI.C
new file mode 100644 (file)
index 0000000..d604b06
--- /dev/null
@@ -0,0 +1,685 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+
+///////////////////////////////////////////////////////////////////////////////
+//                                                                           //
+//  Time Projection Chamber                                                  //
+//  Comparison macro for TPC                                                 //
+//  responsible: 
+//  marian.ivanov@cern.ch                                                    //
+//
+//
+/*
+//
+
+
+.L $ALICE_ROOT/STEER/AliGenInfo.C+
+.L $ALICE_ROOT/STEER/AliTPCComparisonMI.C+
+
+
+TPCCmpTr *t2 = new TPCCmpTr("genTracks.root","cmpTracks.root","galice.root",-1,1,0);
+t2->Exec();
+
+
+TCut cprim("cprim","MC.fVDist[3]<1");
+TCut csec("cprim","MC.fVDist[3]>1");
+TCut crec("crec","fReconstructed==1");
+TCut cteta1("cteta1","abs(MC.fTrackRef.Theta()/3.1415-0.5)<0.25");
+TCut cpos1("cpos1","abs(MC.fParticle.fVz/sqrt(MC.fParticle.fVx*MC.fParticle.fVx+MC.fParticle.fVy*MC.fParticle.fVy))<1");
+TCut csens("csens","abs(sqrt(fVDist[0]**2+fVDist[1]**2)-170)<50");
+TCut cmuon("cmuon","abs(MC.fParticle.fPdgCode==-13)");
+AliTPCComparisonDraw comp; 
+comp.SetIO();
+
+
+//example
+comp.DrawXY("fTPCinP0[3]","fTPCDelta[4]/fTPCinP1[3]","fReconstructed==1&&fPdg==-211"+cprim,"1",4,0.2,1.5,-0.06,0.06)
+comp.fRes->Draw();
+comp.fMean->Draw();  
+comp.Eff("fTPCinP0[3]","fRowsWithDigits>120"+cteta1+cpos1,"fReconstructed>0",10,0,1.5)
+comp.fRes->Draw();
+*/
+
+
+
+
+
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include <stdio.h>
+#include <string.h>
+//ROOT includes
+#include "Rtypes.h"
+#include "TFile.h"
+#include "TTree.h"
+#include "TChain.h"
+#include "TCut.h"
+#include "TString.h"
+#include "TBenchmark.h"
+#include "TStopwatch.h"
+#include "TParticle.h"
+#include "TSystem.h"
+#include "TTimer.h"
+#include "TVector3.h"
+#include "TPad.h"
+#include "TCanvas.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TF1.h"
+#include "TText.h"
+#include "Getline.h"
+#include "TStyle.h"
+
+//ALIROOT includes
+#include "AliRun.h"
+#include "AliStack.h"
+#include "AliTPCtrack.h"
+#include "AliSimDigits.h"
+#include "AliTPCParam.h"
+#include "AliTPC.h"
+#include "AliTPCLoader.h"
+#include "AliDetector.h"
+#include "AliTrackReference.h"
+#include "AliRun.h"
+#include "AliTPCParamSR.h"
+#include "AliTracker.h"
+#include "AliComplexCluster.h"
+#include "AliMagF.h"
+#endif
+#include "AliGenInfo.h"
+#include "AliTPCComparisonMI.h"
+
+//
+//
+  
+
+
+void AliTPCRecInfo::Update(AliMCInfo* info,AliTPCParam * param, Bool_t reconstructed, Int_t direction)
+{
+  //
+  //
+  //calculates derived variables
+  //
+  
+  if (reconstructed==kFALSE){
+    fReconstructed = kFALSE;
+    return;
+  }
+  fReconstructed = kTRUE;
+  // Find the nearest track reference
+  //
+  Double_t radius = TMath::Sqrt(fTrack.GetX()*fTrack.GetX()+fTrack.GetY()*fTrack.GetY());
+  TClonesArray * tpcreferences = info->fTPCReferences;
+  Int_t nentries = tpcreferences->GetEntriesFast();
+  AliTrackReference * ref = 0;
+  Double_t dr = 1000000;
+  for (Int_t i=0;i<nentries;i++){
+    AliTrackReference * refi = (AliTrackReference*)tpcreferences->UncheckedAt(i);
+    if (!refi) continue;
+    //if  ( (direction<0) && (refi->R()<radius-10)) continue;
+    //if  ( (direction>0) && (refi->R()>radius+10)) continue;    
+    if  (TMath::Abs(refi->R()-radius)<dr){
+      dr = TMath::Abs(refi->R()-radius);
+      ref = refi;
+    }
+  }  
+  if (ref==0) {
+    fReconstructed = kFALSE;
+    return;
+  }
+  //
+  //
+  fdEdx = fTrack.GetdEdx();
+  //
+  // propagate track to nearest reference in direction
+  TVector3 local = TPCCmpTr::TR2Local(ref,param);
+  local.GetXYZ(fTRLocalCoord);
+  //
+  // rotate if neccessary
+  Float_t ymax = local.X()*TMath::Tan(0.5*param->GetInnerAngle());  
+  Float_t y = fTrack.GetYat(local.X());
+  if (y > ymax) {
+    fTrack.Rotate(param->GetInnerAngle());
+  } else if (y <-ymax) {
+    fTrack.Rotate(-param->GetInnerAngle());
+  }   
+  Double_t xtrack = fTrack.GetX();
+  Double_t delta = (local.X()-xtrack)*0.1;
+  for (Int_t i=0;i<9;i++){
+    fTrack.PropagateTo(xtrack+delta*float(i));
+  }
+  fTrack.PropagateTo(local.X());
+
+  Double_t par[5], cov[15], localX;
+  fTrack.GetExternalParameters(localX,par);
+  fTrack.GetExternalCovariance(cov);
+  //
+  //
+  fRecPhi=TMath::ASin(par[2]) + fTrack.GetAlpha();
+  Double_t phi2 =TMath::ATan2(ref->Py(),ref->Px());
+  if (phi2<0) phi2+=2*TMath::Pi();
+  fGenPhi =phi2;
+  //
+  if (fRecPhi<0) fRecPhi+=2*TMath::Pi();
+  if (fRecPhi>=2*TMath::Pi()) fRecPhi-=2*TMath::Pi();
+  fLambda = TMath::ATan(par[3]);
+  fRecPt_1 = TMath::Abs(par[4]);  
+  //TPC
+  //
+  //
+  Double_t phi=TMath::ATan2(par[0],localX) + fTrack.GetAlpha();
+  Double_t r=TMath::Sqrt(localX*localX + par[0]*par[0]);
+  fTPCinR1[0]= r*TMath::Cos(phi); 
+  fTPCinR1[1]= r*TMath::Sin(phi); 
+  fTPCinR1[2]= par[1];
+  fTPCinR1[3] = TMath::Sqrt(fTPCinR1[0]*fTPCinR1[0]+fTPCinR1[1]*fTPCinR1[1]);
+  fTPCinR1[4] = TMath::ATan2(fTPCinR1[1],fTPCinR1[0]); 
+  //
+  fTPCinP1[0] = fTrack.Px();
+  fTPCinP1[1] = fTrack.Py();
+  fTPCinP1[2] = fTrack.Pz();
+  fTPCinP1[3] = TMath::Sqrt(fTPCinP1[0]*fTPCinP1[0]+fTPCinP1[1]*fTPCinP1[1]);
+  //
+  fTPCinR0[0] = ref->X();      
+  fTPCinR0[1] = ref->Y();      
+  fTPCinR0[2] = ref->Z();
+  fTPCinR0[3] = TMath::Sqrt(fTPCinR0[0]*fTPCinR0[0]+fTPCinR0[1]*fTPCinR0[1]);
+  fTPCinR0[4] = TMath::ATan2(fTPCinR0[1],fTPCinR0[0]);
+  //
+  fTPCinP0[0]=ref->Px();
+  fTPCinP0[1]=ref->Py();
+  fTPCinP0[2]=ref->Pz();
+  fTPCinP0[3]=ref->Pt();
+  //
+  //
+  if (fTPCinP1[3]>0.0000001){
+    fTPCAngle1[0] = TMath::ATan2(fTPCinP1[1],fTPCinP1[0]);
+    fTPCAngle1[1] = TMath::ATan(fTPCinP1[2]/fTPCinP1[3]);  
+    //
+    fTPCAngle0[0] = TMath::ATan2(fTPCinP0[1],fTPCinP0[0]);
+    fTPCAngle0[1] = TMath::ATan(fTPCinP0[2]/fTPCinP0[3]);
+  }    
+  //
+  fTPCDelta[0] = (fTPCinR0[4]-fTPCinR1[4])*fTPCinR1[3];  //delta rfi
+  fTPCPools[0] = fTPCDelta[0]/TMath::Sqrt(cov[0]);
+  fTPCDelta[1] = (fTPCinR0[2]-fTPCinR1[2]);              //delta z
+  fTPCPools[1] = fTPCDelta[1]/TMath::Sqrt(cov[2]);
+  //
+  fTPCDelta[2] = (fTPCAngle0[0]-fTPCAngle1[0]);
+  fTPCPools[2] = fTPCDelta[2]/TMath::Sqrt(cov[5]);
+  //
+  fTPCDelta[3] = TMath::Tan(fTPCAngle0[1])-TMath::Tan(fTPCAngle1[1]);
+  fTPCPools[3] = fTPCDelta[3]/TMath::Sqrt(cov[9]);
+  fTPCDelta[4] = (fTPCinP0[3]-fTPCinP1[3]);
+  Double_t sign = (fTrack.GetC()>0) ? 1:-1;
+  fTPCPools[4] = sign*(1./fTPCinP0[3]-1./fTPCinP1[3])/TMath::Sqrt(cov[14]);      
+  
+}
+
+
+
+
+////////////////////////////////////////////////////////////////////////
+TPCCmpTr::TPCCmpTr()
+{
+  Reset();
+}
+
+////////////////////////////////////////////////////////////////////////
+TPCCmpTr::TPCCmpTr(const char* fnGenTracks,
+                  const char* fnCmp,
+                  const char* fnGalice, Int_t direction,
+                  Int_t nEvents, Int_t firstEvent)
+{
+  Reset();
+  //  fFnGenTracks = fnGenTracks;
+  //  fFnCmp = fnCmp;
+  sprintf(fFnGenTracks,"%s",fnGenTracks);
+  sprintf(fFnCmp,"%s",fnCmp);
+
+  fFirstEventNr = firstEvent;
+  fEventNr = firstEvent;
+  fNEvents = nEvents;
+  fDirection = direction;
+  //
+  fLoader = AliRunLoader::Open(fnGalice);
+  if (gAlice){
+    //delete gAlice->GetRunLoader();
+    delete gAlice;
+    gAlice = 0x0;
+  }
+  if (fLoader->LoadgAlice()){
+    cerr<<"Error occured while l"<<endl;
+  }
+  Int_t nall = fLoader->GetNumberOfEvents();
+  if (nEvents==0) {
+    nEvents =nall;
+    fNEvents=nall;
+    fFirstEventNr=0;
+  }    
+
+  if (nall<=0){
+    cerr<<"no events available"<<endl;
+    fEventNr = 0;
+    return;
+  }
+  if (firstEvent+nEvents>nall) {
+    fEventNr = nall-firstEvent;
+    cerr<<"restricted number of events availaible"<<endl;
+  }
+  AliMagF * magf = gAlice->Field();
+  AliTracker::SetFieldMap(magf);
+
+}
+
+
+////////////////////////////////////////////////////////////////////////
+TPCCmpTr::~TPCCmpTr()
+{
+  //if (fLoader) {
+  //  delete fLoader;
+  //}
+}
+
+//////////////////////////////////////////////////////////////
+Int_t TPCCmpTr::SetIO()
+{
+  //
+  // 
+  CreateTreeCmp();
+  if (!fTreeCmp) return 1;
+  fParamTPC = GetTPCParam();
+  //
+  if (!ConnectGenTree()) {
+    cerr<<"Cannot connect tree with generated tracks"<<endl;
+    return 1;
+  }
+  return 0;
+}
+
+//////////////////////////////////////////////////////////////
+
+Int_t TPCCmpTr::SetIO(Int_t eventNr)
+{
+  //
+  // 
+  // SET INPUT
+  //gAlice->GetEvent(eventNr);
+  //  fLoader->SetEventNumber(eventNr);  
+  //fLoader = AliRunLoader::Open("galice.root");
+  fLoader->GetEvent(eventNr);  
+  //
+  AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
+  tpcl->LoadTracks();
+  fTreeRecTracks = tpcl->TreeT();
+  
+  //
+  return 0;
+}
+
+
+
+////////////////////////////////////////////////////////////////////////
+void TPCCmpTr::Reset()
+{
+  fEventNr = 0;
+  fNEvents = 0;
+  fTreeCmp = 0;
+  //
+  fFileGenTracks = 0;
+  fDebug = 0;
+  //
+  fParamTPC = 0;
+  fTreeRecTracks = 0;
+  fTreePoints =0;
+
+  fTPCTrack = 0; 
+  fTracks   = 0;
+  fTrackPoints =0;
+}
+
+////////////////////////////////////////////////////////////////////////
+Int_t TPCCmpTr::Exec(Int_t nEvents, Int_t firstEventNr)
+{
+  fNEvents = nEvents;
+  fFirstEventNr = firstEventNr;
+  return Exec();
+}
+
+////////////////////////////////////////////////////////////////////////
+Int_t TPCCmpTr::Exec()
+{
+  TStopwatch timer;
+  timer.Start();
+
+  if (SetIO()==1) 
+    return 1;
+   
+  cerr<<"fFirstEventNr, fNEvents: "<<fFirstEventNr<<" "<<fNEvents<<endl;
+  for (Int_t eventNr = fFirstEventNr; eventNr < fFirstEventNr+fNEvents;
+       eventNr++) {
+    fEventNr=eventNr;
+    SetIO(fEventNr);
+    fNParticles = gAlice->GetEvent(fEventNr);    
+
+    fIndexRecTracks = new Int_t[fNParticles*4];  //write at maximum 4 tracks corresponding to particle
+    fFakeRecTracks = new Int_t[fNParticles];
+    fMultiRecTracks = new Int_t[fNParticles];
+    for (Int_t i = 0; i<fNParticles; i++) {
+      fIndexRecTracks[4*i] = -1;
+      fIndexRecTracks[4*i+1] = -1;
+      fIndexRecTracks[4*i+2] = -1;
+      fIndexRecTracks[4*i+3] = -1;
+
+      fFakeRecTracks[i] = 0;
+      fMultiRecTracks[i] = 0;
+    }
+  
+    cout<<"Start to process event "<<fEventNr<<endl;
+    cout<<"\tfNParticles = "<<fNParticles<<endl;
+    if (fDebug>2) cout<<"\tStart loop over TreeT"<<endl;
+    if (TreeTLoop(eventNr)>0) return 1;
+
+    if (fDebug>2) cout<<"\tStart loop over tree genTracks"<<endl;
+    if (TreeGenLoop(eventNr)>0) return 1;
+    if (fDebug>2) cout<<"\tEnd loop over tree genTracks"<<endl;
+
+    delete [] fIndexRecTracks;
+    delete [] fFakeRecTracks;
+    delete [] fMultiRecTracks;
+  }
+
+  CloseOutputFile();
+
+  cerr<<"Exec finished"<<endl;
+  timer.Stop();
+  timer.Print();
+  return 0;
+
+}
+////////////////////////////////////////////////////////////////////////
+Bool_t TPCCmpTr::ConnectGenTree()
+{
+//
+// connect all branches from the genTracksTree
+// use the same variables as for the new cmp tree, it may work
+//
+  fFileGenTracks = TFile::Open(fFnGenTracks,"READ");
+  if (!fFileGenTracks) {
+    cerr<<"Error in ConnectGenTree: cannot open file "<<fFnGenTracks<<endl;
+    return kFALSE;
+  }
+  fTreeGenTracks = (TTree*)fFileGenTracks->Get("genTracksTree");
+  if (!fTreeGenTracks) {
+    cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file "
+       <<fFnGenTracks<<endl;
+    return kFALSE;
+  }
+  //
+  fMCInfo = new AliMCInfo;
+  fTreeGenTracks->SetBranchAddress("MC",&fMCInfo);
+
+
+  if (fDebug > 1) {
+    cout<<"Number of gen. tracks with TR: "<<fTreeGenTracks->GetEntries()<<endl;
+  }
+  return kTRUE;
+}
+
+
+////////////////////////////////////////////////////////////////////////
+void TPCCmpTr::CreateTreeCmp() 
+{
+  fFileCmp = TFile::Open(fFnCmp,"RECREATE");
+  if (!fFileCmp) {
+    cerr<<"Error in CreateTreeCmp: cannot open file "<<fFnCmp<<endl;
+    return;
+  }
+
+
+  fTreeCmp    = new TTree("TPCcmpTracks","TPCcmpTracks");
+  //
+  fMCInfo = new AliMCInfo;
+  fRecInfo = new AliTPCRecInfo;
+  //
+  fTPCTrack = new AliTPCtrack;
+   //
+  fTreeCmp->Branch("MC","AliMCInfo",&fMCInfo);
+  fTreeCmp->Branch("RC","AliTPCRecInfo",&fRecInfo);
+  fTreeCmp->Branch("fTPCTrack","AliTPCtrack",&fTPCTrack);
+  //
+  fTreeCmp->AutoSave(); 
+
+}
+////////////////////////////////////////////////////////////////////////
+void TPCCmpTr::CloseOutputFile()  
+{
+  if (!fFileCmp) {
+    cerr<<"File "<<fFnCmp<<" not found as an open file."<<endl;
+    return;
+  }
+  fFileCmp->cd();
+  fTreeCmp->Write();    
+  delete fTreeCmp;
+  
+  fFileCmp->Close();
+  delete fFileCmp;
+  return;
+}
+////////////////////////////////////////////////////////////////////////
+
+TVector3 TPCCmpTr::TR2Local(AliTrackReference *trackRef,
+                           AliTPCParam *paramTPC) {
+
+  Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
+  Int_t index[4];
+  paramTPC->Transform0to1(x,index);
+  paramTPC->Transform1to2(x,index);
+  return TVector3(x);
+}
+////////////////////////////////////////////////////////////////////////
+
+Int_t TPCCmpTr::TreeTLoop(Int_t eventNr)
+{
+  //
+  // loop over all TPC reconstructed tracks and store info in memory
+  //
+  TStopwatch  timer;
+  timer.Start();
+  
+  if (!fTreeRecTracks) {
+    cerr<<"Can't get a tree with TPC rec. tracks  "<<endl;
+    return 1;
+  }
+  //fTreePoints=(TTree*)fFileRecTracks->Get("trackDebug");
+  
+  Int_t nEntries = (Int_t) fTreeRecTracks->GetEntries();
+  if (fDebug > 2) cout<<"Event, rec. tracks: "<<eventNr<<" "
+                     <<nEntries<<endl;
+  TBranch * br= fTreeRecTracks->GetBranch("tracks");
+  br->SetAddress(&fTPCTrack);
+  TBranch *brp = 0;
+  if (fTreePoints) brp = fTreePoints->GetBranch("debug");
+
+  if (fTracks){
+    fTracks->Delete();    
+    delete fTracks;
+  }
+  if (fTrackPoints){
+    fTrackPoints->Delete();
+    delete fTrackPoints;
+    fTrackPoints =0;
+  }
+  fTracks      = new TObjArray(nEntries);
+  if (brp){
+    fTrackPoints = new TObjArray(nEntries);
+  }
+  else fTrackPoints = 0;
+
+  //
+  //load tracks to the memory
+  for (Int_t i=0; i<nEntries;i++){
+    AliTPCtrack * track = new AliTPCtrack;
+    br->SetAddress(&track);
+    br->GetEntry(i);
+    fTracks->AddAt(track,i);
+  }
+  //
+  //load track points to the memory
+  if (brp) for (Int_t i=0; i<nEntries;i++){
+    TClonesArray * arr = new TClonesArray("AliTPCTrackPoint2");
+    brp->SetAddress(&arr);
+    brp->GetEntry(i);
+    if (arr!=0)
+      for (Int_t j=0;j<arr->GetEntriesFast();j++){
+       AliTPCTrackPoint2 * point = (AliTPCTrackPoint2*)arr->UncheckedAt(j);
+       if (point && point->fID>=0){
+         fTrackPoints->AddAt(arr,point->fID);
+         break;
+       }
+      }    
+  }
+  //
+
+  //
+  for (Int_t iEntry=0; iEntry<nEntries;iEntry++){
+    //br->GetEntry(iEntry);
+    fTPCTrack = (AliTPCtrack*)fTracks->UncheckedAt(iEntry);
+    //
+    Int_t label = fTPCTrack->GetLabel();
+    Int_t absLabel = abs(label);
+    if (absLabel < fNParticles) {
+      //      fIndexRecTracks[absLabel] =  iEntry;
+      if (label < 0) fFakeRecTracks[absLabel]++;
+      
+      if (fMultiRecTracks[absLabel]>0){
+       if (fMultiRecTracks[absLabel]<4)
+         fIndexRecTracks[absLabel*4+fMultiRecTracks[absLabel]] =  iEntry;      
+      }
+      else      
+       fIndexRecTracks[absLabel*4] =  iEntry;
+      fMultiRecTracks[absLabel]++;
+    }
+  }  
+  printf("Time spended in TreeTLoop\n");
+  timer.Print();
+  
+  if (fDebug > 2) cerr<<"end of TreeTLoop"<<endl;
+
+  return 0;
+}
+
+////////////////////////////////////////////////////////////////////////
+Int_t TPCCmpTr::TreeGenLoop(Int_t eventNr)
+{
+//
+// loop over all entries for a given event, find corresponding 
+// rec. track and store in the fTreeCmp
+//
+  TStopwatch timer;
+  timer.Start();
+  Int_t entry = 0;
+  Double_t nParticlesTR = fTreeGenTracks->GetEntriesFast();
+  TBranch * branch = fTreeCmp->GetBranch("RC");
+  branch->SetAddress(&fRecInfo); // set all pointers
+
+  while (entry < nParticlesTR) {
+    fTreeGenTracks->GetEntry(entry);
+    entry++;
+    if (eventNr < fMCInfo->fEventNr) continue;
+    if (eventNr > fMCInfo->fEventNr) continue;
+    //
+    if (fDebug > 2 && fMCInfo->fLabel < 10) {
+      cerr<<"Fill track with a label "<<fMCInfo->fLabel<<endl;
+    }
+    fRecInfo->Reset();
+    //
+    if (fIndexRecTracks[fMCInfo->fLabel*4] >= 0) {
+      fTPCTrack= (AliTPCtrack*)fTracks->UncheckedAt(fIndexRecTracks[fMCInfo->fLabel*4]);
+      
+      //      if (nBytes > 0) {
+      if (fTPCTrack) {
+       //
+       //
+       //if multifound find track with biggest pt - returning tracks - loos energy
+       //
+       if (fIndexRecTracks[fMCInfo->fLabel*4]+1){
+         Double_t pt = TMath::Abs(1/fTPCTrack->Get1Pt());
+         for (Int_t i=1;i<4;i++){
+           if (fIndexRecTracks[fMCInfo->fLabel*4+i]>=0){
+             AliTPCtrack * track = (AliTPCtrack*)fTracks->UncheckedAt(fIndexRecTracks[fMCInfo->fLabel*4+i]);
+             Double_t pt2 = TMath::Abs(1/track->Get1Pt());           
+             if  (pt2>pt)
+               fTPCTrack = track;                 
+           }
+         }
+       }
+       //
+       fRecInfo->fTP=0;
+       fRecInfo->fTrack =*fTPCTrack; 
+       fRecInfo->fReconstructed = 1;
+       fRecInfo->fFake = fFakeRecTracks[fMCInfo->fLabel];
+       fRecInfo->fMultiple = fMultiRecTracks[fMCInfo->fLabel];
+       fRecInfo->Update(fMCInfo, fParamTPC, kTRUE, 1);
+      }      
+    } 
+    else{
+      fRecInfo->Update(fMCInfo, fParamTPC, kFALSE, 1);      
+    }    
+    fTreeCmp->Fill();
+  }
+  fTreeCmp->AutoSave();
+  fTracks->Delete();
+  fTPCTrack =0;
+  if (fTrackPoints){
+    fTrackPoints->Delete();
+    delete fTrackPoints;
+    fTrackPoints =0;
+  } 
+  printf("Time spended in TreeGenLoop\n");
+  timer.Print();
+  if (fDebug > 2) cerr<<"end of TreeGenLoop"<<endl;  
+  return 0;
+}
+
+
+////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////
+
+void AliTPCComparisonDraw::SetIO(const char *fname)
+{
+  TFile *file = TFile::Open(fname);
+  if (!file){
+    printf("Comparison file couldn't be generated\n"); 
+    return;
+   }
+  //
+  fTree = (TTree*) file->Get("TPCcmpTracks");
+  if (!fTree) {
+    printf("no track comparison tree found\n");
+    file->Close();
+    delete file;
+  }
+
+}
+
+void AliTPCComparisonDraw::Eff()
+{
+  
+}
+
+void AliTPCComparisonDraw::ResPt()
+{
+}
diff --git a/STEER/AliTPCComparisonMI.h b/STEER/AliTPCComparisonMI.h
new file mode 100644 (file)
index 0000000..1a8a71a
--- /dev/null
@@ -0,0 +1,333 @@
+
+
+////////////////////////////////////////////////////////////////////////
+//
+// Start of implementation of the class digitRow
+//
+////////////////////////////////////////////////////////////////////////
+const Int_t kgRowBytes = 32;
+
+class digitRow: public TObject {
+
+public:
+  digitRow();
+  virtual ~digitRow(){;}
+  void SetRow(Int_t row);
+  Bool_t TestRow(Int_t row);
+  digitRow & operator=(const digitRow &digOld);
+  Int_t RowsOn(Int_t upto=8*kgRowBytes);
+  Int_t Last();
+  Int_t First();
+  void Reset();
+
+//private:
+  UChar_t fDig[kgRowBytes];
+
+  ClassDef(digitRow,1)  // container for digit pattern
+};
+ClassImp(digitRow)
+
+
+////////////////////////////////////////////////////////////////////////
+//
+// Start of implementation of the class AliTPCGenInfo
+//
+////////////////////////////////////////////////////////////////////////
+
+class AliTPCGenInfo: public TObject {
+
+public:
+  AliTPCGenInfo();
+  ~AliTPCGenInfo();
+
+  AliTrackReference fTrackRef;      // track reference saved in the output tree
+  AliTrackReference fTrackRefOut;   // decay track reference saved in the output tree
+  TParticle fParticle;           // generated particle 
+  Int_t fLabel;                   // track label
+  Int_t fEventNr;                 // event number
+
+  Float_t fDecayCoord[3];         // position of particle decay
+  Double_t fVDist[4];             //distance of the particle vertex from primary vertex
+  
+  Int_t fRowsWithDigitsInn;    // number of rows with digits in the inner sectors
+  Int_t fRowsWithDigits;       // number of rows with digits in the outer sectors
+  Int_t fRowsTrackLength;      // last - first row with digit
+  Int_t fDigitsInSeed;         // digits in the default seed rows
+  Float_t fPrim;               // theoretical dedx in tpc according particle momenta and mass
+  digitRow fRow;               // information about digits row pattern
+  TClonesArray * fReferences;  //containner with all track references
+  ClassDef(AliTPCGenInfo,1)  // container for 
+};
+ClassImp(AliTPCGenInfo)
+
+
+
+class AliTPCGenV0Info: public TObject {
+public:
+  AliTPCGenInfo fMCd;      //info about daughter particle
+  AliTPCGenInfo fMCm;      //info about mother particle
+  void Update();        // put some derived info to special field 
+  Double_t    fDist1;    //info about closest distance according closest MC - linear DCA
+  Double_t    fDist2;    //info about closest distance parabolic DCA
+  //
+  Double_t     fPdr[3];    //momentum at vertex daughter  - according approx at DCA
+  Double_t     fPd[4];     //exact momentum from MC info
+  Double_t     fX[3];      //exact position of the vertex
+  Double_t     fXr[3];     //rec. position according helix
+  //
+  Double_t     fPm[3];    //momentum at the vertex mother
+  Double_t     fAngle[3]; //three angels
+  Double_t     fRr;       // rec position of the vertex 
+  Double_t     fR;        //exact r position of the vertex
+  Int_t        fPdg[2];   //pdg code of mother and daugter particles
+  Int_t        fLab[2];   //MC label of the partecle
+  ClassDef(AliTPCGenV0Info,1)  // container for  
+};
+ClassImp(AliTPCGenV0Info)
+
+
+
+void AliTPCGenV0Info::Update()
+{
+  fPd[0] = fMCd.fParticle.Px();
+  fPd[1] = fMCd.fParticle.Py();
+  fPd[2] = fMCd.fParticle.Pz();
+  fPd[3] = fMCd.fParticle.P();
+  fX[0]  = fMCd.fParticle.Vx();
+  fX[1]  = fMCd.fParticle.Vy();
+  fX[2]  = fMCd.fParticle.Vz();
+  fR     = TMath::Sqrt( fX[0]*fX[0]+
+                                fX[1]*fX[1]);
+  fPdg[0]    = fMCd.fParticle.GetPdgCode();
+  fPdg[1]    = fMCm.fParticle.GetPdgCode();
+  //
+  fLab[0]    = fMCd.fParticle.GetUniqueID();
+  fLab[1]    = fMCm.fParticle.GetUniqueID();
+
+}
+
+
+
+
+     /////////////////////////////////////////////////////////////////////////
+class AliTPCRecInfo: public TObject {
+
+public:
+  AliTPCRecInfo(){fTP = new TClonesArray("AliTPCTrackPoint2",0);}
+  ~AliTPCRecInfo(){if (fTP) {fTP->Delete();delete fTP;}}
+  //
+  AliTPCtrack fTPCTrack;          // tpc track
+  Float_t fTRLocalCoord[3];       //local coordinates of the track ref.
+  Int_t   fReconstructed;         //flag if track was reconstructed
+  Double_t fRecPhi;         // reconstructed phi angle (0;2*kPI)
+  Double_t fLambda;         // reconstructed 
+  Double_t fRecPt_1;        // reconstructed 
+  Float_t fdEdx;           // reconstructed  dEdx      
+  Int_t fFake;             // fake track
+  Int_t fMultiple;         // number of reconstructions
+  TClonesArray *fTP;        //container with track  points 
+  void Reset();
+  //
+  ClassDef(AliTPCRecInfo,1)  // container for 
+};
+ClassImp(AliTPCRecInfo)
+
+void AliTPCRecInfo::Reset()
+{
+  fMultiple =0; 
+  fFake     =0;
+  fReconstructed=0;
+  fRecPhi =0;
+  fLambda =0;
+}
+
+
+/////////////////////////////////////////////////////////
+/////////////////////////////////////////////////////////
+/////////////////////////////////////////////////////////
+
+
+class AliTPCRecV0Info: public TObject {
+public:
+  AliTPCRecInfo  fT1;      //track1
+  AliTPCRecInfo  fT2;      //track2  
+  Double_t    fDist1;    //info about closest distance according closest MC - linear DCA
+  Double_t    fDist2;    //info about closest distance parabolic DCA
+  //
+  Double_t     fPdr[3];    //momentum at vertex daughter  - according approx at DCA
+  Double_t     fXr[3];     //rec. position according helix
+  //
+  Double_t     fPm[3];    //momentum at the vertex mother
+  Double_t     fAngle[3]; //three angles
+  Double_t     fRr;       // rec position of the vertex 
+  Int_t        fLab[2];   //MC label of the partecle
+  ClassDef(AliTPCRecV0Info,1)  // container for  
+};
+
+ClassImp(AliTPCRecV0Info)
+
+
+
+  
+
+
+////////////////////////////////////////////////////////////////////////
+// 
+// Start of implementation of the class TPCFindGenTracks
+//
+////////////////////////////////////////////////////////////////////////
+
+class TPCFindGenTracks {
+
+public:
+  TPCFindGenTracks();
+  TPCFindGenTracks(char* fnHits,
+                  char* fnDigits ="tpc.digits.root",
+                  char* fnRes    ="genTracks.root",
+                  Int_t nEvents=1, Int_t firstEvent=0);
+  virtual ~TPCFindGenTracks();
+  void Reset();
+  Int_t Exec();
+  Int_t Exec(Int_t nEvents, Int_t firstEventNr);
+  void CreateTreeGenTracks();
+  void CloseOutputFile();
+  Int_t TreeKLoop();
+  Int_t TreeTRLoop();
+  Int_t TreeDLoop();
+  //void  FillInfo(Int_t iParticle);
+  void SetFirstEventNr(Int_t i) {fFirstEventNr = i;}
+  void SetNEvents(Int_t i) {fNEvents = i;}
+  void SetDebug(Int_t level) {fDebug = level;}
+  Int_t SetIO(Int_t eventNr);
+  Int_t SetIO();
+  Float_t TR2LocalX(AliTrackReference *trackRef,
+                   AliTPCParam *paramTPC);
+
+public:
+  AliTPCGenInfo*  fMCInfo;           //! information writen per particle
+  Int_t fDebug;                   //! debug flag  
+  Int_t fEventNr;                 //! current event number
+  Int_t fLabel;                   //! track label
+  Int_t fNEvents;                 //! number of events to process
+  Int_t fFirstEventNr;            //! first event to process
+  Int_t fNParticles;              //! number of particles in TreeK
+  TTree *fTreeGenTracks;          //! output tree with generated tracks
+  char *fFnRes;                   //! output file name with stored tracks
+  char *fFnHits;                  //! input file name with hits
+  char *fFnDigits;                //! input file name with digits
+  TFile *fFileGenTracks;          //! output file with stored fTreeGenTracks
+  TFile *fFileHits;               //! input file with hits
+  TFile *fFileTreeD;              //! input file with digits
+  //
+  TTree * fTreeD;                 //! current tree with digits
+  TTree * fTreeTR;                //! current tree with TR
+  AliStack *fStack;               //! current stack
+  //
+  digitRow *fContainerDigitRow;   //! big container for partial information
+  //
+  AliTrackReference *fReferences; //! container with track references
+  Int_t *fReferenceIndex0;        //! first index for given track
+  Int_t *fReferenceIndex1;        //! last  index for given track
+  //
+  AliTPCParam* fParamTPC;         //! AliTPCParam
+  Double_t fVPrim[3];             //! primary vertex position
+                                  // the fVDist[3] contains size of the 3-vector
+
+private:
+
+// some constants for the original non-pareller tracking (by Y.Belikov)
+  static const Int_t seedRow11 = 158;  // nRowUp - 1
+  static const Int_t seedRow12 = 139;  // nRowUp - 1 - (Int_t) 0.125*nRowUp
+  static const Int_t seedRow21 = 149;  // seedRow11 - shift
+  static const Int_t seedRow22 = 130;  // seedRow12 - shift
+  static const Double_t kRaddeg = 180./kPI;
+
+  static const Int_t fgMaxIndexTR = 50000; // maximum number of tracks with a track ref
+  static const Int_t fgMaxTR = 1000000; // maximum number of  track refs
+
+  static const Int_t fgMaxParticles = 2000000; // maximum number of generated particles
+  static const Double_t fgPtCut = .1; // do not store particles with generated pT less than this
+  static const Float_t fgTrackRefLocalXMax = 82.95;
+  static const Float_t fgTrackRefLocalXMaxDelta = 5.;
+
+  ClassDef(TPCFindGenTracks,1)    // class which creates and fills tree with TPCGenTrack objects
+};
+ClassImp(TPCFindGenTracks)
+
+
+
+////////////////////////////////////////////////////////////////////////
+// 
+// Start of implementation of the class TPCCmpTr
+//
+////////////////////////////////////////////////////////////////////////
+
+class TPCCmpTr {
+
+public:
+  TPCCmpTr();
+  TPCCmpTr(char* fnRecTracks,
+          char* fnGenTracks   ="genTracks.root",
+          char* fnCmpRes      ="cmpTracks.root", 
+          char* fnGalice      ="galice.root",
+          Int_t nEvents=1, Int_t firstEvent=0);
+  virtual ~TPCCmpTr();
+  void Reset();
+  Int_t Exec();
+  Int_t Exec(Int_t nEvents, Int_t firstEventNr);
+  void CreateTreeCmp();
+  void CloseOutputFile();
+  Bool_t ConnectGenTree();
+  Int_t TreeGenLoop(Int_t eventNr);
+  Int_t TreeTLoop(Int_t eventNr);
+  void SetFirstEventNr(Int_t i) {fFirstEventNr = i;}
+  void SetNEvents(Int_t i) {fNEvents = i;}
+  void SetDebug(Int_t level) {fDebug = level;}
+
+// tmp method, should go to TrackReferenceTPC
+  TVector3 TR2Local(AliTrackReference *trackRef,
+                   AliTPCParam *paramTPC);
+
+private:
+
+  Int_t fEventNr;                 //! current event number
+  Int_t fNEvents;                 //! number of events to process
+  Int_t fFirstEventNr;            //! first event to process
+  //
+  char *fFnCmp;                   //! output file name with cmp tracks
+  TFile *fFileCmp;                //! output file with cmp tracks
+  TTree *fTreeCmp;                //! output tree with cmp tracks
+  //
+  char *fFnGenTracks;             //! input file name with gen tracks
+  TFile *fFileGenTracks;
+  TTree *fTreeGenTracks;
+  //
+  char *fFnHits;                  //! input file name with gAlice object (needed for B)
+  TFile *fFileHits;               //! input file with gAlice
+  //
+  char *fFnRecTracks;             //! input file name with tpc rec. tracks
+  TFile *fFileRecTracks;          //! input file with reconstructed tracks
+  TTree *fTreeRecTracks;          //! tree with reconstructed tracks
+  TTree *fTreePoints;             //! tree with reconstructed points
+  //
+  Int_t *fIndexRecTracks;         //! index of particle label in the TreeT_TPC
+  Int_t *fFakeRecTracks;          //! number of fake tracks
+  Int_t *fMultiRecTracks;         //! number of multiple reconstructions
+  //
+  TObjArray * fTracks;            //!container with tracks 
+  TObjArray * fTrackPoints;       //! container with track points
+  //
+  AliTPCParam* fParamTPC;         //! AliTPCParam
+  Int_t fNParticles;              //! number of particles in the input tree genTracks
+  Int_t fDebug;                   //! debug flag  
+  Int_t fNextTreeGenEntryToRead;    //! last entry already read from genTracks tree
+  //
+  AliTPCGenInfo*  fMCInfo;           //! MC information writen per particle
+  AliTPCRecInfo*  fRecInfo;          //! Rec. information writen per particle
+  //
+  AliTPCtrack *fTPCTrack;         //! pointer to TPC track to connect branch
+
+  ClassDef(TPCCmpTr,1)    // class which creates and fills tree with TPCGenTrack objects
+};
+ClassImp(TPCCmpTr)
+