PWG3/muon/AliMuonTrackCuts.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 #include "AliMuonTrackCuts.h"
  17
  18 #include "TMath.h"
  19 #include "TList.h"
  20 #include "TArrayD.h"
  21 #include "TObjArray.h"
  22 #include "TObjString.h"
  23 #include "TFile.h"
  24 #include "TParameter.h"
  25 #include "TKey.h"
  26 #include "TVector3.h"
  27
  28 #include "AliLog.h"
  29 #include "AliVParticle.h"
  30 #include "AliESDMuonTrack.h"
  31 #include "AliAODTrack.h"
  32
  33 /// \cond CLASSIMP
  34 ClassImp(AliMuonTrackCuts) // Class implementation in ROOT context
  35 /// \endcond
  36
  37
  38 //________________________________________________________________________
  39 AliMuonTrackCuts::AliMuonTrackCuts() :
  40   AliAnalysisCuts(),
  41   fIsESD(kFALSE),
  42   fIsMC(kFALSE),
  43   fUseCustomParam(kFALSE),
  44   fParameters(TArrayD(kNParameters))
  45 {
  46   /// Default ctor.
  47   fParameters.Reset();
  48 }
  49
  50 //________________________________________________________________________
  51 AliMuonTrackCuts::AliMuonTrackCuts(const char* name, const char* title, Bool_t isESD ) :
  52   AliAnalysisCuts(name, title),
  53   fIsESD(isESD),
  54   fIsMC(kFALSE),
  55   fUseCustomParam(kFALSE),
  56   fParameters(TArrayD(kNParameters))
  57 {
  58   /// Constructor
  59   fParameters.Reset();
  60   SetDefaultFilterMask();
  61 }
  62
  63
  64 //________________________________________________________________________
  65 AliMuonTrackCuts::AliMuonTrackCuts(const AliMuonTrackCuts& obj) :
  66   AliAnalysisCuts(obj),
  67   fIsESD(obj.fIsESD),
  68   fIsMC(obj.fIsMC),
  69   fUseCustomParam(obj.fUseCustomParam),
  70   fParameters(obj.fParameters)
  71 {
  72   /// Copy constructor
  73 }
  74
  75
  76 //________________________________________________________________________
  77 AliMuonTrackCuts::~AliMuonTrackCuts()
  78 {
  79   /// Destructor
  80 }
  81
  82 //________________________________________________________________________
  83 Bool_t AliMuonTrackCuts::RunMatchesRange( Int_t runNumber, const Char_t* objName )
  84 {
  85   /// Check if the object contains the run
  86   TString sname(objName);
  87   TObjArray* array = sname.Tokenize("_");
  88   array->SetOwner();
  89   Int_t runRange[2] = { -1, -1 };
  90   if ( array->GetEntries() >= 3 ) {
  91     for ( Int_t irun=0; irun<2; ++irun ) {
  92       TString currRun = array->At(irun+1)->GetName();
  93       if ( currRun.IsDigit() ) runRange[irun] = currRun.Atoi();
  94     }
  95   }
  96   delete array;
  97   return ( runNumber >= runRange[0] && runNumber <= runRange[1]);
  98 }
  99
 100 //________________________________________________________________________
 101 void AliMuonTrackCuts::SetUseCustomParam( Bool_t useCustomParam, Int_t runNumber  )
 102 {
 103   /// Flag to select custom parameters
 104   /// It first searches the default parameters in OADB
 105   /// then disables the access to the OADB
 106   /// and allows to manually modify parameters
 107
 108   if ( ! fUseCustomParam && useCustomParam ) SetRun(runNumber);
 109   fUseCustomParam = useCustomParam;
 110 }
 111
 112 //________________________________________________________________________
 113 Bool_t AliMuonTrackCuts::SetRun( Int_t runNumber )
 114 {
 115   /// Get parameters from OADB for runNumber
 116
 117   if ( fUseCustomParam ) return kFALSE;
 118   return StreamParameters(runNumber, -1);
 119 }
 120
 121
 122 //________________________________________________________________________
 123 Bool_t AliMuonTrackCuts::StreamParameters( Int_t runNumber,  Int_t runMax )
 124 {
 125   if ( runMax > 0 ) { // Stream to OADB
 126     if ( ! fUseCustomParam ) {
 127       AliError("Users are not allowed to update OADB. Use SetUseCustomParam() instead");
 128       return kFALSE;
 129     }
 130   }
 131
 132   TString filename = "$ALICE_ROOT/OADB/PWG3/MuonTrackCuts.root";
 133   if ( fIsMC ) filename.ReplaceAll(".root", "_MC.root");
 134
 135   TString parNames[kNParameters];
 136   parNames[kMeanDcaX]       = "MeanDcaX";
 137   parNames[kMeanDcaY]       = "MeanDcaY";
 138   parNames[kMeanDcaZ]       = "MeanDcaZ";
 139   parNames[kMeanPCorr23]    = "MeanPCorr23";
 140   parNames[kMeanPCorr310]   = "MeanPCorr310";
 141   parNames[kSigmaPdca23]    = "SigmaPdca23";
 142   parNames[kSigmaPdca310]   = "SigmaPdca310";
 143   parNames[kNSigmaPdcaCut]  = "NSigmaPdcaCut";
 144   parNames[kChi2NormCut]    = "Chi2NormCut";
 145   parNames[kRelPResolution] = "RelPResolution";
 146
 147   TObjArray* paramList = 0x0;
 148
 149   if ( runMax < 0 ) { // Get from OADB
 150     TFile* file = TFile::Open(filename.Data(), "READ");
 151     if ( ! file ) {
 152       AliError(Form("OADB file %s not found!", filename.Data()));
 153       return kFALSE;
 154     }
 155
 156     TList* listOfKeys = file->GetListOfKeys();
 157     TIter next(listOfKeys);
 158     TObject* key = 0x0;
 159     Bool_t foundMatch = kFALSE;
 160     TObject* defaultObj = 0x0;
 161     while ( ( key = next() ) ) {
 162       TString objName = key->GetName();
 163       objName.ToUpper();
 164       if ( RunMatchesRange(runNumber, objName.Data()) ) {
 165         paramList = static_cast<TObjArray*>(file->Get(key->GetName()));
 166         foundMatch = kTRUE;
 167         break;
 168       }
 169       if ( objName.Contains("DEFAULT") ) defaultObj = file->Get(key->GetName());
 170     }
 171
 172     if ( ! foundMatch ) {
 173       AliWarning("Run number not found in OADB: using default");
 174       if ( defaultObj ) paramList = static_cast<TObjArray*>(defaultObj);
 175       else {
 176         file->Close();
 177         AliError("Default parameters not found in OADB!");
 178         return kFALSE;
 179       }
 180     }
 181
 182     AliInfo(Form("Required run %i. Param. set: %s", runNumber, paramList->GetName()));
 183
 184     for ( Int_t ipar=0; ipar<kNParameters; ++ipar ) {
 185       TParameter<Double_t>* param = static_cast<TParameter<Double_t>*>(paramList->FindObject(parNames[ipar].Data()));
 186       if ( ! param ) {
 187         AliWarning(Form("Parameter %s not set", parNames[ipar].Data()));
 188         continue;
 189       }
 190       fParameters[ipar] = param->GetVal();
 191     }
 192
 193     file->Close();
 194   }
 195   else {
 196     if ( ! paramList ) {
 197       paramList = new TObjArray(kNParameters);
 198       paramList->SetOwner();
 199     }
 200     for ( Int_t ipar=0; ipar<kNParameters; ++ipar ) {
 201       TParameter<Double_t>* param= new TParameter<Double_t>(parNames[ipar].Data(), fParameters[ipar]);
 202       paramList->AddAt(param, ipar);
 203     }
 204
 205     TString paramListName = "MuonCuts_";
 206     paramListName += ( runNumber < 0 ) ? "Default" : Form("%i_%i",runNumber, runMax);
 207     AliInfo(Form("Adding %s to file %s", paramListName.Data(), filename.Data()));
 208     paramList->SetName(paramListName.Data());
 209     TFile* file = TFile::Open(filename.Data(), "UPDATE");
 210     paramList->Write(paramListName.Data(), TObject::kSingleKey);
 211     file->Close();
 212     delete paramList;
 213   }
 214   return kTRUE;
 215 }
 216
 217 //________________________________________________________________________
 218 Bool_t AliMuonTrackCuts::SetParameter(Int_t iparam, Float_t value)
 219 {
 220   /// Set parameter
 221   if ( fUseCustomParam ) {
 222     fParameters.SetAt(value, iparam);
 223     return kTRUE;
 224   }
 225
 226   AliWarning("Parameters automatically taken from OADB. If you want to use with custom parameters, use SetUseCustomParam()");
 227   return kFALSE;
 228 }
 229
 230
 231 //________________________________________________________________________
 232 Bool_t AliMuonTrackCuts::IsSelected( TObject* obj )
 233 {
 234   /// Track is selected
 235   UInt_t filterMask = GetFilterMask();
 236   UInt_t selectionMask = GetSelectionMask(obj);
 237
 238   return ( ( selectionMask & filterMask ) == filterMask );
 239 }
 240
 241
 242 //________________________________________________________________________
 243 UInt_t AliMuonTrackCuts::GetSelectionMask( const TObject* obj )
 244 {
 245   /// Get selection mask
 246
 247   const AliVParticle* track = static_cast<const AliVParticle*> ( obj );
 248
 249   UInt_t selectionMask = 0;
 250
 251   Bool_t isMuon = ( fIsESD ) ? ((AliESDMuonTrack*)track)->ContainTrackerData() :  ((AliAODTrack*)track)->IsMuonTrack();
 252
 253   if ( ! isMuon ) return selectionMask;
 254
 255   Double_t eta = track->Eta();
 256   if ( eta > -4. && eta < -2.5 ) selectionMask |= kMuEta;
 257
 258   Double_t rAbsEnd =  ( fIsESD ) ? ((AliESDMuonTrack*)track)->GetRAtAbsorberEnd() : ((AliAODTrack*)track)->GetRAtAbsorberEnd();
 259   Double_t thetaAbsEndDeg = TMath::ATan( rAbsEnd / 505. ) * TMath::RadToDeg();
 260
 261   if ( thetaAbsEndDeg > 2. && thetaAbsEndDeg < 10. ) selectionMask |= kMuThetaAbs;
 262
 263   Int_t matchTrig = ( fIsESD ) ? ((AliESDMuonTrack*)track)->GetMatchTrigger() : ((AliAODTrack*)track)->GetMatchTrigger();
 264   if ( matchTrig >= 1 ) selectionMask |= kMuMatchApt;
 265   if ( matchTrig >= 2 ) selectionMask |= kMuMatchLpt;
 266   if ( matchTrig >= 3 ) selectionMask |= kMuMatchHpt;
 267
 268   Double_t chi2norm = ( fIsESD ) ? ((AliESDMuonTrack*)track)->GetNormalizedChi2() : ((AliAODTrack*)track)->Chi2perNDF();
 269   if ( chi2norm < GetChi2NormCut() ) selectionMask |= kMuTrackChiSquare;
 270
 271   TVector3 dcaAtVz = GetCorrectedDCA(track);
 272   Double_t pTotMean = GetAverageMomentum(track);
 273
 274   Double_t pDca = pTotMean * dcaAtVz.Mag();
 275
 276   Double_t pTot = track->P();
 277
 278   // Momentum resolution only
 279   // The cut depends on the momentum resolution. In particular:
 280   // Sigma_pDCA = Sqrt( Sigma_pDCA_measured^2 + Sigma_p/p * pDCA)
 281   // The relative momentum distribution Sigma_p/p is estimated from the track Delta_p,
 282   // and from the error on sagitta Delta_s:
 283   // Delta_p/p = p*Delta_s/(1+p*Delta_s)
 284   // A cut at N sigmas requres a cut in N Delta_s, so
 285   // Sigma_p/p(N) = p*N*Delta_s/(1+p*N*Delta_s)
 286   //Double_t pResolutionEffect = pDca * pTot * GetRelPResolution() / ( 1. + GetNSigmaPdca() * GetRelPResolution()*pTot );
 287
 288   //Double_t pResolutionEffect = 0.4 * pTot;  // Values used in 2010 data
 289   //Double_t pResolutionEffect = 0.32 * pTot; // Values in 2011
 290   //Double_t sigmaPdca = GetSigmaPdca(rAbsEnd);
 291   //Double_t sigmaPdcaWithRes = TMath::Sqrt( sigmaPdca*sigmaPdca + pResolutionEffect*pResolutionEffect );
 292
 293   // Momentum resolution and slope resolution
 294   // Due to the momentum resolution, the measured momentum is biased
 295   // Since we want to keep as much signal as possible, we want to avoid
 296   // that a measured pxDCA is rejected since the momentum is overestimated
 297   // p_true = p_meas - Delta_p
 298   // p_true = p_meas - N*Delta_s*p_meas / (1+n*Delta_s*p_meas)
 299   // Hence:
 300   // p_true x DCA < N * Sigma_pDCA_meas =>
 301   // p_meas x DCA < N * Sigma_pDCA_meas / ( 1 - N*Delta_s*p_meas / (1+n*Delta_s*p_meas))
 302   // Finally the cut value has to be summed in quadrature with the error on DCA,
 303   // which is given by the slope resolution
 304   // p_meas x DCA < N * Sqrt( ( Sigma_pDCA_meas / ( 1 - N*Delta_s*p_meas / (1+n*Delta_s*p_meas)) )^2 + (distance * sigma_slope * p_meas )^2)
 305   Double_t nrp = GetNSigmaPdca() * GetRelPResolution() * pTot;
 306   Double_t pResolutionEffect = GetSigmaPdca(rAbsEnd) / ( 1. - nrp / ( 1. + nrp ) );
 307   Double_t slopeResolutionEffect = 535. * GetSlopeResolution() * pTot;
 308
 309   Double_t sigmaPdcaWithRes = TMath::Sqrt( pResolutionEffect*pResolutionEffect + slopeResolutionEffect*slopeResolutionEffect );
 310
 311   if ( pDca < GetNSigmaPdca() * sigmaPdcaWithRes ) selectionMask |= kMuPdca;
 312
 313   return selectionMask;
 314 }
 315
 316
 317 //________________________________________________________________________
 318 Bool_t AliMuonTrackCuts::IsSelected( TList* /* list */)
 319 {
 320   /// Not implemented
 321   AliError("Function not implemented: Use IsSelected(TObject*)");
 322   return kFALSE;
 323 }
 324
 325
 326 //________________________________________________________________________
 327 Int_t AliMuonTrackCuts::GetThetaAbsBin ( Double_t rAtAbsEnd ) const
 328 {
 329   /// Get theta abs bin
 330   Double_t thetaAbsEndDeg = TMath::ATan( rAtAbsEnd / 505. ) * TMath::RadToDeg();
 331   Int_t thetaAbsBin = ( thetaAbsEndDeg < 3. ) ? kThetaAbs23 : kThetaAbs310;
 332   return thetaAbsBin;
 333 }
 334
 335
 336 //________________________________________________________________________
 337 TVector3 AliMuonTrackCuts::GetCorrectedDCA ( const AliVParticle* track ) const
 338 {
 339   /// Get corrected DCA
 340
 341   Double_t vtxPos[3];
 342   if ( fIsESD ) {
 343     vtxPos[0] = ((AliESDMuonTrack*)track)->GetNonBendingCoor();
 344     vtxPos[1] = ((AliESDMuonTrack*)track)->GetBendingCoor();
 345     vtxPos[2] = ((AliESDMuonTrack*)track)->GetZ();
 346   }
 347   else ((AliAODTrack*)track)->GetXYZ(vtxPos);
 348
 349   TVector3 vertex(vtxPos);
 350
 351   TVector3 dcaTrack(0.,0., vtxPos[2]);
 352   dcaTrack.SetX( ( fIsESD ) ? ((AliESDMuonTrack*)track)->GetNonBendingCoorAtDCA() : ((AliAODTrack*)track)->XAtDCA() );
 353   dcaTrack.SetY( ( fIsESD ) ? ((AliESDMuonTrack*)track)->GetBendingCoorAtDCA() : ((AliAODTrack*)track)->YAtDCA() );
 354
 355   TVector3 dcaAtVz = dcaTrack - vertex - GetMeanDCA();
 356
 357   return dcaAtVz;
 358 }
 359
 360 //________________________________________________________________________
 361 Double_t AliMuonTrackCuts::GetAverageMomentum ( const AliVParticle* track ) const
 362 {
 363   /// Get average momentum before and after the absorber
 364
 365   Double_t pTotMean = 0.;
 366
 367   Double_t pTot = track->P();
 368   //if ( fIsESD ) pTotMean = 0.5 * ( pTot + ((AliESDMuonTrack*)track)->PUncorrected() );
 369   if ( fIsESD ) pTotMean = ((AliESDMuonTrack*)track)->PUncorrected(); // Increased stability if using uncorrected value
 370   else {
 371     pTotMean = pTot - GetMeanPCorr(((AliAODTrack*)track)->GetRAtAbsorberEnd());
 372   }
 373
 374   return pTotMean;
 375 }
 376
 377
 378 //________________________________________________________________________
 379 void AliMuonTrackCuts::SetMeanDCA ( Double_t xAtDca, Double_t yAtDca, Double_t zAtDca )
 380 {
 381     /// Set mean DCA from track
 382   SetParameter(kMeanDcaX, xAtDca);
 383   SetParameter(kMeanDcaY, yAtDca);
 384   SetParameter(kMeanDcaZ, zAtDca);
 385 }
 386
 387 //________________________________________________________________________
 388 TVector3 AliMuonTrackCuts::GetMeanDCA () const
 389 {
 390     /// Get mean DCA from track
 391   return TVector3(fParameters[kMeanDcaX], fParameters[kMeanDcaY], fParameters[kMeanDcaZ]);
 392 }
 393
 394 //________________________________________________________________________
 395 void AliMuonTrackCuts::SetMeanPCorr ( Double_t pCorrThetaAbs23, Double_t pCorrThetaAbs310 )
 396 {
 397   /// Set mean p correction
 398   SetParameter(kMeanPCorr23, pCorrThetaAbs23);
 399   SetParameter(kMeanPCorr310, pCorrThetaAbs310);
 400 }
 401
 402 //________________________________________________________________________
 403 Double_t AliMuonTrackCuts::GetMeanPCorr ( Double_t rAtAbsEnd ) const
 404 {
 405   /// Get mean p correction
 406   return fParameters[kMeanPCorr23+GetThetaAbsBin(rAtAbsEnd)];
 407 }
 408
 409 //________________________________________________________________________
 410 void AliMuonTrackCuts::SetSigmaPdca ( Double_t sigmaThetaAbs23, Double_t sigmaThetaAbs310 )
 411 {
 412   /// Set sigma pdca
 413   SetParameter(kSigmaPdca23, sigmaThetaAbs23);
 414   SetParameter(kSigmaPdca310, sigmaThetaAbs310);
 415 }
 416
 417 //________________________________________________________________________
 418 Double_t AliMuonTrackCuts::GetSigmaPdca ( Double_t rAtAbsEnd ) const
 419 {
 420   /// Get mean pdca
 421   return fParameters[kSigmaPdca23+GetThetaAbsBin(rAtAbsEnd)];
 422 }
 423
 424 //________________________________________________________________________
 425 void AliMuonTrackCuts::SetNSigmaPdca ( Double_t nSigmas )
 426 {
 427   /// Set N sigma pdca cut
 428   SetParameter(kNSigmaPdcaCut, nSigmas);
 429 }
 430
 431 //________________________________________________________________________
 432 Double_t AliMuonTrackCuts::GetNSigmaPdca () const
 433 {
 434   /// Get N sigma pdca cut
 435   return fParameters[kNSigmaPdcaCut];
 436 }
 437
 438 //________________________________________________________________________
 439 void AliMuonTrackCuts::SetChi2NormCut ( Double_t chi2normCut )
 440 {
 441   /// Set cut on normalized chi2 of tracks
 442   SetParameter(kChi2NormCut,chi2normCut);
 443 }
 444
 445 //________________________________________________________________________
 446 Double_t AliMuonTrackCuts::GetChi2NormCut () const
 447 {
 448   /// Get cut on normalized chi2 of tracks
 449   return fParameters[kChi2NormCut];
 450 }
 451
 452 //________________________________________________________________________
 453 void AliMuonTrackCuts::SetRelPResolution ( Double_t relPResolution )
 454 {
 455   /// Set relative momentum resolution
 456   SetParameter(kRelPResolution,relPResolution);
 457 }
 458
 459 //________________________________________________________________________
 460 Double_t AliMuonTrackCuts::GetRelPResolution () const
 461 {
 462   /// Get relative momentum resolution
 463   return fParameters[kRelPResolution];
 464 }
 465
 466
 467 //________________________________________________________________________
 468 void AliMuonTrackCuts::SetSlopeResolution ( Double_t slopeResolution )
 469 {
 470   /// Set slope resolution
 471   SetParameter(kSlopeResolution,slopeResolution);
 472 }
 473
 474 //________________________________________________________________________
 475 Double_t AliMuonTrackCuts::GetSlopeResolution () const
 476 {
 477   /// Get slope resolution
 478   return fParameters[kSlopeResolution];
 479 }
 480
 481 //________________________________________________________________________
 482 void AliMuonTrackCuts::SetDefaultFilterMask ()
 483 {
 484   /// Standard cuts for single muon
 485   SetFilterMask ( kMuEta | kMuThetaAbs | kMuPdca | kMuMatchApt );
 486 }
 487
 488 //________________________________________________________________________
 489 void AliMuonTrackCuts::Print(Option_t* option) const
 490 {
 491   //
 492   /// Print info
 493   //
 494   TString sopt(option);
 495   sopt.ToLower();
 496   if ( sopt.IsNull() || sopt.Contains("*") || sopt.Contains("all") ) sopt = "mask param";
 497   UInt_t filterMask = GetFilterMask();
 498   if ( sopt.Contains("mask") ) {
 499     printf(" *** Muon track filter mask: *** \n");
 500     printf("  0x%x\n", filterMask);
 501     if ( filterMask & kMuEta ) printf("  -4 < eta < -2.5\n");
 502     if ( filterMask & kMuThetaAbs ) printf("  2 < theta_abs < 10 deg\n");
 503     if ( filterMask & kMuPdca ) printf("  pxDCA cut\n");
 504     if ( filterMask & kMuMatchApt ) printf("  match Apt\n");
 505     if ( filterMask & kMuMatchLpt ) printf("  match Lpt\n");
 506     if ( filterMask & kMuMatchHpt ) printf("  match Hpt\n");
 507     if ( filterMask & kMuTrackChiSquare ) printf("  Chi2 cut on track\n");
 508     printf(" ******************** \n");
 509   }
 510   if ( sopt.Contains("param") ) {
 511     printf(" *** Muon track parameter summary: ***\n");
 512     printf("  Mean vertex DCA: (%g, %g, %g)\n", fParameters[kMeanDcaX], fParameters[kMeanDcaY], fParameters[kMeanDcaZ]);
 513     printf("  Mean p correction (GeV/c): theta2-3 = %g  theta3-10 = %g\n", fParameters[kMeanPCorr23], fParameters[kMeanPCorr310]);
 514     printf("  Sigma p x DCA (cm x GeV/c): theta2-3 = %g  theta3-10 = %g\n", fParameters[kSigmaPdca23], fParameters[kSigmaPdca310]);
 515     printf("  Cut p x DCA in units of sigma: %g\n", fParameters[kNSigmaPdcaCut]);
 516     printf("  Cut on track chi2/NDF: %g\n", fParameters[kChi2NormCut]);
 517     printf("  Momentum resolution: %g\n", fParameters[kRelPResolution]);
 518     printf("  Slope resolution: %g\n", fParameters[kSlopeResolution]);
 519     printf(" ********************************\n");
 520   }
 521 }