First version of some classes for comparison and performance studies (Marian)
authorhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 18 May 2007 16:45:01 +0000 (16:45 +0000)
committerhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 18 May 2007 16:45:01 +0000 (16:45 +0000)
PWG1/AliESDComparisonMI.C [new file with mode: 0644]
PWG1/AliESDComparisonMI.h [new file with mode: 0644]
PWG1/AliGenInfo.cxx [new file with mode: 0644]
PWG1/AliGenInfo.h [new file with mode: 0644]
PWG1/AliTreeDraw.cxx [new file with mode: 0644]
PWG1/AliTreeDraw.h [new file with mode: 0644]
PWG1/PWG1LinkDef.h [new file with mode: 0644]
PWG1/libPWG1.pkg [new file with mode: 0644]
PWG1/rootlogon.C [new file with mode: 0644]

diff --git a/PWG1/AliESDComparisonMI.C b/PWG1/AliESDComparisonMI.C
new file mode 100644 (file)
index 0000000..39c5b6c
--- /dev/null
@@ -0,0 +1,1699 @@
+/**************************************************************************
+ * 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,0,0);
+t2->Exec();
+
+//
+//some cuts definition
+TCut cprim("cprim","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<0.01&&abs(MC.fVDist[2])<0.01")
+//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 citsin("citsin","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<5");
+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.fParticle.Theta()/3.1415-0.5)<0.25");
+TCut cteta05("cteta05","abs(MC.fParticle.Theta()/3.1415-0.5)<0.1");
+
+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(); 
+TFile f("genHits.root");
+TTree * treel = (TTree*)f.Get("HitLines");
+if (treel) comp->fTree->AddFriend(treel,"L");
+
+//
+//example
+comp.fTree->SetAlias("radius","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)");
+comp.fTree->SetAlias("direction","MC.fParticle.fVx*MC.fParticle.fPx+MC.fParticle.fVy*MC.fParticle.fPy");
+comp.fTree->SetAlias("decaydir","MC.fTRdecay.fX*MC.fTRdecay.fPx+MC.fTRdecay.fY*MC.fTRdecay.fPy");
+comp.fTree->SetAlias("theta","MC.fTrackRef.Theta()");
+comp.fTree->SetAlias("primdca","sqrt(RC.fITStrack.fD[0]**2+RC.fITStrack.fD[1]**2)");
+comp.fTree->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN");
+comp.fTree->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN");
+
+
+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.DrawXY("fTPCinP0[3]","fTPCDelta[4]/fTPCinP1[3]","fReconstructed==1"+cprim,"1",4,0.2,1.5,-0.06,0.06)
+comp.fRes->Draw();
+comp.fMean->Draw();  
+
+comp.DrawXY("fITSinP0[3]","fITSDelta[4]/fITSinP1[3]","fReconstructed==1&&fITSOn"+cprim,"1",4,0.2,1.5,-0.06,0.06)
+comp.fRes->Draw();
+
+comp.Eff("fTPCinP0[3]","fRowsWithDigits>120"+cteta1+cpos1+cprim,"fTPCOn",20,0.2,1.5)
+comp.fRes->Draw();
+
+comp.Eff("fTPCinP0[3]","fRowsWithDigits>120"+cteta1+cpos1+cprim,"fTPCOn&&fITSOn&&fESDTrack.fITSFakeRatio<0.1",10,0.2,1.5)
+comp.fRes->Draw();
+comp.Eff("fTPCinP0[3]","fRowsWithDigits>120"+cteta1+cpos1+cprim,"fTPCOn&&fITSOn&&fESDTrack.fITSFakeRatio>0.1",10,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,550); hdedx0->SetMarkerColor(1);
+TH2F * hdedx1 = new TH2F("dEdx1","dEdx1",100, 0,2,200,0,550); hdedx1->SetMarkerColor(4);
+TH2F * hdedx2 = new TH2F("dEdx2","dEdx2",100, 0,2,200,0,550); hdedx2->SetMarkerColor(3);
+TH2F * hdedx3 = new TH2F("dEdx3","dEdx3",100, 0,2,200,0,550); hdedx3->SetMarkerColor(2);
+
+comp.fTree->Draw("fESDTrack.fITSsignal:MC.fParticle.P()>>dEdx0","fITSOn&&abs(fPdg)==211&&fITStrack.fN==6"+cprim) 
+comp.fTree->Draw("fESDTrack.fITSsignal:MC.fParticle.P()>>dEdx1","fITSOn&&abs(fPdg)==2212&&fITStrack.fN==6"+cprim) 
+comp.fTree->Draw("fESDTrack.fITSsignal:MC.fParticle.P()>>dEdx2","fITSOn&&abs(fPdg)==321&&fITStrack.fN==6"+cprim) 
+comp.fTree->Draw("fESDTrack.fITSsignal:MC.fParticle.P()>>dEdx3","fITSOn&&abs(fPdg)==11&&fITStrack.fN==6"+cprim) 
+
+
+comp.fTree->Draw("fESDTrack.fTRDsignal:MC.fParticle.P()>>dEdx0","fTRDOn&&abs(fPdg)==211&&fTRDtrack.fN>40&&fStatus[2]>1") 
+comp.fTree->Draw("fESDTrack.fTRDsignal:MC.fParticle.P()>>dEdx1","fTRDOn&&abs(fPdg)==2212&&fTRDtrack.fN>40&&fStatus[2]>1") 
+comp.fTree->Draw("fESDTrack.fTRDsignal:MC.fParticle.P()>>dEdx2","fTRDOn&&abs(fPdg)==321&&fTRDtrack.fN>40&&fStatus[2]>1") 
+comp.fTree->Draw("fESDTrack.fTRDsignal:MC.fParticle.P()>>dEdx3","fTRDOn&&abs(fPdg)==11&&fTRDtrack.fN>40&&fStatus[2]>1") 
+
+comp.fTree->Draw("fESDTrack.fTPCsignal:fTPCinP0[4]>>dEdx0","fTPCOn&&abs(fPdg)==211&&fESDTrack.fTPCncls>180&&fESDTrack.fTPCsignal>10"+cteta1); 
+comp.fTree->Draw("fESDTrack.fTPCsignal:fTPCinP0[4]>>dEdx1","fTPCOn&&abs(fPdg)==2212&&fESDTrack.fTPCncls>180&&fESDTrack.fTPCsignal>10"+cteta1); 
+comp.fTree->Draw("fESDTrack.fTPCsignal:fTPCinP0[4]>>dEdx2","fTPCOn&&abs(fPdg)==321&&fESDTrack.fTPCncls>180&&fESDTrack.fTPCsignal>10"+cteta1); 
+comp.fTree->Draw("fESDTrack.fTPCsignal:fTPCinP0[4]>>dEdx3","fTPCOn&&abs(fPdg)==11&&fESDTrack.fTPCncls>180&&fESDTrack.fTPCsignal>10"+cteta1); 
+
+hdedx3->SetXTitle("P(GeV/c)");
+hdedx3->SetYTitle("dEdx(unit)");
+hdedx3->Draw(); hdedx1->Draw("same"); hdedx2->Draw("same"); hdedx0->Draw("same");
+
+comp.DrawXY("fITSinP0[3]","fITSPools[4]","fReconstructed==1&&fPdg==-211&&fITSOn"+cprim,"1",4,0.2,1.0,-8,8)
+
+TProfile prof("prof","prof",10,0.5,5);
+
+
+
+
+*/
+
+
+#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 "AliESDfriend.h"
+#include "AliESDtrack.h"
+#include "AliITStrackMI.h"
+#include "AliTRDtrack.h"
+#include "AliHelix.h"
+#include "AliESDVertex.h"
+#include "AliExternalTrackParam.h"
+#include "AliESDkink.h"
+#include "AliESDv0.h"
+#include "AliV0.h"
+
+#endif
+#include "AliGenInfo.h"
+#include "AliESDComparisonMI.h"
+
+
+
+
+
+void MakeAliases(AliESDComparisonDraw&comp)
+{
+  //
+  // aliases definition
+  //
+  comp.fTree->SetAlias("radius","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)");
+  comp.fTree->SetAlias("direction","MC.fParticle.fVx*MC.fParticle.fPx+MC.fParticle.fVy*MC.fParticle.fPy");
+  comp.fTree->SetAlias("decaydir","MC.fTRdecay.fX*MC.fTRdecay.fPx+MC.fTRdecay.fY*MC.fTRdecay.fPy");
+  comp.fTree->SetAlias("theta","MC.fTrackRef.Theta()");
+  comp.fTree->SetAlias("primdca","sqrt(RC.fITStrack.fD[0]**2+RC.fITStrack.fD[1]**2)");
+  comp.fTree->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN");
+  comp.fTree->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN");
+  
+  comp.fTree->SetAlias("trddedx","(RC.fESDTrack.fTRDsignals[0]+RC.fESDTrack.fTRDsignals[1]+RC.fESDTrack.fTRDsignals[2]+RC.fESDTrack.fTRDsignals[3]+RC.fESDTrack.fTRDsignals[4]+RC.fESDTrack.fTRDsignals[5])/6.");
+  
+  comp.fTree->SetAlias("dtofmc2","fESDTrack.fTrackTime[2]-(10^12*MC.fTOFReferences[0].fTime)");
+  comp.fTree->SetAlias("dtofrc2","(fESDTrack.fTrackTime[2]-fESDTrack.fTOFsignal)");
+
+  comp.fTree->SetAlias("psum","fESDTrack.fTOFr[4]+fESDTrack.fTOFr[3]+fESDTrack.fTOFr[2]+fESDTrack.fTOFr[1]+fESDTrack.fTOFr[0]");
+  comp.fTree->SetAlias("P0","fESDTrack.fTOFr[0]/psum");
+  comp.fTree->SetAlias("P1","fESDTrack.fTOFr[1]/psum");
+  comp.fTree->SetAlias("P2","fESDTrack.fTOFr[2]/psum");
+  comp.fTree->SetAlias("P3","fESDTrack.fTOFr[3]/psum");
+  comp.fTree->SetAlias("P4","fESDTrack.fTOFr[4]/psum");
+  comp.fTree->SetAlias("MaxP","max(max(max(P0,P1),max(P2,P3)),P4)");
+}
+
+
+
+
+void  AliESDRecInfo::UpdatePoints(AliESDtrack*track)
+{
+  //
+  //
+  Int_t iclusters[200];
+  Float_t density[160];
+  for (Int_t i=0;i<160;i++) density[i]=-1.;
+  fTPCPoints[0]= 160;
+  fTPCPoints[1] = -1;
+  //
+  if (fTPCPoints[0]<fTPCPoints[1]) return;
+  //  Int_t nclusters=track->GetTPCclusters(iclusters);
+
+  Int_t ngood=0;
+  Int_t undeff=0;
+  Int_t nall =0;
+  Int_t range=20;
+  for (Int_t i=0;i<160;i++){
+    Int_t last = i-range;
+    if (nall<range) nall++;
+    if (last>=0){
+      if (iclusters[last]>0&& (iclusters[last]&0x8000)==0) ngood--;
+      if (iclusters[last]==-1) undeff--;
+    }
+    if (iclusters[i]>0&& (iclusters[i]&0x8000)==0)   ngood++;
+    if (iclusters[i]==-1) undeff++;
+    if (nall==range &&undeff<range/2) density[i-range/2] = Float_t(ngood)/Float_t(nall-undeff);
+  }
+  Float_t maxdens=0;
+  Int_t indexmax =0;
+  for (Int_t i=0;i<160;i++){
+    if (density[i]<0) continue;
+    if (density[i]>maxdens){
+      maxdens=density[i];
+      indexmax=i;
+    }
+  }
+  //
+  //max dens point
+  fTPCPoints[3] = maxdens;
+  fTPCPoints[1] = indexmax;
+  //
+  // last point
+  for (Int_t i=indexmax;i<160;i++){
+    if (density[i]<0) continue;
+    if (density[i]<maxdens/2.) {
+      break;
+    }
+    fTPCPoints[2]=i;
+  }
+  //
+  // first point
+  for (Int_t i=indexmax;i>0;i--){
+    if (density[i]<0) continue;
+    if (density[i]<maxdens/2.) {
+      break;
+    }
+    fTPCPoints[0]=i;
+  }
+  //
+  // Density at the last 30 padrows
+  //
+  // 
+  nall  = 0;
+  ngood = 0;
+  for (Int_t i=159;i>0;i--){
+    if (iclusters[i]==-1) continue; //dead zone
+    nall++;
+    if (iclusters[i]>0)   ngood++;
+    if (nall>20) break;
+  }
+  fTPCPoints[4] = Float_t(ngood)/Float_t(nall);
+  //
+  if ((track->GetStatus()&AliESDtrack::kITSrefit)>0) fTPCPoints[0]=-1;
+  //
+  //
+  // check TRDPoints
+  /*
+  nclusters=track->GetTRDclusters(iclusters);
+  for (Int_t i=nclusters;i>0;i--){
+    
+  }
+  */
+
+
+}
+
+//
+//
+void AliESDRecInfo::Update(AliMCInfo* info,AliTPCParam * /*par*/, Bool_t reconstructed, AliESD */*event*/)
+{
+  //
+  //
+  //calculates derived variables
+  //  
+  //
+  UpdatePoints(&fESDTrack);
+  fBestTOFmatch=1000;
+  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]);
+  }
+  //
+  for (Int_t i=0;i<4;i++) fStatus[i] =0;
+  fReconstructed = kFALSE;
+  fTPCOn = kFALSE;
+  fITSOn = kFALSE;
+  fTRDOn = kFALSE;  
+  if (reconstructed==kFALSE) return;
+
+  fLabels[0] = info->fLabel;
+  fLabels[1] = info->fPrimPart;
+  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 ((fESDTrack.GetStatus()&AliESDtrack::kTPCrefit)>0){
+    fStatus[1] =3;
+  }
+  else{
+    if ((fESDTrack.GetStatus()&AliESDtrack::kTPCout)>0){
+      fStatus[1] =2;
+    }
+    else{
+      if ((fESDTrack.GetStatus()&AliESDtrack::kTPCin)>0)
+       fStatus[1]=1;
+    }      
+  }
+  //
+  if ((fESDTrack.GetStatus()&AliESDtrack::kITSout)>0){
+    fStatus[0] =2;
+  }
+  else{
+    if ((fESDTrack.GetStatus()&AliESDtrack::kITSrefit)>0){
+      fStatus[0] =1;
+    }
+    else{
+      fStatus[0]=0;
+    }      
+  }
+
+  //
+  //
+  if ((fESDTrack.GetStatus()&AliESDtrack::kTRDrefit)>0){
+    fStatus[2] =2;
+  }
+  else{
+    if ((fESDTrack.GetStatus()&AliESDtrack::kTRDout)>0){
+      fStatus[2] =1;
+    }
+  }
+  if ((fESDTrack.GetStatus()&AliESDtrack::kTRDStop)>0){
+    fStatus[2] =10;
+  }
+
+  //
+  //TOF 
+  // 
+  if (((fESDTrack.GetStatus()&AliESDtrack::kTOFout)>0)){
+    //
+    // best tof match
+    Double_t times[5];
+    fESDTrack.GetIntegratedTimes(times);    
+    for (Int_t i=0;i<5;i++){
+      if ( TMath::Abs(fESDTrack.GetTOFsignal()-times[i]) <TMath::Abs(fBestTOFmatch) ){
+       fBestTOFmatch = fESDTrack.GetTOFsignal()-times[i];
+      }
+    }
+    Int_t toflabel[3];
+    fESDTrack.GetTOFLabel(toflabel);
+    Bool_t toffake=kTRUE;
+    Bool_t tofdaughter=kFALSE;
+    for (Int_t i=0;i<3;i++){
+      if (toflabel[i]<0) continue;      
+      if (toflabel[i]== TMath::Abs(fESDTrack.GetLabel()))  toffake=kFALSE;     
+      if (toflabel[i]==info->fParticle.GetDaughter(0) || (toflabel[i]==info->fParticle.GetDaughter(1))) tofdaughter=kTRUE;  // decay product of original particle
+      fStatus[3]=1;
+    }
+    if (toffake) fStatus[3] =3;       //total fake
+    if (tofdaughter) fStatus[3]=2;    //fake because of decay
+  }else{
+    fStatus[3]=0;
+  }
+
+
+  if (fStatus[1]>0 &&info->fNTPCRef>0&&TMath::Abs(fTPCinP0[3])>0.0001){
+    //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.000000000000001){
+      fTPCAngle1[0] = TMath::ATan2(fTPCinP1[1],fTPCinP1[0]);
+      fTPCAngle1[1] = TMath::ATan(fTPCinP1[2]/fTPCinP1[3]);  
+    }    
+    Double_t cov[15], param[5],x, alpha;
+    fESDTrack.GetInnerExternalCovariance(cov);
+    fESDTrack.GetInnerExternalParameters(alpha, x,param);
+    if (x<50) return ;
+    //
+    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; 
+    fSign =sign;
+    fTPCPools[4] = sign*(1./fTPCinP0[3]-1./fTPCinP1[3])/TMath::Sqrt(TMath::Abs(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; 
+    fSign = sign;
+    fITSPools[4] = sign*(1./fITSinP0[3]-1./fITSinP1[3])/TMath::Sqrt(cov[14]);    
+  }
+  
+}
+
+
+void  AliESDRecV0Info::Update(Float_t vertex[3])
+{ 
+
+  if ( (fT1.fStatus[1]>0)&& (fT2.fStatus[1]>0)){
+    Float_t distance1,distance2;
+    Double_t xx[3],pp[3];
+    //
+    Double_t xd[3],pd[3],signd;
+    Double_t xm[3],pm[3],signm;
+    //
+    //
+    if (fT1.fITSOn&&fT2.fITSOn){
+      for (Int_t i=0;i<3;i++){
+       xd[i] = fT2.fITSinR1[i];
+       pd[i] = fT2.fITSinP1[i];
+       xm[i] = fT1.fITSinR1[i];
+       pm[i] = fT1.fITSinP1[i];
+      }
+    }
+    else{
+      
+      for (Int_t i=0;i<3;i++){
+       xd[i] = fT2.fTPCinR1[i];
+       pd[i] = fT2.fTPCinP1[i];
+       xm[i] = fT1.fTPCinR1[i];
+       pm[i] = fT1.fTPCinP1[i];
+      }
+    }
+    //
+    //
+    signd =  fT2.fSign<0 ? -1:1;
+    signm =  fT1.fSign<0 ? -1:1;
+
+    AliHelix dhelix1(xd,pd,signd);
+    dhelix1.GetMomentum(0,pp,0);
+    dhelix1.Evaluate(0,xx);      
+    // 
+    //  Double_t x2[3],p2[3];
+    //            
+    AliHelix mhelix(xm,pm,signm);    
+    //
+    //find intersection linear
+    //
+    Double_t phase[2][2],radius[2];
+    Int_t  points = dhelix1.GetRPHIintersections(mhelix, phase, radius,200);
+    Double_t delta1=10000,delta2=10000;  
+
+    if (points==1){
+      fRs[0] = TMath::Sqrt(radius[0]);
+      fRs[1] = TMath::Sqrt(radius[0]);
+    }
+    if (points==2){
+      fRs[0] =TMath::Min(TMath::Sqrt(radius[0]),TMath::Sqrt(radius[1]));
+      fRs[1] =TMath::Max(TMath::Sqrt(radius[0]),TMath::Sqrt(radius[1]));
+    }
+    
+    if (points>0){
+      dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
+      dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
+      dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
+    }
+    if (points==2){    
+      dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
+      dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
+      dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
+    }
+    if (points==1){
+      fRs[0] = TMath::Sqrt(radius[0]);
+      fRs[1] = TMath::Sqrt(radius[0]);
+      fDistMinR = delta1;
+    }
+    if (points==2){
+      if (radius[0]<radius[1]){
+       fRs[0] = TMath::Sqrt(radius[0]);
+       fRs[1] = TMath::Sqrt(radius[1]);
+       fDistMinR = delta1;
+      }
+      else{
+       fRs[0] = TMath::Sqrt(radius[1]);
+       fRs[1] = TMath::Sqrt(radius[0]);
+       fDistMinR = delta2;
+      }
+    }
+    //
+    //
+    distance1 = TMath::Min(delta1,delta2);
+    //
+    //find intersection parabolic
+    //
+    points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
+    delta1=10000,delta2=10000;  
+    
+    if (points>0){
+      dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
+    }
+    if (points==2){    
+      dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
+    }
+    
+    distance2 = TMath::Min(delta1,delta2);
+    if (distance2>100) fDist2 =100;
+    return;
+    if (delta1<delta2){
+      //get V0 info
+      dhelix1.Evaluate(phase[0][0],fXr);
+      dhelix1.GetMomentum(phase[0][0],fPdr);
+      mhelix.GetMomentum(phase[0][1],fPm);
+      dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],fAngle);
+      fRr = TMath::Sqrt(radius[0]);
+    }
+    else{
+      dhelix1.Evaluate(phase[1][0],fXr);
+      dhelix1.GetMomentum(phase[1][0], fPdr);
+      mhelix.GetMomentum(phase[1][1], fPm);
+      dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],fAngle);
+      fRr = TMath::Sqrt(radius[1]);
+    }
+    fDist1 = TMath::Sqrt(distance1);
+    fDist2 = TMath::Sqrt(distance2);      
+    
+    if (fDist2<10.5){
+      Double_t x,alpha,param[5],cov[15];
+      //
+      fT1.fESDTrack.GetInnerExternalParameters(alpha,x,param);
+      fT1.fESDTrack.GetInnerExternalCovariance(cov);
+      AliExternalTrackParam paramm(x,alpha,param,cov);
+      //
+      fT2.fESDTrack.GetInnerExternalParameters(alpha,x,param);
+      fT2.fESDTrack.GetInnerExternalCovariance(cov);
+      AliExternalTrackParam paramd(x,alpha,param,cov);
+    }    
+    //            
+    //   
+    
+    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);
+  }
+}
+
+////
+void  AliESDRecKinkInfo::Update()
+{
+
+  if ( (fT1.fTPCOn)&& (fT2.fTPCOn)){
+    //
+    // IF BOTH RECONSTRUCTED
+    Float_t distance1,distance2;
+    Double_t xx[3],pp[3];
+    //
+    Double_t xd[3],pd[3],signd;
+    Double_t xm[3],pm[3],signm;
+    for (Int_t i=0;i<3;i++){
+      xd[i] = fT2.fTPCinR1[i];
+      pd[i] = fT2.fTPCinP1[i];
+      xm[i] = fT1.fTPCinR1[i];
+      pm[i] = fT1.fTPCinP1[i];
+    }
+    signd =  fT2.fSign<0 ? -1:1;
+    signm =  fT1.fSign<0 ? -1:1;
+
+    AliHelix dhelix1(xd,pd,signd);
+    dhelix1.GetMomentum(0,pp,0);
+    dhelix1.Evaluate(0,xx);      
+    // 
+    //  Double_t x2[3],p2[3];
+    //            
+    AliHelix mhelix(xm,pm,signm);    
+    //
+    //find intersection linear
+    //
+    Double_t phase[2][2],radius[2];
+    Int_t  points = dhelix1.GetRPHIintersections(mhelix, phase, radius,200);
+    Double_t delta1=10000,delta2=10000;  
+
+    if (points==1){
+      fMinR = TMath::Sqrt(radius[0]);
+    }
+    if (points==2){
+      fMinR =TMath::Min(TMath::Sqrt(radius[0]),TMath::Sqrt(radius[1]));
+    }
+    
+    if (points>0){
+      dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
+      dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
+      dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
+    }
+    if (points==2){    
+      dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
+      dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
+      dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
+    }
+    if (points==1){
+      fMinR = TMath::Sqrt(radius[0]);
+      fDistMinR = delta1;
+    }
+    if (points==2){
+      if (radius[0]<radius[1]){
+       fMinR = TMath::Sqrt(radius[0]);
+       fDistMinR = delta1;
+      }
+      else{
+       fMinR = TMath::Sqrt(radius[1]);
+       fDistMinR = delta2;
+      }
+    }
+    //
+    //
+    distance1 = TMath::Min(delta1,delta2);
+    //
+    //find intersection parabolic
+    //
+    points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
+    delta1=10000,delta2=10000;  
+    
+    if (points>0){
+      dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
+    }
+    if (points==2){    
+      dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
+    }
+    
+    distance2 = TMath::Min(delta1,delta2);
+    if (delta1<delta2){
+      //get V0 info
+      dhelix1.Evaluate(phase[0][0],fXr);
+      dhelix1.GetMomentum(phase[0][0],fPdr);
+      mhelix.GetMomentum(phase[0][1],fPm);
+      dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],fAngle);
+      fRr = TMath::Sqrt(radius[0]);
+    }
+    else{
+      dhelix1.Evaluate(phase[1][0],fXr);
+      dhelix1.GetMomentum(phase[1][0], fPdr);
+      mhelix.GetMomentum(phase[1][1], fPm);
+      dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],fAngle);
+      fRr = TMath::Sqrt(radius[1]);
+    }
+    fDist1 = TMath::Sqrt(distance1);
+    fDist2 = TMath::Sqrt(distance2);      
+    
+    if (fDist2<10.5){
+      Double_t x,alpha,param[5],cov[15];
+      //
+      fT1.fESDTrack.GetInnerExternalParameters(alpha,x,param);
+      fT1.fESDTrack.GetInnerExternalCovariance(cov);
+      AliExternalTrackParam paramm(x,alpha,param,cov);
+      //
+      fT2.fESDTrack.GetInnerExternalParameters(alpha,x,param);
+      fT2.fESDTrack.GetInnerExternalCovariance(cov);
+      AliExternalTrackParam paramd(x,alpha,param,cov);
+      /*
+      AliESDkink kink;
+      kink.Update(&paramm,&paramd);
+      //      kink.Dump();
+      Double_t diff  = kink.fRr-fRr;
+      Double_t diff2 = kink.fDist2-fDist2;
+      printf("Diff\t%f\t%f\n",diff,diff2);
+      */
+    }
+    
+    //            
+    //
+  }
+
+}
+
+
+////////////////////////////////////////////////////////////////////////
+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,0);
+
+}
+
+
+////////////////////////////////////////////////////////////////////////
+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->SetBranchStatus("*",1);
+    tree->SetBranchAddress("ESD", &fEvent);
+    tree->SetBranchAddress("ESDfriend.",&fESDfriend);
+    tree->GetEntry(eventNr);
+    fEvent->SetESDfriend(fESDfriend);    
+  }
+
+
+  /*
+  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;
+  fTreeCmpKinks =0;
+  fTreeCmpV0 =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;
+  fNextKinkToRead = 0;
+  fNextV0ToRead   =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 Short_t[fNParticles*20];  //write at maximum 4 tracks corresponding to particle
+    fIndexRecKinks  = new Short_t[fNParticles*20];  //write at maximum 20 tracks corresponding to particle
+    fIndexRecV0  = new Short_t[fNParticles*20];  //write at maximum 20 tracks corresponding to particle
+
+    fFakeRecTracks = new Short_t[fNParticles];
+    fMultiRecTracks = new Short_t[fNParticles];
+    fMultiRecKinks = new Short_t[fNParticles];
+    fMultiRecV0 = new Short_t[fNParticles];
+
+    for (Int_t i = 0; i<fNParticles; i++) {
+      for (Int_t j=0;j<20;j++){
+       fIndexRecTracks[20*i+j] = -1;
+       fIndexRecKinks[20*i+j]  = -1;
+       fIndexRecV0[20*i+j]  = -1;
+      }
+      fFakeRecTracks[i] = 0;
+      fMultiRecTracks[i] = 0;
+      fMultiRecKinks[i] = 0;
+      fMultiRecV0[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;
+    BuildKinkInfo0(eventNr);
+    BuildV0Info(eventNr);
+    fRecArray->Delete();
+
+    if (fDebug>2) cout<<"\tEnd loop over tree genTracks"<<endl;
+
+    delete [] fIndexRecTracks;
+    delete [] fIndexRecKinks;
+    delete [] fIndexRecV0;
+    delete [] fFakeRecTracks;
+    delete [] fMultiRecTracks;
+    delete [] fMultiRecKinks;
+    delete [] fMultiRecV0;
+  }
+
+  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);
+  //
+  //
+  fTreeGenKinks = (TTree*)fFileGenTracks->Get("genKinksTree");
+  if (!fTreeGenKinks) {
+    cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file "
+       <<fFnGenTracks<<endl;
+    //return kFALSE;
+  }
+  else{
+    fGenKinkInfo = new AliGenKinkInfo;
+    fTreeGenKinks->SetBranchAddress("MC",&fGenKinkInfo);
+  }
+
+  fTreeGenV0 = (TTree*)fFileGenTracks->Get("genV0Tree");
+  if (!fTreeGenV0) {
+    cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file "
+       <<fFnGenTracks<<endl;
+    //return kFALSE;
+  }
+  else{
+    fGenV0Info = new AliGenV0Info;
+    fTreeGenV0->SetBranchAddress("MC",&fGenV0Info);
+  }
+  //
+  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; 
+  //  AliITStrackMI * itsTrack = new AliITStrackMI;  
+  fTreeCmp->Branch("MC","AliMCInfo",&fMCInfo,256000);
+  fTreeCmp->Branch("RC","AliESDRecInfo",&fRecInfo,256000);
+  //  fTreeCmp->Branch("ITS","AliITStrackMI",&itsTrack);
+  delete esdTrack;
+  //
+  //
+  fTreeCmpKinks    = new TTree("ESDcmpKinks","ESDcmpKinks"); 
+  fGenKinkInfo     = new AliGenKinkInfo;
+  fRecKinkInfo     = new AliESDRecKinkInfo;
+  fTreeCmpKinks->Branch("MC.","AliGenKinkInfo",&fGenKinkInfo,256000);
+  fTreeCmpKinks->Branch("RC.","AliESDRecKinkInfo",&fRecKinkInfo,256000);
+  //
+  //
+  fTreeCmpV0       = new TTree("ESDcmpV0","ESDcmpV0"); 
+  fGenV0Info     = new AliGenV0Info;
+  fRecV0Info     = new AliESDRecV0Info;
+  fTreeCmpV0->Branch("MC.","AliGenV0Info",   &fGenV0Info,256000);
+  fTreeCmpV0->Branch("RC.","AliESDRecV0Info",&fRecV0Info,256000);
+  //
+  fTreeCmp->AutoSave(); 
+  fTreeCmpKinks->AutoSave(); 
+  fTreeCmpV0->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->Transform1to2Ideal(x,index);
+  return TVector3(x);
+}
+////////////////////////////////////////////////////////////////////////
+
+Int_t ESDCmpTr::TreeTLoop()
+{
+  //
+  // loop over all ESD reconstructed tracks and store info in memory
+  //
+  // + loop over all reconstructed kinks
+  TStopwatch  timer;
+  timer.Start();
+  //  
+  Int_t nEntries = (Int_t)fEvent->GetNumberOfTracks();  
+  Int_t nKinks = (Int_t) fEvent->GetNumberOfKinks();
+  Int_t nV0MIs = (Int_t) fEvent->GetNumberOfV0s();
+  fSignedKinks = new Short_t[nKinks];
+  fSignedV0    = new Short_t[nV0MIs];
+  //
+  // load kinks to the memory
+  for (Int_t i=0; i<nKinks;i++){
+    AliESDkink * kink =fEvent->GetKink(i);
+    fSignedKinks[i]=0;
+    if (kink->fStatus<0) continue;
+  }
+  //
+  for (Int_t i=0; i<nV0MIs;i++){
+    AliV0 * v0MI = (AliV0*)fEvent->GetV0(i);
+    fSignedV0[i]=0;
+    if (v0MI->fStatus<0) continue;
+  }
+  
+  //
+  //
+  AliESDtrack * track=0;
+  for (Int_t iEntry=0; iEntry<nEntries;iEntry++){
+    //track = (AliESDtrack*)fTracks->UncheckedAt(iEntry);
+    track = (AliESDtrack*)fEvent->GetTrack(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]<20)
+         fIndexRecTracks[absLabel*20+fMultiRecTracks[absLabel]] =  iEntry;     
+      }
+      else      
+       fIndexRecTracks[absLabel*20] =  iEntry;
+      fMultiRecTracks[absLabel]++;
+    }
+  }
+  // sort reconstructed kinks  
+  //
+  AliESDkink * kink=0;
+  for (Int_t iEntry=0; iEntry<nKinks;iEntry++){
+    kink = (AliESDkink*)fEvent->GetKink(iEntry);
+    if (!kink) continue;
+    //
+    Int_t label0 = TMath::Abs(kink->GetLabel(0));
+    Int_t label1 = TMath::Abs(kink->GetLabel(1));
+    Int_t absLabel = TMath::Min(label0,label1);
+    if (absLabel < fNParticles) {
+      if (fMultiRecKinks[absLabel]>0){
+       if (fMultiRecKinks[absLabel]<20)
+         fIndexRecKinks[absLabel*20+fMultiRecKinks[absLabel]] =  iEntry;       
+      }
+      else      
+       fIndexRecKinks[absLabel*20] =  iEntry;
+      fMultiRecKinks[absLabel]++;
+    }
+  }  
+  // --sort reconstructed V0
+  //
+  AliV0 * v0MI=0;
+  for (Int_t iEntry=0; iEntry<nV0MIs;iEntry++){
+    v0MI = (AliV0*)fEvent->GetV0(iEntry);
+    if (!v0MI) continue;
+    //
+    //    Int_t label0 = TMath::Abs(v0MI->fLab[0]);
+    //Int_t label1 = TMath::Abs(v0MI->fLab[1]);
+    //
+    for (Int_t i=0;i<2;i++){
+      // Int_t absLabel = TMath::Min(label0,label1);      
+      //      Int_t absLabel =  TMath::Abs(v0MI->fLab[i]);
+      Int_t absLabel =  TMath::Abs(0);
+      if (absLabel < fNParticles) {
+       if (fMultiRecV0[absLabel]>0){
+         if (fMultiRecV0[absLabel]<20)
+           fIndexRecV0[absLabel*20+fMultiRecV0[absLabel]] =  iEntry;   
+       }
+       else      
+         fIndexRecV0[absLabel*20] =  iEntry;
+       fMultiRecV0[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");
+  TBranch * branchF = fTreeCmp->GetBranch("F");
+  
+  branch->SetAddress(&fRecInfo); // set all pointers
+  fRecArray = new TObjArray(fNParticles);
+  AliESDtrack dummytrack;  //
+  AliESDfriendTrack dummytrackF;  //
+
+  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;
+    }
+    //    if (fMCInfo->fNTPCRef<1) continue; // not TPCref
+    //
+    fRecInfo->Reset();
+    AliESDtrack * track=0;
+    fRecInfo->fReconstructed =0;
+    TVector3 local = TR2Local(&(fMCInfo->fTrackRef),fParamTPC);
+    local.GetXYZ(fRecInfo->fTRLocalCoord);     
+    //
+    if (fIndexRecTracks[fMCInfo->fLabel*20] >= 0) {
+      //track= (AliESDtrack*)fTracks->UncheckedAt(fIndexRecTracks[fMCInfo->fLabel*4]);
+      track= (AliESDtrack*)fEvent->GetTrack(fIndexRecTracks[fMCInfo->fLabel*20]);
+      //
+      //
+      // find nearest track if multifound
+      //Int_t sign = Int_t(track->GetSign()*fMCInfo->fCharge);
+      //
+      Int_t status = 0;
+      if  ((track->GetStatus()&AliESDtrack::kITSrefit)>0) status++;
+      if  ((track->GetStatus()&AliESDtrack::kTPCrefit)>0) status++;
+      if  ((track->GetStatus()&AliESDtrack::kTRDrefit)>0) status++;
+
+      //
+      if (fIndexRecTracks[fMCInfo->fLabel*20+1]>0){
+       //
+       Double_t p[3];
+       track->GetInnerPxPyPz(p);
+       Float_t maxp = p[0]*p[0]+p[1]*p[1]+p[2]*p[2];
+       //
+       for (Int_t i=1;i<20;i++){
+         if (fIndexRecTracks[fMCInfo->fLabel*20+i]>=0){
+           AliESDtrack * track2 = (AliESDtrack*)fEvent->GetTrack(fIndexRecTracks[fMCInfo->fLabel*20+i]);
+           if (!track2) continue;
+           //Int_t sign2 = track->GetSign()*fMCInfo->fCharge; //           
+           //if (sign2<0) continue;
+           track2->GetInnerPxPyPz(p);
+           Float_t mom = p[0]*p[0]+p[1]*p[1]+p[2]*p[2];
+           /*
+           if (sign<0){
+             sign = sign2;
+             track = track2;
+             maxp = mom;
+             continue;
+           }
+           */
+           //
+           Int_t status2 = 0;
+           if  ((track2->GetStatus()&AliESDtrack::kITSrefit)>0) status2++;
+           if  ((track2->GetStatus()&AliESDtrack::kTPCrefit)>0) status2++;
+           if  ((track2->GetStatus()&AliESDtrack::kTRDrefit)>0) status2++;
+           if (status2<status) continue;
+           //
+           if (mom<maxp) continue;
+           maxp = mom;
+           track = track2;
+           //
+         }
+       }
+      }        
+      //
+      if (track) {
+       new (&(fRecInfo->fESDTrack)) AliESDtrack(*track);
+       if (track->GetFriendTrack()){
+         new( &(fRecInfo->fTrackF)) AliESDfriendTrack(*(track->GetFriendTrack()));
+       }
+      }else{
+       fRecInfo->fESDTrack = dummytrack;
+       fRecInfo->fTrackF = dummytrackF;
+      }
+
+
+
+      if (track->GetITStrack())
+       //fRecInfo->fITStrack = *((AliITStrackMI*)track->GetITStrack());
+       new (&(fRecInfo->fITStrack)) AliITStrackMI(*((AliITStrackMI*)track->GetITStrack()));
+      else{
+       fRecInfo->fITStrack = *track;
+      }
+      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, fEvent);          
+    }
+    else{
+      fRecInfo->fESDTrack = dummytrack;
+      fRecInfo->Update(fMCInfo,fParamTPC,kFALSE, fEvent);
+    }
+    fRecArray->AddAt(new AliESDRecInfo(*fRecInfo),fMCInfo->fLabel);
+    fTreeCmp->Fill();
+  }
+  fTreeCmp->AutoSave();
+  //fTracks->Delete();
+  printf("Time spended in TreeGenLoop\n");
+  timer.Print();
+  if (fDebug > 2) cerr<<"end of TreeGenLoop"<<endl;
+
+  return 0;
+}
+
+
+
+////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////
+Int_t ESDCmpTr::BuildKinkInfo0(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 = fNextKinkToRead;
+  Double_t nParticlesTR = fTreeGenKinks->GetEntriesFast();
+  cerr<<"fNParticles, nParticlesTR, fNextKinkToRead: "<<fNParticles<<" "
+      <<nParticlesTR<<" "<<fNextKinkToRead<<endl;
+  //
+  TBranch * branch = fTreeCmpKinks->GetBranch("RC.");
+  branch->SetAddress(&fRecKinkInfo); // set all pointers
+  
+  //
+  while (entry < nParticlesTR) {
+    fTreeGenKinks->GetEntry(entry);
+    entry++;
+    if (eventNr < fGenKinkInfo->fMCm.fEventNr) continue;
+    if (eventNr > fGenKinkInfo->fMCm.fEventNr) continue;;
+    //
+    fNextKinkToRead = entry-1;
+    //
+    //
+    AliESDRecInfo*  fRecInfo1 = (AliESDRecInfo*)fRecArray->At(fGenKinkInfo->fMCm.fLabel);
+    AliESDRecInfo*  fRecInfo2 = (AliESDRecInfo*)fRecArray->At(fGenKinkInfo->fMCd.fLabel);
+    fRecKinkInfo->fT1 = (*fRecInfo1);
+    fRecKinkInfo->fT2 = (*fRecInfo2);
+    fRecKinkInfo->fStatus =0;
+    if (fRecInfo1 && fRecInfo1->fTPCOn) fRecKinkInfo->fStatus+=1;
+    if (fRecInfo2 && fRecInfo2->fTPCOn) fRecKinkInfo->fStatus+=2;
+    if (fRecKinkInfo->fStatus==3&&fRecInfo1->fSign!=fRecInfo2->fSign) fRecKinkInfo->fStatus*=-1;
+    
+    if (fRecKinkInfo->fStatus==3){
+      fRecKinkInfo->Update();    
+    }
+    Int_t label =  TMath::Min(fGenKinkInfo->fMCm.fLabel,fGenKinkInfo->fMCd.fLabel);
+    Int_t label2 = TMath::Max(fGenKinkInfo->fMCm.fLabel,fGenKinkInfo->fMCd.fLabel);
+    
+    AliESDkink *kink=0;
+    fRecKinkInfo->fRecStatus   =0;
+    fRecKinkInfo->fMultiple    = fMultiRecKinks[label];
+    fRecKinkInfo->fKinkMultiple=0;
+    //
+    if (fMultiRecKinks[label]>0){
+
+      //      for (Int_t j=0;j<TMath::Min(fMultiRecKinks[label],100);j++){
+      for (Int_t j=TMath::Min(fMultiRecKinks[label],Short_t(20))-1;j>=0;j--){
+       Int_t index = fIndexRecKinks[label*20+j];
+       //AliESDkink *kink2  = (AliESDkink*)fKinks->At(index);
+       AliESDkink *kink2  = (AliESDkink*)fEvent->GetKink(index);
+       if (TMath::Abs(kink2->fLab[0])==label &&TMath::Abs(kink2->fLab[1])==label2) {
+         fRecKinkInfo->fKinkMultiple++;
+         fSignedKinks[index]=1;
+         Int_t c0=0;
+         if (kink){
+           //      if (kink->fTRDOn) c0++;
+           //if (kink->fITSOn) c0++;
+           if (kink->GetStatus(2)>0) c0++;
+           if (kink->GetStatus(0)>0) c0++;
+         }
+         Int_t c2=0;
+         //      if (kink2->fTRDOn) c2++;
+         //if (kink2->fITSOn) c2++;
+         if (kink2->GetStatus(2)>0) c2++;
+         if (kink2->GetStatus(0)>0) c2++;
+
+         if (c2<c0) continue;
+         kink =kink2;
+       }
+       if (TMath::Abs(kink2->fLab[1])==label &&TMath::Abs(kink2->fLab[0])==label2) {
+         fRecKinkInfo->fKinkMultiple++;
+         fSignedKinks[index]=1;
+         Int_t c0=0;
+         if (kink){
+           //if (kink->fTRDOn) c0++;
+           //if (kink->fITSOn) c0++;
+           if (kink->GetStatus(2)>0) c0++;
+           if (kink->GetStatus(0)>0) c0++;
+
+         }
+         Int_t c2=0;
+         //      if (kink2->fTRDOn) c2++;
+         //if (kink2->fITSOn) c2++;
+         if (kink2->GetStatus(2)>0) c2++;
+         if (kink2->GetStatus(0)>0) c2++;
+
+         if (c2<c0) continue;
+         kink =kink2;
+       }
+      }
+    }
+    if (kink){
+      fRecKinkInfo->fKink = *kink;
+      fRecKinkInfo->fRecStatus=1;
+    }
+    fTreeCmpKinks->Fill();
+  }
+  //  Int_t nkinks = fKinks->GetEntriesFast();
+  Int_t nkinks = fEvent->GetNumberOfKinks();
+  for (Int_t i=0;i<nkinks;i++){
+    if (fSignedKinks[i]==0){
+      //      AliESDkink *kink  = (AliESDkink*)fKinks->At(i);
+      AliESDkink *kink  = (AliESDkink*)fEvent->GetKink(i);
+      if (!kink) continue;
+      //
+      fRecKinkInfo->fKink = *kink;
+      fRecKinkInfo->fRecStatus =-2;
+      //
+      AliESDRecInfo*  fRecInfo1 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(kink->fLab[0]));
+      AliESDRecInfo*  fRecInfo2 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(kink->fLab[1]));
+      if (fRecInfo1 && fRecInfo2){
+       fRecKinkInfo->fT1 = (*fRecInfo1);
+       fRecKinkInfo->fT2 = (*fRecInfo2);
+       fRecKinkInfo->fRecStatus =-1;
+      }
+      fTreeCmpKinks->Fill();
+    }
+  }
+
+
+  fTreeCmpKinks->AutoSave();
+  printf("Time spended in BuilKinkInfo Loop\n");
+  timer.Print();
+  if (fDebug > 2) cerr<<"end of BuildKinkInfo Loop"<<endl;
+  return 0;
+}
+
+
+
+////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////
+
+
+
+Int_t ESDCmpTr::BuildV0Info(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 = fNextV0ToRead;
+  Double_t nParticlesTR = fTreeGenV0->GetEntriesFast();
+  cerr<<"fNParticles, nParticlesTR, fNextV0ToRead: "<<fNParticles<<" "
+      <<nParticlesTR<<" "<<fNextV0ToRead<<endl;
+  //
+  TBranch * branch = fTreeCmpV0->GetBranch("RC.");
+  branch->SetAddress(&fRecV0Info); // set all pointers
+  const AliESDVertex * esdvertex = fEvent->GetVertex();
+  Float_t vertex[3]= {esdvertex->GetXv(), esdvertex->GetYv(),esdvertex->GetZv()};
+  
+  //
+  while (entry < nParticlesTR) {
+    fTreeGenV0->GetEntry(entry);
+    entry++;
+    if (eventNr < fGenV0Info->fMCm.fEventNr) continue;
+    if (eventNr > fGenV0Info->fMCm.fEventNr) continue;;
+    //
+    fNextV0ToRead = entry-1;
+    //
+    //
+    AliESDRecInfo*  fRecInfo1 = (AliESDRecInfo*)fRecArray->At(fGenV0Info->fMCm.fLabel);
+    AliESDRecInfo*  fRecInfo2 = (AliESDRecInfo*)fRecArray->At(fGenV0Info->fMCd.fLabel);
+    if (fGenV0Info->fMCm.fCharge*fGenV0Info->fMCd.fCharge>0) continue;  // interactions
+    if (!fRecInfo1 || !fRecInfo2) continue;
+    fRecV0Info->fT1 = (*fRecInfo1);
+    fRecV0Info->fT2 = (*fRecInfo2);
+    fRecV0Info->fV0Status =0;
+    if (fRecInfo1 && fRecInfo1->fStatus[1]>0) fRecV0Info->fV0Status+=1;
+    if (fRecInfo2 && fRecInfo2->fStatus[1]>0) fRecV0Info->fV0Status+=2;
+
+    if (fRecV0Info->fV0Status==3&&fRecInfo1->fSign==fRecInfo2->fSign) fRecV0Info->fV0Status*=-1;
+
+
+    if (abs(fRecV0Info->fV0Status)==3){
+      fRecV0Info->Update(vertex);
+      {
+       //
+       // TPC V0 Info
+       Double_t x,alpha, param[5],cov[15];
+       if ( fRecV0Info->fT1.fESDTrack.GetInnerParam() && fRecV0Info->fT2.fESDTrack.GetInnerParam()){
+         fRecV0Info->fT1.fESDTrack.GetInnerExternalParameters(alpha,x,param);
+         fRecV0Info->fT1.fESDTrack.GetInnerExternalCovariance(cov);
+         AliExternalTrackParam paramP(x,alpha,param,cov);
+         //
+         fRecV0Info->fT2.fESDTrack.GetInnerExternalParameters(alpha,x,param);
+         fRecV0Info->fT2.fESDTrack.GetInnerExternalCovariance(cov);
+         AliExternalTrackParam paramM(x,alpha,param,cov);
+         //
+         fRecV0Info->fV0tpc.SetParamN(paramM);
+         fRecV0Info->fV0tpc.SetParamP(paramP);
+         Double_t pid1[5],pid2[5];
+         fRecV0Info->fT1.fESDTrack.GetESDpid(pid1);
+         fRecV0Info->fT1.fESDTrack.GetESDpid(pid2);
+         //
+         //fRecV0Info->fV0tpc.UpdatePID(pid1,pid2);
+         fRecV0Info->fV0tpc.Update(vertex);
+       
+         //
+         //
+         fRecV0Info->fT1.fESDTrack.GetExternalParameters(x,param);
+         fRecV0Info->fT1.fESDTrack.GetExternalCovariance(cov);
+         alpha = fRecV0Info->fT1.fESDTrack.GetAlpha();
+         new (&paramP) AliExternalTrackParam(x,alpha,param,cov);
+         //
+         fRecV0Info->fT2.fESDTrack.GetExternalParameters(x,param);
+         fRecV0Info->fT2.fESDTrack.GetExternalCovariance(cov);
+         alpha = fRecV0Info->fT2.fESDTrack.GetAlpha();
+         new (&paramM) AliExternalTrackParam(x,alpha,param,cov);
+         //
+         fRecV0Info->fV0its.SetParamN(paramM);
+         fRecV0Info->fV0its.SetParamP(paramP);
+         //      fRecV0Info->fV0its.UpdatePID(pid1,pid2);
+         fRecV0Info->fV0its.Update(vertex);
+       }
+      }
+      if (TMath::Abs(fGenV0Info->fMCm.fPdg)==11 &&TMath::Abs(fGenV0Info->fMCd.fPdg)==11){
+       if (fRecV0Info->fDist2>10){
+         fRecV0Info->Update(vertex);
+       }
+       if (fRecV0Info->fDist2>10){
+         fRecV0Info->Update(vertex);
+       }
+      }
+    }   
+    //
+    // take the V0 from reconstruction
+    Int_t label =  TMath::Min(fGenV0Info->fMCm.fLabel,fGenV0Info->fMCd.fLabel);
+    Int_t label2 = TMath::Max(fGenV0Info->fMCm.fLabel,fGenV0Info->fMCd.fLabel);    
+    AliV0 *v0MI=0;
+    fRecV0Info->fRecStatus   =0;
+    fRecV0Info->fMultiple    = fMultiRecV0[label];
+    fRecV0Info->fV0Multiple=0;
+    //
+    if (fMultiRecV0[label]>0 || fMultiRecV0[label2]>0){
+
+      //      for (Int_t j=0;j<TMath::Min(fMultiRecV0s[label],100);j++){
+      for (Int_t j=TMath::Min(fMultiRecV0[label],Short_t(20))-1;j>=0;j--){
+       Int_t index = fIndexRecV0[label*20+j];
+       if (index<0) continue;
+       AliV0 *v0MI2  = (AliV0*)fEvent->GetV0(index);
+//     if (TMath::Abs(v0MI2->fLab[0])==label &&TMath::Abs(v0MI2->fLab[1])==label2) {
+//       v0MI =v0MI2;
+//       fRecV0Info->fV0Multiple++;
+//       fSignedV0[index]=1;
+//     }
+//     if (TMath::Abs(v0MI2->fLab[1])==label &&TMath::Abs(v0MI2->fLab[0])==label2) {
+//       v0MI =v0MI2;
+//       fRecV0Info->fV0Multiple++;
+//       fSignedV0[index]=1;
+//     }
+      }
+    }
+    if (v0MI){
+      fRecV0Info->fV0rec = *v0MI;
+      fRecV0Info->fRecStatus=1;
+    }
+
+    fTreeCmpV0->Fill();
+  }
+  //
+  // write fake v0s
+
+  Int_t nV0MIs = fEvent->GetNumberOfV0s();
+  for (Int_t i=0;i<nV0MIs;i++){
+    if (fSignedV0[i]==0){
+      AliV0 *v0MI  = (AliV0*)fEvent->GetV0(i);
+      if (!v0MI) continue;
+      //
+      fRecV0Info->fV0rec = *v0MI;
+      fRecV0Info->fV0Status  =-10;
+      fRecV0Info->fRecStatus =-2;
+      //
+ //      AliESDRecInfo*  fRecInfo1 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(v0MI->fLab[0]));
+//       AliESDRecInfo*  fRecInfo2 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(v0MI->fLab[1]));
+//       if (fRecInfo1 && fRecInfo2){
+//     fRecV0Info->fT1 = (*fRecInfo1);
+//     fRecV0Info->fT2 = (*fRecInfo2);
+//     fRecV0Info->fRecStatus =-1;
+//       }
+      fRecV0Info->Update(vertex);
+      fTreeCmpV0->Fill();
+    }
+  }
+
+
+
+  fTreeCmpV0->AutoSave();
+  printf("Time spended in BuilV0Info Loop\n");
+  timer.Print();
+  if (fDebug > 2) cerr<<"end of BuildV0Info Loop"<<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/PWG1/AliESDComparisonMI.h b/PWG1/AliESDComparisonMI.h
new file mode 100644 (file)
index 0000000..bf8f797
--- /dev/null
@@ -0,0 +1,233 @@
+
+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 UpdatePoints(AliESDtrack* track);
+  void Update(AliMCInfo* info,AliTPCParam * par, Bool_t reconstructed, AliESD *event);
+  void Reset();
+  Float_t  fTPCPoints[10]; //start , biggest end points,max density .. density at the last 30 pad-rows
+  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
+  AliESDfriendTrack fTrackF;      // friend track
+  AliITStrackMI fITStrack;        //its track
+  AliTRDtrack fTRDtrack;        //its track
+  Float_t fBestTOFmatch;        //best matching between times
+  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
+  Int_t  fStatus[4];        // status -0 not found - 1 -only in - 2 -in-out -3 -in -out-refit
+  Bool_t fITSOn;           // ITS refitted inward
+  Bool_t fTRDOn;           // ITS refitted inward
+  Float_t fDeltaP;          //delta of momenta
+  Double_t fSign;           // sign
+  Int_t fLabels[2];         // labels
+  ClassDef(AliESDRecInfo,2)  // 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       fRs[2];     // minimum radius in rphi intersection
+  Double_t       fDistMinR; // distance at minimal radius
+  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
+  Int_t          fV0Status;       // status of the kink
+  AliV0          fV0tpc;           // Vo information from reconsturction according TPC
+  AliV0          fV0its;           // Vo information from reconsturction according ITS
+  AliV0          fV0rec;           // V0 information form the reconstruction
+  Int_t          fMultiple;
+  Int_t          fV0Multiple;
+  Int_t          fRecStatus;    // status form the reconstuction
+  ClassDef(AliESDRecV0Info,2)   // container for  
+};
+
+ClassImp(AliESDRecV0Info)
+
+
+class AliESDRecKinkInfo: public TObject {
+public:
+  void Update();
+  AliESDRecInfo  fT1;      //track1
+  AliESDRecInfo  fT2;      //track2  
+  AliESDkink     fKink;    //kink
+  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 
+  Double_t       fMinR;     // minimum radius in rphi intersection
+  Double_t       fDistMinR; // distance at minimal radius
+  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
+  Int_t          fStatus;       //status -tracks 
+  Int_t          fRecStatus;    //kink -status- 0 - not found  1-good -  fake
+  Int_t          fMultiple;
+  Int_t          fKinkMultiple;
+  ClassDef(AliESDRecKinkInfo,1)   // container for  
+};
+
+ClassImp(AliESDRecKinkInfo)
+
+
+
+////////////////////////////////////////////////////////////////////////
+// 
+// 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();
+  Int_t BuildKinkInfo0(Int_t eventNr); // build kink info 0
+  Int_t BuildV0Info(Int_t eventNr); // build kink info 0
+  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
+  TTree *fTreeCmpKinks;                //! output tree with cmp Kinks
+  TTree *fTreeCmpV0;                //! output tree with cmp V0
+  //
+  char  fFnGenTracks[1000];             //! input file name with gen tracks
+  TFile *fFileGenTracks;
+  TTree *fTreeGenTracks;
+  TTree *fTreeGenKinks;            // tree with gen kinks
+  TTree *fTreeGenV0;            // tree with gen V0
+  //
+  //
+  Int_t  fDirection;
+  //
+  AliRunLoader * fLoader;         //! pointer to the run loader
+  //TTree *fTreeRecTracks;          //! tree with reconstructed tracks
+  //
+  Short_t *fIndexRecTracks;         //! index of particle label in the TreeT_ESD
+  Short_t *fFakeRecTracks;          //! number of fake tracks
+  Short_t *fMultiRecTracks;         //! number of multiple reconstructions
+  //
+  Short_t *fIndexRecKinks;         //! index of particle label in treeesd
+  Short_t *fMultiRecKinks;         //! number of multiple reconstructions
+  Short_t *fSignedKinks;           //! indicator that kink was not fake
+  //
+  Short_t *fIndexRecV0;         //! index of particle label in treeesd
+  Short_t *fMultiRecV0;         //! number of multiple reconstructions
+  Short_t *fSignedV0;                //! indicator that kink was not fake
+  //
+  TObjArray *fRecArray;           // container with rec infos
+  AliESD *fEvent;                 //!event
+  AliESDfriend *fESDfriend;              //!event friend
+  //
+  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
+  Int_t fNextKinkToRead;            //! last entry already read from genKinks tree
+  Int_t fNextV0ToRead;            //! last entry already read from genV0 tree
+  //
+  AliMCInfo*  fMCInfo;           //! MC information writen per particle
+  AliGenKinkInfo* fGenKinkInfo;      //! MC information writen per Kink
+  AliGenV0Info* fGenV0Info;      //! MC information writen per Kink
+  AliESDRecInfo*  fRecInfo;          //! Rec. information writen per particle
+  AliESDfriendTrack*  fFriend;          //! friend track
+  AliESDRecKinkInfo* fRecKinkInfo;    //! reconstructed kink info
+  AliESDRecV0Info* fRecV0Info;    //! reconstructed kink info
+  //
+
+  ClassDef(ESDCmpTr,1)    // class which creates and fills tree with ESDGenTrack objects
+};
+ClassImp(ESDCmpTr)
+
+
+
diff --git a/PWG1/AliGenInfo.cxx b/PWG1/AliGenInfo.cxx
new file mode 100644 (file)
index 0000000..0214cf0
--- /dev/null
@@ -0,0 +1,1360 @@
+/**************************************************************************
+ * 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
+
+Generate complex MC information - used for Comparison later on
+How to use it?
+
+.L 
+AliGenInfoMaker *t = new AliGenInfoMaker("galice.root","genTracks.root",0,0)
+t->Exec();
+
+*/
+
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include <stdio.h>
+#include <string.h>
+//ROOT includes
+#include "TROOT.h"
+#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"
+#include "TView.h"
+#include "TView3D.h"
+#include "TGeometry.h"
+#include "TPolyLine3D.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"
+#include "AliHelix.h"
+#include "AliPoints.h"
+#include "AliTrackPointArray.h"
+
+#endif
+#include "AliGenInfo.h" 
+//
+// 
+
+ClassImp(AliTPCdigitRow);
+ClassImp(AliMCInfo);
+ClassImp(AliGenV0Info)
+ClassImp(AliGenKinkInfo)
+ClassImp(AliGenInfoMaker)
+
+
+
+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);
+  fTOFReferences  = new TClonesArray("AliTrackReference",10);
+  fTRdecay.SetTrack(-1);
+  fCharge = 0;
+}
+
+AliMCInfo::~AliMCInfo()
+{
+  if (fTPCReferences) {
+    delete fTPCReferences;
+  }
+  if (fITSReferences){
+    delete fITSReferences;
+  }
+  if (fTRDReferences){    
+    delete fTRDReferences;  
+  }
+  if (fTOFReferences){    
+    delete fTOFReferences;  
+  }
+
+}
+
+
+
+void AliMCInfo::Update()
+{
+  //
+  //
+  fMCtracks =1;
+  if (!fTPCReferences) {
+    fNTPCRef =0;
+    return;
+  }
+  Float_t direction=1;
+  //Float_t rlast=0;
+  fNTPCRef = fTPCReferences->GetEntriesFast();
+  fNITSRef = fITSReferences->GetEntriesFast();
+  fNTRDRef = fTRDReferences->GetEntriesFast();
+  fNTOFRef = fTOFReferences->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();
+  //  
+}
+*/
+
+void AliGenV0Info::Update(Float_t vertex[3])
+{
+  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();
+  //
+  //
+  //
+  Double_t x1[3],p1[3];
+  TParticle & pdaughter = fMCd.fParticle;
+  x1[0] = pdaughter.Vx();      
+  x1[1] = pdaughter.Vy();
+  x1[2] = pdaughter.Vz();
+  p1[0] = pdaughter.Px();
+  p1[1] = pdaughter.Py();
+  p1[2] = pdaughter.Pz();
+  Double_t sign = (pdaughter.GetPDG()->Charge()>0)? -1:1;
+  AliHelix dhelix1(x1,p1,sign);
+  //
+  //
+  Double_t x2[3],p2[3];            
+  //
+  TParticle & pmother = fMCm.fParticle;
+  x2[0] = pmother.Vx();      
+  x2[1] = pmother.Vy();      
+  x2[2] = pmother.Vz();      
+  p2[0] = pmother.Px();
+  p2[1] = pmother.Py();
+  p2[2] = pmother.Pz();
+  //
+  //
+  sign = (pmother.GetPDG()->Charge()>0) ? -1:1;
+  AliHelix mhelix(x2,p2,sign);
+  
+  //
+  //
+  //
+  //find intersection linear
+  //
+  Double_t distance1, distance2;
+  Double_t phase[2][2],radius[2];
+  Int_t  points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
+  Double_t delta1=10000,delta2=10000;    
+  if (points>0){
+    dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
+    dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
+    dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
+  }
+  else{
+    fInvMass=-1;
+    return;
+  }
+  if (points==2){    
+    dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
+    dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
+    dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
+  }
+  distance1 = TMath::Min(delta1,delta2);
+  //
+  //find intersection parabolic
+  //
+  points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
+  delta1=10000,delta2=10000;  
+  
+  if (points>0){
+    dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
+  }
+  if (points==2){    
+    dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
+  }
+  
+  distance2 = TMath::Min(delta1,delta2);
+  //
+  if (delta1<delta2){
+    //get V0 info
+    dhelix1.Evaluate(phase[0][0],fMCXr);
+    dhelix1.GetMomentum(phase[0][0],fMCPdr);
+    mhelix.GetMomentum(phase[0][1],fMCPm);
+    dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],fMCAngle);
+    fMCRr = TMath::Sqrt(radius[0]);
+  }
+  else{
+    dhelix1.Evaluate(phase[1][0],fMCXr);
+    dhelix1.GetMomentum(phase[1][0],fMCPdr);
+    mhelix.GetMomentum(phase[1][1],fMCPm);
+    dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],fMCAngle);
+    fMCRr = TMath::Sqrt(radius[1]);
+  }     
+  //            
+  //
+  fMCDist1 = TMath::Sqrt(distance1);
+  fMCDist2 = TMath::Sqrt(distance2);      
+  //
+  //
+  // 
+  Float_t v[3] = {fMCXr[0]-vertex[0],fMCXr[1]-vertex[1],fMCXr[2]-vertex[2]};
+  //Float_t v[3] = {fMCXr[0],fMCXr[1],fMCXr[2]};
+  Float_t p[3] = {fMCPdr[0]+fMCPm[0], fMCPdr[1]+fMCPm[1],fMCPdr[2]+fMCPm[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);
+  Double_t mass1 = fMCd.fMass;
+  Double_t mass2 = fMCm.fMass;
+
+  
+  Double_t e1    = TMath::Sqrt(mass1*mass1+
+                             fMCPd[0]*fMCPd[0]+
+                             fMCPd[1]*fMCPd[1]+
+                             fMCPd[2]*fMCPd[2]);
+  Double_t e2    = TMath::Sqrt(mass2*mass2+
+                             fMCPm[0]*fMCPm[0]+
+                             fMCPm[1]*fMCPm[1]+
+                             fMCPm[2]*fMCPm[2]);
+  
+  fInvMass =  
+    (fMCPm[0]+fMCPd[0])*(fMCPm[0]+fMCPd[0])+
+    (fMCPm[1]+fMCPd[1])*(fMCPm[1]+fMCPd[1])+
+    (fMCPm[2]+fMCPd[2])*(fMCPm[2]+fMCPd[2]);
+  
+  //  fInvMass = TMath::Sqrt((e1+e2)*(e1+e2)-fInvMass);
+  fInvMass = (e1+e2)*(e1+e2)-fInvMass;
+  if (fInvMass>0) fInvMass = TMath::Sqrt(fInvMass);
+
+    
+}
+
+
+
+/////////////////////////////////////////////////////////////////////////////////
+void AliGenKinkInfo::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();
+  //
+  //
+  //
+  Double_t x1[3],p1[3];
+  TParticle & pdaughter = fMCd.fParticle;
+  x1[0] = pdaughter.Vx();      
+  x1[1] = pdaughter.Vy();
+  x1[2] = pdaughter.Vz();
+  p1[0] = pdaughter.Px();
+  p1[1] = pdaughter.Py();
+  p1[2] = pdaughter.Pz();
+  Double_t sign = (pdaughter.GetPDG()->Charge()>0)? -1:1;
+  AliHelix dhelix1(x1,p1,sign);
+  //
+  //
+  Double_t x2[3],p2[3];            
+  //
+  TParticle & pmother = fMCm.fParticle;
+  x2[0] = pmother.Vx();      
+  x2[1] = pmother.Vy();      
+  x2[2] = pmother.Vz();      
+  p2[0] = pmother.Px();
+  p2[1] = pmother.Py();
+  p2[2] = pmother.Pz();
+  //
+  AliTrackReference & pdecay = fMCm.fTRdecay;
+  x2[0] = pdecay.X();      
+  x2[1] = pdecay.Y();      
+  x2[2] = pdecay.Z();      
+  p2[0] = pdecay.Px();
+  p2[1] = pdecay.Py();
+  p2[2] = pdecay.Pz();
+  //
+  sign = (pmother.GetPDG()->Charge()>0) ? -1:1;
+  AliHelix mhelix(x2,p2,sign);
+  
+  //
+  //
+  //
+  //find intersection linear
+  //
+  Double_t distance1, distance2;
+  Double_t phase[2][2],radius[2];
+  Int_t  points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
+  Double_t delta1=10000,delta2=10000;    
+  if (points>0){
+    dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
+    dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
+    dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
+  }
+  if (points==2){    
+    dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
+    dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
+    dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
+  }
+  distance1 = TMath::Min(delta1,delta2);
+  //
+  //find intersection parabolic
+  //
+  points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
+  delta1=10000,delta2=10000;  
+  
+  if (points>0){
+    dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
+  }
+  if (points==2){    
+    dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
+  }
+  
+  distance2 = TMath::Min(delta1,delta2);
+  //
+  if (delta1<delta2){
+    //get V0 info
+    dhelix1.Evaluate(phase[0][0],fMCXr);
+    dhelix1.GetMomentum(phase[0][0],fMCPdr);
+    mhelix.GetMomentum(phase[0][1],fMCPm);
+    dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],fMCAngle);
+    fMCRr = TMath::Sqrt(radius[0]);
+  }
+  else{
+    dhelix1.Evaluate(phase[1][0],fMCXr);
+    dhelix1.GetMomentum(phase[1][0],fMCPdr);
+    mhelix.GetMomentum(phase[1][1],fMCPm);
+    dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],fMCAngle);
+    fMCRr = TMath::Sqrt(radius[1]);
+  }     
+  //            
+  //
+  fMCDist1 = TMath::Sqrt(distance1);
+  fMCDist2 = TMath::Sqrt(distance2);      
+    
+}
+
+
+Float_t AliGenKinkInfo::GetQt(){
+  //
+  //
+  Float_t momentumd = TMath::Sqrt(fMCPd[0]*fMCPd[0]+fMCPd[1]*fMCPd[1]+fMCPd[2]*fMCPd[2]);
+  return TMath::Sin(fMCAngle[2])*momentumd;
+}
+
+
+
+  
+////////////////////////////////////////////////////////////////////////
+
+////////////////////////////////////////////////////////////////////////
+//
+// End of implementation of the class AliMCInfo
+//
+////////////////////////////////////////////////////////////////////////
+
+
+
+////////////////////////////////////////////////////////////////////////
+AliTPCdigitRow::AliTPCdigitRow()
+{
+  Reset();
+}
+////////////////////////////////////////////////////////////////////////
+AliTPCdigitRow & AliTPCdigitRow::operator=(const AliTPCdigitRow &digOld)
+{
+  for (Int_t i = 0; i<kgRowBytes; i++) fDig[i] = digOld.fDig[i];
+  return (*this);
+}
+////////////////////////////////////////////////////////////////////////
+void AliTPCdigitRow::SetRow(Int_t row) 
+{
+  if (row >= 8*kgRowBytes) {
+    cerr<<"AliTPCdigitRow::SetRow: index "<<row<<" out of bounds."<<endl;
+    return;
+  }
+  Int_t iC = row/8;
+  Int_t iB = row%8;
+  SETBIT(fDig[iC],iB);
+}
+
+////////////////////////////////////////////////////////////////////////
+Bool_t AliTPCdigitRow::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 AliTPCdigitRow::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 AliTPCdigitRow::Reset()
+{
+//
+// resets all rows to zero
+//
+  for (Int_t i = 0; i<kgRowBytes; i++) {
+    fDig[i] <<= 8;
+  }
+}
+////////////////////////////////////////////////////////////////////////
+Int_t AliTPCdigitRow::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 AliTPCdigitRow::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 AliTPCdigitRow
+//
+////////////////////////////////////////////////////////////////////////
+  
+////////////////////////////////////////////////////////////////////////
+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,0);
+}
+
+
+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();
+  fLoader->LoadHits();
+  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;
+    //
+    if (BuildKinkInfo()>0) return 1;
+    if (BuildV0Info()>0) return 1;
+    //if (BuildHitLines()>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; 
+  //
+  AliGenKinkInfo *kinkinfo = new AliGenKinkInfo;
+  fTreeKinks = new TTree("genKinksTree","genKinksTree");  
+  fTreeKinks->Branch("MC","AliGenKinkInfo",&kinkinfo);
+  delete kinkinfo;
+  //
+  AliGenV0Info *v0info = new AliGenV0Info;
+  fTreeV0 = new TTree("genV0Tree","genV0Tree");  
+  fTreeV0->Branch("MC","AliGenV0Info",&v0info);
+  delete v0info;
+  //
+  //
+  AliTrackPointArray * points0 = new AliTrackPointArray;  
+  AliTrackPointArray * points1 = new AliTrackPointArray;  
+  AliTrackPointArray * points2 = new AliTrackPointArray;  
+  fTreeHitLines = new TTree("HitLines","HitLines");
+  fTreeHitLines->Branch("TPC.","AliTrackPointArray",&points0,32000,10);
+  fTreeHitLines->Branch("TRD.","AliTrackPointArray",&points1,32000,10);
+  fTreeHitLines->Branch("ITS.","AliTrackPointArray",&points2,32000,10);
+  Int_t label=0;
+  fTreeHitLines->Branch("Label",&label,"label/I");
+  //
+  fTreeGenTracks->AutoSave();
+  fTreeKinks->AutoSave();
+  fTreeV0->AutoSave();
+  fTreeHitLines->AutoSave();
+}
+////////////////////////////////////////////////////////////////////////
+void AliGenInfoMaker::CloseOutputFile() 
+{
+  if (!fFileGenTracks) {
+    cerr<<"File "<<fFnRes<<" not found as an open file."<<endl;
+    return;
+  }
+  fFileGenTracks->cd();
+  fTreeGenTracks->Write();  
+  delete fTreeGenTracks;
+  fTreeKinks->Write();
+  delete fTreeKinks;
+  fTreeV0->Write();
+  delete fTreeV0;
+
+  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::BuildKinkInfo()
+{
+  //
+  // Build tree with Kink Information
+  //
+  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 = fTreeKinks->GetBranch("MC");
+  //
+  AliGenKinkInfo * kinkinfo = new AliGenKinkInfo;
+  //
+  for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
+    // load only particles with TR
+    AliMCInfo * info = GetInfo(iParticle);
+    if (!info) continue;
+    if (info->fCharge==0) continue;  
+    if (info->fTRdecay.P()<0.13) continue;  //momenta cut 
+    if (info->fTRdecay.R()>500)  continue;  //R cut - decay outside barrel
+    TParticle & particle = info->fParticle;
+    Int_t first = particle.GetDaughter(0);
+    if (first ==0) continue;
+
+    Int_t last = particle.GetDaughter(1);
+    if (last ==0) last=first;
+    AliMCInfo * info2 =  0;
+    AliMCInfo * dinfo =  0;
+    
+    for (Int_t id2=first;id2<=last;id2++){
+      info2 = GetInfo(id2);
+      if (!info2) continue;
+      TParticle &p2 = info2->fParticle;
+      Double_t r = TMath::Sqrt(p2.Vx()*p2.Vx()+p2.Vy()*p2.Vy());
+      if ( TMath::Abs(info->fTRdecay.R()-r)>0.01) continue;
+      if (!(p2.GetPDG())) continue;
+      dinfo =info2;
+     
+      kinkinfo->fMCm = (*info);
+      kinkinfo->fMCd = (*dinfo);
+      if (kinkinfo->fMCm.fParticle.GetPDG()==0 ||  kinkinfo->fMCd.fParticle.GetPDG()==0) continue;
+      kinkinfo->Update();
+      br->SetAddress(&kinkinfo);    
+      fTreeKinks->Fill();
+    }
+    /*
+    if (dinfo){
+      kinkinfo->fMCm = (*info);
+      kinkinfo->fMCd = (*dinfo);
+      kinkinfo->Update();
+      br->SetAddress(&kinkinfo);    
+      fTreeKinks->Fill();
+    }
+    */
+  }
+  fTreeGenTracks->AutoSave();
+  if (fDebug > 2) cerr<<"end of Kink Loop"<<endl;
+  return 0;
+}
+
+
+////////////////////////////////////////////////////////////////////////////////////
+Int_t  AliGenInfoMaker::BuildV0Info()
+{
+  //
+  // Build tree with V0 Information
+  //
+  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 = fTreeV0->GetBranch("MC");
+  //
+  AliGenV0Info * v0info = new AliGenV0Info;
+  //
+  //
+  for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
+    // load only particles with TR
+    AliMCInfo * info = GetInfo(iParticle);
+    if (!info) continue;
+    if (info->fCharge==0) continue;  
+    //
+    //
+    TParticle & particle = info->fParticle;
+    Int_t mother = particle.GetMother(0);
+    if (mother <=0) continue;
+    //
+    TParticle * motherparticle = stack->Particle(mother);
+    if (!motherparticle) continue;
+    //
+    Int_t last = motherparticle->GetDaughter(1);
+    if (last==(int)iParticle) continue;
+    AliMCInfo * info2 =  info;
+    AliMCInfo * dinfo =  GetInfo(last);
+    if (!dinfo) continue;
+    if (!dinfo->fParticle.GetPDG()) continue;
+    if (!info2->fParticle.GetPDG()) continue;
+    if (dinfo){
+      v0info->fMCm = (*info);
+      v0info->fMCd = (*dinfo);
+      v0info->fMotherP = (*motherparticle);
+      v0info->Update(fVPrim);
+      br->SetAddress(&v0info);    
+      fTreeV0->Fill();
+    }
+  }
+  fTreeV0->AutoSave();
+  if (fDebug > 2) cerr<<"end of V0 Loop"<<endl;
+  return 0;
+}
+
+
+
+
+////////////////////////////////////////////////////////////////////////
+Int_t AliGenInfoMaker::BuildHitLines()
+{
+
+//
+// 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
+//   AliTrackPointArray *tpcp = new AliTrackPointArray;
+//   AliTrackPointArray *trdp = new AliTrackPointArray;
+//   AliTrackPointArray *itsp = new AliTrackPointArray;
+//   Int_t label =0;
+//   //
+//   TBranch * brtpc = fTreeHitLines->GetBranch("TPC.");
+//   TBranch * brtrd = fTreeHitLines->GetBranch("TRD.");  
+//   TBranch * brits = fTreeHitLines->GetBranch("ITS.");
+//   TBranch * brlabel = fTreeHitLines->GetBranch("Label");
+//   brlabel->SetAddress(&label);
+//   brtpc->SetAddress(&tpcp);
+//   brtrd->SetAddress(&trdp);
+//   brits->SetAddress(&itsp);
+//   //
+//   AliDetector *dtpc = gAlice->GetDetector("TPC");
+//   AliDetector *dtrd = gAlice->GetDetector("TRD");
+//   AliDetector *dits = gAlice->GetDetector("ITS");
+//   for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
+//     // load only particles with TR
+//     AliMCInfo * info = GetInfo(iParticle);
+//     if (!info) continue;
+//     Int_t primpart = info->fPrimPart;
+//     ipdg = info->fParticle.GetPdgCode();
+//     label = iParticle;
+//     //
+//     gAlice->ResetHits();
+//     tpcp->Reset();
+//     itsp->Reset();
+//     trdp->Reset();
+//     tpcp->fLabel1 = ipdg;
+//     trdp->fLabel1 = ipdg;
+//     itsp->fLabel1 = ipdg;
+//     if (dtpc->TreeH()->GetEvent(primpart)){
+//       dtpc->LoadPoints(primpart);
+//       tpcp->Reset(dtpc,iParticle);
+//     }
+//     if (dtrd->TreeH()->GetEvent(primpart)){
+//       dtrd->LoadPoints(primpart);
+//       trdp->Reset(dtrd,iParticle);
+//     }
+//     if (dits->TreeH()->GetEvent(primpart)){
+//       dits->LoadPoints(primpart);
+//       itsp->Reset(dits,iParticle);
+//     }    
+//     //    
+//     fTreeHitLines->Fill();
+//     dtpc->ResetPoints();
+//     dtrd->ResetPoints();
+//     dits->ResetPoints();
+//   }
+//   delete tpcp;
+//   delete trdp;
+//   delete itsp;
+//   fTreeHitLines->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* TOFArrayTR = 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("TOF"))    treeTR->GetBranch("TOF")->SetAddress(&TOFArrayTR);
+  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->TestBit(BIT(2))){  
+       //if decay 
+       if (trackRef->P()<fgTPCPtCut) continue;
+       Int_t label = trackRef->GetTrack(); 
+       AliMCInfo * info = GetInfo(label);
+       if (!info) info = MakeInfo(label);
+       info->fTRdecay = *trackRef;      
+      }
+      //
+      if (trackRef->P()<fgTPCPtCut) continue;
+      Int_t label = trackRef->GetTrack();      
+      AliMCInfo * info = GetInfo(label);
+      if (!info) info = MakeInfo(label);
+      if (!info) continue;
+      info->fPrimPart =  iPrimPart;
+      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->P()<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;
+      info->fPrimPart =  iPrimPart;
+      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->P()<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;
+      info->fPrimPart =  iPrimPart;
+      TClonesArray & arr = *(info->fTRDReferences);
+      new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);     
+    }
+    //
+    // Loop over TOF references
+    //
+    for (Int_t iTrackRef = 0; iTrackRef < TOFArrayTR->GetEntriesFast(); iTrackRef++) {
+      AliTrackReference *trackRef = (AliTrackReference*)TOFArrayTR->At(iTrackRef);            
+      Int_t label = trackRef->GetTrack();      
+      AliMCInfo * info = GetInfo(label);
+      if (!info){
+       if (trackRef->Pt()<fgTPCPtCut) continue;
+       if ( (!info) && trackRef->Pt()<fgTOFPtCut) continue;
+      }
+      if (!info) info = MakeInfo(label);
+      if (!info) continue;
+      info->fPrimPart =  iPrimPart;
+      TClonesArray & arr = *(info->fTOFReferences);
+      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;
+  TOFArrayTR->Delete();
+  delete  TOFArrayTR;
+
+  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];
+}
+////////////////////////////////////////////////////////////////////////
+
+
+
diff --git a/PWG1/AliGenInfo.h b/PWG1/AliGenInfo.h
new file mode 100644 (file)
index 0000000..70d9dbc
--- /dev/null
@@ -0,0 +1,237 @@
+#ifndef ALIGENINFO_H
+#define ALIGENINFO_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+
+
+//////////////////////////////////////////////////////////////////////////////
+//                          Class AliGenInfo                               //
+//   collect together MC info for comparison purposes - effieciency studies and so on//                                                                 //
+//   marian.ivanov@cern.ch                                                  //
+//////////////////////////////////////////////////////////////////////////////
+
+
+
+////////////////////////////////////////////////////////////////////////
+//
+// Start of implementation of the class AliTPCdigitRow
+//
+////////////////////////////////////////////////////////////////////////
+
+#include <TParticle.h>
+#include "AliTrackReference.h"
+
+class TFile;
+class AliRunLoader;
+class AliStack;
+class AliTPCParam;
+
+const Int_t kgRowBytes = 32;
+
+class AliTPCdigitRow: public TObject {
+public:
+  AliTPCdigitRow();
+  virtual ~AliTPCdigitRow(){;}
+  void SetRow(Int_t row);
+  Bool_t TestRow(Int_t row);
+  AliTPCdigitRow & operator=(const AliTPCdigitRow &digOld);
+  Int_t RowsOn(Int_t upto=8*kgRowBytes);
+  Int_t Last();
+  Int_t First();
+  void Reset();
+
+//private:
+  UChar_t fDig[kgRowBytes];
+  ClassDef(AliTPCdigitRow,1)  // container for digit pattern
+};
+
+
+////////////////////////////////////////////////////////////////////////
+//
+// 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
+  //
+  Int_t     fPrimPart;               // index of primary particle in TreeH
+  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
+  AliTPCdigitRow fTPCRow;                  // information about digits row pattern
+  Int_t fNTPCRef;                    // tpc references counter
+  Int_t fNITSRef;                    // ITS references counter
+  Int_t fNTRDRef;                    // TRD references counter
+  Int_t fNTOFRef;                    // TOF references counter
+  TClonesArray * fTPCReferences;     //containner with all track references -in the TPC
+  TClonesArray * fITSReferences;     //container with ITS references
+  TClonesArray * fTRDReferences;     //container with TRD references  
+  TClonesArray * fTOFReferences;     //container with TRD references  
+  //
+  ClassDef(AliMCInfo,1)  // container for 
+};
+
+
+
+class AliGenV0Info: public TObject {
+public:
+  AliMCInfo fMCd;                   //info about daughter particle - second particle for V0
+  AliMCInfo fMCm;      //info about mother particle   - first particle for V0
+  TParticle fMotherP;   //particle info about mother particle
+  void Update(Float_t vertex[3]);        // 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  
+  //
+  Double_t       fInvMass;  //reconstructed invariant mass -
+  Float_t        fPointAngleFi; //point angle fi
+  Float_t        fPointAngleTh; //point angle theta
+  Float_t        fPointAngle;   //point angle full
+  //
+  ClassDef(AliGenV0Info,1)  // container for  
+};
+
+
+
+class AliGenKinkInfo: public TObject {
+public:
+  AliMCInfo fMCd;            //info about daughter particle - second particle for V0
+  AliMCInfo fMCm;            //info about mother particle   - first particle for V0
+  void Update();             // put some derived info to special field 
+  Float_t GetQt();           //
+  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(AliGenKinkInfo,1) // container for  
+};
+
+
+
+////////////////////////////////////////////////////////////////////////
+// 
+// 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();
+  Int_t BuildKinkInfo();  // build information about MC kinks
+  Int_t BuildV0Info();  // build information about MC kinks
+  Int_t BuildHitLines();  // build information about MC kinks
+  //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
+  TTree * fTreeKinks;             //!  output tree with Kinks
+  TTree * fTreeV0;                //!  output tree with V0
+  TTree * fTreeHitLines;          //! tree with hit lines
+  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
+  Float_t fVPrim[3];             //! primary vertex position
+                                  // the fVDist[3] contains size of the 3-vector
+
+private:
+  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.03; // 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
+  static const Double_t fgTOFPtCut = 0.15; // do not store particles with generated pT less than this
+  ClassDef(AliGenInfoMaker,1)    // class which creates and fills tree with TPCGenTrack objects
+};
+
+
+
+
+
+AliTPCParam * GetTPCParam();
+Float_t TPCBetheBloch(Float_t bg);
+
+#endif
diff --git a/PWG1/AliTreeDraw.cxx b/PWG1/AliTreeDraw.cxx
new file mode 100644 (file)
index 0000000..9592e76
--- /dev/null
@@ -0,0 +1,387 @@
+/**************************************************************************
+ * 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
+Frequenlty used function for visualization 
+marian.ivanov@cern.ch
+*/
+
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include <stdio.h>
+#include <string.h>
+//ROOT includes
+#include "TROOT.h"
+#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"
+#include "TView.h"
+#include "TView3D.h"
+#include "TGeometry.h"
+#include "TPolyLine3D.h"
+#include "TPolyMarker3D.h"
+
+//ALIROOT includes
+#include "AliTrackPointArray.h"
+#include "AliTreeDraw.h" 
+
+#endif
+
+//
+//     Class for visualization and some statistacal analysis using tree
+//     To be used in comparisons
+//                and calib viewers based on tree    
+
+
+ClassImp(AliTreeDraw)
+
+
+AliTreeDraw::AliTreeDraw():
+  fTree(0),
+  fRes(0),
+  fMean(0),
+  fPoints(0){
+  //
+  // default constructor
+  //
+}
+
+void  AliTreeDraw::ClearHisto(){
+  //
+  //
+  delete fRes; 
+  delete fMean;
+  fRes=0;
+  fMean=0;
+}
+
+
+
+TH1F * AliTreeDraw::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, Int_t nBinsRes)
+{
+  //
+  Double_t* bins = CreateLogBins(nbins, minx, maxx);
+  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;
+  ClearHisto();
+  fRes  = hRes;
+  fMean = hMean;
+  return hRes;
+}
+
+
+
+TH1F * AliTreeDraw::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, Int_t nBinsRes)
+{
+  //
+  // 
+  //
+  Double_t* bins = CreateLogBins(nbins, minx, maxx);
+  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;
+  ClearHisto();
+  fRes  = hRes;
+  fMean = hMean;
+  return hRes;
+}
+
+///////////////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////////////
+TH1F * AliTreeDraw::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 * AliTreeDraw::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* AliTreeDraw::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 AliTreeDraw::AliLabelAxes(TH1* histo, const char* xAxisTitle, const char* yAxisTitle)
+{
+  //
+  histo->GetXaxis()->SetTitle(xAxisTitle);
+  histo->GetYaxis()->SetTitle(yAxisTitle);
+}
+
+
+TH1F* AliTreeDraw::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* AliTreeDraw::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
+  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;
+}
+
+
+
+
+void   AliTreeDraw::GetPoints3D(const char * label, const char * chpoints, 
+                               const char* selection, TTree * tpoints, Int_t color,Float_t rmin){
+  //
+  // load selected points from tree
+  //
+   if (!fPoints) fPoints = new TObjArray;
+   if (tpoints->GetIndex()==0) tpoints->BuildIndex("fLabel","Label");
+   TBranch * br = tpoints->GetBranch(chpoints);
+   if (!br) return;
+   AliTrackPointArray * points = new AliTrackPointArray;
+   br->SetAddress(&points);
+
+   Int_t npoints = fTree->Draw(label,selection);
+   Float_t xyz[30000];
+   rmin*=rmin;
+   for (Int_t i=0;i<npoints;i++){
+     Int_t index = (Int_t)fTree->GetV1()[i];
+     tpoints->GetEntryWithIndex(index,index);
+     if (points->GetNPoints()<2) continue;
+     Int_t goodpoints=0;
+     for (Int_t i=0;i<points->GetNPoints(); i++){
+       xyz[goodpoints*3]   = points->GetX()[i];
+       xyz[goodpoints*3+1] = points->GetY()[i];
+       xyz[goodpoints*3+2] = points->GetZ()[i];
+       if ( points->GetX()[i]*points->GetX()[i]+points->GetY()[i]*points->GetY()[i]>rmin) goodpoints++;
+     } 
+     TPolyMarker3D * marker = new TPolyMarker3D(goodpoints,xyz); 
+     marker->SetMarkerColor(color);
+     marker->SetMarkerStyle(1);
+     fPoints->AddLast(marker);
+   }
+}
+
diff --git a/PWG1/AliTreeDraw.h b/PWG1/AliTreeDraw.h
new file mode 100644 (file)
index 0000000..9ea3b81
--- /dev/null
@@ -0,0 +1,62 @@
+#ifndef ALITREEDRAW_H
+#define ALITREEDRAW_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+
+
+//////////////////////////////////////////////////////////////////////////////
+//                          Class AliTreeDraw                               //
+//                                                                          //
+//   marian.ivanov@cern.ch                                                  //
+//////////////////////////////////////////////////////////////////////////////
+
+
+
+#include <TObject.h>
+#include <TObjArray.h>
+
+class TH1;
+class TH1F;
+class TH2F;
+class TTree;
+
+class AliTreeDraw: public TObject{
+public:
+  AliTreeDraw();
+  ~AliTreeDraw(){;}
+  void SetTree(TTree *tree){fTree=tree;}
+  const TH1 * GetRes() const{ return (TH1*)fRes;}
+  const TH1 * GetMean() const{ return (TH1*)fMean;}
+  const TObjArray* GetPoints() const {return fPoints;}
+  void  ClearHisto();
+  void  ClearPoints(){if (fPoints) fPoints->Clear();}
+  //
+  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, Int_t nBinsRes=30);
+  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, Int_t nBinsRes=30); 
+  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);
+  //
+  void   GetPoints3D(const char * label, const char * chpoints, const char* selection, TTree * tpoints, Int_t color=6, Float_t rmin=4.);
+
+ private: 
+  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);
+
+  TTree * fTree;    //the tree for visualization - NOT OWNER
+  TH1F  * fRes;     //temporary histogram        - OWNER  
+  TH1F  * fMean;    //temporary histogram        - OWNER
+  TObjArray *fPoints;//                          - OWNER
+  ClassDef(AliTreeDraw,0)
+};
+
+#endif
diff --git a/PWG1/PWG1LinkDef.h b/PWG1/PWG1LinkDef.h
new file mode 100644 (file)
index 0000000..212156b
--- /dev/null
@@ -0,0 +1,16 @@
+#ifdef __CINT__
+
+#pragma link off all glols;
+#pragma link off all classes;
+#pragma link off all functions;
+
+#pragma link C++ class AliTreeDraw+;
+#pragma link C++ class AliTPCdigitRow+;
+#pragma link C++ class AliGenV0Info+;
+#pragma link C++ class AliGenKinkInfo+;
+#pragma link C++ class AliGenInfoMaker+;
+#pragma link C++ class AliMCInfo+;
+
+
+
+#endif
diff --git a/PWG1/libPWG1.pkg b/PWG1/libPWG1.pkg
new file mode 100644 (file)
index 0000000..0626586
--- /dev/null
@@ -0,0 +1,9 @@
+
+SRCS:= AliTreeDraw.cxx \
+       AliGenInfo.cxx
+
+HDRS:= $(SRCS:.cxx=.h)
+
+DHDR:= PWG1LinkDef.h
+
+EINCLUDE:= TPC
diff --git a/PWG1/rootlogon.C b/PWG1/rootlogon.C
new file mode 100644 (file)
index 0000000..21171f1
--- /dev/null
@@ -0,0 +1,20 @@
+
+{
+  printf("\nWELCOME to ALICE\n\n"); 
+  gSystem->SetIncludePath("-I$ROOTSYS/include -I$ALICE_ROOT/ -I$ALICE_ROOT/include -I$ALICE_ROOT/STEER\
+   -I$ALICE_ROOT/TPC -I$ALICE_ROOT/ITS -I$ALICE_ROOT/TRD -I$ALICE_ROOT/TOF -I$ALICE_ROOT/RAW -I$ALICE_ROOT/kdtree/include");
+  //load library neccesary for TPC alone 
+  gSystem->Load("$(ROOTSYS)/lib/libEG");
+  gSystem->Load("$(ROOTSYS)/lib/libGeom");  
+  gSystem->Load("$(ROOTSYS)/lib/libVMC");
+  gSystem->Load("$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET)/libSTEER");
+  gSystem->Load("$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET)/libTPCrec");
+  gSystem->Load("$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET)/libTPCsim");
+  gSystem->Load("$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET)/libTRDrec");
+  gSystem->Load("$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET)/libTRDsim");
+  //gSystem->Load("$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET)/libTRD");
+  //gSystem->Load("$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET)/libTRACKING");
+  
+  
+}
+