[u/mrichter/AliRoot.git] / HLT / exa / evalrowhough.C

// $Id$

struct GoodTrack 
{
  Int_t label;
  Double_t eta;
  Int_t code;
  Double_t px,py,pz;
  Double_t x,y,z;
  Int_t nhits;
  Int_t sector;
};
typedef struct GoodTrack GoodTrack;

//Histograms
TNtuple *fNtuppel;
TH1F *fPtRes;
TH1F *fNGoodTracksPt;
TH1F *fNFoundTracksPt;
TH1F *fNFakeTracksPt;
TH1F *fTrackEffPt;
TH1F *fFakeTrackEffPt;
TH1F *fNGoodTracksEta;
TH1F *fNFoundTracksEta;
TH1F *fNFakeTracksEta;
TH1F *fTrackEffEta;
TH1F *fFakeTrackEffEta;
TH1F *fNGoodTracksPhi;
TH1F *fNFoundTracksPhi;
TH1F *fNFakeTracksPhi;
TH1F *fTrackEffPhi;
TH1F *fFakeTrackEffPhi;
TProfile *fFakesPt;
TProfile *fFakesPhi;
TProfile *fSlicesPt;
TProfile *fSlicesPhi;
TProfile2D *fFakesPtvsPhi;
TNtuple *fNtupleRes;
TH1F *fPtRes2;
TH1F *fEtaRes;
TH1F *fPsiRes;

TH1F *fNGoodTracksRad;
TH1F *fNFoundTracksRad;
TH1F *fNGoodTracksZ;
TH1F *fNFoundTracksZ;
TH1F *fRadEff;
TH1F *fZEff;

void evalrowhough(Char_t *path="./fitter",Char_t *offlinepath="./",Int_t ievent=0, Float_t ptmin=0.4)  
{
  CreateHistos(25,0.1,4.1);

  Char_t fname[1024];
  sprintf(fname,"%s/tracks_%d.raw",path,ievent);
  AliHLTFileHandler *tfile = new AliHLTFileHandler();
  if(!tfile->SetBinaryInput(fname)){
    LOG(AliHLTLog::kError,"AliHLTEvaluation::Setup","File Open")
      <<"Inputfile "<<fname<<" does not exist"<<ENDLOG; 
    return;
  }

  AliHLTTrackArray *fTracks = new AliHLTTrackArray("AliHLTHoughTrack");
  tfile->Binary2TrackArray(fTracks);
  //fTracks->QSort();
  tfile->CloseBinaryInput();
  delete tfile;

  for(Int_t i=0; i<fTracks->GetNTracks(); i++)
    {
      AliHLTHoughTrack *track = fTracks->GetCheckedTrack(i);
      if(!track) continue; 
      
      track->SetEta(fTracks->GetCheckedTrack(i)->GetPseudoRapidity());
      UInt_t *ids = track->GetHitNumbers();
      Int_t slice = (ids[0]>>25)&0x7f;
      track->SetSlice(slice);
      cout<<"Track with label "<<ids<<" "<<track->GetNHits()<<" "<<track->GetMCid()<<" "<<ids[0]<<" "<<track->GetWeight()<<" "<<track->GetSlice()<<endl;
    }

  Char_t filename[1024];
  sprintf(filename,"%s/good_tracks_tpc_%d",offlinepath,ievent);
  ifstream in(filename);
  if(!in)
    {
      cerr<<"AliHLTEvaluate::GetGoodParticles : Problems opening file :"<<filename<<endl;
      return;
    }

  Int_t MaxTracks=20000;
  Int_t fGoodGen=0;
  fGoodTracks = new GoodTrack[MaxTracks];
  
  while (in>>fGoodTracks[fGoodGen].label>>fGoodTracks[fGoodGen].code>>
	 fGoodTracks[fGoodGen].px>>fGoodTracks[fGoodGen].py>>fGoodTracks[fGoodGen].pz>>
	 fGoodTracks[fGoodGen].x>>fGoodTracks[fGoodGen].y>>fGoodTracks[fGoodGen].z)
    {
      fGoodTracks[fGoodGen].nhits=-1;
      fGoodTracks[fGoodGen].sector=-1; 
      fGoodGen++;
      if (fGoodGen==MaxTracks) 
	{
	  cerr<<"AliHLTEvaluate::GetGoodParticles : Too many good tracks !\n";
	  break;
	}
    }

  Int_t fGoodGen2 = 0;
  Int_t fGoodFound = 0;
  Int_t fGoodFound2 = 0;
  Float_t fMinGoodPt = ptmin;
  Float_t fMaxGoodPt = 99999.;

  if(!fGoodTracks)
    {
      cerr<<"AliHLTEvaluate::FillEffHistos : No good tracks"<<endl;
      return;
    }
  cout<<"Comparing "<<fGoodGen<<" good tracks ..."<<endl;
  for(Int_t i=0; i<fGoodGen; i++)
    {
      //cout<<"Checking particle "<<i<<endl;
      Float_t ptg = TMath::Sqrt(fGoodTracks[i].px*fGoodTracks[i].px + fGoodTracks[i].py*fGoodTracks[i].py);
      Float_t pg = TMath::Sqrt(fGoodTracks[i].px*fGoodTracks[i].px + fGoodTracks[i].py*fGoodTracks[i].py + fGoodTracks[i].pz*fGoodTracks[i].pz);
      if(ptg < fMinGoodPt || ptg > fMaxGoodPt) continue;
      Float_t pzg=fGoodTracks[i].pz;
      Float_t dipangle=TMath::ATan2(pzg,ptg)*180./TMath::Pi();
      Float_t etag=0.5 * TMath::Log((pg+pzg)/(pg-pzg));
      Float_t phig=TMath::ATan2(fGoodTracks[i].py,fGoodTracks[i].px);
      Float_t rg = TMath::Sqrt(fGoodTracks[i].x*fGoodTracks[i].x + fGoodTracks[i].y*fGoodTracks[i].y);
      Float_t zg = TMath::Abs(fGoodTracks[i].z);
      if (phig<0) phig+=2*TMath::Pi();

      fNGoodTracksPt->Fill(ptg);
      fNGoodTracksEta->Fill(etag);
      fNGoodTracksPhi->Fill(phig);
      if((zg != 0) || (rg != 0)) { 
	fNGoodTracksRad->Fill(rg);
	fNGoodTracksZ->Fill(zg);
      }
      fGoodGen2++;
      Int_t found = 0;
      Int_t nslices = 1;
      Int_t org_slice;
      for(Int_t k=0; k<fTracks->GetNTracks(); k++)
	{
	  AliHLTHoughTrack *track = (AliHLTHoughTrack *)fTracks->GetCheckedTrack(k);
	  if(!track) continue;
	  //	  Int_t nHits = track->GetNumberOfPoints();
	  //	  if(nHits < fMinPointsOnTrack) break;
	  Int_t nHits = 1;
	  Int_t tracklabel;
	  tracklabel = track->GetMCid();
	  
	  if(TMath::Abs(tracklabel) != fGoodTracks[i].label) continue;
	  if(found == 0) {
	    found=1;
	    Float_t pt=track->GetPt();
	    org_slice = track->GetSlice();
	    Float_t eta = track->GetEta();
	    Float_t phi = track->GetPsi();
	    //cout<<eta<<" "<<phi<<" "<<etag<<" "<<phig<<" "<<track->GetEtaIndex()<<endl;
	    if(tracklabel == fGoodTracks[i].label) 
	      {
		fGoodFound++;
		
		fNFoundTracksPt->Fill(ptg); 
		fNFoundTracksEta->Fill(etag);
		fNFoundTracksPhi->Fill(phig);
		if((zg != 0) || (rg != 0)) { 
		  fNFoundTracksRad->Fill(rg);
		  fNFoundTracksZ->Fill(zg);
		}
		fNtuppel->Fill(ptg,pt,nHits);
		fPtRes->Fill((pt-ptg)/ptg*100.);
		fEtaRes->Fill(eta-etag);
		fPsiRes->Fill(phi-phig);
	      }
	    else 
	      {
		fNFakeTracksPt->Fill(ptg); 
		fNFakeTracksEta->Fill(etag);
		fNFakeTracksPhi->Fill(phig);
	      }
	    //fPtRes->Fill((pt-ptg)/ptg*100.);
	    //fNtuppel->Fill(ptg,pt,nHits);
	  }
	  else {
	    found++;
	    if(track->GetSlice() != org_slice)
	      nslices = 2;
	  }
	}
      if(!found)
	//cout<<"Track "<<fGoodTracks[i].label<<" was not found"<<endl;
	continue;
      else {
	fFakesPt->Fill(ptg,(Double_t)found-1.0);
	fFakesPhi->Fill(phig,(Double_t)found-1.0);
	fSlicesPt->Fill(ptg,(Double_t)nslices);
	fSlicesPhi->Fill(phig,(Double_t)nslices);
	//Float_t local_phi = phi-(phi/)
	//fFakesPtvsPhi->Fill(1.0/ptg,
      }
      fGoodFound2 += found;
    }

  cout<<fGoodFound<<"("<<fGoodFound2<<") tracks found out of "<<fGoodGen2<<" generated tracks"<<endl;
  cout<<" The overall efficiency is "<<(Float_t)fGoodFound/(Float_t)fGoodGen2<<endl;
  
  CalcEffHistos();
  plotptres(ptmin);
  Write2File();
}

void CreateHistos(Int_t nbin,Float_t xlow,Float_t xup)
{
  //Create the histograms 
  
  fNtuppel = new TNtuple("fNtuppel","Pt resolution","pt_gen:pt_found:nHits");
  fNtuppel->SetDirectory(0);
  fPtRes = new TH1F("fPtRes","Relative Pt resolution",30,-10.,10.); 
  fNGoodTracksPt = new TH1F("fNGoodTracksPt","Good tracks vs pt",nbin,xlow,xup);    
  fNFoundTracksPt = new TH1F("fNFoundTracksPt","Found tracks vs pt",nbin,xlow,xup);
  fNFakeTracksPt = new TH1F("fNFakeTracksPt","Fake tracks vs pt",nbin,xlow,xup);
  fTrackEffPt = new TH1F("fTrackEffPt","Tracking efficiency vs pt",nbin,xlow,xup);
  fFakeTrackEffPt = new TH1F("fFakeTrackEffPt","Efficiency for fake tracks vs pt",nbin,xlow,xup);
  
  fNGoodTracksEta = new TH1F("fNGoodTracksEta","Good tracks vs eta",nbin,-1,1);
  fNFoundTracksEta = new TH1F("fNFoundTracksEta","Found tracks vs eta",nbin,-1,1);
  fNFakeTracksEta = new TH1F("fNFakeTracksEta","Fake tracks vs eta",nbin,-1,1);
  fTrackEffEta = new TH1F("fTrackEffEta","Tracking efficienct vs eta",nbin,-1,1);
  fFakeTrackEffEta = new TH1F("fFakeTrackEffEta","Efficiency for fake tracks vs eta",nbin,-1,1);

  
  fNGoodTracksPhi = new TH1F("fNGoodTracksPhi","Good tracks vs phi",nbin,0,6.28);
  fNFoundTracksPhi = new TH1F("fNFoundTracksPhi","Found tracks vs phi",nbin,0,6.28);
  fNFakeTracksPhi = new TH1F("fNFakeTracksPhi","Fake tracks vs phi",nbin,0,6.28);
  fTrackEffPhi = new TH1F("fTrackEffPhi","Tracking efficienct vs phi",nbin,0,6.28);
  fFakeTrackEffPhi = new TH1F("fFakeTrackEffPhi","Efficiency for fake tracks vs phi",nbin,0,6.28);
  fFakesPt = new TProfile("fFakesPt","Number of ghosts vs pt",nbin,xlow,xup,0,100);
  fFakesPhi = new TProfile("fFakesPhi","Number of ghosts vs phi",nbin,0,6.28,0,100);
  fSlicesPt = new TProfile("fSlicesPt","Number of slices vs pt",nbin,xlow,xup,0,100);
  fSlicesPhi = new TProfile("fSlicesPhi","Number of slices vs phi",nbin,0,6.28,0,100);
  fPtRes2 = new TH1F("fPtRes2","Relative Pt resolution (Pt>2 GeV)",20,-100.,100.);
  fEtaRes = new TH1F("fEtaRes","Eta resolution",30,-0.1,0.1);
  fPsiRes = new TH1F("fPsiRes","Psi resolution",30,-0.1,0.1);

  fNGoodTracksRad = new TH1F("fNGoodTracksRad","Good tracks vs distance(x,y)",1,0,0.2);    
  fNFoundTracksRad = new TH1F("fNFoundTracksRad","Found tracks vs distance(x,y)",1,0,0.2);
  fRadEff = new TH1F("fRadEff","Tracking efficiency vs distance(x,y)",1,0,0.2);
  fNGoodTracksZ = new TH1F("fNGoodTracksZ","Good tracks vs distance(z)",1,0,0.2);    
  fNFoundTracksZ = new TH1F("fNFoundTracksZ","Found tracks vs distance(z)",1,0,0.2);
  fZEff = new TH1F("fZEff","Tracking efficiency vs distance(z)",1,0,0.2);
}

void Write2File()
{
  //Write histograms to file:
  
  TFile *of = TFile::Open("hough_efficiencies.root","RECREATE");
  if(!of->IsOpen())
    {
      cout<<"Problems opening rootfile"<<endl;
      return;
    }
  
  of->cd();
  fNtuppel->Write();
  fPtRes->Write();
  fNGoodTracksPt->Write();
  fNFoundTracksPt->Write();
  fNFakeTracksPt->Write();
  fTrackEffPt->Write();
  fFakeTrackEffPt->Write();
  fNGoodTracksEta->Write();
  fNFoundTracksEta->Write();
  fNFakeTracksEta->Write();
  fTrackEffEta->Write();
  fFakeTrackEffEta->Write();

  fNGoodTracksPhi->Write();
  fNFoundTracksPhi->Write();
  fNFakeTracksPhi->Write();
  fTrackEffPhi->Write();
  fFakeTrackEffPhi->Write();
  
  fFakesPt->Write();
  fFakesPhi->Write();
  fSlicesPt->Write();
  fSlicesPhi->Write();

  fPtRes2->Write();
  fEtaRes->Write();
  fPsiRes->Write();

  fNGoodTracksRad->Write();
  fNFoundTracksRad->Write(); 
  fRadEff->Write();
  fNGoodTracksZ->Write();
  fNFoundTracksZ->Write();
  fZEff->Write();

  of->Close();
}

void CalcEffHistos()
{  
  fNFoundTracksPt->Sumw2(); fNGoodTracksPt->Sumw2();
  fTrackEffPt->Divide(fNFoundTracksPt,fNGoodTracksPt,1,1.,"b");
  fFakeTrackEffPt->Divide(fNFakeTracksPt,fNGoodTracksPt,1,1.,"b");
  fTrackEffPt->SetMaximum(1.4);
  fTrackEffPt->SetXTitle("P_{T} [GeV]");
  fTrackEffPt->SetLineWidth(2);
  fFakeTrackEffPt->SetFillStyle(3013);
  fTrackEffPt->SetLineColor(4);
  fFakeTrackEffPt->SetFillColor(2);

  fNFoundTracksEta->Sumw2(); fNGoodTracksEta->Sumw2();
  fTrackEffEta->Divide(fNFoundTracksEta,fNGoodTracksEta,1,1.,"b");
  fFakeTrackEffEta->Divide(fNFakeTracksEta,fNGoodTracksEta,1,1.,"b");
  fTrackEffEta->SetMaximum(1.4);
  fTrackEffEta->SetXTitle("#eta");
  fTrackEffEta->SetLineWidth(2);
  fFakeTrackEffEta->SetFillStyle(3013);
  fTrackEffEta->SetLineColor(4);
  fFakeTrackEffEta->SetFillColor(2);
       
  fNFoundTracksPhi->Sumw2(); fNGoodTracksPhi->Sumw2();
  fTrackEffPhi->Divide(fNFoundTracksPhi,fNGoodTracksPhi,1,1.,"b");
  fFakeTrackEffPhi->Divide(fNFakeTracksPhi,fNGoodTracksPhi,1,1.,"b");
  fTrackEffPhi->SetMaximum(1.4);
  fTrackEffPhi->SetXTitle("#psi [rad]");
  fTrackEffPhi->SetLineWidth(2);
  fFakeTrackEffPhi->SetFillStyle(3013);
  fTrackEffPhi->SetLineColor(4);
  fFakeTrackEffPhi->SetFillColor(2);

  fEtaRes->SetXTitle("#eta");
  fPsiRes->SetXTitle("#psi [rad]");

  fRadEff->Divide(fNFoundTracksRad,fNGoodTracksRad,1,1.,"b");
  fRadEff->SetXTitle("Distance [cm]");
  fZEff->Divide(fNFoundTracksZ,fNGoodTracksZ,1,1.,"b");
  fZEff->SetXTitle("Distance [cm]");
}

void plotptres(Float_t ptmin)
{
  //Plot the relative pt resolution vs pt.
  
  const Int_t n=8;
  
  Double_t hltx[9] = {0.2,0.4,0.6,0.8,1.0,1.2,1.4,1.6,2.0};
  Double_t hltxerr[n];
  Double_t hltavx[n];
  Double_t hlty[n];
  Double_t hltyerr[n];
  Char_t string[1024];

  for(Int_t i=0; i<n; i++)
    {
      TH1F *hist = new TH1F("hist","",100,-50.*(float)(i+1),50.*(float)(i+1));
      sprintf(string,"pt_gen > %f && pt_gen <= %f",hltx[i],hltx[i+1]);
      fNtuppel->Draw("(pt_found-pt_gen)/pt_gen*100>>hist",string,"goff");
      TF1 *f1 = new TF1("f1","gaus");
      hist->Fit("f1","E");
      hlty[i] = f1->GetParameter(2);
      hltyerr[i] = f1->GetParError(2);
      //hlty[i] = hist->GetRMS();
      //hltyerr[i] = 0;
      hltavx[i] = (hltx[i]+hltx[i+1])/2.0;
      hltxerr[i] = (hltx[i+1]-hltx[i])/2.0;
      cout<<string<<" "<<i<<" "<<hlty[i]<<" "<<hltyerr[i]<<endl;
      delete f1;
      delete hist;
    }

  sprintf(string,"pt_gen > %f",2.0);
  fNtuppel->Draw("(pt_found-pt_gen)/pt_gen*100>>fPtRes2",string,"goff");
  TF1 *f1 = new TF1("f1","gaus");
  fPtRes2->Fit("f1","E");
  fPtRes2->SetXTitle("#Delta P_{t} / P_{t} [%]");
  fPtRes2->SetYTitle("#Delta N");
  delete f1;

  TGraphErrors *gr1 = new TGraphErrors(n,hltavx,hlty,hltxerr,hltyerr);
  TCanvas *c1 = new TCanvas("c1","",1);
  fPtRes = c1->DrawFrame(0.1,0,2,15);
  gr1->SetTitle("");
  gr1->Draw("P");
  gr1->SetMarkerStyle(20);
  gr1->SetLineWidth(2);
  gr1->SetMarkerSize(1.3);
  f1 = new TF1("f1","pol1",ptmin,2.0);
  gr1->Fit("f1","ER");
  fPtRes->SetXTitle("p_{t} [GeV]");
  fPtRes->SetYTitle("#Delta P_{t} / P_{t} [%]");
  delete f1;
}
Commit	Line	Data
0bd0c1ef	1	// $Id$
	2
	3	struct GoodTrack
	4	{
	5	Int_t label;
	6	Double_t eta;
	7	Int_t code;
	8	Double_t px,py,pz;
	9	Double_t x,y,z;
	10	Int_t nhits;
	11	Int_t sector;
	12	};
	13	typedef struct GoodTrack GoodTrack;
	14
	15	//Histograms
	16	TNtuple *fNtuppel;
	17	TH1F *fPtRes;
	18	TH1F *fNGoodTracksPt;
	19	TH1F *fNFoundTracksPt;
	20	TH1F *fNFakeTracksPt;
	21	TH1F *fTrackEffPt;
	22	TH1F *fFakeTrackEffPt;
	23	TH1F *fNGoodTracksEta;
	24	TH1F *fNFoundTracksEta;
	25	TH1F *fNFakeTracksEta;
	26	TH1F *fTrackEffEta;
	27	TH1F *fFakeTrackEffEta;
	28	TH1F *fNGoodTracksPhi;
	29	TH1F *fNFoundTracksPhi;
	30	TH1F *fNFakeTracksPhi;
	31	TH1F *fTrackEffPhi;
	32	TH1F *fFakeTrackEffPhi;
	33	TProfile *fFakesPt;
	34	TProfile *fFakesPhi;
	35	TProfile *fSlicesPt;
	36	TProfile *fSlicesPhi;
	37	TProfile2D *fFakesPtvsPhi;
	38	TNtuple *fNtupleRes;
	39	TH1F *fPtRes2;
	40	TH1F *fEtaRes;
	41	TH1F *fPsiRes;
	42
	43	TH1F *fNGoodTracksRad;
	44	TH1F *fNFoundTracksRad;
	45	TH1F *fNGoodTracksZ;
	46	TH1F *fNFoundTracksZ;
	47	TH1F *fRadEff;
	48	TH1F *fZEff;
	49
	50	void evalrowhough(Char_t path="./fitter",Char_t offlinepath="./",Int_t ievent=0, Float_t ptmin=0.4)
	51	{
	52	CreateHistos(25,0.1,4.1);
	53
	54	Char_t fname[1024];
	55	sprintf(fname,"%s/tracks_%d.raw",path,ievent);
4aa41877	56	AliHLTFileHandler *tfile = new AliHLTFileHandler();
0bd0c1ef	57	if(!tfile->SetBinaryInput(fname)){
4aa41877	58	LOG(AliHLTLog::kError,"AliHLTEvaluation::Setup","File Open")
0bd0c1ef	59	<<"Inputfile "<<fname<<" does not exist"<<ENDLOG;
	60	return;
	61	}
	62
4aa41877	63	AliHLTTrackArray *fTracks = new AliHLTTrackArray("AliHLTHoughTrack");
0bd0c1ef	64	tfile->Binary2TrackArray(fTracks);
	65	//fTracks->QSort();
	66	tfile->CloseBinaryInput();
	67	delete tfile;
	68
	69	for(Int_t i=0; i<fTracks->GetNTracks(); i++)
	70	{
4aa41877	71	AliHLTHoughTrack *track = fTracks->GetCheckedTrack(i);
0bd0c1ef	72	if(!track) continue;
	73
	74	track->SetEta(fTracks->GetCheckedTrack(i)->GetPseudoRapidity());
	75	UInt_t *ids = track->GetHitNumbers();
	76	Int_t slice = (ids[0]>>25)&0x7f;
	77	track->SetSlice(slice);
	78	cout<<"Track with label "<<ids<<" "<<track->GetNHits()<<" "<<track->GetMCid()<<" "<<ids[0]<<" "<<track->GetWeight()<<" "<<track->GetSlice()<<endl;
	79	}
	80
	81	Char_t filename[1024];
	82	sprintf(filename,"%s/good_tracks_tpc_%d",offlinepath,ievent);
	83	ifstream in(filename);
	84	if(!in)
	85	{
4aa41877	86	cerr<<"AliHLTEvaluate::GetGoodParticles : Problems opening file :"<<filename<<endl;
0bd0c1ef	87	return;
	88	}
	89
	90	Int_t MaxTracks=20000;
	91	Int_t fGoodGen=0;
	92	fGoodTracks = new GoodTrack[MaxTracks];
	93
	94	while (in>>fGoodTracks[fGoodGen].label>>fGoodTracks[fGoodGen].code>>
	95	fGoodTracks[fGoodGen].px>>fGoodTracks[fGoodGen].py>>fGoodTracks[fGoodGen].pz>>
	96	fGoodTracks[fGoodGen].x>>fGoodTracks[fGoodGen].y>>fGoodTracks[fGoodGen].z)
	97	{
	98	fGoodTracks[fGoodGen].nhits=-1;
	99	fGoodTracks[fGoodGen].sector=-1;
	100	fGoodGen++;
	101	if (fGoodGen==MaxTracks)
	102	{
4aa41877	103	cerr<<"AliHLTEvaluate::GetGoodParticles : Too many good tracks !\n";
0bd0c1ef	104	break;
	105	}
	106	}
	107
	108	Int_t fGoodGen2 = 0;
	109	Int_t fGoodFound = 0;
	110	Int_t fGoodFound2 = 0;
	111	Float_t fMinGoodPt = ptmin;
	112	Float_t fMaxGoodPt = 99999.;
	113
	114	if(!fGoodTracks)
	115	{
4aa41877	116	cerr<<"AliHLTEvaluate::FillEffHistos : No good tracks"<<endl;
0bd0c1ef	117	return;
	118	}
	119	cout<<"Comparing "<<fGoodGen<<" good tracks ..."<<endl;
	120	for(Int_t i=0; i<fGoodGen; i++)
	121	{
	122	//cout<<"Checking particle "<<i<<endl;
	123	Float_t ptg = TMath::Sqrt(fGoodTracks[i].pxfGoodTracks[i].px + fGoodTracks[i].pyfGoodTracks[i].py);
	124	Float_t pg = TMath::Sqrt(fGoodTracks[i].pxfGoodTracks[i].px + fGoodTracks[i].pyfGoodTracks[i].py + fGoodTracks[i].pz*fGoodTracks[i].pz);
	125	if(ptg < fMinGoodPt \|\| ptg > fMaxGoodPt) continue;
	126	Float_t pzg=fGoodTracks[i].pz;
	127	Float_t dipangle=TMath::ATan2(pzg,ptg)*180./TMath::Pi();
	128	Float_t etag=0.5 * TMath::Log((pg+pzg)/(pg-pzg));
	129	Float_t phig=TMath::ATan2(fGoodTracks[i].py,fGoodTracks[i].px);
	130	Float_t rg = TMath::Sqrt(fGoodTracks[i].xfGoodTracks[i].x + fGoodTracks[i].yfGoodTracks[i].y);
	131	Float_t zg = TMath::Abs(fGoodTracks[i].z);
	132	if (phig<0) phig+=2*TMath::Pi();
	133
	134	fNGoodTracksPt->Fill(ptg);
	135	fNGoodTracksEta->Fill(etag);
	136	fNGoodTracksPhi->Fill(phig);
	137	if((zg != 0) \|\| (rg != 0)) {
	138	fNGoodTracksRad->Fill(rg);
	139	fNGoodTracksZ->Fill(zg);
	140	}
	141	fGoodGen2++;
	142	Int_t found = 0;
	143	Int_t nslices = 1;
	144	Int_t org_slice;
	145	for(Int_t k=0; k<fTracks->GetNTracks(); k++)
	146	{
4aa41877	147	AliHLTHoughTrack track = (AliHLTHoughTrack )fTracks->GetCheckedTrack(k);
0bd0c1ef	148	if(!track) continue;
	149	// Int_t nHits = track->GetNumberOfPoints();
	150	// if(nHits < fMinPointsOnTrack) break;
	151	Int_t nHits = 1;
	152	Int_t tracklabel;
	153	tracklabel = track->GetMCid();
	154
	155	if(TMath::Abs(tracklabel) != fGoodTracks[i].label) continue;
	156	if(found == 0) {
	157	found=1;
	158	Float_t pt=track->GetPt();
	159	org_slice = track->GetSlice();
	160	Float_t eta = track->GetEta();
	161	Float_t phi = track->GetPsi();
	162	//cout<<eta<<" "<<phi<<" "<<etag<<" "<<phig<<" "<<track->GetEtaIndex()<<endl;
	163	if(tracklabel == fGoodTracks[i].label)
	164	{
	165	fGoodFound++;
	166
	167	fNFoundTracksPt->Fill(ptg);
	168	fNFoundTracksEta->Fill(etag);
	169	fNFoundTracksPhi->Fill(phig);
	170	if((zg != 0) \|\| (rg != 0)) {
	171	fNFoundTracksRad->Fill(rg);
	172	fNFoundTracksZ->Fill(zg);
	173	}
	174	fNtuppel->Fill(ptg,pt,nHits);
	175	fPtRes->Fill((pt-ptg)/ptg*100.);
	176	fEtaRes->Fill(eta-etag);
	177	fPsiRes->Fill(phi-phig);
	178	}
	179	else
	180	{
	181	fNFakeTracksPt->Fill(ptg);
	182	fNFakeTracksEta->Fill(etag);
	183	fNFakeTracksPhi->Fill(phig);
	184	}
	185	//fPtRes->Fill((pt-ptg)/ptg*100.);
	186	//fNtuppel->Fill(ptg,pt,nHits);
	187	}
	188	else {
	189	found++;
	190	if(track->GetSlice() != org_slice)
	191	nslices = 2;
	192	}
	193	}
	194	if(!found)
	195	//cout<<"Track "<<fGoodTracks[i].label<<" was not found"<<endl;
	196	continue;
	197	else {
	198	fFakesPt->Fill(ptg,(Double_t)found-1.0);
	199	fFakesPhi->Fill(phig,(Double_t)found-1.0);
	200	fSlicesPt->Fill(ptg,(Double_t)nslices);
	201	fSlicesPhi->Fill(phig,(Double_t)nslices);
	202	//Float_t local_phi = phi-(phi/)
	203	//fFakesPtvsPhi->Fill(1.0/ptg,
	204	}
	205	fGoodFound2 += found;
	206	}
	207
	208	cout<<fGoodFound<<"("<<fGoodFound2<<") tracks found out of "<<fGoodGen2<<" generated tracks"<<endl;
	209	cout<<" The overall efficiency is "<<(Float_t)fGoodFound/(Float_t)fGoodGen2<<endl;
	210
	211	CalcEffHistos();
212	plotptres(ptmin);
213	Write2File();
214	}
215
216	void CreateHistos(Int_t nbin,Float_t xlow,Float_t xup)
217	{
218	//Create the histograms
219
220	fNtuppel = new TNtuple("fNtuppel","Pt resolution","pt_gen:pt_found:nHits");
221	fNtuppel->SetDirectory(0);
222	fPtRes = new TH1F("fPtRes","Relative Pt resolution",30,-10.,10.);
223	fNGoodTracksPt = new TH1F("fNGoodTracksPt","Good tracks vs pt",nbin,xlow,xup);
224	fNFoundTracksPt = new TH1F("fNFoundTracksPt","Found tracks vs pt",nbin,xlow,xup);
225	fNFakeTracksPt = new TH1F("fNFakeTracksPt","Fake tracks vs pt",nbin,xlow,xup);
226	fTrackEffPt = new TH1F("fTrackEffPt","Tracking efficiency vs pt",nbin,xlow,xup);
227	fFakeTrackEffPt = new TH1F("fFakeTrackEffPt","Efficiency for fake tracks vs pt",nbin,xlow,xup);
228
229	fNGoodTracksEta = new TH1F("fNGoodTracksEta","Good tracks vs eta",nbin,-1,1);
230	fNFoundTracksEta = new TH1F("fNFoundTracksEta","Found tracks vs eta",nbin,-1,1);
231	fNFakeTracksEta = new TH1F("fNFakeTracksEta","Fake tracks vs eta",nbin,-1,1);
232	fTrackEffEta = new TH1F("fTrackEffEta","Tracking efficienct vs eta",nbin,-1,1);
233	fFakeTrackEffEta = new TH1F("fFakeTrackEffEta","Efficiency for fake tracks vs eta",nbin,-1,1);
234
235
236	fNGoodTracksPhi = new TH1F("fNGoodTracksPhi","Good tracks vs phi",nbin,0,6.28);
237	fNFoundTracksPhi = new TH1F("fNFoundTracksPhi","Found tracks vs phi",nbin,0,6.28);
238	fNFakeTracksPhi = new TH1F("fNFakeTracksPhi","Fake tracks vs phi",nbin,0,6.28);
239	fTrackEffPhi = new TH1F("fTrackEffPhi","Tracking efficienct vs phi",nbin,0,6.28);
240	fFakeTrackEffPhi = new TH1F("fFakeTrackEffPhi","Efficiency for fake tracks vs phi",nbin,0,6.28);
241	fFakesPt = new TProfile("fFakesPt","Number of ghosts vs pt",nbin,xlow,xup,0,100);
242	fFakesPhi = new TProfile("fFakesPhi","Number of ghosts vs phi",nbin,0,6.28,0,100);
243	fSlicesPt = new TProfile("fSlicesPt","Number of slices vs pt",nbin,xlow,xup,0,100);
244	fSlicesPhi = new TProfile("fSlicesPhi","Number of slices vs phi",nbin,0,6.28,0,100);
245	fPtRes2 = new TH1F("fPtRes2","Relative Pt resolution (Pt>2 GeV)",20,-100.,100.);
246	fEtaRes = new TH1F("fEtaRes","Eta resolution",30,-0.1,0.1);
247	fPsiRes = new TH1F("fPsiRes","Psi resolution",30,-0.1,0.1);
248
249	fNGoodTracksRad = new TH1F("fNGoodTracksRad","Good tracks vs distance(x,y)",1,0,0.2);
250	fNFoundTracksRad = new TH1F("fNFoundTracksRad","Found tracks vs distance(x,y)",1,0,0.2);
251	fRadEff = new TH1F("fRadEff","Tracking efficiency vs distance(x,y)",1,0,0.2);
252	fNGoodTracksZ = new TH1F("fNGoodTracksZ","Good tracks vs distance(z)",1,0,0.2);
253	fNFoundTracksZ = new TH1F("fNFoundTracksZ","Found tracks vs distance(z)",1,0,0.2);
254	fZEff = new TH1F("fZEff","Tracking efficiency vs distance(z)",1,0,0.2);
255	}
256
257	void Write2File()
258	{
259	//Write histograms to file:
260
261	TFile *of = TFile::Open("hough_efficiencies.root","RECREATE");
262	if(!of->IsOpen())
263	{
264	cout<<"Problems opening rootfile"<<endl;
265	return;
266	}
267
268	of->cd();
269	fNtuppel->Write();
270	fPtRes->Write();
271	fNGoodTracksPt->Write();
272	fNFoundTracksPt->Write();
273	fNFakeTracksPt->Write();
274	fTrackEffPt->Write();
275	fFakeTrackEffPt->Write();
276	fNGoodTracksEta->Write();
277	fNFoundTracksEta->Write();
278	fNFakeTracksEta->Write();
279	fTrackEffEta->Write();
280	fFakeTrackEffEta->Write();
281
282	fNGoodTracksPhi->Write();
283	fNFoundTracksPhi->Write();
284	fNFakeTracksPhi->Write();
285	fTrackEffPhi->Write();
286	fFakeTrackEffPhi->Write();
287
288	fFakesPt->Write();
289	fFakesPhi->Write();
290	fSlicesPt->Write();
291	fSlicesPhi->Write();
292
293	fPtRes2->Write();
294	fEtaRes->Write();
295	fPsiRes->Write();
296
297	fNGoodTracksRad->Write();
298	fNFoundTracksRad->Write();
299	fRadEff->Write();
300	fNGoodTracksZ->Write();
301	fNFoundTracksZ->Write();
302	fZEff->Write();
303
304	of->Close();
305	}
306
307	void CalcEffHistos()
308	{
309	fNFoundTracksPt->Sumw2(); fNGoodTracksPt->Sumw2();
310	fTrackEffPt->Divide(fNFoundTracksPt,fNGoodTracksPt,1,1.,"b");
311	fFakeTrackEffPt->Divide(fNFakeTracksPt,fNGoodTracksPt,1,1.,"b");
312	fTrackEffPt->SetMaximum(1.4);
313	fTrackEffPt->SetXTitle("P_{T} [GeV]");
314	fTrackEffPt->SetLineWidth(2);
315	fFakeTrackEffPt->SetFillStyle(3013);
316	fTrackEffPt->SetLineColor(4);
317	fFakeTrackEffPt->SetFillColor(2);
318
319	fNFoundTracksEta->Sumw2(); fNGoodTracksEta->Sumw2();
320	fTrackEffEta->Divide(fNFoundTracksEta,fNGoodTracksEta,1,1.,"b");
321	fFakeTrackEffEta->Divide(fNFakeTracksEta,fNGoodTracksEta,1,1.,"b");
322	fTrackEffEta->SetMaximum(1.4);
323	fTrackEffEta->SetXTitle("#eta");
324	fTrackEffEta->SetLineWidth(2);
325	fFakeTrackEffEta->SetFillStyle(3013);
326	fTrackEffEta->SetLineColor(4);
327	fFakeTrackEffEta->SetFillColor(2);
328
329	fNFoundTracksPhi->Sumw2(); fNGoodTracksPhi->Sumw2();
330	fTrackEffPhi->Divide(fNFoundTracksPhi,fNGoodTracksPhi,1,1.,"b");
331	fFakeTrackEffPhi->Divide(fNFakeTracksPhi,fNGoodTracksPhi,1,1.,"b");
332	fTrackEffPhi->SetMaximum(1.4);
333	fTrackEffPhi->SetXTitle("#psi [rad]");
334	fTrackEffPhi->SetLineWidth(2);
335	fFakeTrackEffPhi->SetFillStyle(3013);
336	fTrackEffPhi->SetLineColor(4);
337	fFakeTrackEffPhi->SetFillColor(2);
338
339	fEtaRes->SetXTitle("#eta");
340	fPsiRes->SetXTitle("#psi [rad]");
341
342	fRadEff->Divide(fNFoundTracksRad,fNGoodTracksRad,1,1.,"b");
343	fRadEff->SetXTitle("Distance [cm]");
344	fZEff->Divide(fNFoundTracksZ,fNGoodTracksZ,1,1.,"b");
345	fZEff->SetXTitle("Distance [cm]");
346	}
347
348	void plotptres(Float_t ptmin)
349	{
350	//Plot the relative pt resolution vs pt.
351
352	const Int_t n=8;
353
354	Double_t hltx[9] = {0.2,0.4,0.6,0.8,1.0,1.2,1.4,1.6,2.0};
355	Double_t hltxerr[n];
356	Double_t hltavx[n];
357	Double_t hlty[n];
358	Double_t hltyerr[n];
359	Char_t string[1024];
360
361	for(Int_t i=0; i<n; i++)
362	{
363	TH1F hist = new TH1F("hist","",100,-50.(float)(i+1),50.*(float)(i+1));
364	sprintf(string,"pt_gen > %f && pt_gen <= %f",hltx[i],hltx[i+1]);
365	fNtuppel->Draw("(pt_found-pt_gen)/pt_gen*100>>hist",string,"goff");
366	TF1 *f1 = new TF1("f1","gaus");
367	hist->Fit("f1","E");
368	hlty[i] = f1->GetParameter(2);
369	hltyerr[i] = f1->GetParError(2);
370	//hlty[i] = hist->GetRMS();
371	//hltyerr[i] = 0;
372	hltavx[i] = (hltx[i]+hltx[i+1])/2.0;
373	hltxerr[i] = (hltx[i+1]-hltx[i])/2.0;
374	cout<<string<<" "<<i<<" "<<hlty[i]<<" "<<hltyerr[i]<<endl;
375	delete f1;
376	delete hist;
377	}
378
379	sprintf(string,"pt_gen > %f",2.0);
380	fNtuppel->Draw("(pt_found-pt_gen)/pt_gen*100>>fPtRes2",string,"goff");
381	TF1 *f1 = new TF1("f1","gaus");
382	fPtRes2->Fit("f1","E");
383	fPtRes2->SetXTitle("#Delta P_{t} / P_{t} [%]");
384	fPtRes2->SetYTitle("#Delta N");
385	delete f1;
386
387	TGraphErrors *gr1 = new TGraphErrors(n,hltavx,hlty,hltxerr,hltyerr);
388	TCanvas *c1 = new TCanvas("c1","",1);
389	fPtRes = c1->DrawFrame(0.1,0,2,15);
390	gr1->SetTitle("");
391	gr1->Draw("P");
392	gr1->SetMarkerStyle(20);
393	gr1->SetLineWidth(2);
394	gr1->SetMarkerSize(1.3);
395	f1 = new TF1("f1","pol1",ptmin,2.0);
396	gr1->Fit("f1","ER");
397	fPtRes->SetXTitle("p_{t} [GeV]");
398	fPtRes->SetYTitle("#Delta P_{t} / P_{t} [%]");
399	delete f1;
400	}