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 /*
  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 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 }
  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
 100   delete fCutMatrix;            // range matrix
 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   //
 151   if (!AcceptPair(param0,param1)) return;
 152   //
 153   if (!fResolHistos) MakeHistos();
 154   //
 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
 297   mat(6,0)= -TMath::Pi();   mat(6,1)=TMath::Pi();    mat(6,2)=18*4;   // fi range
 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
 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
 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
 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;
 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
 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 }