[u/mrichter/AliRoot.git] / TOF / AliTOFComparison.C

/****************************************************************************
 *      This macro estimates efficiency of matching with the TOF.           *
 *      TOF "Good" tracks are those originating from the primary vertex,    *
 *      being "good" in the ITS and having at least one digit in the TOF.   * 
 *         (To get the list of "good" tracks one should first run           *
 *          AliTPCComparison.C and AliITSComparisonV2.C macros)             *
 *           Origin: I.Belikov, CERN, Jouri.Belikov@cern.ch                 *
 ****************************************************************************/

#if !defined(__CINT__) || defined(__MAKECINT__)
  #include <Riostream.h>
  #include <fstream.h>

  #include "AliRun.h"
  #include "AliHeader.h"
  #include "AliStack.h"
  #include "AliRunLoader.h"
  #include "AliLoader.h"

  #include "AliESD.h"
  #include "AliTOFdigit.h"

  #include "TKey.h"
  #include "TFile.h"
  #include "TTree.h"
  #include "TH1.h"
  #include "TClonesArray.h"
  #include "TStyle.h"
  #include "TCanvas.h"
  #include "TLine.h"
  #include "TText.h"
  #include "TParticle.h"
#endif

struct GoodTrackTOF {
  Int_t lab;
  Int_t code;
  Float_t px,py,pz;
  Float_t x,y,z;
};

extern AliRun *gAlice;

Int_t AliTOFComparison() {
   Int_t good_tracks_tof(GoodTrackTOF *gt, const Int_t max);

   cerr<<"Doing comparison...\n";

   if (gAlice) { 
     delete gAlice->GetRunLoader();
     delete gAlice;//if everything was OK here it is already NULL
     gAlice = 0x0;
   }

   TClonesArray dummy("AliTOFdigit",10000), *digits=&dummy;   

   {
   AliRunLoader *rl = AliRunLoader::Open("galice.root");
   if (rl == 0x0) {
      cerr<<"Can not open session"<<endl;
      return 1;
   }
   AliLoader* tofl = rl->GetLoader("TOFLoader");
   if (tofl == 0x0) {
      cerr<<"Can not get the TOF loader"<<endl;
      return 2;
   }
   tofl->LoadDigits("read");

   rl->GetEvent(0);

   TTree *tofTree=tofl->TreeD();
   if (!tofTree) {
      cerr<<"Can't get the TOF cluster tree !\n";
      return 3;
   } 
   TBranch *branch=tofTree->GetBranch("TOF");
   if (!branch) { 
      cerr<<"Can't get the branch with the TOF digits !\n";
      return 4;
   }
   branch->SetAddress(&digits);

   tofTree->GetEvent(0);

   delete rl;
   }

   Int_t nd=digits->GetEntriesFast();
   cerr<<"Number of digits: "<<nd<<endl;

   const Int_t MAX=15000;
   GoodTrackTOF gt[MAX];
   Int_t ngood=0;
   ifstream in("good_tracks_tof");
   if (in) {
      cerr<<"Reading good tracks...\n";
      while (in>>gt[ngood].lab>>gt[ngood].code>>
                 gt[ngood].px>>gt[ngood].py>>gt[ngood].pz>>
                 gt[ngood].x >>gt[ngood].y >>gt[ngood].z) {
         ngood++;
         cerr<<ngood<<'\r';
         if (ngood==MAX) {
            cerr<<"Too many good tracks !\n";
            break;
         }
      }
      if (!in.eof()) cerr<<"Read error (good_tracks_its) !\n";
   } else {
      cerr<<"Marking good tracks (this will take a while)...\n";
      ngood=good_tracks_tof(gt,MAX);
      ofstream out("good_tracks_tof");
      if (out) {
        for (Int_t ngd=0; ngd<ngood; ngd++)
          out<<gt[ngd].lab<<' '<<gt[ngd].code<<' '
             <<gt[ngd].px<<' '<<gt[ngd].py<<' '<<gt[ngd].pz<<' '
             <<gt[ngd].x <<' '<<gt[ngd].y <<' '<<gt[ngd].z <<endl;
      } else cerr<<"Can not open file (good_tracks_tof) !\n";
      out.close();
   }

   Double_t pmin=0.2;
   Double_t pmax=4.0;

   TH1F *hgood=new TH1F("hgood","Good tracks",30,pmin,pmax);    
   TH1F *hfound=new TH1F("hfound","Matched tracks",30,pmin,pmax);
   TH1F *hfake=new TH1F("hfake","Mismatched tracks",30,pmin,pmax);
   TH1F *hgp=new TH1F("hgp","",30,pmin,pmax); //efficiency for good tracks
   hgp->SetLineColor(4); hgp->SetLineWidth(2);
   TH1F *hfp=new TH1F("hfp","Probability of mismatching",30,pmin,pmax);
   hfp->SetFillColor(1); hfp->SetFillStyle(3013); hfp->SetLineWidth(2);

   TH1F *hgoo=new TH1F("hgoo","Good tracks",30,-1,1);    
   TH1F *hfoun=new TH1F("hfoun","Matched tracks",30,-1,1);
   TH1F *hfak=new TH1F("hfak","Mismatched tracks",30,-1,1);
   TH1F *hgl=new TH1F("hgl","",30,-1,1); //efficiency for good tracks
   hgl->SetLineColor(4); hgl->SetLineWidth(2);
   TH1F *hfl=new TH1F("hfl","Probability of mismatching",30,-1,1);
   hfl->SetFillColor(1); hfl->SetFillStyle(3013); hfl->SetLineWidth(2);

   TCanvas *c1=new TCanvas("c1","",0,0,600,900);
   c1->Divide(1,2);

   TFile *ef=TFile::Open("AliESDs.root");
   if (!ef->IsOpen()) {cerr<<"Can't AliESDs.root !\n"; return 1;}

   TIter next(ef->GetListOfKeys());
   TKey *key=0;
   Int_t nev=0;
   while ((key=(TKey*)next())!=0) {
     cerr<<"Processing event number : "<<nev++<<endl;

     AliESD *event=(AliESD*)key->ReadObj();
     Int_t ntrk=event->GetNumberOfTracks();
     cerr<<"Number of ESD tracks : "<<ntrk<<endl; 

     Int_t matched=0;
     Int_t mismatched=0;
     for (Int_t i=0; i<ngood; i++) {
         Int_t lab=gt[i].lab;
         Double_t pxg=gt[i].px, pyg=gt[i].py, pzg=gt[i].pz;
         Double_t ptg=TMath::Sqrt(pxg*pxg+pyg*pyg);

	 if (ptg<0.1) continue;

         Double_t tgl=pzg/ptg; //tan(lambda)

         if (ptg>pmin) { hgood->Fill(ptg); hgoo->Fill(tgl); }

         Int_t j;
	 AliESDtrack *t=0;
         for (j=0; j<ntrk; j++) {
             AliESDtrack *tt=event->GetTrack(j);
             if (lab!=TMath::Abs(tt->GetLabel())) continue;
             t=tt;
             //if ((tt->GetStatus()&AliESDtrack::kTOFpid) == 0) continue;
             if (tt->GetTOFsignal() < 0) continue;
             UInt_t idx=tt->GetTOFcluster();
             if ((Int_t)idx>=nd) {
	       cerr<<"Wrong digit index ! "<<idx<<endl;
               return 5;
             }
             AliTOFdigit *dig=(AliTOFdigit*)digits->UncheckedAt(idx);
             Int_t *label=dig->GetTracks();
             if (label[0]!=lab)
             if (label[1]!=lab)
	       if (label[2]!=lab) {
                 mismatched++; 
                 if (ptg>pmin) { hfake->Fill(ptg); hfak->Fill(tgl); } 
                 break;
               }
             if (ptg>pmin) { hfound->Fill(ptg); hfoun->Fill(tgl); }
             matched++;
             break;
         }
         if (j==ntrk) {
	    cerr<<"Not matched: "<<lab<<"   ";
            if (t) {
               cerr<<(t->GetStatus()&AliESDtrack::kITSout)<<' '
	           <<(t->GetStatus()&AliESDtrack::kTPCout)<<' '
	           <<(t->GetStatus()&AliESDtrack::kTRDout)<<' '
	           <<(t->GetStatus()&AliESDtrack::kTIME);
	    } else cerr<<"No ESD track !";
            cerr<<endl;
         }
     }

     cerr<<"Number of good tracks: "<<ngood<<endl;
     cerr<<"Number of matched tracks: "<<matched<<endl;
     cerr<<"Number of mismatched tracks: "<<mismatched<<endl;
     if (ngood!=0) cerr<<"Efficiency: "<<Float_t(matched)/ngood<<endl;

     hfound->Sumw2(); hgood->Sumw2(); hfake->Sumw2();
     hgp->Divide(hfound,hgood,1,1.,"b");
     hfp->Divide(hfake,hgood,1,1.,"b");
     hgp->SetMaximum(1.4);
     hgp->SetYTitle("Matching efficiency");
     hgp->SetXTitle("Pt (GeV/c)");

     hfoun->Sumw2(); hgoo->Sumw2(); hfak->Sumw2();
     hgl->Divide(hfoun,hgoo,1,1.,"b");
     hfl->Divide(hfak,hgoo,1,1.,"b");
     hgl->SetMaximum(1.4);
     hgl->SetYTitle("Matching efficiency");
     hgl->SetXTitle("Tan(lambda)");

     c1->cd(1);

     hgp->Draw();
     hfp->Draw("histsame");
     TLine *line1 = new TLine(pmin,1.0,pmax,1.0); line1->SetLineStyle(4);
     line1->Draw("same");
     TLine *line2 = new TLine(pmin,0.9,pmax,0.9); line2->SetLineStyle(4);
     line2->Draw("same");

     c1->cd(2);

     hgl->Draw();
     hfl->Draw("histsame");
     TLine *line3 = new TLine(-1,1.0,1,1.0); line3->SetLineStyle(4);
     line3->Draw("same");
     TLine *line4 = new TLine(-1,0.9,1,0.9); line4->SetLineStyle(4);
     line4->Draw("same");

     break;
   }

   TFile fc("AliTOFComparison.root","RECREATE");
   c1->Write();
   fc.Close();

   return 0;
}

Int_t good_tracks_tof(GoodTrackTOF *gt, const Int_t max) {
   ifstream in("good_tracks_its");
   if (!in) {
     cerr<<"Can't get good_tracks_its !\n"; exit(11);
   }   

   AliRunLoader *rl = AliRunLoader::Open("galice.root","COMPARISON");
   if (!rl) {
       cerr<<"Can't start session !\n";
       exit(1);
   }

   rl->GetEvent(0);


   rl->LoadgAlice();
   rl->LoadHeader();
   Int_t np = rl->GetHeader()->GetNtrack();

   Int_t *good=new Int_t[np];
   Int_t k;
   for (k=0; k<np; k++) good[k]=0;

   AliLoader* tofl = rl->GetLoader("TOFLoader");
   if (tofl == 0x0) {
      cerr<<"Can not get the TOF loader"<<endl;
      exit(2);
   }
   tofl->LoadDigits("read");

   TTree *dTree=tofl->TreeD();
   if (!dTree) {
      cerr<<"Can't get the TOF cluster tree !\n";
      exit(3);
   } 

   TBranch *branch=dTree->GetBranch("TOF");
   if (!branch) { 
     cerr<<"Can't get the branch with the TOF digits !\n";
     exit(4);
   }
   TClonesArray dummy("AliTOFdigit",10000), *digits=&dummy;
   branch->SetAddress(&digits);
   
  dTree->GetEvent(0);
  Int_t nd=digits->GetEntriesFast();
  cerr<<"Number of digits: "<<nd<<endl;

  for (Int_t i=0; i<nd; i++) {
    AliTOFdigit *d=(AliTOFdigit*)digits->UncheckedAt(i);
    Int_t l0=d->GetTrack(0);
       if (l0>=np) {cerr<<"Wrong label: "<<l0<<endl; continue;}
    Int_t l1=d->GetTrack(1);
       if (l1>=np) {cerr<<"Wrong label: "<<l1<<endl; continue;}
    Int_t l2=d->GetTrack(2);
       if (l2>=np) {cerr<<"Wrong label: "<<l2<<endl; continue;}
    if (l0>=0) good[l0]++; 
    if (l1>=0) good[l1]++; 
    if (l2>=0) good[l2]++;
  }

   
  rl->LoadKinematics();
  AliStack* stack = rl->Stack();
  Int_t nt=0;
  Double_t px,py,pz,x,y,z;
  Int_t code,lab;
  while (in>>lab>>code>>px>>py>>pz>>x>>y>>z) {
      if (good[lab] == 0) continue;
      TParticle *p = (TParticle*)stack->Particle(lab);
      if (p == 0x0) {
         cerr<<"Can not get particle "<<lab<<endl;
         exit(1);
      }
      if (TMath::Abs(p->Vx())>0.1) continue;
      if (TMath::Abs(p->Vy())>0.1) continue;
      //      if (TMath::Abs(p->Vz())>0.1) continue;

      gt[nt].lab=lab;
      gt[nt].code=p->GetPdgCode();
//**** px py pz - in global coordinate system
      gt[nt].px=p->Px(); gt[nt].py=p->Py(); gt[nt].pz=p->Pz();
      gt[nt].x=p->Vx(); gt[nt].y=p->Vy(); gt[nt].z=p->Vz();
      nt++;
      if (nt==max) {cerr<<"Too many good tracks !\n"; break;}
  }

  delete[] good;

  rl->UnloadKinematics();
  rl->UnloadHeader();
  rl->UnloadgAlice();
  delete rl;

  return nt;
}
Commit	Line	Data
3f83f224	1	/****************************************************************************
	2	* This macro estimates efficiency of matching with the TOF. *
	3	* TOF "Good" tracks are those originating from the primary vertex, *
	4	* being "good" in the ITS and having at least one digit in the TOF. *
	5	* (To get the list of "good" tracks one should first run *
	6	* AliTPCComparison.C and AliITSComparisonV2.C macros) *
	7	* Origin: I.Belikov, CERN, Jouri.Belikov@cern.ch *
	8	****************************************************************************/
	9
	10	#if !defined(__CINT__) \|\| defined(__MAKECINT__)
	11	#include <Riostream.h>
	12	#include <fstream.h>
	13
	14	#include "AliRun.h"
	15	#include "AliHeader.h"
	16	#include "AliStack.h"
	17	#include "AliRunLoader.h"
	18	#include "AliLoader.h"
	19
	20	#include "AliESD.h"
	21	#include "AliTOFdigit.h"
	22
	23	#include "TKey.h"
	24	#include "TFile.h"
	25	#include "TTree.h"
	26	#include "TH1.h"
	27	#include "TClonesArray.h"
	28	#include "TStyle.h"
	29	#include "TCanvas.h"
	30	#include "TLine.h"
	31	#include "TText.h"
	32	#include "TParticle.h"
	33	#endif
	34
	35	struct GoodTrackTOF {
	36	Int_t lab;
	37	Int_t code;
	38	Float_t px,py,pz;
	39	Float_t x,y,z;
	40	};
	41
	42	extern AliRun *gAlice;
	43
	44	Int_t AliTOFComparison() {
	45	Int_t good_tracks_tof(GoodTrackTOF *gt, const Int_t max);
	46
	47	cerr<<"Doing comparison...\n";
	48
	49	if (gAlice) {
	50	delete gAlice->GetRunLoader();
	51	delete gAlice;//if everything was OK here it is already NULL
	52	gAlice = 0x0;
	53	}
	54
	55	TClonesArray dummy("AliTOFdigit",10000), *digits=&dummy;
	56
	57	{
	58	AliRunLoader *rl = AliRunLoader::Open("galice.root");
	59	if (rl == 0x0) {
	60	cerr<<"Can not open session"<<endl;
	61	return 1;
	62	}
	63	AliLoader* tofl = rl->GetLoader("TOFLoader");
	64	if (tofl == 0x0) {
65	cerr<<"Can not get the TOF loader"<<endl;
66	return 2;
67	}
68	tofl->LoadDigits("read");
69
70	rl->GetEvent(0);
71
72	TTree *tofTree=tofl->TreeD();
73	if (!tofTree) {
74	cerr<<"Can't get the TOF cluster tree !\n";
75	return 3;
76	}
77	TBranch *branch=tofTree->GetBranch("TOF");
78	if (!branch) {
79	cerr<<"Can't get the branch with the TOF digits !\n";
80	return 4;
81	}
82	branch->SetAddress(&digits);
83
84	tofTree->GetEvent(0);
85
86	delete rl;
87	}
88
89	Int_t nd=digits->GetEntriesFast();
90	cerr<<"Number of digits: "<<nd<<endl;
91
92	const Int_t MAX=15000;
93	GoodTrackTOF gt[MAX];
94	Int_t ngood=0;
95	ifstream in("good_tracks_tof");
96	if (in) {
97	cerr<<"Reading good tracks...\n";
98	while (in>>gt[ngood].lab>>gt[ngood].code>>
99	gt[ngood].px>>gt[ngood].py>>gt[ngood].pz>>
100	gt[ngood].x >>gt[ngood].y >>gt[ngood].z) {
101	ngood++;
102	cerr<<ngood<<'\r';
103	if (ngood==MAX) {
104	cerr<<"Too many good tracks !\n";
105	break;
106	}
107	}
108	if (!in.eof()) cerr<<"Read error (good_tracks_its) !\n";
109	} else {
110	cerr<<"Marking good tracks (this will take a while)...\n";
111	ngood=good_tracks_tof(gt,MAX);
112	ofstream out("good_tracks_tof");
113	if (out) {
114	for (Int_t ngd=0; ngd<ngood; ngd++)
115	out<<gt[ngd].lab<<' '<<gt[ngd].code<<' '
116	<<gt[ngd].px<<' '<<gt[ngd].py<<' '<<gt[ngd].pz<<' '
117	<<gt[ngd].x <<' '<<gt[ngd].y <<' '<<gt[ngd].z <<endl;
118	} else cerr<<"Can not open file (good_tracks_tof) !\n";
119	out.close();
120	}
121
122	Double_t pmin=0.2;
123	Double_t pmax=4.0;
124
125	TH1F *hgood=new TH1F("hgood","Good tracks",30,pmin,pmax);
126	TH1F *hfound=new TH1F("hfound","Matched tracks",30,pmin,pmax);
127	TH1F *hfake=new TH1F("hfake","Mismatched tracks",30,pmin,pmax);
128	TH1F *hgp=new TH1F("hgp","",30,pmin,pmax); //efficiency for good tracks
129	hgp->SetLineColor(4); hgp->SetLineWidth(2);
130	TH1F *hfp=new TH1F("hfp","Probability of mismatching",30,pmin,pmax);
131	hfp->SetFillColor(1); hfp->SetFillStyle(3013); hfp->SetLineWidth(2);
132
133	TH1F *hgoo=new TH1F("hgoo","Good tracks",30,-1,1);
134	TH1F *hfoun=new TH1F("hfoun","Matched tracks",30,-1,1);
135	TH1F *hfak=new TH1F("hfak","Mismatched tracks",30,-1,1);
136	TH1F *hgl=new TH1F("hgl","",30,-1,1); //efficiency for good tracks
137	hgl->SetLineColor(4); hgl->SetLineWidth(2);
138	TH1F *hfl=new TH1F("hfl","Probability of mismatching",30,-1,1);
139	hfl->SetFillColor(1); hfl->SetFillStyle(3013); hfl->SetLineWidth(2);
140
1c3e95d8	141	TCanvas *c1=new TCanvas("c1","",0,0,600,900);
	142	c1->Divide(1,2);
	143
3f83f224	144	TFile *ef=TFile::Open("AliESDs.root");
	145	if (!ef->IsOpen()) {cerr<<"Can't AliESDs.root !\n"; return 1;}
	146
	147	TIter next(ef->GetListOfKeys());
	148	TKey *key=0;
	149	Int_t nev=0;
	150	while ((key=(TKey*)next())!=0) {
	151	cerr<<"Processing event number : "<<nev++<<endl;
	152
	153	AliESD event=(AliESD)key->ReadObj();
	154	Int_t ntrk=event->GetNumberOfTracks();
	155	cerr<<"Number of ESD tracks : "<<ntrk<<endl;
	156
	157	Int_t matched=0;
	158	Int_t mismatched=0;
	159	for (Int_t i=0; i<ngood; i++) {
	160	Int_t lab=gt[i].lab;
	161	Double_t pxg=gt[i].px, pyg=gt[i].py, pzg=gt[i].pz;
	162	Double_t ptg=TMath::Sqrt(pxgpxg+pygpyg);
	163
	164	if (ptg<0.1) continue;
	165
	166	Double_t tgl=pzg/ptg; //tan(lambda)
	167
	168	if (ptg>pmin) { hgood->Fill(ptg); hgoo->Fill(tgl); }
	169
	170	Int_t j;
	171	AliESDtrack *t=0;
	172	for (j=0; j<ntrk; j++) {
	173	AliESDtrack *tt=event->GetTrack(j);
	174	if (lab!=TMath::Abs(tt->GetLabel())) continue;
	175	t=tt;
	176	//if ((tt->GetStatus()&AliESDtrack::kTOFpid) == 0) continue;
	177	if (tt->GetTOFsignal() < 0) continue;
	178	UInt_t idx=tt->GetTOFcluster();
	179	if ((Int_t)idx>=nd) {
	180	cerr<<"Wrong digit index ! "<<idx<<endl;
	181	return 5;
	182	}
	183	AliTOFdigit dig=(AliTOFdigit)digits->UncheckedAt(idx);
	184	Int_t *label=dig->GetTracks();
	185	if (label[0]!=lab)
	186	if (label[1]!=lab)
	187	if (label[2]!=lab) {
	188	mismatched++;
	189	if (ptg>pmin) { hfake->Fill(ptg); hfak->Fill(tgl); }
	190	break;
	191	}
	192	if (ptg>pmin) { hfound->Fill(ptg); hfoun->Fill(tgl); }
	193	matched++;
	194	break;
	195	}
	196	if (j==ntrk) {
	197	cerr<<"Not matched: "<<lab<<" ";
	198	if (t) {
	199	cerr<<(t->GetStatus()&AliESDtrack::kITSout)<<' '
	200	<<(t->GetStatus()&AliESDtrack::kTPCout)<<' '
	201	<<(t->GetStatus()&AliESDtrack::kTRDout)<<' '
	202	<<(t->GetStatus()&AliESDtrack::kTIME);
	203	} else cerr<<"No ESD track !";
	204	cerr<<endl;
	205	}
	206	}
	207
208	cerr<<"Number of good tracks: "<<ngood<<endl;
209	cerr<<"Number of matched tracks: "<<matched<<endl;
210	cerr<<"Number of mismatched tracks: "<<mismatched<<endl;
211	if (ngood!=0) cerr<<"Efficiency: "<<Float_t(matched)/ngood<<endl;
212
213	hfound->Sumw2(); hgood->Sumw2(); hfake->Sumw2();
214	hgp->Divide(hfound,hgood,1,1.,"b");
215	hfp->Divide(hfake,hgood,1,1.,"b");
216	hgp->SetMaximum(1.4);
217	hgp->SetYTitle("Matching efficiency");
218	hgp->SetXTitle("Pt (GeV/c)");
219
220	hfoun->Sumw2(); hgoo->Sumw2(); hfak->Sumw2();
221	hgl->Divide(hfoun,hgoo,1,1.,"b");
222	hfl->Divide(hfak,hgoo,1,1.,"b");
223	hgl->SetMaximum(1.4);
224	hgl->SetYTitle("Matching efficiency");
225	hgl->SetXTitle("Tan(lambda)");
226
3f83f224	227	c1->cd(1);
	228
	229	hgp->Draw();
	230	hfp->Draw("histsame");
	231	TLine *line1 = new TLine(pmin,1.0,pmax,1.0); line1->SetLineStyle(4);
	232	line1->Draw("same");
	233	TLine *line2 = new TLine(pmin,0.9,pmax,0.9); line2->SetLineStyle(4);
	234	line2->Draw("same");
	235
	236	c1->cd(2);
	237
	238	hgl->Draw();
	239	hfl->Draw("histsame");
	240	TLine *line3 = new TLine(-1,1.0,1,1.0); line3->SetLineStyle(4);
	241	line3->Draw("same");
	242	TLine *line4 = new TLine(-1,0.9,1,0.9); line4->SetLineStyle(4);
	243	line4->Draw("same");
	244
	245	break;
	246	}
	247
3f7a17bb	248	TFile fc("AliTOFComparison.root","RECREATE");
	249	c1->Write();
	250	fc.Close();
	251
3f83f224	252	return 0;
	253	}
	254
	255	Int_t good_tracks_tof(GoodTrackTOF *gt, const Int_t max) {
	256	ifstream in("good_tracks_its");
	257	if (!in) {
	258	cerr<<"Can't get good_tracks_its !\n"; exit(11);
	259	}
	260
	261	AliRunLoader *rl = AliRunLoader::Open("galice.root","COMPARISON");
	262	if (!rl) {
	263	cerr<<"Can't start session !\n";
	264	exit(1);
	265	}
	266
	267	rl->GetEvent(0);
	268
	269
	270	rl->LoadgAlice();
	271	rl->LoadHeader();
	272	Int_t np = rl->GetHeader()->GetNtrack();
	273
	274	Int_t *good=new Int_t[np];
	275	Int_t k;
	276	for (k=0; k<np; k++) good[k]=0;
	277
	278	AliLoader* tofl = rl->GetLoader("TOFLoader");
	279	if (tofl == 0x0) {
	280	cerr<<"Can not get the TOF loader"<<endl;
	281	exit(2);
	282	}
	283	tofl->LoadDigits("read");
	284
	285	TTree *dTree=tofl->TreeD();
	286	if (!dTree) {
	287	cerr<<"Can't get the TOF cluster tree !\n";
	288	exit(3);
	289	}
	290
	291	TBranch *branch=dTree->GetBranch("TOF");
	292	if (!branch) {
	293	cerr<<"Can't get the branch with the TOF digits !\n";
	294	exit(4);
	295	}
	296	TClonesArray dummy("AliTOFdigit",10000), *digits=&dummy;
	297	branch->SetAddress(&digits);
	298
	299	dTree->GetEvent(0);
	300	Int_t nd=digits->GetEntriesFast();
	301	cerr<<"Number of digits: "<<nd<<endl;
	302
	303	for (Int_t i=0; i<nd; i++) {
	304	AliTOFdigit d=(AliTOFdigit)digits->UncheckedAt(i);
	305	Int_t l0=d->GetTrack(0);
	306	if (l0>=np) {cerr<<"Wrong label: "<<l0<<endl; continue;}
	307	Int_t l1=d->GetTrack(1);
	308	if (l1>=np) {cerr<<"Wrong label: "<<l1<<endl; continue;}
	309	Int_t l2=d->GetTrack(2);
	310	if (l2>=np) {cerr<<"Wrong label: "<<l2<<endl; continue;}
	311	if (l0>=0) good[l0]++;
	312	if (l1>=0) good[l1]++;
	313	if (l2>=0) good[l2]++;
	314	}
	315
316
317	rl->LoadKinematics();
318	AliStack* stack = rl->Stack();
319	Int_t nt=0;
320	Double_t px,py,pz,x,y,z;
321	Int_t code,lab;
322	while (in>>lab>>code>>px>>py>>pz>>x>>y>>z) {
323	if (good[lab] == 0) continue;
324	TParticle p = (TParticle)stack->Particle(lab);
325	if (p == 0x0) {
326	cerr<<"Can not get particle "<<lab<<endl;
327	exit(1);
328	}
329	if (TMath::Abs(p->Vx())>0.1) continue;
330	if (TMath::Abs(p->Vy())>0.1) continue;
2b85d556	331	// if (TMath::Abs(p->Vz())>0.1) continue;
3f83f224	332
	333	gt[nt].lab=lab;
	334	gt[nt].code=p->GetPdgCode();
	335	//**** px py pz - in global coordinate system
	336	gt[nt].px=p->Px(); gt[nt].py=p->Py(); gt[nt].pz=p->Pz();
	337	gt[nt].x=p->Vx(); gt[nt].y=p->Vy(); gt[nt].z=p->Vz();
	338	nt++;
	339	if (nt==max) {cerr<<"Too many good tracks !\n"; break;}
	340	}
	341
	342	delete[] good;
	343
	344	rl->UnloadKinematics();
	345	rl->UnloadHeader();
	346	rl->UnloadgAlice();
	347	delete rl;
	348
	349	return nt;
	350	}