PWG/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 "TArrayI.h"
  22 #include "TObjArray.h"
  23 #include "TObjString.h"
  24 #include "TFile.h"
  25 #include "TParameter.h"
  26 #include "TKey.h"
  27 #include "TVector3.h"
  28 #include "TSystem.h"
  29
  30 #include "AliLog.h"
  31 #include "AliVParticle.h"
  32 #include "AliESDMuonTrack.h"
  33 #include "AliAnalysisManager.h"
  34 #include "AliInputEventHandler.h"
  35 #include "AliVEvent.h"
  36
  37 #include "AliProdInfo.h"
  38 #include "AliOADBContainer.h"
  39
  40 #include "AliAnalysisMuonUtility.h"
  41
  42 /// \cond CLASSIMP
  43 ClassImp(AliMuonTrackCuts) // Class implementation in ROOT context
  44 /// \endcond
  45
  46
  47 //________________________________________________________________________
  48 AliMuonTrackCuts::AliMuonTrackCuts() :
  49   AliAnalysisCuts(),
  50   fIsMC(kFALSE),
  51   fUseCustomParam(kFALSE),
  52   fSharpPtCut(kFALSE),
  53   fAllowDefaultParams(kFALSE),
  54   fPassName(""),
  55   fOADBParam()
  56 {
  57   /// Default ctor.
  58 }
  59
  60 //________________________________________________________________________
  61 AliMuonTrackCuts::AliMuonTrackCuts(const char* name, const char* title ) :
  62 AliAnalysisCuts(name, title),
  63   fIsMC(kFALSE),
  64   fUseCustomParam(kFALSE),
  65   fSharpPtCut(kFALSE),
  66   fAllowDefaultParams(kFALSE),
  67   fPassName(""),
  68   fOADBParam("muonTrackCutsParam")
  69 {
  70   /// Constructor
  71   SetDefaultFilterMask();
  72 }
  73
  74
  75 //________________________________________________________________________
  76 AliMuonTrackCuts::AliMuonTrackCuts(const AliMuonTrackCuts& obj) :
  77   AliAnalysisCuts(obj),
  78   fIsMC(obj.fIsMC),
  79   fUseCustomParam(obj.fUseCustomParam),
  80   fSharpPtCut(obj.fSharpPtCut),
  81   fAllowDefaultParams(obj.fAllowDefaultParams),
  82   fPassName(obj.fPassName),
  83   fOADBParam(obj.fOADBParam)
  84 {
  85   /// Copy constructor
  86 }
  87
  88
  89 //________________________________________________________________________
  90 AliMuonTrackCuts& AliMuonTrackCuts::operator=(const AliMuonTrackCuts& obj)
  91 {
  92   /// Assignment operator
  93   if ( this != &obj ) {
  94     AliAnalysisCuts::operator=(obj);
  95     fIsMC = obj.fIsMC;
  96     fUseCustomParam = obj.fUseCustomParam;
  97     fSharpPtCut = obj.fSharpPtCut;
  98     fAllowDefaultParams = obj.fAllowDefaultParams;
  99     fPassName = obj.fPassName;
 100     fOADBParam = obj.fOADBParam;
 101   }
 102   return *this;
 103 }
 104
 105
 106 //________________________________________________________________________
 107 AliMuonTrackCuts::~AliMuonTrackCuts()
 108 {
 109   /// Destructor
 110 }
 111
 112 //________________________________________________________________________
 113 void AliMuonTrackCuts::SetPassNumber ( Int_t passNumber )
 114 {
 115   /// Set pass number (for backward compatibility)
 116   AliWarning("Obsolete: please use SetPassName instead");
 117   fPassName = Form("pass%i",passNumber);
 118 }
 119
 120 //________________________________________________________________________
 121 void AliMuonTrackCuts::SetCustomParamFromRun( Int_t runNumber, TString passName )
 122 {
 123   /// It first searches the default parameters in OADB
 124   /// then disables the access to the OADB
 125   /// and allows to manually modify parameters
 126
 127   if ( passName.IsNull() ) passName = fPassName;
 128   else fPassName = passName;
 129   ReadParamFromOADB ( runNumber, passName );
 130   fUseCustomParam = kTRUE;
 131   AliWarning ("From now on SetRun does NOTHING!!");
 132 }
 133
 134 //________________________________________________________________________
 135 AliOADBMuonTrackCutsParam* AliMuonTrackCuts::CustomParam ( )
 136 {
 137   /// Returning the muon track cuts parameters (not const, so you can change the parameters)
 138   /// CAVEAT: if you only want to Get the parameters, please use GetMuonTrackCutsParam()
 139   /// If you want to modify the parameters, you need to call SetCustomParamFromRun at least once,
 140   /// otherwise CustomParam returns a null pointer.
 141
 142   if ( ! fUseCustomParam ) {
 143     AliError("This method allows you to modify the paramaters.\nIf you only want to get them, please use GetMuonTrackCutsParam instead.\nOtherwise, please call at least once SetCustomParamFromRun.");
 144     return 0x0;
 145   }
 146   return &fOADBParam;
 147 }
 148
 149 //________________________________________________________________________
 150 void AliMuonTrackCuts::SetCustomParam ( const AliInputEventHandler* eventHandler )
 151 {
 152   /// It tries to sets the parameters from the OADB
 153   /// for the current run, then gives the possiblity to the user
 154   /// to change them in his task
 155
 156   Int_t runNumber = eventHandler->GetEvent()->GetRunNumber();
 157   TString passName = GuessPass(eventHandler);
 158   SetCustomParamFromRun( runNumber, passName );
 159 }
 160
 161 //________________________________________________________________________
 162 TString AliMuonTrackCuts::GuessPass ( const AliInputEventHandler* eventHandler )
 163 {
 164   /// If pass not defined guess it from event handler
 165
 166   TString passName = fPassName;
 167   if ( passName.IsNull() ) {
 168     // Pass number not set by user
 169     // First try to get it from data
 170     //    AliProdInfo prodInfo(eventHandler->GetUserInfo());
 171     //    if ( prodInfo.GetRecoPass() >= 0 ) {
 172     //      passNumber = prodInfo.GetRecoPass();
 173     //      AliInfo(Form("Getting pass number from prodInfo: pass%i", passNumber));
 174     //    }
 175     //    else {
 176     // If not found, try to guess it from data path
 177     passName = AliAnalysisMuonUtility::GetPassName(eventHandler);
 178     AliInfo(Form("Guessing pass name from path: %s", passName.Data()));
 179     //    }
 180   }
 181   return passName;
 182 }
 183
 184 //________________________________________________________________________
 185 Bool_t AliMuonTrackCuts::SetRun ( const AliInputEventHandler* eventHandler )
 186 {
 187   /// Get parameters from OADB for current run
 188
 189   if ( fUseCustomParam ) return kFALSE;
 190   Int_t runNumber = eventHandler->GetEvent()->GetRunNumber();
 191   TString passName = GuessPass(eventHandler);
 192   return ReadParamFromOADB ( runNumber, passName );
 193 }
 194
 195
 196 //________________________________________________________________________
 197 Bool_t AliMuonTrackCuts::ReadParamFromOADB ( Int_t runNumber, TString passName )
 198 {
 199
 200   /// Read parameters from OADB
 201
 202   if ( passName.IsNull() && ! fAllowDefaultParams ) AliFatal("Pass name not specified! Please provide one or allow for default parameters");
 203
 204   TString filename = Form("%s/PWG/MUON/MuonTrackCuts.root",AliAnalysisManager::GetOADBPath());
 205   if ( fIsMC ) filename.ReplaceAll(".root", "_MC.root");
 206
 207   TFile* file = TFile::Open(filename.Data(), "READ");
 208   if ( ! file ) {
 209     AliFatal(Form("OADB file %s not found!", filename.Data()));
 210     return kFALSE;
 211   }
 212
 213   // Search the container name to find the correct pass
 214   AliOADBContainer* matchContainer = 0x0, *defaultContainer = 0x0;
 215   AliOADBMuonTrackCutsParam* matchParams = 0x0, *defaultParams = 0x0;
 216
 217   TList* listOfKeys = file->GetListOfKeys();
 218   TIter next(listOfKeys);
 219   TObject* key = 0x0;
 220   // loop on keys
 221   while ( ( key = next() ) ) {
 222
 223     TString checkName(key->GetName());
 224     checkName.ToUpper();
 225     Bool_t isDefault = checkName.Contains("DEFAULT");
 226     // if user selects a specific pass name, check for it
 227     // otherwise use default
 228     if ( isDefault ) {
 229       if ( ! fAllowDefaultParams ) continue;
 230     }
 231     else if ( passName.CompareTo(key->GetName()) ) continue;
 232
 233     AliOADBContainer* oadbContainer = static_cast<AliOADBContainer*> (file->Get(key->GetName()));
 234     // Check if the found parameters are default or match the requested run
 235     AliOADBMuonTrackCutsParam* currParams = static_cast<AliOADBMuonTrackCutsParam*> (oadbContainer->GetObject(runNumber, "default"));
 236     if ( ! currParams ) continue;
 237     if ( isDefault ) {
 238       defaultContainer = oadbContainer;
 239       defaultParams = currParams;
 240     }
 241     else {
 242       matchContainer = oadbContainer;
 243       matchParams = currParams;
 244       break;
 245     }
 246   } // loop on keys
 247
 248   AliOADBContainer* selectedContainer = 0x0;
 249   if ( matchParams ) {
 250     selectedContainer = matchContainer;
 251     fOADBParam = *matchParams;
 252   }
 253   else if ( defaultParams ) {
 254     selectedContainer = defaultContainer;
 255     fOADBParam = *defaultParams;
 256     AliWarning(Form("Requested run %i not found in %s: using %s (%s)", runNumber, passName.Data(), fOADBParam.GetName(), selectedContainer->GetName()));
 257   }
 258   else {
 259     AliFatal(Form("Requested run %i not found in %s! Please check your pass name or allow default parameters", runNumber, passName.Data()));
 260     return kFALSE; // Coverity fix
 261   }
 262
 263   file->Close();
 264
 265   AliInfo(Form("Requested run %i in pass %s. Param. set: %s (%s)", runNumber, passName.Data(), fOADBParam.GetName(), selectedContainer->GetName()));
 266
 267   return kTRUE;
 268 }
 269
 270 //________________________________________________________________________
 271 Bool_t AliMuonTrackCuts::IsSelected( TObject* obj )
 272 {
 273   /// Track is selected
 274   UInt_t filterMask = GetFilterMask();
 275   UInt_t selectionMask = GetSelectionMask(obj);
 276
 277   AliDebug(1, Form("IsMuon %i  selected %i  mask 0x%x", AliAnalysisMuonUtility::IsMuonTrack(static_cast<AliVParticle*>(obj)), ( selectionMask & filterMask ) == filterMask, selectionMask ));
 278
 279   return ( ( selectionMask & filterMask ) == filterMask );
 280 }
 281
 282
 283 //________________________________________________________________________
 284 UInt_t AliMuonTrackCuts::GetSelectionMask( const TObject* obj )
 285 {
 286   /// Get selection mask
 287
 288   const AliVParticle* track = static_cast<const AliVParticle*> ( obj );
 289
 290   UInt_t selectionMask = 0;
 291
 292   if ( ! AliAnalysisMuonUtility::IsMuonTrack(track) ) return selectionMask;
 293
 294   Double_t eta = track->Eta();
 295   if ( eta > -4. && eta < -2.5 ) selectionMask |= kMuEta;
 296
 297   Double_t thetaAbsEndDeg = AliAnalysisMuonUtility::GetThetaAbsDeg(track);
 298
 299   if ( thetaAbsEndDeg > 2. && thetaAbsEndDeg < 10. ) selectionMask |= kMuThetaAbs;
 300
 301   Int_t matchTrig = AliAnalysisMuonUtility::GetMatchTrigger(track);
 302   Int_t cutLevel[3] = {kMuMatchApt, kMuMatchLpt, kMuMatchHpt};
 303   Double_t pt = track->Pt();
 304   for ( Int_t ilevel=0; ilevel<3; ilevel++ ) {
 305     if ( matchTrig < ilevel+1 ) break;
 306     if ( fSharpPtCut && pt < fOADBParam.GetSharpPtCut(ilevel) ) break;
 307     selectionMask |= cutLevel[ilevel];
 308   }
 309
 310   if ( AliAnalysisMuonUtility::GetChi2perNDFtracker(track) < fOADBParam.GetChi2NormCut() ) selectionMask |= kMuTrackChiSquare;
 311
 312   TVector3 dcaAtVz = GetCorrectedDCA(track);
 313   Double_t pTotMean = GetAverageMomentum(track);
 314
 315   Double_t pDca = pTotMean * dcaAtVz.Mag();
 316
 317   Double_t pTot = track->P();
 318
 319   Double_t sigmaPdca = IsThetaAbs23(track) ? fOADBParam.GetSigmaPdca23() : fOADBParam.GetSigmaPdca310();
 320
 321   // Momentum resolution only
 322   // The cut depends on the momentum resolution. In particular:
 323   // Sigma_pDCA = Sqrt( Sigma_pDCA_measured^2 + Sigma_p/p * pDCA)
 324   // The relative momentum distribution Sigma_p/p is estimated from the track Delta_p,
 325   // and from the error on sagitta Delta_s:
 326   // Delta_p/p = p*Delta_s/(1+p*Delta_s)
 327   // A cut at N sigmas requres a cut in N Delta_s, so
 328   // Sigma_p/p(N) = p*N*Delta_s/(1+p*N*Delta_s)
 329   //Double_t pResolutionEffect = pDca * pTot * GetRelPResolution() / ( 1. + GetNSigmaPdca() * GetRelPResolution()*pTot );
 330
 331   //Double_t pResolutionEffect = 0.4 * pTot;  // Values used in 2010 data
 332   //Double_t pResolutionEffect = 0.32 * pTot; // Values in 2011
 333   //Double_t sigmaPdcaWithRes = TMath::Sqrt( sigmaPdca*sigmaPdca + pResolutionEffect*pResolutionEffect );
 334
 335   // Momentum resolution and slope resolution
 336   // Due to the momentum resolution, the measured momentum is biased
 337   // Since we want to keep as much signal as possible, we want to avoid
 338   // that a measured pxDCA is rejected since the momentum is overestimated
 339   // p_true = p_meas - Delta_p
 340   // p_true = p_meas - N*Delta_s*p_meas / (1+n*Delta_s*p_meas)
 341   // Hence:
 342   // p_true x DCA < N * Sigma_pDCA_meas =>
 343   // p_meas x DCA < N * Sigma_pDCA_meas / ( 1 - N*Delta_s*p_meas / (1+n*Delta_s*p_meas))
 344   // Finally the cut value has to be summed in quadrature with the error on DCA,
 345   // which is given by the slope resolution
 346   // 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)
 347   Double_t nrp = fOADBParam.GetNSigmaPdca() * fOADBParam.GetRelPResolution() * pTot;
 348   Double_t pResolutionEffect = sigmaPdca / ( 1. - nrp / ( 1. + nrp ) );
 349   Double_t slopeResolutionEffect = 535. * fOADBParam.GetSlopeResolution() * pTot;
 350
 351   Double_t sigmaPdcaWithRes = TMath::Sqrt( pResolutionEffect*pResolutionEffect + slopeResolutionEffect*slopeResolutionEffect );
 352
 353   if ( pDca < fOADBParam.GetNSigmaPdca() * sigmaPdcaWithRes ) selectionMask |= kMuPdca;
 354
 355   AliDebug(1, Form("Selection mask 0x%x\n", selectionMask));
 356
 357   return selectionMask;
 358 }
 359
 360
 361 //________________________________________________________________________
 362 Bool_t AliMuonTrackCuts::IsSelected( TList* /* list */)
 363 {
 364   /// Not implemented
 365   AliError("Function not implemented: Use IsSelected(TObject*)");
 366   return kFALSE;
 367 }
 368
 369
 370 //________________________________________________________________________
 371 Bool_t AliMuonTrackCuts::IsThetaAbs23 ( const AliVParticle* track ) const
 372 {
 373   /// Check if theta_abs is smaller than 3 degrees
 374   return ( AliAnalysisMuonUtility::GetThetaAbsDeg(track) < 3. );
 375 }
 376
 377
 378 //________________________________________________________________________
 379 TVector3 AliMuonTrackCuts::GetCorrectedDCA ( const AliVParticle* track ) const
 380 {
 381   /// Get corrected DCA
 382
 383   TVector3 vertex(AliAnalysisMuonUtility::GetXatVertex(track), AliAnalysisMuonUtility::GetYatVertex(track), AliAnalysisMuonUtility::GetZatVertex(track));
 384
 385   TVector3 dcaTrack(AliAnalysisMuonUtility::GetXatDCA(track), AliAnalysisMuonUtility::GetYatDCA(track), AliAnalysisMuonUtility::GetZatDCA(track));
 386
 387   TVector3 dcaAtVz = dcaTrack - vertex - fOADBParam.GetMeanDCA();
 388
 389   return dcaAtVz;
 390 }
 391
 392 //________________________________________________________________________
 393 Double_t AliMuonTrackCuts::GetAverageMomentum ( const AliVParticle* track ) const
 394 {
 395   /// Get average momentum before and after the absorber
 396
 397   Double_t pTotMean = 0.;
 398   Double_t pTot = track->P();
 399   //if ( isESD ) pTotMean = 0.5 * ( pTot + ((AliESDMuonTrack*)track)->PUncorrected() );
 400   if ( ! AliAnalysisMuonUtility::IsAODTrack(track) ) pTotMean = ((AliESDMuonTrack*)track)->PUncorrected(); // Increased stability if using uncorrected value
 401   else {
 402     Double_t meanPcorr = IsThetaAbs23(track) ? fOADBParam.GetMeanPCorr23() : fOADBParam.GetMeanPCorr310();
 403     pTotMean = pTot - meanPcorr;
 404   }
 405
 406   return pTotMean;
 407 }
 408
 409
 410 //________________________________________________________________________
 411 void AliMuonTrackCuts::SetDefaultFilterMask ()
 412 {
 413   /// Standard cuts for single muon
 414   SetFilterMask ( kMuEta | kMuThetaAbs | kMuPdca | kMuMatchApt );
 415 }
 416
 417 //________________________________________________________________________
 418 Bool_t AliMuonTrackCuts::TrackPtCutMatchTrigClass ( const AliVParticle* track, const TArrayI ptCutFromClass) const
 419 {
 420   /// Check if track passes the trigger pt cut level used in the trigger class
 421   Int_t matchTrig = AliAnalysisMuonUtility::GetMatchTrigger(track);
 422   Bool_t matchTrackerPt = kTRUE;
 423   if ( IsApplySharpPtCutInMatching() ) {
 424     matchTrackerPt = ( track->Pt() >= fOADBParam.GetSharpPtCut(ptCutFromClass[0]-1,kFALSE) );
 425   }
 426   Bool_t passCut = ( ( matchTrig >= ptCutFromClass[0] ) && matchTrackerPt );
 427   AliDebug(1,Form("Class matchTrig %i %i  trackMatchTrig %i  trackPt %g (required %i)  passCut %i", ptCutFromClass[0], ptCutFromClass[1], matchTrig, track->Pt(), IsApplySharpPtCutInMatching(),passCut));
 428   return passCut;
 429 }
 430
 431
 432 //________________________________________________________________________
 433 void AliMuonTrackCuts::Print(Option_t* option) const
 434 {
 435   //
 436   /// Print info
 437   //
 438   TString sopt(option);
 439   sopt.ToLower();
 440   if ( sopt.IsNull() || sopt.Contains("*") || sopt.Contains("all") ) sopt = "mask param";
 441   UInt_t filterMask = GetFilterMask();
 442   Int_t cutLevel[3] = {kMuMatchApt, kMuMatchLpt, kMuMatchHpt};
 443   TString cutLevelName[3] = {"Apt", "Lpt", "Hpt"};
 444   if ( sopt.Contains("mask") ) {
 445     printf(" *** Muon track filter mask: *** \n");
 446     printf("  0x%x\n", filterMask);
 447     if ( filterMask & kMuEta ) printf("  -4 < eta < -2.5\n");
 448     if ( filterMask & kMuThetaAbs ) printf("  2 < theta_abs < 10 deg\n");
 449     if ( filterMask & kMuPdca ) printf("  pxDCA cut\n");
 450     for ( Int_t ilevel=0; ilevel<3; ilevel++ ) {
 451       if ( filterMask & cutLevel[ilevel] ) {
 452         printf("  match %s", cutLevelName[ilevel].Data());
 453         if ( fSharpPtCut ) printf(" && sharp pt from tracker");
 454         printf("\n");
 455       }
 456     }
 457     if ( filterMask & kMuTrackChiSquare ) printf("  Chi2 cut on track\n");
 458     printf(" ******************** \n");
 459   }
 460   if ( sopt.Contains("param") ) fOADBParam.Print();
 461 }