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