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
  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   //
 131   //if (TMath::Abs(param0->GetX()-param1->GetX())>kDeltaX) return;
 132   //if (TMath::Abs(param0->GetAlpha()-param1->GetAlpha())>kDeltaA) return;
 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
 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
 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;
 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
 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 }