[u/mrichter/AliRoot.git] / TRD / Cal / AliTRDCalTrkAttach.cxx

/**************************************************************************
* 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.                  *
**************************************************************************/

//////////////////////////////////////////////////////////////////////////
//                                                                      //
// Container for calibration parameters for AliTRDseedV1::AttachClusters()
// .... Longer description of content ...
//
// For calibration procedure check AliTRDtrackleOflHelper::CalibrateAttach()
// .... some reference to the calibration procedure .........
//                                                                      //
// Authors:                                                             //
//   Alex Bercuci <a.bercuci@gsi.de>                                    //
//                                                                      //
//////////////////////////////////////////////////////////////////////////

#include <TFile.h>
#include <TROOT.h>
#include <TMath.h>
#include "TH1.h"
#include "TGraphErrors.h"
#include "TObjArray.h"

#include "AliLog.h"

#include "AliTRDcalibDB.h"
#include "AliTRDgeometry.h"
#include "AliTRDCalTrkAttach.h"

ClassImp(AliTRDCalTrkAttach)

//______________________________________________________________________________
AliTRDCalTrkAttach::AliTRDCalTrkAttach()
  :TNamed("AliTRDCalTrkAttach", "Calibration of AliTRDseedV1::AttachClusters")
  ,fRClikeLimit(0.65)
  ,fScaleCov(2.)
  ,fLike(NULL)
{
// Default constructor

  fNsgmDy[0] = 5; fNsgmDy[1] = 7;
  fLikeMinRelDecrease[0] = .2; fLikeMinRelDecrease[1] = .3;
}

//______________________________________________________________________________
AliTRDCalTrkAttach::~AliTRDCalTrkAttach()
{
// Destructor
  if(fLike) delete fLike;
}


//______________________________________________________________________________
Double_t AliTRDCalTrkAttach::CookLikelihood(Bool_t chg, Int_t ly, Float_t pt, Float_t phiTrk, Int_t n, Double_t dyr, Double_t dphi, Double_t sr) const
{
// Calculate likelihood for a segment to belong to a tracklet
// Based on calibrated values

  if(n<4){
    AliDebug(4, Form("Failed basic cut[s] : n[%d] ...", n));
    return 0.;
  }
  //check likelihood array  
  if (!fLike || fLike->GetEntries() != AliTRDgeometry::kNlayer*kNcharge*kNcalib) {
    AliError("No usable AttachClusters calib object.");
    return 0.;
  }
  Int_t offset(ly*kNcalib*kNcharge); // offset for layer
  TGraphErrors *g(NULL);
  if(!(g = (TGraphErrors*)fLike->At(offset+Int_t(kResPos)*kNcharge+Int_t(chg)))){
    AliError("Failed retrieving AttachClusters graph.");
    return 0.;
  }
  // Interpolate p_t
  Int_t npts(g->GetN()), ip(0), jp(-1);
  Double_t x0, y0, x1, y1, dd(0.), invdx(0.), f[4]={0., 0., 0., 0.};
  for(Int_t kp(0); kp<npts; kp++){
    g->GetPoint(kp, x1, y1);
    if(x1>=pt){jp=kp; break;}
  }
  Bool_t boundary(kFALSE);
  if(jp<0){
    jp = npts-1; 
    g->GetPoint(jp, x1, y1);
    ip = npts-1;
    boundary = kTRUE;
  }else if(jp==0){ 
    ip = 0;
    boundary = kTRUE;
  }else{ 
    ip = jp-1;
    g->GetPoint(ip, x0, y0);
    invdx = 1./(x0-x1);
  }
  // process pt dependences
  // +++ process dy
  Double_t dym = boundary?y1:((pt*(y0-y1) + (x0*y1-x1*y0))*invdx),
           sym = 0.5*(g->GetErrorY(ip)+g->GetErrorY(jp));
  dd      = (dyr - dym)/sym; dd*=dd;
  f[0] = TMath::Exp(-0.5*dd);
  // +++ process dphi
  if(!(g = (TGraphErrors*)fLike->At(offset+Int_t(kResAng)*kNcharge+Int_t(chg)))){
    AliError("Failed retrieving AttachClusters graph.");
    return 0.;
  }
  g->GetPoint(ip, x0, y0);g->GetPoint(jp, x1, y1);
  Double_t dpm = boundary?y1:((pt*(y0-y1) + (x0*y1-x1*y0))*invdx),
           spm = 0.5*(g->GetErrorY(ip)+g->GetErrorY(jp));
  dd      = (dphi - dpm)/spm; dd*=dd;
  f[1] = TMath::Exp(-0.5*dd);
  // +++ process no of clusters
  if(!(g = (TGraphErrors*)fLike->At(offset+Int_t(kNclMean)*kNcharge+Int_t(chg)))){
    AliError("Failed retrieving AttachClusters graph.");
    return 0.;
  }
  g->GetPoint(ip, x0, y0);g->GetPoint(jp, x1, y1);
  Double_t nm = boundary?y1:((pt*(y0-y1) + (x0*y1-x1*y0))*invdx);
  f[2] = (nm-TMath::Abs(n-nm))/nm;
 
  // process phi dependences
  // +++ process <s>/s
  if(!(g = (TGraphErrors*)fLike->At(offset+Int_t(kSigma)*kNcharge+Int_t(chg)))){
    AliError("Failed retrieving AttachClusters graph.");
    return 0.;
  }
  // Interpolate phi [deg]
  npts=g->GetN(); jp=-1;
  for(Int_t kp(0); kp<npts; kp++){
    g->GetPoint(kp, x1, y1);
    if(x1>=phiTrk){jp=kp; break;}
  }
  if(jp<0){
    jp = npts-1; 
    g->GetPoint(jp, x1, y1);
    ip = jp;
    boundary = kTRUE;
  }else if(jp==0){ 
    ip = jp;
    boundary = kTRUE;
  }else{ 
    ip = jp-1;
    g->GetPoint(ip, x0, y0);
    invdx = 1./(x0-x1);
    boundary = kFALSE;
  }
  Double_t sm = boundary?y1:((phiTrk*(y0-y1) + (x0*y1-x1*y0))*invdx),
           ssm = 0.5*(g->GetErrorY(ip)+g->GetErrorY(jp));
  dd      = (sr - sm)/ssm; dd*=dd;
  f[3] = 1.;//TMath::Exp(-0.5*dd);

  // Calculate likelihood
  Double_t length = f[0]*f[0]+f[1]*f[1]+f[2]*f[2]+f[3]*f[3];
  length = TMath::Sqrt(length);
  Double_t cosTht = f[0]+f[1]+f[2]+f[3];
  cosTht /= (4.*length);
  AliDebug(2, Form("Like[%5.3f] ThtLike[%6.2f](deg)\n"
    "    F_dy (%+6.2f)=%4.2f\n" 
    "    F_phi(%+6.2f)=%4.2f\n"
    "    F_ncl(%+6d)=%4.2f\n"
    "    F_<s>(%+6.2f)=%4.2f", 
    length, TMath::ACos(cosTht)*TMath::RadToDeg(),
    dyr, f[0], dphi, f[1], n, f[2], sr, f[3]));

  return length;
}

//______________________________________________________________________________
void AliTRDCalTrkAttach::Help()
{
// Display help message
  printf(
    "Draw likelihood distribution. Possible options are of the form \"lt\"\n" 
    "where \"l\" is the layer number [0-5] and\n"
    "      \"t\" is the type. This is one of the following\n"
    "            \"y\" - r-phi roads\n"
    "            \"p\" - angular (deflection) roads\n"
    "            \"n\" - number of clusters\n"
    "            \"s\" - cluster spread\n");
}

//______________________________________________________________________________
void AliTRDCalTrkAttach::Draw(Option_t* opt)
{
// Draw likelihood distribution. Possible options are of the form "lt" 
// where "l" is the layer number [0-5] and
//       "t" is the type. This is one of the following
//             "y" - r-phi roads
//             "p" - angular (deflection) roads
//             "n" - number of clusters
//             "s" - cluster spread
  if (!fLike || fLike->GetEntries() != AliTRDgeometry::kNlayer*kNcharge*kNcalib) {
    AliError("No likelihood distributions stored");
    return;
  }
  if(!opt){
    Help();
    return;
  }
  Int_t ly(-1); Char_t cly[] = {'0','1','2','3','4','5'}; 
  for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
    if(opt[0] == cly[ily]){
      ly = ily; break;
    }
  }
  if(ly<0){
    Help();
    return;
  }
  Int_t typ(-1); const Char_t ctyp[] = {'y', 'p', 's', 'n'};
  for(Int_t it(0); it<kNcalib; it++){
    if(opt[1] == ctyp[it]){
      typ = it; break;
    }
  }
  if(typ<0){
    Help();
    return;
  }
  Int_t offset(ly*kNcalib*kNcharge+typ*kNcharge); // offset for layer and type
  TGraphErrors *g[2]={(TGraphErrors*)fLike->At(offset), (TGraphErrors*)fLike->At(offset+1)};
  if(!g[0] || !g[1]){
    AliError(Form("Failed retrieving graphs for Ly[%d] Typ[%c]", ly, ctyp[typ]));
    return;
  }

  // draw !!
  const Char_t *ttyp[] = {"#Deltay [cm]", "#Delta#phi [deg]", "#Delta#sigma/<#sigma>", "n_{cl}^{tracklet}"};
  g[0]->Draw("apl"); g[0]->SetLineColor(kBlue); g[0]->SetMarkerColor(kBlue);
  g[1]->Draw("pl");  g[1]->SetLineColor(kRed); g[1]->SetMarkerColor(kRed);
  TH1 *h=g[0]->GetHistogram();
  h->GetYaxis()->SetTitle(ttyp[typ]);
  h->GetYaxis()->SetLabelFont(52);
  h->GetYaxis()->SetTitleFont(62);
  h->GetXaxis()->SetTitle("p_{t} [GeV/c]");
  h->GetXaxis()->SetLabelFont(52);
  h->GetXaxis()->SetTitleFont(62);
}

//______________________________________________________________________________
Bool_t AliTRDCalTrkAttach::LoadReferences(const Char_t *file)
{
// Load calibration data from file

  if(!file || !TFile::Open(file)){
    AliError("Parametrization file missing or unreadable.");
    return kFALSE;
  }
  TGraphErrors *g(NULL);
  Char_t co[kNcalib] = {'y', 'p', 's', 'n'},
         cs[2] = {'n', 'p'};

  if(fLike) fLike->Clear();
  else fLike = new TObjArray(AliTRDgeometry::kNlayer*kNcharge*kNcalib);
  for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
    for(Int_t icalib(0); icalib<kNcalib; icalib++){
      for(Int_t isgn(0); isgn<kNcharge; isgn++){
        if(!(g = (TGraphErrors*)gFile->Get(Form("%c%c%d", co[icalib], cs[isgn], ily)))) return kFALSE;
        fLike->AddAt(g, kNcharge*(ily*kNcalib+icalib)+isgn);
      }
    }
  }
  return kTRUE;
}
Commit	Line	Data
4b44a8e1	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	//////////////////////////////////////////////////////////////////////////
	17	// //
	18	// Container for calibration parameters for AliTRDseedV1::AttachClusters()
	19	// .... Longer description of content ...
	20	//
	21	// For calibration procedure check AliTRDtrackleOflHelper::CalibrateAttach()
	22	// .... some reference to the calibration procedure .........
	23	// //
	24	// Authors: //
	25	// Alex Bercuci <a.bercuci@gsi.de> //
	26	// //
	27	//////////////////////////////////////////////////////////////////////////
	28
	29	#include <TFile.h>
	30	#include <TROOT.h>
	31	#include <TMath.h>
3da83347	32	#include "TH1.h"
4b44a8e1	33	#include "TGraphErrors.h"
	34	#include "TObjArray.h"
	35
	36	#include "AliLog.h"
	37
4b44a8e1	38	#include "AliTRDcalibDB.h"
803dc399	39	#include "AliTRDgeometry.h"
803dc399	40	#include "AliTRDCalTrkAttach.h"
4b44a8e1	41
	42	ClassImp(AliTRDCalTrkAttach)
	43
	44	//______________________________________________________________________________
	45	AliTRDCalTrkAttach::AliTRDCalTrkAttach()
	46	:TNamed("AliTRDCalTrkAttach", "Calibration of AliTRDseedV1::AttachClusters")
	47	,fRClikeLimit(0.65)
	48	,fScaleCov(2.)
	49	,fLike(NULL)
	50	{
	51	// Default constructor
	52
	53	fNsgmDy[0] = 5; fNsgmDy[1] = 7;
	54	fLikeMinRelDecrease[0] = .2; fLikeMinRelDecrease[1] = .3;
	55	}
	56
	57	//______________________________________________________________________________
	58	AliTRDCalTrkAttach::~AliTRDCalTrkAttach()
	59	{
	60	// Destructor
	61	if(fLike) delete fLike;
	62	}
	63
	64
	65	//______________________________________________________________________________
	66	Double_t AliTRDCalTrkAttach::CookLikelihood(Bool_t chg, Int_t ly, Float_t pt, Float_t phiTrk, Int_t n, Double_t dyr, Double_t dphi, Double_t sr) const
	67	{
	68	// Calculate likelihood for a segment to belong to a tracklet
	69	// Based on calibrated values
	70
	71	if(n<4){
	72	AliDebug(4, Form("Failed basic cut[s] : n[%d] ...", n));
	73	return 0.;
	74	}
	75	//check likelihood array
803dc399	76	if (!fLike \|\| fLike->GetEntries() != AliTRDgeometry::kNlayerkNchargekNcalib) {
4b44a8e1	77	AliError("No usable AttachClusters calib object.");
	78	return 0.;
	79	}
803dc399	80	Int_t offset(lykNcalibkNcharge); // offset for layer
4b44a8e1	81	TGraphErrors *g(NULL);
803dc399	82	if(!(g = (TGraphErrors)fLike->At(offset+Int_t(kResPos)kNcharge+Int_t(chg)))){
4b44a8e1	83	AliError("Failed retrieving AttachClusters graph.");
	84	return 0.;
	85	}
	86	// Interpolate p_t
	87	Int_t npts(g->GetN()), ip(0), jp(-1);
	88	Double_t x0, y0, x1, y1, dd(0.), invdx(0.), f[4]={0., 0., 0., 0.};
	89	for(Int_t kp(0); kp<npts; kp++){
	90	g->GetPoint(kp, x1, y1);
	91	if(x1>=pt){jp=kp; break;}
	92	}
	93	Bool_t boundary(kFALSE);
	94	if(jp<0){
	95	jp = npts-1;
	96	g->GetPoint(jp, x1, y1);
	97	ip = npts-1;
	98	boundary = kTRUE;
	99	}else if(jp==0){
	100	ip = 0;
	101	boundary = kTRUE;
	102	}else{
	103	ip = jp-1;
	104	g->GetPoint(ip, x0, y0);
	105	invdx = 1./(x0-x1);
	106	}
	107	// process pt dependences
	108	// +++ process dy
	109	Double_t dym = boundary?y1:((pt(y0-y1) + (x0y1-x1y0))invdx),
	110	sym = 0.5*(g->GetErrorY(ip)+g->GetErrorY(jp));
	111	dd = (dyr - dym)/sym; dd*=dd;
	112	f[0] = TMath::Exp(-0.5*dd);
	113	// +++ process dphi
803dc399	114	if(!(g = (TGraphErrors)fLike->At(offset+Int_t(kResAng)kNcharge+Int_t(chg)))){
4b44a8e1	115	AliError("Failed retrieving AttachClusters graph.");
	116	return 0.;
	117	}
	118	g->GetPoint(ip, x0, y0);g->GetPoint(jp, x1, y1);
	119	Double_t dpm = boundary?y1:((pt(y0-y1) + (x0y1-x1y0))invdx),
	120	spm = 0.5*(g->GetErrorY(ip)+g->GetErrorY(jp));
	121	dd = (dphi - dpm)/spm; dd*=dd;
	122	f[1] = TMath::Exp(-0.5*dd);
	123	// +++ process no of clusters
803dc399	124	if(!(g = (TGraphErrors)fLike->At(offset+Int_t(kNclMean)kNcharge+Int_t(chg)))){
4b44a8e1	125	AliError("Failed retrieving AttachClusters graph.");
	126	return 0.;
	127	}
	128	g->GetPoint(ip, x0, y0);g->GetPoint(jp, x1, y1);
	129	Double_t nm = boundary?y1:((pt(y0-y1) + (x0y1-x1y0))invdx);
	130	f[2] = (nm-TMath::Abs(n-nm))/nm;
	131
	132	// process phi dependences
	133	// +++ process <s>/s
803dc399	134	if(!(g = (TGraphErrors)fLike->At(offset+Int_t(kSigma)kNcharge+Int_t(chg)))){
4b44a8e1	135	AliError("Failed retrieving AttachClusters graph.");
	136	return 0.;
	137	}
	138	// Interpolate phi [deg]
	139	npts=g->GetN(); jp=-1;
	140	for(Int_t kp(0); kp<npts; kp++){
	141	g->GetPoint(kp, x1, y1);
	142	if(x1>=phiTrk){jp=kp; break;}
	143	}
	144	if(jp<0){
	145	jp = npts-1;
	146	g->GetPoint(jp, x1, y1);
	147	ip = jp;
	148	boundary = kTRUE;
	149	}else if(jp==0){
	150	ip = jp;
	151	boundary = kTRUE;
	152	}else{
	153	ip = jp-1;
	154	g->GetPoint(ip, x0, y0);
	155	invdx = 1./(x0-x1);
	156	boundary = kFALSE;
	157	}
	158	Double_t sm = boundary?y1:((phiTrk(y0-y1) + (x0y1-x1y0))invdx),
	159	ssm = 0.5*(g->GetErrorY(ip)+g->GetErrorY(jp));
	160	dd = (sr - sm)/ssm; dd*=dd;
803dc399	161	f[3] = 1.;//TMath::Exp(-0.5*dd);
4b44a8e1	162
	163	// Calculate likelihood
	164	Double_t length = f[0]f[0]+f[1]f[1]+f[2]f[2]+f[3]f[3];
	165	length = TMath::Sqrt(length);
	166	Double_t cosTht = f[0]+f[1]+f[2]+f[3];
	167	cosTht /= (4.*length);
	168	AliDebug(2, Form("Like[%5.3f] ThtLike[%6.2f](deg)\n"
	169	" F_dy (%+6.2f)=%4.2f\n"
	170	" F_phi(%+6.2f)=%4.2f\n"
	171	" F_ncl(%+6d)=%4.2f\n"
	172	" F_<s>(%+6.2f)=%4.2f",
	173	length, TMath::ACos(cosTht)*TMath::RadToDeg(),
	174	dyr, f[0], dphi, f[1], n, f[2], sr, f[3]));
	175
	176	return length;
	177	}
	178
3da83347	179	//______________________________________________________________________________
	180	void AliTRDCalTrkAttach::Help()
	181	{
	182	// Display help message
	183	printf(
	184	"Draw likelihood distribution. Possible options are of the form \"lt\"\n"
	185	"where \"l\" is the layer number [0-5] and\n"
	186	" \"t\" is the type. This is one of the following\n"
	187	" \"y\" - r-phi roads\n"
	188	" \"p\" - angular (deflection) roads\n"
	189	" \"n\" - number of clusters\n"
	190	" \"s\" - cluster spread\n");
	191	}
	192
	193	//______________________________________________________________________________
	194	void AliTRDCalTrkAttach::Draw(Option_t* opt)
	195	{
	196	// Draw likelihood distribution. Possible options are of the form "lt"
	197	// where "l" is the layer number [0-5] and
	198	// "t" is the type. This is one of the following
	199	// "y" - r-phi roads
	200	// "p" - angular (deflection) roads
	201	// "n" - number of clusters
	202	// "s" - cluster spread
	203	if (!fLike \|\| fLike->GetEntries() != AliTRDgeometry::kNlayerkNchargekNcalib) {
	204	AliError("No likelihood distributions stored");
	205	return;
	206	}
	207	if(!opt){
	208	Help();
	209	return;
	210	}
	211	Int_t ly(-1); Char_t cly[] = {'0','1','2','3','4','5'};
	212	for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
	213	if(opt[0] == cly[ily]){
	214	ly = ily; break;
	215	}
	216	}
	217	if(ly<0){
	218	Help();
	219	return;
	220	}
	221	Int_t typ(-1); const Char_t ctyp[] = {'y', 'p', 's', 'n'};
	222	for(Int_t it(0); it<kNcalib; it++){
	223	if(opt[1] == ctyp[it]){
	224	typ = it; break;
	225	}
	226	}
	227	if(typ<0){
	228	Help();
	229	return;
	230	}
	231	Int_t offset(lykNcalibkNcharge+typ*kNcharge); // offset for layer and type
	232	TGraphErrors g[2]={(TGraphErrors)fLike->At(offset), (TGraphErrors*)fLike->At(offset+1)};
	233	if(!g[0] \|\| !g[1]){
	234	AliError(Form("Failed retrieving graphs for Ly[%d] Typ[%c]", ly, ctyp[typ]));
	235	return;
	236	}
	237
	238	// draw !!
	239	const Char_t *ttyp[] = {"#Deltay [cm]", "#Delta#phi [deg]", "#Delta#sigma/<#sigma>", "n_{cl}^{tracklet}"};
	240	g[0]->Draw("apl"); g[0]->SetLineColor(kBlue); g[0]->SetMarkerColor(kBlue);
	241	g[1]->Draw("pl"); g[1]->SetLineColor(kRed); g[1]->SetMarkerColor(kRed);
	242	TH1 *h=g[0]->GetHistogram();
243	h->GetYaxis()->SetTitle(ttyp[typ]);
244	h->GetYaxis()->SetLabelFont(52);
245	h->GetYaxis()->SetTitleFont(62);
246	h->GetXaxis()->SetTitle("p_{t} [GeV/c]");
247	h->GetXaxis()->SetLabelFont(52);
248	h->GetXaxis()->SetTitleFont(62);
249	}
250
4b44a8e1	251	//______________________________________________________________________________
	252	Bool_t AliTRDCalTrkAttach::LoadReferences(const Char_t *file)
	253	{
	254	// Load calibration data from file
	255
	256	if(!file \|\| !TFile::Open(file)){
	257	AliError("Parametrization file missing or unreadable.");
	258	return kFALSE;
	259	}
	260	TGraphErrors *g(NULL);
	261	Char_t co[kNcalib] = {'y', 'p', 's', 'n'},
	262	cs[2] = {'n', 'p'};
	263
	264	if(fLike) fLike->Clear();
803dc399	265	else fLike = new TObjArray(AliTRDgeometry::kNlayerkNchargekNcalib);
803dc399	266	for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
4b44a8e1	267	for(Int_t icalib(0); icalib<kNcalib; icalib++){
803dc399	268	for(Int_t isgn(0); isgn<kNcharge; isgn++){
4b44a8e1	269	if(!(g = (TGraphErrors*)gFile->Get(Form("%c%c%d", co[icalib], cs[isgn], ily)))) return kFALSE;
803dc399	270	fLike->AddAt(g, kNcharge(ilykNcalib+icalib)+isgn);
4b44a8e1	271	}
	272	}
	273	}
	274	return kTRUE;
	275	}