TPC/AliExternalComparison.cxx

   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 /* $Id:$ */
  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
  57 AliExternalComparison::AliExternalComparison(const AliExternalComparison& comp)
  58   :TNamed(comp),
  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
  69 AliExternalComparison& AliExternalComparison::operator=(const AliExternalComparison& comp)
  70 {
  71   //
  72   //
  73   //
  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   }
  81   return *this;
  82 }
  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
  94   delete fCutMatrix;            // range matrix
  95 }
  96
  97
  98 void AliExternalComparison::Analyze() {
  99   //
 100   //
 101   //
 102 }
 103
 104
 105 Long64_t AliExternalComparison::Merge(TCollection *const li) {
 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
 123 void  AliExternalComparison::Add(AliExternalComparison *const comp){
 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   //
 145   if (!AcceptPair(param0,param1)) return;
 146   //
 147   if (!fResolHistos) MakeHistos();
 148   //
 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
 291   mat(6,0)= -TMath::Pi();   mat(6,1)=TMath::Pi();    mat(6,2)=18*4;   // fi range
 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
 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
 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
 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;
 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
 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 }