PWGPP/EvTrkSelection/AliAnalysisTrackingUncertainties.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, proviyaded 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 purapose. It is         *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 //////////////////////////////////////////////////////////////////////////////
  17 //                                                                          //
  18 // Analysis task for the systematic study of the uncertainties related to   //
  19 // the tracking and ITS-TPC matching efficiency for different particle      //
  20 // species.                                                                 //
  21 //                                                                          //
  22 //////////////////////////////////////////////////////////////////////////////
  23
  24
  25 #include "Riostream.h"
  26 #include "TH1F.h"
  27 #include "TH2D.h"
  28 #include "TH3F.h"
  29 #include "THn.h"
  30 #include "TList.h"
  31 #include "TMath.h"
  32 #include "TParticlePDG.h"
  33 //
  34 #include "AliAnalysisTaskSE.h"
  35 #include "AliAnalysisManager.h"
  36 #include "AliPID.h"
  37 #include "AliESDtrackCuts.h"
  38 #include "AliESDVertex.h"
  39 #include "AliESDEvent.h"
  40 #include "AliESDInputHandler.h"
  41 #include "AliESDtrack.h"
  42 #include "AliESDpid.h"
  43 #include "AliESDUtils.h"
  44 #include "AliMCEventHandler.h"
  45 #include "AliMCEvent.h"
  46 #include "AliStack.h"
  47 #include "AliLog.h"
  48 //
  49 #include "AliAnalysisTrackingUncertainties.h"
  50
  51
  52 ClassImp(AliAnalysisTrackingUncertainties)
  53
  54 // histogram constants
  55 const Int_t kNumberOfAxes = 5;
  56
  57 //________________________________________________________________________
  58 AliAnalysisTrackingUncertainties::AliAnalysisTrackingUncertainties()
  59   : AliAnalysisTaskSE("TaskTestPA"),
  60     fESD(0),
  61     fESDpid(0),
  62     fUtils(0),
  63     fMCtrue(0),
  64     fListHist(0),
  65     fESDtrackCuts(0),
  66     fMatchTr(),
  67     fMatchChi()
  68 {
  69   // default Constructor
  70   /* fast compilation test
  71
  72      gSystem->Load("libANALYSIS");
  73      gSystem->Load("libANALYSISalice");
  74      .L AliAnalysisTrackingUncertainties.cxx++
  75    */
  76 }
  77
  78
  79 //________________________________________________________________________
  80 AliAnalysisTrackingUncertainties::AliAnalysisTrackingUncertainties(const char *name)
  81   : AliAnalysisTaskSE(name),
  82     fESD(0),
  83     fESDpid(0),
  84     fUtils(0),
  85     fMCtrue(0),
  86     fListHist(0),
  87     fESDtrackCuts(0),
  88     fMatchTr(),
  89     fMatchChi()
  90 {
  91   //
  92   // standard constructur which should be used
  93   //
  94   fMCtrue = kTRUE;
  95   //
  96   // create track cuts
  97   //
  98   fESDtrackCuts = new AliESDtrackCuts("AliESDtrackCuts","AliESDtrackCuts");
  99   fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
 100   fESDtrackCuts->SetEtaRange(-1., 1.);
 101   //
 102   // analysis utils if needed
 103   //
 104   fUtils = new AliAnalysisUtils();
 105   //
 106
 107   // Output slot #0 writes into a TList container
 108   DefineOutput(1, TList::Class());
 109
 110 }
 111
 112
 113 //________________________________________________________________________
 114 void AliAnalysisTrackingUncertainties::UserCreateOutputObjects()
 115 {
 116   //
 117   // Create histograms
 118   // Called once
 119   //
 120   fListHist = new TList();
 121   fListHist->SetOwner(kTRUE);
 122   //
 123   // (1.) basic QA and statistics histograms
 124   //
 125   TH2F * histVertexSelection = new TH2F("histVertexSelection", "vertex selection; vertex z (cm); accepted/rejected", 100, -50., 50., 2, -0.5, 1.5);
 126   fListHist->Add(histVertexSelection);
 127   //
 128   // (2.) track cut variation histograms
 129   //
 130   InitializeTrackCutHistograms();
 131   //
 132   // (3.) ITS -> TPC matching histograms
 133   //
 134   Int_t    binsMatch[kNumberOfAxes] = { 10,   50,    20,            18,  6};
 135   Double_t minMatch[kNumberOfAxes]  = {  0,  0.1,    -1,             0, -0.5};
 136   Double_t maxMatch[kNumberOfAxes]  = {200,   20,    +1, 2*TMath::Pi(),  5.5};
 137   //
 138   TString axisNameMatch[kNumberOfAxes]  = {"matchChi2","pT","eta","phi","pid"};
 139   TString axisTitleMatch[kNumberOfAxes] = {"matchChi2","pT","eta","phi","pid"};
 140   //
 141   THnF * hBestMatch = new THnF("hBestMatch","ITS -> TPC matching ",kNumberOfAxes, binsMatch, minMatch, maxMatch);
 142   for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
 143     hBestMatch->GetAxis(iaxis)->SetName(axisNameMatch[iaxis]);
 144     hBestMatch->GetAxis(iaxis)->SetTitle(axisTitleMatch[iaxis]);
 145   }
 146   BinLogAxis(hBestMatch, 1);
 147   fListHist->Add(hBestMatch);
 148   //
 149   //
 150   //
 151   const int nbPt=40;
 152   const double ptMax=5;
 153   //
 154   TH2F * hNMatch    = new TH2F("hNMatch","N Matches",nbPt,0,ptMax,kMaxMatch+1,-0.5,kMaxMatch+0.5);
 155   TH2F * hAllMatch  = new TH2F("hAllMatch","All Matches Chi2",nbPt,0,ptMax,2*int(TMath::Max(1.1,kMaxChi2)),0,kMaxChi2);
 156   TH2F * hAllMatchGlo  = new TH2F("hAllMatchGlo","All Matches Chi2",nbPt,0,ptMax,2*int(TMath::Max(1.1,kMaxChi2)),0,kMaxChi2);
 157   fListHist->Add(hNMatch);
 158   fListHist->Add(hAllMatch);
 159   fListHist->Add(hAllMatchGlo);
 160   //
 161   TH2F * hNMatchBg    = new TH2F("hNMatchBg","N Matches",nbPt,0,ptMax,kMaxMatch+1,-0.5,kMaxMatch+0.5);
 162   TH2F * hBestMatchBg = new TH2F("hBestMatchBg","Best Match Chi2",nbPt,0,ptMax,2*int(TMath::Max(1.1,kMaxChi2)),0,kMaxChi2);
 163   TH2F * hAllMatchBg  = new TH2F("hAllMatchBg","All Matches Chi2",nbPt,0,ptMax,2*int(TMath::Max(1.1,kMaxChi2)),0,kMaxChi2);
 164   TH2F * hAllMatchGloBg  = new TH2F("hAllMatchGloBg","All Matches Chi2",nbPt,0,ptMax,2*int(TMath::Max(1.1,kMaxChi2)),0,kMaxChi2);
 165   fListHist->Add(hNMatchBg);
 166   fListHist->Add(hBestMatchBg);
 167   fListHist->Add(hAllMatchBg);
 168   fListHist->Add(hAllMatchGloBg);
 169   //add default track cuts in the output list
 170   fListHist->Add(fESDtrackCuts);
 171   //
 172   // post data
 173   //
 174   PostData(1, fListHist);
 175
 176
 177 }
 178
 179
 180
 181 //________________________________________________________________________
 182 void AliAnalysisTrackingUncertainties::UserExec(Option_t *)
 183 {
 184   //
 185   // main event loop
 186   //
 187   fESD = dynamic_cast<AliESDEvent*>( InputEvent() );
 188   if (!fESD) {
 189     PostData(1, fListHist);
 190     return;
 191   }
 192   //
 193   if (!fESDpid) fESDpid = ((AliESDInputHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetESDpid();
 194   //
 195   // Check Monte Carlo information and other access first:
 196   //
 197   AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 198   if (!eventHandler) {
 199     fMCtrue = kFALSE;
 200   }
 201   //
 202   // extract generated particles information
 203   //
 204   AliMCEvent* mcEvent = 0x0;
 205   AliStack* stack = 0x0;
 206   if (eventHandler) mcEvent = eventHandler->MCEvent();
 207   if (!mcEvent) {
 208     if (fMCtrue) return;
 209   }
 210   if (fMCtrue) {
 211     //
 212     stack = mcEvent->Stack();
 213     if (!stack) return;
 214     //
 215     for(Int_t i = 0; i < stack->GetNtrack(); i++) {
 216       /* at the moment nothing is needed here
 217       TParticle * trackMC = stack->Particle(i);
 218       Double_t rap = trackMC->Eta();
 219       Double_t y = trackMC->Y();
 220       Double_t pT  = trackMC->Pt();
 221       */
 222       //
 223     }
 224
 225   }
 226   //
 227   if (!fESDtrackCuts) {
 228     PostData(1, fListHist);
 229     return;
 230   }
 231   //
 232   // monitor vertex position and selection
 233   //
 234   TH2F * histVertexSelection = (TH2F *) fListHist->FindObject("histVertexSelection");
 235   //
 236   Float_t vertexZ = 0.;
 237   if (IsVertexAccepted(fESD, vertexZ)) {
 238     histVertexSelection->Fill(vertexZ, 0);
 239   } else {
 240     histVertexSelection->Fill(vertexZ, 1);
 241     return;
 242   }
 243   //
 244   // fill track cut variation histograms
 245   //
 246   ProcessTrackCutVariation();
 247   //
 248   // fill ITS->TPC matching histograms
 249   //
 250   ProcessItsTpcMatching();
 251
 252   // Post output data
 253   PostData(1, fListHist);
 254
 255 }
 256
 257
 258 //________________________________________________________________________
 259 void  AliAnalysisTrackingUncertainties::ProcessItsTpcMatching(){
 260   //
 261   // check how many its-sa tracks get matched to TPC
 262   //
 263   int ntr = fESD->GetNumberOfTracks();
 264   //
 265   // initialize histograms
 266   //
 267   TH2F * hNMatch      = (TH2F*) fListHist->FindObject("hNMatch");
 268   THnF * hBestMatch   = (THnF*) fListHist->FindObject("hBestMatch");
 269   TH2F * hAllMatch    = (TH2F*) fListHist->FindObject("hAllMatch");
 270   TH2F * hAllMatchGlo = (TH2F*) fListHist->FindObject("hAllMatchGlo");
 271   //
 272   TH2F * hNMatchBg      = (TH2F*) fListHist->FindObject("hNMatchBg");
 273   TH2F * hBestMatchBg   = (TH2F*) fListHist->FindObject("hBestMatchBg");
 274   TH2F * hAllMatchBg    = (TH2F*) fListHist->FindObject("hAllMatchBg");
 275   TH2F * hAllMatchGloBg = (TH2F*) fListHist->FindObject("hAllMatchGloBg");
 276   //
 277   for (int it=0;it<ntr;it++) {
 278     AliESDtrack* trSA = fESD->GetTrack(it);
 279     if (!trSA->IsOn(AliESDtrack::kITSpureSA) || !trSA->IsOn(AliESDtrack::kITSrefit)) continue;
 280     double pt = trSA->Pt();
 281     //
 282     Int_t nmatch = 0;
 283     for (int i=kMaxMatch;i--;) {fMatchChi[i]=0; fMatchTr[i]=0;} // reset array
 284     for (int it1=0;it1<ntr;it1++) {
 285       if (it1==it) continue;
 286       AliESDtrack* trESD = fESD->GetTrack(it1);
 287       if (!trESD->IsOn(AliESDtrack::kTPCrefit)) continue;
 288       Match(trSA,trESD, nmatch);
 289     }
 290     //
 291     hNMatch->Fill(pt,nmatch);
 292     if (nmatch>0) { // matched tracks
 293       Double_t vecHistMatch[kNumberOfAxes] = {fMatchChi[0], pt, trSA->Eta(), trSA->Phi(), 0};
 294       hBestMatch->Fill(vecHistMatch);
 295     } else {        // un-matched tracks -> should be in overflow bin
 296       Double_t vecHistMatch[kNumberOfAxes] = {9999999., pt, trSA->Eta(), trSA->Phi(), 0};
 297       hBestMatch->Fill(vecHistMatch);
 298     }
 299     //
 300     for (int imt=nmatch;imt--;) {
 301       hAllMatch->Fill(pt,fMatchChi[imt]);
 302       if (fMatchTr[imt]->IsOn(AliESDtrack::kITSrefit)) hAllMatchGlo->Fill(pt,fMatchChi[imt]);
 303     }
 304     //
 305     nmatch = 0;
 306     for (int i=kMaxMatch;i--;) {fMatchChi[i]=0; fMatchTr[i]=0;}
 307     for (int it1=0;it1<ntr;it1++) {
 308       if (it1==it) continue;
 309       AliESDtrack* trESD = fESD->GetTrack(it1);
 310       if (!trESD->IsOn(AliESDtrack::kTPCrefit)) continue;
 311       Match(trSA,trESD, nmatch, TMath::Pi());
 312     }
 313     //
 314     hNMatchBg->Fill(pt,nmatch);
 315     if (nmatch>0) hBestMatchBg->Fill(pt,fMatchChi[0]);
 316     for (int imt=nmatch;imt--;) {
 317       hAllMatchBg->Fill(pt,fMatchChi[imt]);
 318       if (fMatchTr[imt]->IsOn(AliESDtrack::kITSrefit)) hAllMatchGloBg->Fill(pt,fMatchChi[imt]);
 319     }
 320     //
 321   }
 322
 323
 324 }
 325
 326
 327 void AliAnalysisTrackingUncertainties::Match(const AliESDtrack* tr0, const AliESDtrack* tr1, Int_t &nmatch, Double_t rotate) {
 328   //
 329   // check if two tracks are matching, possible rotation for combinatoric backgr.
 330   //
 331   Float_t bField = fESD->GetMagneticField();
 332   //
 333   const AliExternalTrackParam* trtpc0 = tr1->GetInnerParam();
 334   if (!trtpc0) return;
 335   AliExternalTrackParam trtpc(*trtpc0);
 336   //
 337   if (TMath::Abs(rotate)>1e-5) {
 338     const double *par = trtpc.GetParameter();
 339     const double *cov = trtpc.GetCovariance();
 340     double alp = trtpc.GetAlpha() + rotate;
 341     trtpc.Set(trtpc.GetX(),alp,par,cov);
 342   }
 343   //
 344   if (!trtpc.Rotate(tr0->GetAlpha())) return;
 345   if (!trtpc.PropagateTo(tr0->GetX(),bField)) return;
 346   double chi2 = tr0->GetPredictedChi2(&trtpc);
 347   if (chi2>kMaxChi2) return;
 348   //
 349   int ins;
 350   for (ins=0;ins<nmatch;ins++) if (chi2<fMatchChi[ins]) break;
 351   if (ins>=kMaxMatch) return;
 352
 353   for (int imv=nmatch;imv>ins;imv--) {
 354     if (imv>=kMaxMatch) continue;
 355     fMatchTr[imv]  = fMatchTr[imv-1];
 356     fMatchChi[imv] = fMatchChi[imv-1];
 357   }
 358   fMatchTr[ins] = tr1;
 359   fMatchChi[ins] = chi2;
 360   nmatch++;
 361   if (nmatch>=kMaxMatch) nmatch = kMaxMatch;
 362   //
 363 }
 364
 365
 366 //________________________________________________________________________
 367 void AliAnalysisTrackingUncertainties::ProcessTrackCutVariation() {
 368   //
 369   // fill track cut variation histograms - undo cuts step-by-step and fill histograms
 370   //
 371   //
 372   // initialize histograms
 373   //
 374   THnF * histNcl         = (THnF *) fListHist->FindObject("histNcl");
 375   THnF * histChi2Tpc     = (THnF *) fListHist->FindObject("histChi2Tpc");
 376   THnF * histDcaZ        = (THnF *) fListHist->FindObject("histDcaZ");
 377   THnF * histSpd         = (THnF *) fListHist->FindObject("histSpd");
 378   THnF * histNcr         = (THnF *) fListHist->FindObject("histNcr");
 379   THnF * histCRoverFC    = (THnF *) fListHist->FindObject("histCRoverFC");
 380   THnF * histChi2Its     = (THnF *) fListHist->FindObject("histChi2Its");
 381   THnF * histTpcLength   = (THnF *) fListHist->FindObject("histTpcLength");
 382   THnF * histTpcItsMatch = (THnF *) fListHist->FindObject("histTpcItsMatch");
 383   //
 384   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z for the vertex cut
 385   //
 386   for (Int_t i=0;i<fESD->GetNumberOfTracks();++i) {
 387     //
 388     AliESDtrack *track =fESD->GetTrack(i);
 389     //
 390     // relevant variables
 391     //
 392     //Double_t pid        = Double_t(GetPid(track));
 393     //
 394     Int_t nclsTPC       = track->GetTPCncls();
 395     Float_t pT          = track->Pt();
 396     Float_t eta         = track->Eta();
 397     Float_t phi         = track->Phi();
 398     Float_t chi2TPC     = track->GetTPCchi2();
 399     Float_t ncrTPC      = track->GetTPCCrossedRows();
 400     Int_t nclsTPCF      = track->GetTPCNclsF();
 401     Float_t nCRoverFC   = track->GetTPCCrossedRows();
 402     Double_t chi2ITS    = track->GetITSchi2();
 403     Int_t nclsITS       = track->GetITSclusters(0);
 404     Float_t tpcLength   = 0.;
 405
 406     if (track->GetInnerParam() && track->GetESDEvent()) {
 407       tpcLength = track->GetLengthInActiveZone(1, 1.8, 220, track->GetESDEvent()->GetMagneticField());
 408     }
 409
 410     if (nclsTPC != 0) {
 411       chi2TPC /= nclsTPC;
 412     } else {
 413       chi2TPC = 999.;
 414     }
 415
 416     if (nclsTPCF !=0) {
 417       nCRoverFC /= nclsTPCF;
 418     } else {
 419       nCRoverFC = 999.;
 420     }
 421
 422     if (nclsITS != 0){
 423       chi2ITS /= nclsITS;
 424     }else {
 425       chi2ITS = 999.;
 426     }
 427     //
 428     track->GetImpactParameters(dca, cov);
 429     //
 430     // (1.) fill number of clusters histogram
 431     //
 432     Int_t minNclsTPC = fESDtrackCuts->GetMinNClusterTPC();
 433     fESDtrackCuts->SetMinNClustersTPC(0);
 434     if (fESDtrackCuts->AcceptTrack(track)) {
 435       for(Int_t iPid = 0; iPid < 6; iPid++) {
 436         Double_t vecHistNcl[kNumberOfAxes] = {nclsTPC, pT, eta, phi, iPid};
 437         if (IsConsistentWithPid(iPid, track)) histNcl->Fill(vecHistNcl);
 438       }
 439     }
 440     fESDtrackCuts->SetMinNClustersTPC(minNclsTPC);
 441     //
 442     // (2.) fill chi2 TPC histogram
 443     //
 444     Float_t maxChi2 = fESDtrackCuts->GetMaxChi2PerClusterTPC();
 445     fESDtrackCuts->SetMaxChi2PerClusterTPC(999.);
 446     if (fESDtrackCuts->AcceptTrack(track)) {
 447       for(Int_t iPid = 0; iPid < 6; iPid++) {
 448         Double_t vecHistChi2Tpc[kNumberOfAxes] = {chi2TPC, pT, eta, phi, iPid};
 449         if (IsConsistentWithPid(iPid, track)) histChi2Tpc->Fill(vecHistChi2Tpc);
 450       }
 451     }
 452     fESDtrackCuts->SetMaxChi2PerClusterTPC(maxChi2);
 453     //
 454     // (3.) fill dca_z histogram
 455     //
 456     Float_t maxDcaZ = fESDtrackCuts->GetMaxDCAToVertexZ();
 457     fESDtrackCuts->SetMaxDCAToVertexZ(999.);
 458     if (fESDtrackCuts->AcceptTrack(track)) {
 459       for(Int_t iPid = 0; iPid < 6; iPid++) {
 460         Double_t vecHistDcaZ[kNumberOfAxes] = {TMath::Abs(dca[1]), pT, eta, phi, iPid};
 461         if (IsConsistentWithPid(iPid, track)) histDcaZ->Fill(vecHistDcaZ);
 462       }
 463     }
 464     fESDtrackCuts->SetMaxDCAToVertexZ(maxDcaZ);
 465     //
 466     // (4.) fill hit in SPD histogram
 467     //
 468     fESDtrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kOff);
 469     if (fESDtrackCuts->AcceptTrack(track)) {
 470       Int_t hasPoint = 0;
 471       if (track->HasPointOnITSLayer(0) || track->HasPointOnITSLayer(1)) hasPoint = 1;
 472       for(Int_t iPid = 0; iPid < 6; iPid++) {
 473         Double_t vecHistSpd[kNumberOfAxes] = {hasPoint, pT, eta, phi, iPid};
 474         if (IsConsistentWithPid(iPid, track)) histSpd->Fill(vecHistSpd);
 475       }
 476     }
 477     fESDtrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
 478     //
 479     // (5.) fill number of crossed rows histogram
 480     //
 481     //Int_t minNcrTPC = fESDtrackCuts->GetMinNCrossedRowsTPC(); //wait for getter in ESDtrackCuts
 482     Int_t minNcrTPC = 0;  //for now use standard value from 2010 !!
 483     fESDtrackCuts->SetMinNCrossedRowsTPC(0);
 484     if (fESDtrackCuts->AcceptTrack(track)) {
 485       for(Int_t iPid = 0; iPid < 6; iPid++) {
 486         Double_t vecHistNcr[kNumberOfAxes] = {ncrTPC, pT, eta, phi, iPid};
 487         if (IsConsistentWithPid(iPid, track)) histNcr->Fill(vecHistNcr);
 488       }
 489     }
 490     fESDtrackCuts->SetMinNCrossedRowsTPC(minNcrTPC);
 491     //
 492     // (6.) fill crossed rows over findable clusters histogram
 493     //
 494     //Int_t minCRoverFC = fESDtrackCuts->GetMinRatioCrossedRowsOverFindableClustersTPC(); //wait for getter in ESDtrackCuts
 495     Int_t minCRoverFC = 0.; //for now use standard value from 2010 !!
 496     fESDtrackCuts->SetMinRatioCrossedRowsOverFindableClustersTPC(0.);
 497     if (fESDtrackCuts->AcceptTrack(track)) {
 498        for(Int_t iPid = 0; iPid < 6; iPid++) {
 499          Double_t vecHistCRoverFC[kNumberOfAxes] = {nCRoverFC, pT, eta, phi, iPid};
 500          if (IsConsistentWithPid(iPid, track)) histCRoverFC->Fill(vecHistCRoverFC);
 501        }
 502     }
 503     fESDtrackCuts->SetMinRatioCrossedRowsOverFindableClustersTPC(minCRoverFC);
 504     //
 505     // (7.) fill chi2 ITS histogram
 506     //
 507     Float_t maxChi2ITS = fESDtrackCuts->GetMaxChi2PerClusterITS();
 508     fESDtrackCuts->SetMaxChi2PerClusterITS(999.);
 509     if (fESDtrackCuts->AcceptTrack(track)) {
 510       for(Int_t iPid = 0; iPid < 6; iPid++) {
 511         Double_t vecHistChi2ITS[kNumberOfAxes] = {chi2ITS, pT, eta, phi, iPid};
 512         if (IsConsistentWithPid(iPid, track)) histChi2Its->Fill(vecHistChi2ITS);
 513       }
 514     }
 515     fESDtrackCuts->SetMaxChi2PerClusterITS(maxChi2ITS);
 516     //
 517     // (8.) fill active length in TPC histogram
 518     //
 519     Int_t minTpcLength = fESDtrackCuts->GetMinLengthActiveVolumeTPC();
 520     fESDtrackCuts->SetMinLengthActiveVolumeTPC(0);
 521     if (fESDtrackCuts->AcceptTrack(track)) {
 522       for(Int_t iPid = 0; iPid < 6; iPid++) {
 523         Double_t vecHistTpcLength[kNumberOfAxes] = {tpcLength, pT, eta, phi, iPid};
 524         if (IsConsistentWithPid(iPid, track)) histTpcLength->Fill(vecHistTpcLength);
 525       }
 526     }
 527     fESDtrackCuts->SetMinLengthActiveVolumeTPC(minTpcLength);
 528     //
 529     // (9.) fill TPC->ITS matching efficiency histogram
 530     //
 531     Bool_t isMatched = kFALSE;
 532     // remove all ITS requirements
 533     //
 534     // Leonardo and Emilia:
 535     //  -> if MC is available: fill it only for true primaries,
 536     //        --to be done for every cut?
 537     //  -> Postprocessing: plot histogram with 1 divided by histogram with 0 as a function of pT/eta/phi
 538     //  -> Do we want to remove the DCA cut?
 539     Bool_t refit=fESDtrackCuts->GetRequireITSRefit();
 540     Float_t chi2tpc= fESDtrackCuts->GetMaxChi2TPCConstrainedGlobal();
 541     Float_t chi2its= fESDtrackCuts->GetMaxChi2PerClusterITS();
 542     //TString str = fESDtrackCuts->GetMaxDCAToVertexXYPtDep();
 543
 544     fESDtrackCuts->SetRequireITSRefit(kFALSE);
 545     fESDtrackCuts->SetMaxChi2TPCConstrainedGlobal(99999.);
 546     fESDtrackCuts->SetMaxChi2PerClusterITS(999999.);
 547         //TString str = fESDtrackCuts->GetMaxDCAToVertexXYPtDep();
 548     //fESDtrackCuts->SetMaxDCAToVertexXYPtDep();
 549     fESDtrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kOff);
 550
 551     if (fESDtrackCuts->AcceptTrack(track)) {
 552       for(Int_t iPid = 0; iPid < 6; iPid++) {
 553         Double_t vecHistTpcItsMatch[kNumberOfAxes] = {isMatched, pT, eta, phi, iPid};
 554         if (IsConsistentWithPid(iPid, track)) histTpcItsMatch->Fill(vecHistTpcItsMatch); // fill with 1 here
 555       }
 556     }
 557     //apply back the cuts
 558     fESDtrackCuts->SetRequireITSRefit(refit);
 559     fESDtrackCuts->SetMaxChi2TPCConstrainedGlobal(chi2tpc);
 560     fESDtrackCuts->SetMaxChi2PerClusterITS(chi2its);
 561     //fESDtrackCuts->SetMaxDCAToVertexXYPtDep(str.Data());
 562     fESDtrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
 563     //set is matched
 564     isMatched=kTRUE;
 565     if (fESDtrackCuts->AcceptTrack(track)) {
 566       for(Int_t iPid = 0; iPid < 6; iPid++) {
 567         Double_t vecHistTpcItsMatch[kNumberOfAxes] = {isMatched, pT, eta, phi, iPid};
 568         if (IsConsistentWithPid(iPid, track)) histTpcItsMatch->Fill(vecHistTpcItsMatch); // fill with 0 here
 569       }
 570     }
 571
 572   } // end of track loop
 573
 574
 575 }
 576
 577
 578
 579 //________________________________________________________________________
 580 void AliAnalysisTrackingUncertainties::Terminate(Option_t *)
 581 {
 582   // Draw result to the screen
 583   // Called once at the end of the query
 584
 585
 586 }
 587
 588
 589 //________________________________________________________________________
 590 void AliAnalysisTrackingUncertainties::BinLogAxis(const THn *h, Int_t axisNumber) {
 591   //
 592   // Method for the correct logarithmic binning of histograms
 593   //
 594   TAxis *axis = h->GetAxis(axisNumber);
 595   int bins = axis->GetNbins();
 596
 597   Double_t from = axis->GetXmin();
 598   Double_t to = axis->GetXmax();
 599   Double_t *newBins = new Double_t[bins + 1];
 600
 601   newBins[0] = from;
 602   Double_t factor = pow(to/from, 1./bins);
 603
 604   for (int i = 1; i <= bins; i++) {
 605    newBins[i] = factor * newBins[i-1];
 606   }
 607   axis->Set(bins, newBins);
 608   delete [] newBins;
 609
 610 }
 611
 612
 613 //________________________________________________________________________
 614 Bool_t AliAnalysisTrackingUncertainties::IsVertexAccepted(AliESDEvent * esd, Float_t &vertexZ) {
 615   //
 616   // function to check if a proper vertex is reconstructed and write z-position in vertexZ
 617   //
 618   vertexZ = -999.;
 619   Bool_t vertexOkay = kFALSE;
 620   const AliESDVertex *vertex = esd->GetPrimaryVertexTracks();
 621   if (vertex->GetNContributors() < 1) {
 622     //
 623     vertex = esd->GetPrimaryVertexSPD();
 624     if (vertex->GetNContributors() < 1) {
 625       vertexOkay = kFALSE; }
 626     else {
 627       vertexOkay = kTRUE;
 628     }
 629     //
 630     TString vtxTyp = vertex->GetTitle();
 631     Double_t cov[6]={0};
 632     vertex->GetCovarianceMatrix(cov);
 633     Double_t zRes = TMath::Sqrt(cov[5]);
 634     if (vtxTyp.Contains("vertexer:Z") && (zRes>0.25)) vertexOkay = kFALSE;
 635   }
 636   else {
 637     vertexOkay = kTRUE;
 638   }
 639
 640   vertexZ = vertex->GetZ();
 641   return vertexOkay;
 642
 643 }
 644
 645 //________________________________________________________________________
 646 AliAnalysisTrackingUncertainties::ESpecies_t AliAnalysisTrackingUncertainties::GetPid(const AliESDtrack * const tr, Bool_t useTPCTOF) const {
 647     //
 648     // Determine particle species for a given track
 649     // Two approaches can be used: As default the selection is done using TPC-only, in addition
 650     // the TOF usage is optional. In case of TPC-TOF, a valid TOF signal has to be provided for
 651     // the given track. The identification is delegated to helper function for each species.
 652     // Tracks which are selected as more than one species (ambiguous decision) are rejected.
 653     //
 654     // @Return: Particles species (kUndef in case no identification is possible)
 655     //
 656     if(!fESDpid) return kUndef;
 657     if(useTPCTOF && !(tr->GetStatus() & AliVTrack::kTOFpid)) return kUndef;
 658
 659     Bool_t isElectron(kFALSE), isPion(kFALSE), isKaon(kFALSE), isProton(kFALSE);
 660     Int_t nspec(0);
 661     if((isElectron = IsElectron(tr, useTPCTOF))) nspec++;
 662     if((isPion = IsPion(tr, useTPCTOF))) nspec++;
 663     if((isKaon = IsKaon(tr, useTPCTOF))) nspec++;
 664     if((isProton = IsProton(tr,useTPCTOF))) nspec++;
 665     if(nspec != 1) return kUndef;   // No decision or ambiguous decision;
 666     if(isElectron) return kSpecElectron;
 667     if(isPion) return kSpecPion;
 668     if(isProton) return kSpecProton;
 669     if(isKaon) return kSpecKaon;
 670     return kUndef;
 671 }
 672
 673 //________________________________________________________________________
 674 Bool_t AliAnalysisTrackingUncertainties::IsElectron(const AliESDtrack * const tr, Bool_t useTPCTOF) const {
 675     //
 676     // Selection of electron candidates using the upper half of the TPC sigma band, starting at
 677     // the mean ignoring its shift, and going up to 3 sigma above the mean. In case TOF information
 678     // is available, tracks which are incompatible with electrons within 3 sigma are rejected. If
 679     // no TOF information is used, the momentum regions where the kaon and the proton line cross
 680     // the electron line are cut out using a 3 sigma cut around the kaon or proton line.
 681     //
 682
 683     Float_t nsigmaElectronTPC = fESDpid->NumberOfSigmasTPC(tr, AliPID::kElectron);
 684     if(nsigmaElectronTPC < 0 || nsigmaElectronTPC  > 3) return kFALSE;
 685
 686     if(useTPCTOF){
 687         Float_t nsigmaElectronTOF = fESDpid->NumberOfSigmasTOF(tr, AliPID::kElectron);
 688         if(TMath::Abs(nsigmaElectronTOF) > 3) return kFALSE;
 689         else return kTRUE;
 690     } else {
 691         Float_t nsigmaKaonTPC = fESDpid->NumberOfSigmasTPC(tr, AliPID::kKaon),
 692                 nsigmaProtonTPC =fESDpid->NumberOfSigmasTPC(tr, AliPID::kProton);
 693         if(TMath::Abs(nsigmaKaonTPC < 3) || TMath::Abs(nsigmaProtonTPC < 3)) return kFALSE;
 694         else return kTRUE;
 695     }
 696 }
 697
 698 //________________________________________________________________________
 699 Bool_t AliAnalysisTrackingUncertainties::IsConsistentWithPid(Int_t type, const AliESDtrack * const tr) {
 700   //
 701   // just check if the PID is consistent with a given hypothesis in order to
 702   // investigate effects which are only dependent on the energy loss.
 703   //
 704   if (type == kSpecPion)     return IsPion(tr);
 705   if (type == kSpecKaon)     return IsKaon(tr);
 706   if (type == kSpecProton)   return IsProton(tr);
 707   if (type == kSpecElectron) return IsElectron(tr);
 708   if (type == kAll)          return kTRUE;
 709   return kFALSE;
 710
 711 }
 712
 713 //________________________________________________________________________
 714 Bool_t AliAnalysisTrackingUncertainties::IsPion(const AliESDtrack * const tr, Bool_t /*useTPCPTOF*/) const{
 715   //
 716   // Selectron of pion candidates
 717   // @TODO: To be implemented
 718   //
 719   Float_t nsigmaPionTPC = fESDpid->NumberOfSigmasTPC(tr, AliPID::kPion);
 720   if (TMath::Abs(nsigmaPionTPC) < 3) return kTRUE;
 721   return kFALSE;
 722
 723 }
 724
 725 //________________________________________________________________________
 726 Bool_t AliAnalysisTrackingUncertainties::IsKaon(const AliESDtrack * const tr, Bool_t /*useTPCPTOF*/) const {
 727   //
 728   // Selection of kaon candidates
 729   // @TODO: To be implemented
 730   //
 731   Float_t nsigmaKaonTPC = fESDpid->NumberOfSigmasTPC(tr, AliPID::kKaon);
 732   if (TMath::Abs(nsigmaKaonTPC) < 3) return kTRUE;
 733   return kFALSE;
 734
 735 }
 736
 737 //________________________________________________________________________
 738 Bool_t AliAnalysisTrackingUncertainties::IsProton(const AliESDtrack * const tr, Bool_t /*useTPCPTOF*/) const{
 739   //
 740   // Selection of proton candidates
 741   // @TODO: To be implemented
 742   //
 743   Float_t nsigmaProtonTPC = fESDpid->NumberOfSigmasTPC(tr, AliPID::kProton);
 744   if (TMath::Abs(nsigmaProtonTPC) < 3) return kTRUE;
 745   return kFALSE;
 746
 747 }
 748
 749 //________________________________________________________________________
 750 void AliAnalysisTrackingUncertainties::InitializeTrackCutHistograms() {
 751   //
 752   // create histograms for the track cut studies
 753   //
 754   //
 755   // (1.) number of clusters
 756   //                               0-ncl, 1-pt, 2-eta,         3-phi, 4-pid(0-3 -> electron-proton, 4 -> undef, 5 -> all)
 757   Int_t    binsNcl[kNumberOfAxes] = { 40,   50,    20,            18,  6};
 758   Double_t minNcl[kNumberOfAxes]  = { 20,  0.1,    -1,             0, -0.5};
 759   Double_t maxNcl[kNumberOfAxes]  = {160,   20,    +1, 2*TMath::Pi(),  5.5};
 760   //
 761   TString axisNameNcl[kNumberOfAxes]  = {"ncl","pT","eta","phi","pid"};
 762   TString axisTitleNcl[kNumberOfAxes] = {"ncl","pT","eta","phi","pid"};
 763   //
 764   THnF * histNcl = new THnF("histNcl","number of clusters histogram",kNumberOfAxes, binsNcl, minNcl, maxNcl);
 765   BinLogAxis(histNcl, 1);
 766   fListHist->Add(histNcl);
 767   //
 768   for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
 769     histNcl->GetAxis(iaxis)->SetName(axisNameNcl[iaxis]);
 770     histNcl->GetAxis(iaxis)->SetTitle(axisTitleNcl[iaxis]);
 771   }
 772   //
 773   // (2.) chi2/cls-TPC
 774   //                                  0-chi2, 1-pt, 2-eta,         3-phi, 4-pid(0-3 -> electron-proton, 4 -> undef, 5 -> all)
 775   Int_t    binsChi2Tpc[kNumberOfAxes] = { 40,   50,    20,            18,  6};
 776   Double_t minChi2Tpc[kNumberOfAxes]  = {  0,  0.1,    -1,             0, -0.5};
 777   Double_t maxChi2Tpc[kNumberOfAxes]  = {  8,   20,    +1, 2*TMath::Pi(),  5.5};
 778   //
 779   TString axisNameChi2Tpc[kNumberOfAxes]  = {"chi2tpc","pT","eta","phi","pid"};
 780   TString axisTitleChi2Tpc[kNumberOfAxes] = {"chi2tpc","pT","eta","phi","pid"};
 781   //
 782   THnF * histChi2Tpc = new THnF("histChi2Tpc","chi2 per cls. in TPC",kNumberOfAxes, binsChi2Tpc, minChi2Tpc, maxChi2Tpc);
 783   BinLogAxis(histChi2Tpc, 1);
 784   fListHist->Add(histChi2Tpc);
 785   //
 786   for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
 787     histChi2Tpc->GetAxis(iaxis)->SetName(axisNameChi2Tpc[iaxis]);
 788     histChi2Tpc->GetAxis(iaxis)->SetTitle(axisTitleChi2Tpc[iaxis]);
 789   }
 790   //
 791   // (3.) dca_z
 792   //                               0-dcaZ, 1-pt, 2-eta,         3-phi, 4-pid(0-3 -> electron-proton, 4 -> undef, 5 -> all)
 793   Int_t    binsDcaZ[kNumberOfAxes] = { 20,   50,    20,            18,  6};
 794   Double_t minDcaZ[kNumberOfAxes]  = {  0,  0.1,    -1,             0, -0.5};
 795   Double_t maxDcaZ[kNumberOfAxes]  = {  4,   20,    +1, 2*TMath::Pi(),  5.5};
 796   //
 797   TString axisNameDcaZ[kNumberOfAxes]  = {"dcaZ","pT","eta","phi","pid"};
 798   TString axisTitleDcaZ[kNumberOfAxes] = {"dcaZ","pT","eta","phi","pid"};
 799   //
 800   THnF * histDcaZ = new THnF("histDcaZ","dca_z to prim. vtx.",kNumberOfAxes, binsDcaZ, minDcaZ, maxDcaZ);
 801   BinLogAxis(histDcaZ, 1);
 802   fListHist->Add(histDcaZ);
 803   //
 804   for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
 805     histDcaZ->GetAxis(iaxis)->SetName(axisNameDcaZ[iaxis]);
 806     histDcaZ->GetAxis(iaxis)->SetTitle(axisTitleDcaZ[iaxis]);
 807   }
 808   //
 809   // (4.) hit in SPD layer
 810   //                              0-spdHit, 1-pt, 2-eta,         3-phi, 4-pid(0-3 -> electron-proton, 4 -> undef, 5 -> all)
 811   Int_t    binsSpd[kNumberOfAxes] = {    2,   50,    20,            18,  6};
 812   Double_t minSpd[kNumberOfAxes]  = { -0.5,  0.1,    -1,             0, -0.5};
 813   Double_t maxSpd[kNumberOfAxes]  = {  1.5,   20,    +1, 2*TMath::Pi(),  5.5};
 814   //
 815   TString axisNameSpd[kNumberOfAxes]  = {"spdHit","pT","eta","phi","pid"};
 816   TString axisTitleSpd[kNumberOfAxes] = {"spdHit","pT","eta","phi","pid"};
 817   //
 818   THnF * histSpd = new THnF("histSpd","hit in SPD layer or not",kNumberOfAxes, binsSpd, minSpd, maxSpd);
 819   BinLogAxis(histSpd, 1);
 820   fListHist->Add(histSpd);
 821   //
 822   for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
 823     histSpd->GetAxis(iaxis)->SetName(axisNameSpd[iaxis]);
 824     histSpd->GetAxis(iaxis)->SetTitle(axisTitleSpd[iaxis]);
 825   }
 826   //
 827   // (5.) number of crossed rows
 828   //                               0-ncr, 1-pt, 2-eta,         3-phi, 4-pid(0,unid,etc.)
 829   Int_t    binsNcr[kNumberOfAxes] = { 40,   50,    20,            18,  6};
 830   Double_t minNcr[kNumberOfAxes]  = { 20,  0.1,    -1,             0, -0.5};
 831   Double_t maxNcr[kNumberOfAxes]  = {160,   20,    +1, 2*TMath::Pi(),  5.5};
 832   //
 833   TString axisNameNcr[kNumberOfAxes]  = {"Ncr","pT","eta","phi","pid"};
 834   TString axisTitleNcr[kNumberOfAxes] = {"Ncr","pT","eta","phi","pid"};
 835   //
 836   THnF * histNcr = new THnF("histNcr","number of crossed rows TPC histogram",kNumberOfAxes, binsNcr, minNcr, maxNcr);
 837   BinLogAxis(histNcr, 1);
 838   fListHist->Add(histNcr);
 839   //
 840   for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
 841     histNcr->GetAxis(iaxis)->SetName(axisNameNcr[iaxis]);
 842     histNcr->GetAxis(iaxis)->SetTitle(axisTitleNcr[iaxis]);
 843   }
 844   //
 845   // (6.) ratio crossed rows over findable clusters
 846   //                                0-CRoverFC, 1-pt, 2-eta,         3-phi, 4-pid(0,unid,etc.)
 847   Int_t    binsCRoverFC[kNumberOfAxes] = {  26,   50,    20,            18,  6};
 848   Double_t minCRoverFC[kNumberOfAxes]  = { 0.4,  0.1,    -1,             0, -0.5};
 849   Double_t maxCRoverFC[kNumberOfAxes]  = { 1.8,   20,    +1, 2*TMath::Pi(),  5.5};
 850   //
 851   TString axisNameCRoverFC[kNumberOfAxes]  = {"CRoverFC","pT","eta","phi","pid"};
 852   TString axisTitleCRoverFC[kNumberOfAxes] = {"CRoverFC","pT","eta","phi","pid"};
 853   //
 854   THnF * histCRoverFC = new THnF("histCRoverFC","number of crossed rows over findable clusters histogram",kNumberOfAxes, binsCRoverFC, minCRoverFC, maxCRoverFC);
 855   BinLogAxis(histCRoverFC, 1);
 856   fListHist->Add(histCRoverFC);
 857   //
 858   for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
 859     histCRoverFC->GetAxis(iaxis)->SetName(axisNameCRoverFC[iaxis]);
 860     histCRoverFC->GetAxis(iaxis)->SetTitle(axisTitleCRoverFC[iaxis]);
 861   }
 862   //
 863   // (7.) max chi2 / ITS cluster
 864   //                               0-Chi2Its, 1-pt, 2-eta,         3-phi, 4-pid(0,unid,etc.)
 865   Int_t    binsChi2Its[kNumberOfAxes] = { 25,   50,    20,            18,  6};
 866   Double_t minChi2Its[kNumberOfAxes]  = {  0,  0.1,    -1,             0, -0.5};
 867   Double_t maxChi2Its[kNumberOfAxes]  = { 50,   20,    +1, 2*TMath::Pi(),  5.5};
 868   //
 869   TString axisNameChi2Its[kNumberOfAxes]  = {"Chi2Its","pT","eta","phi","pid"};
 870   TString axisTitleChi2Its[kNumberOfAxes] = {"Chi2Its","pT","eta","phi","pid"};
 871   //
 872   THnF * histChi2Its = new THnF("histChi2Its","number of crossed rows TPC histogram",kNumberOfAxes, binsChi2Its, minChi2Its, maxChi2Its);
 873   BinLogAxis(histChi2Its, 1);
 874   fListHist->Add(histChi2Its);
 875   //
 876   for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
 877     histChi2Its->GetAxis(iaxis)->SetName(axisNameChi2Its[iaxis]);
 878     histChi2Its->GetAxis(iaxis)->SetTitle(axisTitleChi2Its[iaxis]);
 879   }
 880   //
 881   // (8.) tpc active volume length
 882   //                                0-TpcLength, 1-pt, 2-eta,         3-phi, 4-pid(0,unid,etc.)
 883   Int_t    binsTpcLength[kNumberOfAxes] = {  40,   50,    20,            18,  6};
 884   Double_t minTpcLength[kNumberOfAxes]  = {   0,  0.1,    -1,             0, -0.5};
 885   Double_t maxTpcLength[kNumberOfAxes]  = { 170,   20,    +1, 2*TMath::Pi(),  5.5};
 886   //
 887   TString axisNameTpcLength[kNumberOfAxes]  = {"TpcLength","pT","eta","phi","pid"};
 888   TString axisTitleTpcLength[kNumberOfAxes] = {"TpcLength","pT","eta","phi","pid"};
 889   //
 890   THnF * histTpcLength = new THnF("histTpcLength","number of crossed rows TPC histogram",kNumberOfAxes, binsTpcLength, minTpcLength, maxTpcLength);
 891   BinLogAxis(histTpcLength, 1);
 892   fListHist->Add(histTpcLength);
 893   //
 894   for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
 895     histTpcLength->GetAxis(iaxis)->SetName(axisNameTpcLength[iaxis]);
 896     histTpcLength->GetAxis(iaxis)->SetTitle(axisTitleTpcLength[iaxis]);
 897   }
 898   //
 899   // (9.) match TPC->ITS
 900   //                                  0-is matched, 1-pt, 2-eta,         3-phi, 4-pid(0-3 -> electron-proton, 4 -> undef, 5 -> all)
 901   Int_t    binsTpcItsMatch[kNumberOfAxes] = {    2,   50,    20,            18,  6};
 902   Double_t minTpcItsMatch[kNumberOfAxes]  = { -0.5,  0.1,    -1,             0, -0.5};
 903   Double_t maxTpcItsMatch[kNumberOfAxes]  = {  1.5,   20,    +1, 2*TMath::Pi(),  5.5};
 904   //
 905   TString axisNameTpcItsMatch[kNumberOfAxes]  = {"isMatched","pT","eta","phi","pid"};
 906   TString axisTitleTpcItsMatch[kNumberOfAxes] = {"isMatched","pT","eta","phi","pid"};
 907   //
 908   THnF * histTpcItsMatch = new THnF("histTpcItsMatch","TPC -> ITS matching",kNumberOfAxes, binsTpcItsMatch, minTpcItsMatch, maxTpcItsMatch);
 909   BinLogAxis(histTpcItsMatch, 1);
 910   fListHist->Add(histTpcItsMatch);
 911   //
 912   for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
 913     histTpcItsMatch->GetAxis(iaxis)->SetName(axisNameTpcItsMatch[iaxis]);
 914     histTpcItsMatch->GetAxis(iaxis)->SetTitle(axisTitleTpcItsMatch[iaxis]);
 915   }
 916
 917
 918
 919
 920 }
 921