[u/mrichter/AliRoot.git] / TPC / AliExternalComparison.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.                  *
 **************************************************************************/

/*


*/


#include "AliExternalComparison.h"
#include "AliExternalTrackParam.h"
#include "AliTrackReference.h"
#include "THnSparse.h"
#include "TParticle.h"
#include "TMatrix.h"
#include "TParticlePDG.h"
#include "TParticle.h"
#include "TTreeStream.h"
#include "TAxis.h"

ClassImp(AliExternalComparison)


AliExternalComparison::AliExternalComparison() 
  :TNamed(),
   fResolHistos(0),
   fPullHistos(0),
   fRangeMatrix(0),
   fCutMatrix(0)
{  
  //
  // Default constructor
  //
}


AliExternalComparison::AliExternalComparison(const Text_t *name, const Text_t *title)
  :TNamed(name,title),
   fResolHistos(0),
   fPullHistos(0),
   fRangeMatrix(0),
   fCutMatrix(0)
{
  //
  // Non default cosntructor
  //
}


AliExternalComparison::~AliExternalComparison(){
  //
  //
  //
  if (fResolHistos) fResolHistos->Delete();
  if (fPullHistos)  fPullHistos->Delete();
  delete fResolHistos;             // resolution histogram
  delete fPullHistos;              // pull       histogram
  delete fRangeMatrix;            // range matrix
}
 

void AliExternalComparison::Analyze() {
  //
  //
  //
}


Long64_t AliExternalComparison::Merge(TCollection *li) {
  //
  //
  //   
  TIterator* iter = li->MakeIterator();
  AliExternalComparison* comp = 0;
  TString strName(GetName());
  while ((comp = (AliExternalComparison*)iter->Next())) {
    if (!comp->InheritsFrom(AliExternalComparison::Class())) {
      return -1;
    }
    if (strName.CompareTo(comp->GetName())!=0) return -1;
    // add histograms here...
    Add(comp);
  }
  return 0;  
}

void  AliExternalComparison:: Add(AliExternalComparison*comp){
  //
  // Add comparison object
  //
  if (!fResolHistos) return;
  for (Int_t i=0;i<5;i++){
    THnSparse * h0 = (THnSparse*)fResolHistos->At(i);
    THnSparse * h1 = (THnSparse*)comp->fResolHistos->At(i);
    if (h0&&h1) h0->Add(h1);
    h0 = (THnSparse*)fPullHistos->At(i);
    h1 = (THnSparse*)comp->fPullHistos->At(i);
    if (h0&&h1) h0->Add(h1);
  }
}


void   AliExternalComparison::Process(const AliExternalTrackParam *param0, const AliExternalTrackParam *param1){
  //
  // Process- fill histogram with residuals
  //    Tracks has to be in the same local X 
  //    
  const Double_t kDeltaX=0.2;
  const Double_t kDeltaA=0.01;
  if (!AcceptPair(param0,param1)) return;
  //
  if (!fResolHistos) MakeHistos(); 
  //
  //if (TMath::Abs(param0->GetX()-param1->GetX())>kDeltaX) return;
  //if (TMath::Abs(param0->GetAlpha()-param1->GetAlpha())>kDeltaA) return;
  //
  const Double_t *p0 = param0->GetParameter();
  const Double_t *p1 = param1->GetParameter();
  const Double_t *c0 = param0->GetCovariance();
  const Double_t *c1 = param1->GetCovariance();
  Double_t xyz[3];
  param1->GetXYZ(xyz);
  //
  Double_t vec[10];
  vec[0] = p1[0];
  vec[1] = p1[1];
  vec[2] = p1[2];
  vec[3] = p1[3];
  vec[4] = TMath::Sqrt(TMath::Abs(p1[4]))* ((p1[4]>0)? 1:-1);
  vec[5] = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
  vec[6] = TMath::ATan2(xyz[1],xyz[0]);
  //
  //
  THnSparse * his = 0;
  //
  vec[7] = p0[0]-p1[0];
  his = (THnSparse*)fResolHistos->At(0);
  his->Fill(vec);
  vec[7]=(p0[0]-p1[0])/TMath::Sqrt(c0[0]+c1[0]);
  his = (THnSparse*)fPullHistos->At(0);
  his->Fill(vec);
  //
  //
  vec[7] = p0[1]-p1[1];
  his = (THnSparse*)fResolHistos->At(1);
  his->Fill(vec);
  vec[7]=(p0[1]-p1[1])/TMath::Sqrt(c0[2]+c1[2]);
  his = (THnSparse*)fPullHistos->At(1);
  his->Fill(vec);
  //
  vec[7] = p0[2]-p1[2];
  his = (THnSparse*)fResolHistos->At(2);
  his->Fill(vec);
  vec[7]=(p0[2]-p1[2])/TMath::Sqrt(c0[5]+c1[5]);
  his = (THnSparse*)fPullHistos->At(2);
  his->Fill(vec);
  //
  vec[7] = p0[3]-p1[3];
  his = (THnSparse*)fResolHistos->At(3);
  his->Fill(vec);
  vec[7]=(p0[3]-p1[3])/TMath::Sqrt(c0[9]+c1[9]);
  his = (THnSparse*)fPullHistos->At(3);
  his->Fill(vec);
  //
  vec[7] = p0[4]-p1[4];
  his = (THnSparse*)fResolHistos->At(4);
  his->Fill(vec);
  vec[7]=(p0[4]-p1[4])/TMath::Sqrt(c0[14]+c1[14]);
  his = (THnSparse*)fPullHistos->At(4);
  his->Fill(vec);
}

void  AliExternalComparison::Process(const AliExternalTrackParam *param0, TParticle *part){
  //
  //
  //
  AliExternalTrackParam *param1 = MakeExternalParam(part);
  param1->Rotate(param0->GetAlpha());
  if (param0&&param1) Process(param0, param1);
}
void AliExternalComparison::Process(const AliExternalTrackParam *param0, TParticle *part, const AliTrackReference *ref){
  //
  //
  //
  AliExternalTrackParam *param1 = MakeExternalParam(part,ref);
  param1->Rotate(param0->GetAlpha());
  if (param0&&param1) Process(param0, param1);
}

THnSparse * AliExternalComparison::GetHisto(Int_t ivar, Int_t type) {
  //
  //
  //
  if (!fResolHistos) return 0;
  if (ivar<0 || ivar>4)  return 0;
  if (type==0) return (THnSparse*)fResolHistos->At(ivar);
  if (type==1) return (THnSparse*)fPullHistos->At(ivar);
  return 0;
}

void AliExternalComparison::MakeHistos(){
  //
  //
  // 
  if (!fRangeMatrix) SetDefaultRange();
  TMatrixD & mat = *fRangeMatrix;  
  //
  Double_t xmin[8], xmax[8];
  Double_t xpmin[8], xpmax[8];
  Int_t nbins[10];
  for (Int_t i=0;i<7;i++){
    xmin[i]=mat(i,0);
    xmax[i]=mat(i,1);
    xpmin[i]=mat(i,0);
    xpmax[i]=mat(i,1);
    nbins[i]=TMath::Nint(mat(i,2));
  }
  //
  fResolHistos = new TObjArray(5);
  fPullHistos = new TObjArray(5);
  //
  for (Int_t idelta=0;idelta<5;idelta++){
    xmin[7]=mat(idelta+7,0); 
    xmax[7]=mat(idelta+7,1); 
    nbins[7]=TMath::Nint(mat(idelta+7,2));
    fResolHistos->AddAt(new THnSparseF(Form("Resol%s",GetName()),Form("Resol%s",GetTitle()),
				     8,nbins,xmin,xmax),idelta);
    xmin[7]=-10;
    xmax[7]=10;
    fPullHistos->AddAt(new THnSparseF(Form("Pull%s",GetName()),Form("Pull%s",GetTitle()),
				      8,nbins,xmin,xmax),idelta);
  }
  //
}


void  AliExternalComparison::SetDefaultRange(Float_t scale,Float_t arm, Int_t nbins){
  //
  // set default binning
  // scale - define position range
  // arm  - Level arm
  // Binning:
  // 7 variables
  // 5 deltas

  if (!fRangeMatrix) fRangeMatrix = new TMatrixD(12,3);
  TMatrixD & mat = *fRangeMatrix;  
  //
  // variables
  //
  mat(0,0)=-50;    mat(0,1)=50;   mat(0,2)=50;  // P0 -y            range
  mat(1,0)=-250;   mat(1,1)=250;  mat(1,2)=50;  // P1 -Z            range
  mat(2,0)=-0.99;  mat(2,1)=0.99; mat(2,2)=20;  // P2 -snp          range
  mat(3,0)=-1.5;   mat(3,1)=1.5;  mat(3,2)=20;  // P3 -tantheta     range
  mat(4,0)=-3;     mat(4,1)=3;    mat(4,2)=50;   // sqrt(P4) -sqrt(1/pt) range  
  //
  mat(5,0)= 0;     mat(5,1)=250;  mat(5,2)=50;   // R  range  
  mat(6,0)= TMath::Pi();   mat(6,1)=TMath::Pi();    mat(6,2)=18*4;   // fi range  
  //
  // Resolution
  //
  mat(7,0)=-scale;      mat(7,1)=scale;   mat(7,2)=nbins;       // P0 -y            range
  mat(8,0)=-scale;      mat(8,1)=scale;   mat(8,2)=nbins;       // P1 -Z            range
  mat(9,0)=-scale/arm;  mat(9,1)=scale/arm; mat(9,2)=nbins;     // P2 -snp          range
  mat(10,0)=-scale/arm; mat(10,1)=scale/arm;  mat(10,2)=nbins;  // P3 -tantheta     range
  mat(11,0)=-1;         mat(11,1)=1;   mat(11,2)=nbins;         // sqrt(P4) -sqrt(1/pt) range
  //
}

void  AliExternalComparison::SetDefaultCuts(){
  //

  if (!fCutMatrix) fCutMatrix = new TMatrixD(10,1);
  TMatrixD & mat = *fCutMatrix;  
  //
  mat(0,0)=10;     //dP0
  mat(1,0)=20;     //dP1
  mat(2,0)=0.05;   //dP2
  mat(3,0)=0.05;   //dP3
  mat(4,0)=1;      //dP4
  //
  mat(5,0)=10;      //dnP0
  mat(6,0)=50;      //dnP1
  mat(7,0)=10;      //dnP2
  mat(8,0)=10;      //dnP3
  mat(9,0)=10;     //dnP4
  //
}

void AliExternalComparison::SetResolRange(Int_t param, Float_t min, Float_t max, Int_t nbins){
  //
  //
  //
  if (!fRangeMatrix) SetDefaultRange();
  TMatrixD & mat = *fRangeMatrix; 
  if (param<0) return;
  if (param>4) return;
  mat(7+param,0)=min;      mat(7+param,1)=max;   mat(7+param,2)=nbins;
}

void AliExternalComparison::SetDistCut(Float_t dP0, Float_t dP1,Float_t dP2,Float_t dP3, Float_t dP4){
  //
  // Set diff cuts
  //
  if (!fCutMatrix) SetDefaultCuts();
  TMatrixD & mat = *fCutMatrix;  
  mat(0,0)=dP0;     //dP0
  mat(1,0)=dP1;     //dP1
  mat(2,0)=dP2;   //dP2
  mat(3,0)=dP3;   //dP3
  mat(4,0)=dP4;      //dP4

  
}
void AliExternalComparison::SetPullDistCut(Float_t dnP0, Float_t dnP1,Float_t dnP2,Float_t dnP3, Float_t dnP4){
  //
  //
  //
  if (!fCutMatrix) SetDefaultCuts();
  TMatrixD & mat = *fCutMatrix;  
  mat(5,0)=dnP0;     //dnP0
  mat(6,0)=dnP1;     //dnP1
  mat(7,0)=dnP2;     //dnP2
  mat(8,0)=dnP3;     //dnP3
  mat(9,0)=dnP4;     //dnP4
}

Bool_t   AliExternalComparison::AcceptPair(const AliExternalTrackParam *param0, const AliExternalTrackParam *param1){  
  //
  //
  //
  Bool_t isOK =kTRUE;
  if (!fCutMatrix) SetDefaultCuts();
  TMatrixD & mat = *fCutMatrix;  
  //
  const Double_t *p0 = param0->GetParameter();
  const Double_t *p1 = param1->GetParameter();
  const Double_t *c0 = param0->GetCovariance();
  const Double_t *c1 = param1->GetCovariance();
  //
  if (TMath::Abs(p0[0]-p1[0])>mat(0,0)) return kFALSE;
  if (TMath::Abs(p0[1]-p1[1])>mat(1,0)) return kFALSE;
  if (TMath::Abs(p0[2]-p1[2])>mat(2,0)) return kFALSE;
  if (TMath::Abs(p0[3]-p1[3])>mat(3,0)) return kFALSE;
  if (TMath::Abs(p0[4]-p1[4])>mat(4,0)) return kFALSE;
  if ((c0[0]+c1[0])<0) return kFALSE;
  if ((c0[2]+c1[2])<0) return kFALSE;
  if ((c0[5]+c1[5])<0) return kFALSE;
  if ((c0[9]+c1[9])<0) return kFALSE;
  if ((c0[14]+c1[14])<0) return kFALSE;
  
  if (TMath::Abs((p0[0]-p1[0])/TMath::Sqrt(c0[0]+c1[0]))>mat(5,0)) return kFALSE;
  if (TMath::Abs((p0[1]-p1[1])/TMath::Sqrt(c0[2]+c1[2]))>mat(6,0)) return kFALSE;
  if (TMath::Abs((p0[2]-p1[2])/TMath::Sqrt(c0[5]+c1[5]))>mat(7,0)) return kFALSE;
  if (TMath::Abs((p0[3]-p1[3])/TMath::Sqrt(c0[9]+c1[9]))>mat(8,0)) return kFALSE;
  if (TMath::Abs((p0[4]-p1[4])/TMath::Sqrt(c0[14]+c1[14]))>mat(9,0)) return kFALSE;
  return isOK;
}

  
AliExternalTrackParam *AliExternalComparison::MakeExternalParam(TParticle *part, const AliTrackReference *ref){
  //
  // make a Kalman track from track reference
  //
  Double_t xyz[3]={ref->X(),ref->Y(),ref->Z()};
  Double_t pxyz[3]={ref->Px(),ref->Py(),ref->Pz()};
  Double_t cv[21];
  for (Int_t i=0; i<21;i++) cv[i]=0;
  if (!part->GetPDG()) return 0;
  AliExternalTrackParam * param = new AliExternalTrackParam(xyz,pxyz,cv,TMath::Nint(part->GetPDG()->Charge()/3.));
  return param;
}


AliExternalTrackParam * AliExternalComparison::MakeExternalParam(TParticle *part){
  //
  //
  //
  Double_t xyz[3]={part->Vx(),part->Vy(),part->Vz()};
  Double_t pxyz[3]={part->Px(),part->Py(),part->Pz()};
  Double_t cv[21];
  for (Int_t i=0; i<21;i++) cv[i]=0;
  if (!part->GetPDG()) return 0;
  AliExternalTrackParam * param = new AliExternalTrackParam(xyz,pxyz,cv,TMath::Nint(part->GetPDG()->Charge()/3.));
  return param;
}


void AliExternalComparison::MakeComparisonTree(const char * outname){
  //
  // make a comparison tree
  // the tree will be written to the file - outname
  //
  /*
    TFile f("comp.root");
    AliExternalComparison *comp =  (AliExternalComparison*)f.Get("compAlign");
    outname="compTree.root";
  */
  AliExternalComparison *comp = this;
  TTreeSRedirector *pcstream = new TTreeSRedirector(outname);
  //
  THnSparse *his=0;  
  Double_t position[10];
  Double_t value; 
  Int_t *bins = new Int_t[10];
  //
  //
  for (Int_t isPull=0; isPull<2; isPull++){
    for (Int_t ires=0;ires<5; ires++){
      his = comp->GetHisto(ires,isPull);
      if (!his) continue;
      //
      Int_t ndim = his->GetNdimensions();
      //
      for (Long64_t i = 0; i < his->GetNbins(); ++i) {
	value = his->GetBinContent(i, bins);
	for (Int_t idim = 0; idim < ndim; idim++) {
	  position[idim] = his->GetAxis(idim)->GetBinCenter(bins[idim]);
	}      
	(*pcstream)<<"Resol"<<
	  "isPull="<<isPull<<     // normalized error indication
	  "index="<<ires<<        // parameter difference index
	  "bincont="<<value<<     // bin content
	  "val="<<position[7]<<   // parameter difference
	  "p0="<<position[0]<<    //P0
	  "p1="<<position[1]<<    //P1
	  "p2="<<position[2]<<    //P2
	  "p3="<<position[3]<<    //P3
	  "p4="<<position[4]<<    //P4
	  "R="<<position[5]<<     //Radius
	  "phi="<<position[6]<<     //Radius
	  "\n";
      }
    }
  }
  delete pcstream;
}
Commit	Line	Data
d73981ad	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
	19
	20
	21	*/
	22
	23
	24
	25	#include "AliExternalComparison.h"
	26	#include "AliExternalTrackParam.h"
	27	#include "AliTrackReference.h"
	28	#include "THnSparse.h"
	29	#include "TParticle.h"
	30	#include "TMatrix.h"
	31	#include "TParticlePDG.h"
	32	#include "TParticle.h"
	33	#include "TTreeStream.h"
	34	#include "TAxis.h"
	35
	36	ClassImp(AliExternalComparison)
	37
	38
	39	AliExternalComparison::AliExternalComparison()
	40	:TNamed(),
	41	fResolHistos(0),
	42	fPullHistos(0),
	43	fRangeMatrix(0),
	44	fCutMatrix(0)
	45	{
	46	//
	47	// Default constructor
	48	//
	49	}
	50
	51
	52	AliExternalComparison::AliExternalComparison(const Text_t name, const Text_t title)
	53	:TNamed(name,title),
	54	fResolHistos(0),
	55	fPullHistos(0),
	56	fRangeMatrix(0),
	57	fCutMatrix(0)
	58	{
	59	//
	60	// Non default cosntructor
	61	//
	62	}
	63
	64
65
66	AliExternalComparison::~AliExternalComparison(){
67	//
68	//
69	//
70	if (fResolHistos) fResolHistos->Delete();
71	if (fPullHistos) fPullHistos->Delete();
72	delete fResolHistos; // resolution histogram
73	delete fPullHistos; // pull histogram
74	delete fRangeMatrix; // range matrix
75	}
76
77
78	void AliExternalComparison::Analyze() {
79	//
80	//
81	//
82	}
83
84
85	Long64_t AliExternalComparison::Merge(TCollection *li) {
86	//
87	//
88	//
89	TIterator* iter = li->MakeIterator();
90	AliExternalComparison* comp = 0;
91	TString strName(GetName());
92	while ((comp = (AliExternalComparison*)iter->Next())) {
93	if (!comp->InheritsFrom(AliExternalComparison::Class())) {
94	return -1;
95	}
96	if (strName.CompareTo(comp->GetName())!=0) return -1;
97	// add histograms here...
98	Add(comp);
99	}
100	return 0;
101	}
102
103	void AliExternalComparison:: Add(AliExternalComparison*comp){
104	//
105	// Add comparison object
106	//
107	if (!fResolHistos) return;
108	for (Int_t i=0;i<5;i++){
109	THnSparse * h0 = (THnSparse*)fResolHistos->At(i);
110	THnSparse * h1 = (THnSparse*)comp->fResolHistos->At(i);
111	if (h0&&h1) h0->Add(h1);
112	h0 = (THnSparse*)fPullHistos->At(i);
113	h1 = (THnSparse*)comp->fPullHistos->At(i);
114	if (h0&&h1) h0->Add(h1);
115	}
116	}
117
118
119
120	void AliExternalComparison::Process(const AliExternalTrackParam param0, const AliExternalTrackParam param1){
121	//
122	// Process- fill histogram with residuals
123	// Tracks has to be in the same local X
124	//
125	const Double_t kDeltaX=0.2;
126	const Double_t kDeltaA=0.01;
127	if (!AcceptPair(param0,param1)) return;
128	//
129	if (!fResolHistos) MakeHistos();
130	//
8a3aca82	131	//if (TMath::Abs(param0->GetX()-param1->GetX())>kDeltaX) return;
8a3aca82	132	//if (TMath::Abs(param0->GetAlpha()-param1->GetAlpha())>kDeltaA) return;
d73981ad	133	//
	134	const Double_t *p0 = param0->GetParameter();
	135	const Double_t *p1 = param1->GetParameter();
	136	const Double_t *c0 = param0->GetCovariance();
	137	const Double_t *c1 = param1->GetCovariance();
	138	Double_t xyz[3];
	139	param1->GetXYZ(xyz);
	140	//
	141	Double_t vec[10];
	142	vec[0] = p1[0];
	143	vec[1] = p1[1];
	144	vec[2] = p1[2];
	145	vec[3] = p1[3];
	146	vec[4] = TMath::Sqrt(TMath::Abs(p1[4]))* ((p1[4]>0)? 1:-1);
	147	vec[5] = TMath::Sqrt(xyz[0]xyz[0]+xyz[1]xyz[1]);
	148	vec[6] = TMath::ATan2(xyz[1],xyz[0]);
	149	//
	150	//
	151	THnSparse * his = 0;
	152	//
	153	vec[7] = p0[0]-p1[0];
	154	his = (THnSparse*)fResolHistos->At(0);
	155	his->Fill(vec);
	156	vec[7]=(p0[0]-p1[0])/TMath::Sqrt(c0[0]+c1[0]);
	157	his = (THnSparse*)fPullHistos->At(0);
	158	his->Fill(vec);
	159	//
	160	//
	161	vec[7] = p0[1]-p1[1];
	162	his = (THnSparse*)fResolHistos->At(1);
	163	his->Fill(vec);
	164	vec[7]=(p0[1]-p1[1])/TMath::Sqrt(c0[2]+c1[2]);
	165	his = (THnSparse*)fPullHistos->At(1);
	166	his->Fill(vec);
	167	//
	168	vec[7] = p0[2]-p1[2];
	169	his = (THnSparse*)fResolHistos->At(2);
	170	his->Fill(vec);
	171	vec[7]=(p0[2]-p1[2])/TMath::Sqrt(c0[5]+c1[5]);
	172	his = (THnSparse*)fPullHistos->At(2);
	173	his->Fill(vec);
	174	//
	175	vec[7] = p0[3]-p1[3];
	176	his = (THnSparse*)fResolHistos->At(3);
	177	his->Fill(vec);
	178	vec[7]=(p0[3]-p1[3])/TMath::Sqrt(c0[9]+c1[9]);
	179	his = (THnSparse*)fPullHistos->At(3);
	180	his->Fill(vec);
	181	//
	182	vec[7] = p0[4]-p1[4];
	183	his = (THnSparse*)fResolHistos->At(4);
	184	his->Fill(vec);
	185	vec[7]=(p0[4]-p1[4])/TMath::Sqrt(c0[14]+c1[14]);
	186	his = (THnSparse*)fPullHistos->At(4);
	187	his->Fill(vec);
	188	}
	189
	190	void AliExternalComparison::Process(const AliExternalTrackParam param0, TParticle part){
	191	//
	192	//
	193	//
	194	AliExternalTrackParam *param1 = MakeExternalParam(part);
	195	param1->Rotate(param0->GetAlpha());
	196	if (param0&&param1) Process(param0, param1);
197	}
198	void AliExternalComparison::Process(const AliExternalTrackParam param0, TParticle part, const AliTrackReference *ref){
199	//
200	//
201	//
202	AliExternalTrackParam *param1 = MakeExternalParam(part,ref);
203	param1->Rotate(param0->GetAlpha());
204	if (param0&&param1) Process(param0, param1);
205	}
206
207	THnSparse * AliExternalComparison::GetHisto(Int_t ivar, Int_t type) {
208	//
209	//
210	//
211	if (!fResolHistos) return 0;
212	if (ivar<0 \|\| ivar>4) return 0;
213	if (type==0) return (THnSparse*)fResolHistos->At(ivar);
214	if (type==1) return (THnSparse*)fPullHistos->At(ivar);
215	return 0;
216	}
217
218	void AliExternalComparison::MakeHistos(){
219	//
220	//
221	//
222	if (!fRangeMatrix) SetDefaultRange();
223	TMatrixD & mat = *fRangeMatrix;
224	//
225	Double_t xmin[8], xmax[8];
226	Double_t xpmin[8], xpmax[8];
227	Int_t nbins[10];
228	for (Int_t i=0;i<7;i++){
229	xmin[i]=mat(i,0);
230	xmax[i]=mat(i,1);
231	xpmin[i]=mat(i,0);
232	xpmax[i]=mat(i,1);
233	nbins[i]=TMath::Nint(mat(i,2));
234	}
235	//
236	fResolHistos = new TObjArray(5);
237	fPullHistos = new TObjArray(5);
238	//
239	for (Int_t idelta=0;idelta<5;idelta++){
240	xmin[7]=mat(idelta+7,0);
241	xmax[7]=mat(idelta+7,1);
242	nbins[7]=TMath::Nint(mat(idelta+7,2));
243	fResolHistos->AddAt(new THnSparseF(Form("Resol%s",GetName()),Form("Resol%s",GetTitle()),
244	8,nbins,xmin,xmax),idelta);
245	xmin[7]=-10;
246	xmax[7]=10;
247	fPullHistos->AddAt(new THnSparseF(Form("Pull%s",GetName()),Form("Pull%s",GetTitle()),
248	8,nbins,xmin,xmax),idelta);
249	}
250	//
251	}
252
253
254	void AliExternalComparison::SetDefaultRange(Float_t scale,Float_t arm, Int_t nbins){
255	//
256	// set default binning
257	// scale - define position range
258	// arm - Level arm
259	// Binning:
260	// 7 variables
261	// 5 deltas
262
263	if (!fRangeMatrix) fRangeMatrix = new TMatrixD(12,3);
264	TMatrixD & mat = *fRangeMatrix;
265	//
266	// variables
267	//
268	mat(0,0)=-50; mat(0,1)=50; mat(0,2)=50; // P0 -y range
269	mat(1,0)=-250; mat(1,1)=250; mat(1,2)=50; // P1 -Z range
270	mat(2,0)=-0.99; mat(2,1)=0.99; mat(2,2)=20; // P2 -snp range
271	mat(3,0)=-1.5; mat(3,1)=1.5; mat(3,2)=20; // P3 -tantheta range
272	mat(4,0)=-3; mat(4,1)=3; mat(4,2)=50; // sqrt(P4) -sqrt(1/pt) range
273	//
274	mat(5,0)= 0; mat(5,1)=250; mat(5,2)=50; // R range
275	mat(6,0)= TMath::Pi(); mat(6,1)=TMath::Pi(); mat(6,2)=18*4; // fi range
276	//
277	// Resolution
278	//
279	mat(7,0)=-scale; mat(7,1)=scale; mat(7,2)=nbins; // P0 -y range
280	mat(8,0)=-scale; mat(8,1)=scale; mat(8,2)=nbins; // P1 -Z range
281	mat(9,0)=-scale/arm; mat(9,1)=scale/arm; mat(9,2)=nbins; // P2 -snp range
282	mat(10,0)=-scale/arm; mat(10,1)=scale/arm; mat(10,2)=nbins; // P3 -tantheta range
283	mat(11,0)=-1; mat(11,1)=1; mat(11,2)=nbins; // sqrt(P4) -sqrt(1/pt) range
284	//
285	}
286
287	void AliExternalComparison::SetDefaultCuts(){
288	//
289
290	if (!fCutMatrix) fCutMatrix = new TMatrixD(10,1);
291	TMatrixD & mat = *fCutMatrix;
292	//
293	mat(0,0)=10; //dP0
294	mat(1,0)=20; //dP1
295	mat(2,0)=0.05; //dP2
296	mat(3,0)=0.05; //dP3
297	mat(4,0)=1; //dP4
298	//
299	mat(5,0)=10; //dnP0
300	mat(6,0)=50; //dnP1
301	mat(7,0)=10; //dnP2
302	mat(8,0)=10; //dnP3
303	mat(9,0)=10; //dnP4
304	//
305	}
306
8a3aca82	307	void AliExternalComparison::SetResolRange(Int_t param, Float_t min, Float_t max, Int_t nbins){
	308	//
	309	//
	310	//
	311	if (!fRangeMatrix) SetDefaultRange();
	312	TMatrixD & mat = *fRangeMatrix;
	313	if (param<0) return;
	314	if (param>4) return;
	315	mat(7+param,0)=min; mat(7+param,1)=max; mat(7+param,2)=nbins;
	316	}
	317
d73981ad	318	void AliExternalComparison::SetDistCut(Float_t dP0, Float_t dP1,Float_t dP2,Float_t dP3, Float_t dP4){
	319	//
	320	// Set diff cuts
	321	//
	322	if (!fCutMatrix) SetDefaultCuts();
	323	TMatrixD & mat = *fCutMatrix;
	324	mat(0,0)=dP0; //dP0
	325	mat(1,0)=dP1; //dP1
	326	mat(2,0)=dP2; //dP2
	327	mat(3,0)=dP3; //dP3
	328	mat(4,0)=dP4; //dP4
	329
	330
	331	}
	332	void AliExternalComparison::SetPullDistCut(Float_t dnP0, Float_t dnP1,Float_t dnP2,Float_t dnP3, Float_t dnP4){
	333	//
	334	//
	335	//
	336	if (!fCutMatrix) SetDefaultCuts();
	337	TMatrixD & mat = *fCutMatrix;
	338	mat(5,0)=dnP0; //dnP0
	339	mat(6,0)=dnP1; //dnP1
	340	mat(7,0)=dnP2; //dnP2
	341	mat(8,0)=dnP3; //dnP3
	342	mat(9,0)=dnP4; //dnP4
	343	}
	344
	345	Bool_t AliExternalComparison::AcceptPair(const AliExternalTrackParam param0, const AliExternalTrackParam param1){
	346	//
	347	//
	348	//
	349	Bool_t isOK =kTRUE;
	350	if (!fCutMatrix) SetDefaultCuts();
	351	TMatrixD & mat = *fCutMatrix;
	352	//
	353	const Double_t *p0 = param0->GetParameter();
	354	const Double_t *p1 = param1->GetParameter();
	355	const Double_t *c0 = param0->GetCovariance();
	356	const Double_t *c1 = param1->GetCovariance();
	357	//
	358	if (TMath::Abs(p0[0]-p1[0])>mat(0,0)) return kFALSE;
	359	if (TMath::Abs(p0[1]-p1[1])>mat(1,0)) return kFALSE;
	360	if (TMath::Abs(p0[2]-p1[2])>mat(2,0)) return kFALSE;
	361	if (TMath::Abs(p0[3]-p1[3])>mat(3,0)) return kFALSE;
	362	if (TMath::Abs(p0[4]-p1[4])>mat(4,0)) return kFALSE;
8a3aca82	363	if ((c0[0]+c1[0])<0) return kFALSE;
	364	if ((c0[2]+c1[2])<0) return kFALSE;
	365	if ((c0[5]+c1[5])<0) return kFALSE;
	366	if ((c0[9]+c1[9])<0) return kFALSE;
	367	if ((c0[14]+c1[14])<0) return kFALSE;
	368
d73981ad	369	if (TMath::Abs((p0[0]-p1[0])/TMath::Sqrt(c0[0]+c1[0]))>mat(5,0)) return kFALSE;
	370	if (TMath::Abs((p0[1]-p1[1])/TMath::Sqrt(c0[2]+c1[2]))>mat(6,0)) return kFALSE;
	371	if (TMath::Abs((p0[2]-p1[2])/TMath::Sqrt(c0[5]+c1[5]))>mat(7,0)) return kFALSE;
	372	if (TMath::Abs((p0[3]-p1[3])/TMath::Sqrt(c0[9]+c1[9]))>mat(8,0)) return kFALSE;
	373	if (TMath::Abs((p0[4]-p1[4])/TMath::Sqrt(c0[14]+c1[14]))>mat(9,0)) return kFALSE;
	374	return isOK;
	375	}
	376
	377
	378	AliExternalTrackParam AliExternalComparison::MakeExternalParam(TParticle part, const AliTrackReference *ref){
	379	//
	380	// make a Kalman track from track reference
	381	//
	382	Double_t xyz[3]={ref->X(),ref->Y(),ref->Z()};
	383	Double_t pxyz[3]={ref->Px(),ref->Py(),ref->Pz()};
	384	Double_t cv[21];
	385	for (Int_t i=0; i<21;i++) cv[i]=0;
	386	if (!part->GetPDG()) return 0;
	387	AliExternalTrackParam * param = new AliExternalTrackParam(xyz,pxyz,cv,TMath::Nint(part->GetPDG()->Charge()/3.));
	388	return param;
	389	}
	390
	391
	392	AliExternalTrackParam * AliExternalComparison::MakeExternalParam(TParticle *part){
	393	//
	394	//
	395	//
	396	Double_t xyz[3]={part->Vx(),part->Vy(),part->Vz()};
	397	Double_t pxyz[3]={part->Px(),part->Py(),part->Pz()};
	398	Double_t cv[21];
	399	for (Int_t i=0; i<21;i++) cv[i]=0;
	400	if (!part->GetPDG()) return 0;
	401	AliExternalTrackParam * param = new AliExternalTrackParam(xyz,pxyz,cv,TMath::Nint(part->GetPDG()->Charge()/3.));
	402	return param;
	403	}
	404
	405
	406
	407
	408	void AliExternalComparison::MakeComparisonTree(const char * outname){
	409	//
	410	// make a comparison tree
	411	// the tree will be written to the file - outname
	412	//
	413	/*
	414	TFile f("comp.root");
	415	AliExternalComparison comp = (AliExternalComparison)f.Get("compAlign");
	416	outname="compTree.root";
	417	*/
	418	AliExternalComparison *comp = this;
	419	TTreeSRedirector *pcstream = new TTreeSRedirector(outname);
	420	//
	421	THnSparse *his=0;
	422	Double_t position[10];
	423	Double_t value;
	424	Int_t *bins = new Int_t[10];
	425	//
	426	//
	427	for (Int_t isPull=0; isPull<2; isPull++){
	428	for (Int_t ires=0;ires<5; ires++){
	429	his = comp->GetHisto(ires,isPull);
	430	if (!his) continue;
	431	//
	432	Int_t ndim = his->GetNdimensions();
433	//
434	for (Long64_t i = 0; i < his->GetNbins(); ++i) {
435	value = his->GetBinContent(i, bins);
436	for (Int_t idim = 0; idim < ndim; idim++) {
437	position[idim] = his->GetAxis(idim)->GetBinCenter(bins[idim]);
438	}
439	(*pcstream)<<"Resol"<<
440	"isPull="<<isPull<< // normalized error indication
441	"index="<<ires<< // parameter difference index
442	"bincont="<<value<< // bin content
443	"val="<<position[7]<< // parameter difference
444	"p0="<<position[0]<< //P0
445	"p1="<<position[1]<< //P1
446	"p2="<<position[2]<< //P2
447	"p3="<<position[3]<< //P3
448	"p4="<<position[4]<< //P4
449	"R="<<position[5]<< //Radius
450	"phi="<<position[6]<< //Radius
451	"\n";
452	}
453	}
454	}
455	delete pcstream;
456	}