New comparison macro (M.Ivanov)
authorhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 20 Oct 2005 12:42:41 +0000 (12:42 +0000)
committerhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 20 Oct 2005 12:42:41 +0000 (12:42 +0000)
TRD/AliTRDComparison.C [new file with mode: 0644]

diff --git a/TRD/AliTRDComparison.C b/TRD/AliTRDComparison.C
new file mode 100644 (file)
index 0000000..41074b7
--- /dev/null
@@ -0,0 +1,716 @@
+/*
+
+//load mag field
+AliRunLoader *loader = AliRunLoader::Open("galice.root");
+loader->LoadgAlice();
+gAlice = loader->GetAliRun();
+AliMagF * field = gAlice->Field();
+AliKalmanTrack::SetFieldMap(gAlice->Field());
+AliExternalTrackParam::SetFieldMap(gAlice->Field());
+TFile fg("geometry.root");
+gGeoManager = (TGeoManager*)fg.Get("Geo");
+
+.L AliGenInfo.C+
+.L AliESDComparisonMI.C+
+.L AliTRDComparison.C+
+//.L AliESDComparisonPIC.C+
+.x TDR_style.C
+
+MakeCompTr();
+MakeTree();
+
+
+
+MakeComp();
+MakeAlias();
+TFile fff("TRDdebug.root");
+TTree * tree =(TTree*)fff.Get("tracklet");
+TTree * tree2 =(TTree*)fff.Get("Tracks");
+
+tree->Draw("mean0:tracklet.fZ-track.fZ","plane==0&&abs(angleb)<1","");
+
+
+
+*/
+
+//ROOT includes
+#include "Rtypes.h"
+#include "TFile.h"
+#include "TTree.h"
+#include "TChain.h"
+#include "TCut.h"
+#include "TString.h"
+#include "TBenchmark.h"
+#include "TStopwatch.h"
+#include "TParticle.h"
+#include "TSystem.h"
+#include "TTimer.h"
+#include "TVector3.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TCanvas.h"
+#include "TPad.h"
+#include "TF1.h"
+#include "TView.h"
+#include "TGeometry.h"
+#include "TPolyLine3D.h"
+#include "TLegend.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 "AliESDtrack.h"
+#include "AliTRDtrack.h"
+#include "AliITStrackMI.h"
+#include "AliESDV0MI.h"
+#include "AliESDkink.h"
+#include "AliTRDparameter.h"
+#include "AliTRDgeometryFull.h"
+#include "AliTRDcluster.h"
+
+// COMPARISON includes
+#include "../STEER/AliGenInfo.h"
+#include "../STEER/AliESDComparisonMI.h"
+
+
+
+AliComparisonDraw comp;
+AliComparisonDraw comptr;
+AliTRDparameter   par("pica","vyjebana");
+AliTRDgeometryFull    geom;
+AliRunLoader  *fLoader=0; //= AliRunLoader::Open(fnGalice);
+
+TCut cprim("cprim","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<0.01&&abs(MC.fVDist[2])<0.01");
+TCut citsin("cprim","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<0.01&&abs(MC.fVDist[2])<5");
+TCut clength("clength","abs(TRD0.fIntegratedLength-MC.fTOFReferences[0].fLength)<30");
+TCut c1("c1","TRD1.fStopped==0");    // track propagated to refernce
+TCut c0("c0","TRD0.fStopped==0");    // track propagated to reference
+TCut c2("c2","MC.fNTOFRef>0");       //track in tof
+TCut c3("c3","MC.fNTOFRef>0&&MC.fNTRDRef>5");       //track in tof and TRD
+TH1F hpull("hpull","hpull",50,-10,10);
+TH1F hpull6("hpull6","hpull6",50,-6,6);
+
+
+class AliTRDInfoMC: public TObject{
+ public:
+  void Update(AliMCInfo *info);  // update TRD info
+  void Reset();                  // clear TRD info
+  Short_t fLayer[6];         //indicates layer is on
+  Short_t fLayerA[6];        // indicates layer after is ON
+  Short_t fNLayers;          // number of crossed layers
+  Short_t fLast;             // last crossed layer
+  AliTrackReference fRef[6];    // trd reference
+  ClassDef(AliTRDInfoMC,1)   // TRD MC information 
+};
+
+ClassImp(AliTRDInfoMC)
+
+
+class AliTRDInfoCl: public TObject{
+  public:
+  AliTRDcluster fCl[6][20];      // clusters
+  Int_t  fNClusters;           // number of clusters
+  ClassDef(AliTRDInfoCl,1)   // TRD MC information 
+};
+
+ClassImp(AliTRDInfoCl)
+
+
+
+void AliTRDInfoMC::Reset(){
+  //
+  //
+  fNLayers=0;
+  fLast =0;
+  for (Int_t ilayer = 0; ilayer<6; ilayer++){
+    fLayer[ilayer]=0;
+    fLayerA[ilayer]=0;
+  }
+}
+void  AliTRDInfoMC::Update(AliMCInfo *info)
+{
+  // 
+  //
+  Reset();
+  //  Float_t radius[6]={305
+  for (Int_t i=0;i<info->fNTRDRef;i++){
+    AliTrackReference * ref = (AliTrackReference*)info->fTRDReferences->At(i);
+    Float_t r = ref->R();
+    for (Int_t jplane = 0; jplane<6; jplane++){
+      if (TMath::Abs(r - par.GetTime0(jplane))<6){
+      //      if (TMath::Abs(r - 305+jplane*14.)<6){
+       fLayer[jplane]++;
+       for (Int_t kplane=jplane;kplane>=0;kplane--){
+         fLayerA[kplane]++;
+         fRef[jplane]=*ref;
+       }
+      }
+    }
+  }
+  for (Int_t jplane = 0; jplane<6; jplane++){
+    if (fLayer[jplane]>0) fNLayers++;
+    if (fLayer[jplane]>0) fLast = jplane;
+  }
+}
+
+
+
+
+void  MaxBinNorm(const char * expr = "TRD0.fTimBinPlane"){
+  //
+  // max time bin probability distribution
+  //
+  char  var[100];
+  sprintf(var,"%s>>hispi",expr);
+  TH1F * hispi= new TH1F("hispi","hispi",20,1,22);
+  comp.fTree->Draw(var,"abs(MC.fPdg)==211&&TRD0.fN>90&&MC.fParticle.P()>3");
+  hispi->Scale(1./(hispi->GetEntries()));
+  hispi->Draw();
+  sprintf(var,"%s>>hisel",expr);
+  TH1F *hisel = new TH1F("hisel","hisel",20,1,22);
+  comp.fTree->Draw(var,"abs(MC.fPdg)==11&&TRD0.fN>90&&MC.fParticle.P()>3");
+  hisel->Scale(1./(hisel->GetEntries()));
+  hispi->SetMarkerStyle(25); hispi->SetLineStyle(2);
+  hisel->SetMarkerStyle(26); hisel->SetLineStyle(1);
+  hisel->SetXTitle("Maximal Time Bin");
+  hisel->SetYTitle("Probability");
+  hisel->Draw();
+  hispi->Draw("same");
+  TLegend *legend = new TLegend(0.55,0.12,0.85,0.35, "Maximal Time Bin position");
+  legend->SetBorderSize(1);
+  legend->AddEntry(hispi ,"Pions");
+  legend->AddEntry(hisel ,"Electrons");
+  legend->Draw();
+
+}
+
+
+void  dedxNorm(const char * expr = "TRD0.fdEdxPlane/20"){
+  //
+  // max time bin probability distribution
+  //
+  char  var[100];
+  sprintf(var,"%s>>hispi",expr);
+  TH1F * hispi= new TH1F("hispi","hispi",100,1,120);
+  comp.fTree->Draw(var,"abs(MC.fPdg)==211&&TRD0.fN>90&&MC.fParticle.P()>3");
+  hispi->Scale(1./(hispi->GetEntries()));
+  hispi->Draw();
+  sprintf(var,"%s>>hisel",expr);
+  TH1F *hisel = new TH1F("hisel","hisel",100,1,120);
+  comp.fTree->Draw(var,"abs(MC.fPdg)==11&&TRD0.fN>90&&MC.fParticle.P()>3");
+  hisel->Scale(1./(hisel->GetEntries()));
+  hispi->SetMarkerStyle(25); hispi->SetLineStyle(2);
+  hisel->SetMarkerStyle(26); hisel->SetLineStyle(1);
+  hispi->SetXTitle("Normalized dEdx");
+  hispi->SetYTitle("Probability");
+  hispi->Draw();
+  hisel->Draw("same");
+  TLegend *legend = new TLegend(0.55,0.12,0.85,0.35, "Normalized dEdx");
+  legend->SetBorderSize(1);
+  legend->AddEntry(hispi ,"Pions");
+  legend->AddEntry(hisel ,"Electrons");
+  legend->Draw();
+}
+
+
+void EfficiencyPt()
+{
+  //
+  //
+  //
+  comp.Eff("TR.Pt()","InfoMC.fNLayers>3&&abs(MC.fPdg)!=11"+c2+cprim,"TRD0.fN/InfoMC.fNLayers>10&&TRD0.fLab>0",10,0.3,2);
+  TH1F * his = (TH1F*)comp.fRes->Clone();
+  comp.Eff("TR.Pt()","InfoMC.fNLayers>3&&abs(MC.fPdg)!=11"+c2+cprim,"TRD0.fN/InfoMC.fNLayers>10&&TRD0.fLab<0",10,0.3,2);
+  TH1F * hisf = (TH1F*)comp.fRes->Clone();
+  his->SetMarkerStyle(25); his->SetLineStyle(2);
+  hisf->SetMarkerStyle(26); hisf->SetLineStyle(1);
+  his->SetXTitle("P_{t}(GeV/c)");
+  his->SetYTitle("Tracking efficiency (%)");
+  his->Draw();
+  hisf->Draw("same");
+  TLegend *legend = new TLegend(0.55,0.16,0.85,0.4, "TRD efficciency");
+  legend->SetBorderSize(1);
+  legend->AddEntry(his ,"Efficiency");
+  legend->AddEntry(hisf ,"Fake tracks");
+  legend->Draw();
+
+}
+
+void EfficiencyDip()
+{
+  //
+  //
+  //
+  comp.Eff("atan(TR.fPz/TR.Pt())","InfoMC.fNLayers>3&&abs(MC.fPdg)!=11&&TR.Pt()>0.4"+c2+cprim,"TRD0.fN/InfoMC.fNLayers>10&&TRD0.fLab>0",10,-0.9,0.9);
+  TH1F * his = (TH1F*)comp.fRes->Clone();
+  comp.Eff("atan(TR.fPz/TR.Pt())","InfoMC.fNLayers>3&&abs(MC.fPdg)!=11&&TR.Pt()>0.4"+c2+cprim,"TRD0.fN/InfoMC.fNLayers>10&&TRD0.fLab<0",10,-0.9,0.9);
+  TH1F * hisf = (TH1F*)comp.fRes->Clone();
+  his->SetMarkerStyle(25); his->SetLineStyle(2);
+  hisf->SetMarkerStyle(26); hisf->SetLineStyle(1);
+  his->SetXTitle("#lambda(rad)");
+  his->SetYTitle("Tracking efficiency (%)");
+  his->Draw();
+  hisf->Draw("same");
+  TLegend *legend = new TLegend(0.55,0.16,0.85,0.4, "TRD efficciency");
+  legend->SetBorderSize(1);
+  legend->AddEntry(his ,"Efficiency");
+  legend->AddEntry(hisf ,"Fake tracks");
+  legend->Draw();
+}
+
+void ClusterRatio()
+{
+  TH1F *hisr = new TH1F("hisr","Percentage of found clusters",80,0,120);
+  comp.fTree->Draw("100*TRD0.fN/(20*InfoMC.fNLayers)>>hisr","InfoMC.fNLayers>3&&TRD0.fLab>0&&abs(MC.fPdg)!=11&&TR.Pt()>0.4"+cprim+c2);
+  hisr->SetXTitle("(%)");
+  hisr->Draw();
+
+}
+
+
+
+void PtRes(){
+  
+  
+}
+
+
+TH1F * MakeCumul(TH1F *his)
+{
+  TH1F *hcumul = (TH1F*)his->Clone();
+  Float_t all = (Float_t)his->GetSum();
+  for (Int_t i=0;i<=his->GetNbinsX();i++){
+    hcumul->SetBinContent(i,his->Integral(0,i)/all);
+  } 
+  hcumul->SetFillColor(0);
+  return hcumul;
+}
+
+TH1F * PullCumul(const char *var, TCut cut, Int_t div=200, Float_t max=20){
+  TH1F * hpullb = new TH1F("hpullb","hpullb",div,0,max);
+  char v2[1000];
+  sprintf(v2,"%s>>hpullb",var);
+  comp.fTree->Draw(v2, cut);
+  TH1F * res = MakeCumul(hpullb);
+  delete hpullb;
+  return res;
+}
+
+
+
+void MakeComp(){
+  //
+  // Set tree and set aliases
+  //
+  TFile *f = new TFile("trdComp1.root");
+  TTree * tree= (TTree*)f->Get("Comparison");
+  comp.fTree =tree;  
+}
+
+void MakeAlias(){
+
+  comp.fTree->SetAlias("Pt","sqrt(MC.fParticle.fPx**2+MC.fParticle.fPy**2)");
+  comp.fTree->SetAlias("P","sqrt(MC.fParticle.fPx**2+MC.fParticle.fPy**2+MC.fParticle.fPz**2)");
+  comp.fTree->SetAlias("Ptpc","sqrt(fTrackRef.fPx**2+fTrackRef.fPy**2+fTrackRef.fPz**2)");
+  comp.fTree->SetAlias("VRadius","sqrt(MC.fParticle.fVx**2+MC.fParticle.fVy**2)");
+  comp.fTree->SetAlias("DRadius","sqrt(MC.fTRdecay.fX**2+MC.fTRdecay.fY**2)");
+  comp.fTree->SetAlias("ThetaM","MC.fParticle.fPz/Pt");
+  comp.fTree->SetAlias("rel0","(MC.fTOFReferences[0].P()-MC.fTPCReferences[0].P())/MC.fTOFReferences[0].P()");  //relative en. loss
+  comp.fTree->SetAlias("rel1","(TR.P()-MC.fTPCReferences[0].P())/MC.fTPCReferences[0].P()");  //relative en. loss - in the last TRD
+  comp.fTree->SetAlias("pt_tof","MC.fTOFReferences[0].Pt()");
+  comp.fTree->SetAlias("dpt_tof","(abs(1/RC.fRp[4])-MC.fTOFReferences[0].Pt())/MC.fTOFReferences[0].Pt()");
+
+
+  comp.fTree->SetAlias("d1pt_tof","(abs(RC.fRp[4])-1/MC.fTOFReferences[0].Pt())");
+  comp.fTree->SetAlias("d1ptp_tof","(abs(RC.fRp[4])-1/MC.fTOFReferences[0].Pt())/sqrt(RC.fRc[14])");
+
+  comp.fTree->SetAlias("dpt_trd","(abs(1/RC.fRp[4])-TR.Pt())/TR.Pt()");
+
+  comp.fTree->SetAlias("tofth","(MC.fTOFReferences[0].fPz/MC.fTOFReferences[0].Pt())");
+  comp.fTree->SetAlias("dth_tof","TRD0.fT[3]-tofth");
+  comp.fTree->SetAlias("phi_tof","atan2(MC.fTOFReferences[0].fPy,MC.fTOFReferences[0].fPx)");
+  comp.fTree->SetAlias("phi_trd","asin(TRD0.GetSnp())+TRD0.fAlpha");
+  comp.fTree->SetAlias("pullphi0","(TRD0.GetSnp()+TRD0.fAlpha-atan2(MC.fTOFReferences[0].fPy,MC.fTOFReferences[0].fPx))/sqrt(TRD0.fCee)");
+  comp.fTree->SetAlias("dz0","(TRD0.fZ-MC.fTOFReferences[0].fZ)");
+  comp.fTree->SetAlias("pully0","TRD0.fY/sqrt(TRD0.fCyy)");
+  comp.fTree->SetAlias("pullz0","(TRD0.fZ-MC.fTOFReferences[0].fZ)/sqrt(TRD0.fCzz)");
+  //
+  comp.fTree->SetAlias("sigmay","sqrt(TRD0.fCyy)*(1+1/(1+abs(4./(RC.fRp[4]))))");
+  comp.fTree->SetAlias("pully3","TRD0.fY/sigmay");  
+  comp.fTree->SetAlias("deltaz0","(TRD0.fZ-MC.fTOFReferences[0].fZ)");
+  //
+  comp.fTree->SetAlias("pullc1","(abs(RC.fRp[4])-1/TR.Pt())/sqrt(RC.fRc[14])");
+  comp.fTree->SetAlias("deltac1","(abs(RC.fRp[4])-1/TR.Pt())*TR.Pt()");
+  comp.fTree->SetAlias("dist","sqrt(deltaz0**2+TRD0.fY**2)");
+  comp.fTree->SetAlias("pull0","pully0**2+pullz0**2");
+  comp.fTree->SetAlias("pullc0","(abs(TRD0.fC)-(1/(TRD0.fLocalConvConst*MC.fTOFReferences[0].Pt())))/sqrt(TRD0.fCcc)");
+  //
+  comp.fTree->SetAlias("prob5","exp(-(TRD1.fTracklets[5].fChi2)/16.)*max((TRD1.fTracklets[5].fNFound-6)/9.,0)");
+  comp.fTree->SetAlias("prob4","exp(-(TRD1.fTracklets[4].fChi2)/16.)*max((TRD1.fTracklets[4].fNFound-6)/9.,0)");
+  //
+  comp.fTree->SetAlias("TPCE","sqrt(MC.fTPCReferences[4].P()**2+MC.fMass**2)");
+  comp.fTree->SetAlias("TOFE","sqrt(MC.fTOFReferences[0].P()**2+MC.fMass**2)");
+  comp.fTree->SetAlias("TRDE","sqrt(TR.P()**2+MC.fMass**2)");
+  comp.fTree->SetAlias("dtpi","TRD0.fIntegratedTime[2]-10^12*MC.fTOFReferences.fTime");
+  comp.fTree->SetAlias("dtk","TRD0.fIntegratedTime[3]-10^12*MC.fTOFReferences.fTime");
+  comp.fTree->SetAlias("dtp","TRD0.fIntegratedTime[4]-10^12*MC.fTOFReferences.fTime");
+  comp.fTree->SetAlias("dt","(abs(MC.fPdg)==2212)*dtp+(abs(MC.fPdg)==211)*dtpi+(abs(MC.fPdg)==321)*dtk"); //delta time
+  comp.fTree->SetAlias("Beta","sqrt(TR.P()**2/(TR.P()**2+MC.fMass**2))");
+  comp.fTree->SetAlias("dtsigma","10+(1-Beta)*280");
+  comp.fTree->SetAlias("dptrel","(abs(TRD0.GetPt())-MC.fTOFReferences[0].Pt())/MC.fTOFReferences[0].Pt()");
+  //
+  comp.fTree->SetAlias("dphi","atan2(TR.fY,TR..fX)-atan2(TR.fPy,TR..fPx)");
+  comp.fTree->SetAlias("dphi0","acos((TR.fX*TR.fPx+TR.fY*TR.fPy)/(TR.Pt()*TR.R()))");
+
+
+}
+
+
+void MakeCompTr(){
+  //
+  // Set tree and set aliases
+  //
+  TFile *ff = new TFile("TRDdebug.root");
+  TTree * treetr= (TTree*)ff->Get("tracklet");
+  comptr.fTree =treetr;  
+  
+}
+
+
+void DrawDE(TCut cut){
+  comp.DrawXY("(TRD0.fDE)","TPCE-TOFE",c2,cut,5,0.01,0.1,0,0.2);
+  comp.fRes->Draw();
+}
+
+
+void DrawYResol1(TCut cut){
+  comp.fTree->Draw("TRD1.fY","MC.fNTOFRef>0&&MC.fNTRDRef>5"+c1+cut,"");
+}
+
+void DrawYResol0(TCut cut){
+  comp.fTree->Draw("TRD0.fY","MC.fNTOFRef>0&&MC.fNTRDRef>5"+c0+cut,"");
+}
+
+
+void DrawTOFResY(){
+  //
+  comp.DrawXY("TR.Pt()","TRD0.fY",c2,"1",10,0.3,2.,-2,2,30);
+  TH1F * histof = (TH1F*)comp.fRes->Clone();
+  comp.DrawXY("TR.Pt()","TRD1.fY",c2,"1",10,0.3,2.,-2,2,30);
+  TH1F * histrd = (TH1F*)comp.fRes->Clone();
+  histof->SetMarkerStyle(24);
+  histrd->SetMarkerStyle(25);
+  histof->SetXTitle("P_{t}(GeV)");
+  histof->SetYTitle("#Delta_{r-#phi}(cm)");
+  histof->Draw();
+  histrd->Draw("same"); 
+  TLegend *legend = new TLegend(0.55,0.12,0.85,0.35, "TRD resolution r-#phi (cm)");
+  legend->SetBorderSize(1);
+  legend->AddEntry(histof ,"Resolution in TOF plane");
+  legend->AddEntry(histrd ,"Resolution in last TRD plane");
+  legend->Draw();
+}
+
+void DrawTOFResZ(){
+  //
+  comp.DrawXY("TR.Pt()","TRD0.fZ-MC.fTOFReferences.fZ",c2,"1",10,0.3,2.,-4,4,30);
+  TH1F * histof = (TH1F*)comp.fRes->Clone();
+  comp.DrawXY("TR.Pt()","TRD1.fZ-TR.fZ",c2,"1",10,0.3,2.,-4,4,30);
+  TH1F * histrd = (TH1F*)comp.fRes->Clone();
+  histof->SetMarkerStyle(24);
+  histrd->SetMarkerStyle(25);
+  histof->SetXTitle("P_{t}(GeV)");
+  histof->SetYTitle("#Delta_{z}(cm)");
+  histof->Draw();
+  histrd->Draw("same"); 
+  TLegend *legend = new TLegend(0.55,0.12,0.85,0.35, "TRD resolution z (cm)");
+  legend->SetBorderSize(1);
+  legend->AddEntry(histof ,"Resolution in TOF plane");
+  legend->AddEntry(histrd ,"Resolution in last TRD plane");
+  legend->Draw();
+}
+
+
+void MakeCumul(TCut & cut)
+{
+  TH1F * hisy = PullCumul("abs(pully0)",cut+cprim+c0+c2+c3,1000,20);
+  TH1F * hisz = PullCumul("abs(pullz0)",cut+cprim+c0+c2+c3,1000,20);
+  TH1F * hisyz = PullCumul("sqrt(pullz0**2+pully0**2)",cut+cprim+c0+c2+c3,1000,20);
+  hisz->SetXTitle("Pull (unit)");
+  hisz->SetYTitle("Cumulative function");
+  hisy->SetLineColor(1);hisz->SetLineColor(2);hisyz->SetLineColor(3);
+  hisz->Draw(); hisy->Draw("same");hisyz->Draw("same");  
+  TLegend *legend = new TLegend(0.55,0.12,0.85,0.35, "TRD pull (unit)");
+  legend->SetBorderSize(1);
+  legend->AddEntry(hisy ,"Pull y");
+  legend->AddEntry(hisz,"Pull z");
+  legend->AddEntry(hisyz,"Pull yz");
+  legend->Draw();
+}
+
+void MakeCumulAbs(TCut & cut)
+{
+  TH1F * hisy = PullCumul("abs(TRD0.fY)",cut+cprim+c0+c2+c3,1000,6);
+  TH1F * hisz = PullCumul("abs(deltaz0)",cut+cprim+c0+c2+c3,1000,6);  
+  TH1F * hisyz = PullCumul("sqrt(TRD0.fY**2+deltaz0**2)",cut+cprim+c0+c2+c3,1000,6);
+  hisz->SetXTitle("Delta (cm)");
+  hisz->SetYTitle("Cumulative function");
+  //
+  hisy->SetLineColor(1);hisz->SetLineColor(2);hisyz->SetLineColor(3);
+  hisz->Draw(); hisy->Draw("same");hisyz->Draw("same");  
+  TLegend *legend = new TLegend(0.55,0.12,0.85,0.35, "TRD residulas (cm)");
+  legend->SetBorderSize(1);
+  legend->AddEntry(hisy ,"Delta y (cm)");
+  legend->AddEntry(hisz,"Delta z (cm)");
+  legend->AddEntry(hisyz,"Delta  (cm)");
+  legend->Draw();
+}
+
+void LoadClusters(TObjArray *arr, Int_t event)
+{
+  //
+  //
+  // event=0
+  printf("Load cluster - event%d\n",event);
+  fLoader->SetEventNumber(event);
+  AliLoader * trdl = (AliLoader*)fLoader->GetLoader("TRDLoader");
+  trdl->LoadRecPoints();
+
+  TTree * ClusterTree = trdl->TreeR();
+  if (!ClusterTree) return;
+  TBranch * branch = ClusterTree->GetBranch("TRDcluster"); 
+  TObjArray *clusterArray = new TObjArray(1000); 
+  branch->SetAddress(&clusterArray);
+  Int_t nEntries = (Int_t) ClusterTree->GetEntries();
+  //
+  for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {    
+    ClusterTree->GetEvent(iEntry);  
+    for (Int_t icl =0; icl<clusterArray->GetEntriesFast();icl++){
+      AliTRDcluster *cl = (AliTRDcluster*)clusterArray->UncheckedAt(icl);
+      if (!cl) continue;
+      Int_t label = cl->GetLabel(0);
+      Int_t detector=cl->GetDetector();
+      Int_t localTimeBin=cl->GetLocalTimeBin();
+      //      Int_t sector=geom.GetSector(detector);
+      Int_t plane=geom.GetPlane(detector);
+      Int_t pos = label;      
+      if (label<0) continue;
+      if (localTimeBin>19||localTimeBin<0) continue;
+      if (plane<0||plane>5) continue;
+      
+      AliTRDInfoCl * info =( AliTRDInfoCl*)arr->At(pos); 
+      if (!info) {
+       info = new AliTRDInfoCl;
+       arr->AddAt(info,pos);
+      } 
+     
+      info->fCl[plane][localTimeBin]=*cl;
+    }
+  }  
+}
+
+
+void MakeTree()
+{
+  fLoader = AliRunLoader::Open("galice.root");
+  AliESDtrack * esd     = new AliESDtrack;
+  AliTRDtrack * trd     = new AliTRDtrack;
+  AliTRDtrack * trd0    = new AliTRDtrack;
+  AliTRDtrack * trd1    = new AliTRDtrack;
+  AliTRDtrack * trdb    = new AliTRDtrack;
+  AliMCInfo* info       = new AliMCInfo;
+  AliTRDInfoMC *trdinfo = new AliTRDInfoMC;
+  AliTRDInfoCl *infocl = new AliTRDInfoCl;
+  AliTrackReference *ref = new AliTrackReference;
+  TObjArray         * clarray = new TObjArray(6*1000000);
+  Int_t       eventNr   =0;
+  const Int_t kmaxlabel=1000000;
+  //
+  // MC part
+  //
+  TFile f("cmpESDTracks.root");
+  TTree * tree1 = (TTree*) f.Get("ESDcmpTracks");
+  TBranch * brmc  = tree1->GetBranch("MC");
+  brmc->SetAddress(&info);
+  //
+  // 
+  TObjArray *esds = new TObjArray(kmaxlabel);
+  TObjArray *trds = new TObjArray(kmaxlabel);
+  //
+  //  rec part
+  //
+  TFile f2("TRDdebug.root");
+  TTree * tree2 =(TTree*)f2.Get("Tracks");
+  TBranch * br20 = tree2->GetBranch("ESD.");
+  TBranch * br21 = tree2->GetBranch("trd.");
+  TBranch * br22 = tree2->GetBranch("EventNr");
+  br20->SetAddress(&esd);
+  br21->SetAddress(&trd);  
+  br22->SetAddress(&eventNr);
+  //
+  // comparison part
+  //
+  TFile f3("trdComp1.root","new");
+  TTree   * tree3   = new TTree("Comparison","Comparison");
+  tree3->Branch("MC.","AliMCInfo",&info);
+  tree3->Branch("RC.","AliESDtrack",&esd);
+  tree3->Branch("TRD.","AliTRDtrack",&trd);
+  tree3->Branch("TRD0.","AliTRDtrack",&trd0);  // track at tof ref
+  tree3->Branch("TRD1.","AliTRDtrack",&trd1);  // track at last TRD ref
+  tree3->Branch("TRDb.","AliTRDtrack",&trdb);  // track at tof ref
+  tree3->Branch("TR.","AliTrackReference",&ref);  // track at last TRD ref
+  tree3->Branch("InfoMC.","AliTRDInfoMC",&trdinfo);  // info about TRD track MC
+  tree3->Branch("InfoCl.","AliTRDInfoCl",&infocl);  // info about TRD track MC
+  //
+  //
+  Int_t lastevent=-1;
+
+  for (Int_t i=0;i<tree1->GetEntries();i++){
+    brmc->GetEntry(i);
+    if (lastevent!=info->fEventNr){
+      printf("Read Event\t%d\n", info->fEventNr);
+      //
+      // read new event
+      //
+      lastevent = info->fEventNr;
+      //
+      // delete all event
+      for (Int_t ip=0;ip<esds->GetEntriesFast();ip++){
+       if (esds->At(ip)){
+         delete esds->At(ip);
+         delete trds->At(ip);
+       }
+      }
+      esds->Clear();
+      trds->Clear(); 
+      clarray->Clear();
+      LoadClusters(clarray,lastevent);
+      //
+      AliESDtrack * esdp     = new AliESDtrack;
+      AliTRDtrack * trdp     = new AliTRDtrack;
+
+      br20->SetAddress(&esdp);
+      br21->SetAddress(&trdp);
+      //
+      // Load esd and trd tracks to memory
+      Int_t ntracks0 =0;
+      Int_t ntracks1 =0;
+      Int_t ntracks2 =0;
+      
+      for (Int_t ip=0;ip<tree2->GetEntries();ip++){
+       tree2->GetEntry(ip);
+       if (eventNr!=info->fEventNr) continue;
+       //      if (eventNr>info->fEventNr) continue;
+       //
+       Int_t label = TMath::Abs(esdp->GetTPCLabel());
+       if (label==0) continue;
+       if (label>=kmaxlabel) continue;
+       ntracks0++;
+       if (!(esds->At(label))){
+         esds->AddAt(new AliESDtrack(*esdp),label);
+         trds->AddAt(new AliTRDtrack(*trdp),label);
+         ntracks1++;
+       } else{
+         if (esd->GetStatus() & AliESDtrack::kITSin){
+           esds->AddAt(new AliESDtrack(*esdp),label);
+           trds->AddAt(new AliTRDtrack(*trdp),label);
+           ntracks2++;
+         }
+       }
+      }
+      printf(" readed\t%d\t%d\t%d\n",ntracks0, ntracks1,ntracks2);
+    }
+
+    //    if (info->fNTOFRef==0 && info->fNTRDRef<5) continue;
+    Int_t label = info->fLabel;    
+    trdinfo->Update(info);
+    if (clarray->At(label)){
+      infocl = (AliTRDInfoCl*)clarray->At(label);
+    }
+    if (info->fNTPCRef<3) continue;
+    trd0->SetStop(kTRUE);
+    trdb->SetStop(kTRUE);
+    trd1->SetStop(kTRUE);
+    //
+    if (esds->At(label)) {
+      //      if (i%1000==0) printf("\n%d\n",i);
+      trd   = (AliTRDtrack*)(trds->At(label));
+      esd   = (AliESDtrack*)(esds->At(label));
+      //
+      *trd0 = *trd;
+      if (trd->GetBackupTrack()) 
+       *trdb = *trd->GetBackupTrack();
+      else 
+        *trdb = *trd;
+      *trd1 = *trd;
+      if (info->fNTOFRef>0){
+       AliTrackReference * ref = ((AliTrackReference*)(info->fTOFReferences->At(0)));
+       Double_t x   = ((AliTrackReference*)(info->fTOFReferences->At(0)))->X();
+       Double_t y   = ((AliTrackReference*)(info->fTOFReferences->At(0)))->Y();
+       Double_t phi    = TMath::ATan2(y,x);    
+       Double_t radius = TMath::Sqrt(x*x+y*y);
+       //
+       // rotate and propagate to local system
+       //
+       trd0->SetStop(kFALSE);
+       if (!trd0->Rotate(phi-trd0->GetAlpha())) trd0->SetStop(kTRUE);
+       else{
+         Double_t xyz0[3], param[7];
+         trd0->GetGlobalXYZ(xyz0[0],xyz0[1],xyz0[2]);
+         Double_t xyz1[3]={ref->X(),ref->Y(),ref->Z()};
+         AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
+         if (!trd0->PropagateTo(radius,param[1],param[0])) trd0->SetStop(kTRUE);
+       }
+       trdb->SetStop(kFALSE);
+       if (!trdb->Rotate(phi-trdb->GetAlpha())) trdb->SetStop(kTRUE);
+       else
+         if (!trdb->PropagateTo(radius)) trdb->SetStop(kTRUE);
+
+
+      }    
+
+      if (info->fNTRDRef>0){
+       Double_t cradius=0;
+       Double_t cphi=0;
+       for (Int_t itrd =0; itrd<info->fNTRDRef;itrd++){
+         Double_t x   = ((AliTrackReference*)(info->fTRDReferences->At(itrd)))->X();
+         Double_t y   = ((AliTrackReference*)(info->fTRDReferences->At(itrd)))->Y();
+         Double_t phi    = TMath::ATan2(y,x);  
+         //      Double_t snp    =  TMath::ATan2(((AliTrackReference*)(info->fTRDReferences->At(itrd)))->Py(),
+         //                              ((AliTrackReference*)(info->fTRDReferences->At(itrd)))->Px());
+         //      if (TMath::Abs(snp)>0.9 &&cradius>10) break; 
+         Double_t radius = TMath::Sqrt(x*x+y*y);
+         if (radius>cradius){
+           cradius = radius;
+           cphi    = phi;
+           ref     =  ((AliTrackReference*)(info->fTRDReferences->At(itrd)));
+         }else break;   
+       }
+       //
+       // rotate and propagate to local system of the last track reference
+       //
+       trd1->SetStop(kFALSE);
+       if (!trd1->Rotate(cphi-trd1->GetAlpha())) trd1->SetStop(kTRUE);
+       else{
+         if (!trd1->PropagateTo(cradius)) trd1->SetStop(kTRUE);
+       }
+      }    
+      tree3->Fill();
+    }
+  }
+  f3.cd();
+  tree3->Write();
+}
+
+