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

/* $Id:$ */

#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(const AliExternalComparison& comp)
  :TNamed(comp),
   fResolHistos(new TObjArray(*(comp.fResolHistos))),
   fPullHistos(new TObjArray(*(comp.fPullHistos))),
   fRangeMatrix(new TMatrixD(*(comp.fRangeMatrix))),
   fCutMatrix(new TMatrixD(*(comp.fCutMatrix)))
{
  //
  // copy constructor
  //
}

AliExternalComparison& AliExternalComparison::operator=(const AliExternalComparison& comp)
{
  //
  //
  //
  if(this != &comp) {
    TNamed::operator=(comp);
    fResolHistos=new TObjArray(*(comp.fResolHistos));
    fPullHistos=new TObjArray(*(comp.fPullHistos));
    fRangeMatrix=new TMatrixD(*(comp.fRangeMatrix));
    fCutMatrix=new TMatrixD(*(comp.fCutMatrix));
  }
  return *this;
}


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

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


Long64_t AliExternalComparison::Merge(TCollection *const 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 *const 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 
  //    
  if (!AcceptPair(param0,param1)) return;
  //
  if (!fResolHistos) MakeHistos(); 
  //
  //
  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::SetParameterRange(Int_t param, Float_t min, Float_t max, Int_t nbins){
  //
  //
  // 
  if (!fRangeMatrix) SetDefaultRange();
  TMatrixD & mat = *fRangeMatrix; 
  if (param<0) return;
  if (param>6) return;
  mat(param,0)=min;      mat(param,1)=max;   mat(param,2)=nbins; 
}

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