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
 206   TFile* file = TFile::Open(filename.Data(), "READ");
 207   if ( ! file ) {
 208     AliFatal(Form("OADB file %s not found!", filename.Data()));
 209     return kFALSE;
 210   }
 211
 212   TString containerName = "MuonTrackCutsParam_data";
 213   if ( fIsMC ) containerName.ReplaceAll("_data", "_MC");
 214   TString defaultName = "default";
 215
 216   AliOADBContainer* oadbContainer = static_cast<AliOADBContainer*>(file->Get(containerName.Data()));
 217   if ( ! oadbContainer ) {
 218     AliFatal(Form("OADB container %s not found in %s!", containerName.Data(), filename.Data()));
 219     return kFALSE;
 220   }
 221   AliOADBMuonTrackCutsParam* param = static_cast<AliOADBMuonTrackCutsParam*>(oadbContainer->GetObject(runNumber,defaultName.Data(),passName));
 222
 223   TString paramName = param->GetName();
 224   if ( paramName == "default" ) {
 225     if ( ! fAllowDefaultParams ) {
 226       AliFatal(Form("Requested run %i not found in %s! Please check your pass name or allow default parameters", runNumber, passName.Data()));
 227       return kFALSE; // Coverity fix
 228     }
 229
 230     // Search if there is a better tuned default for the considered period
 231     containerName.Append("_def");
 232     oadbContainer = static_cast<AliOADBContainer*>(file->Get(containerName.Data()));
 233     param = static_cast<AliOADBMuonTrackCutsParam*>(oadbContainer->GetObject(runNumber,defaultName.Data(),passName));
 234
 235     AliWarning(Form("Requested run %i not found in %s: using %s", runNumber, passName.Data(), param->GetName()));
 236   }
 237
 238   fOADBParam = *param;
 239   file->Close();
 240
 241   AliInfo(Form("Requested run %i in pass %s. Param. set: %s", runNumber, passName.Data(), fOADBParam.GetName()));
 242
 243   return kTRUE;
 244 }
 245
 246 //________________________________________________________________________
 247 Bool_t AliMuonTrackCuts::IsSelected( TObject* obj )
 248 {
 249   /// Track is selected
 250   UInt_t filterMask = GetFilterMask();
 251   UInt_t selectionMask = GetSelectionMask(obj);
 252
 253   AliDebug(1, Form("IsMuon %i  selected %i  mask 0x%x", AliAnalysisMuonUtility::IsMuonTrack(static_cast<AliVParticle*>(obj)), ( selectionMask & filterMask ) == filterMask, selectionMask ));
 254
 255   return ( ( selectionMask & filterMask ) == filterMask );
 256 }
 257
 258
 259 //________________________________________________________________________
 260 UInt_t AliMuonTrackCuts::GetSelectionMask( const TObject* obj )
 261 {
 262   /// Get selection mask
 263
 264   const AliVParticle* track = static_cast<const AliVParticle*> ( obj );
 265
 266   UInt_t selectionMask = 0;
 267
 268   if ( ! AliAnalysisMuonUtility::IsMuonTrack(track) ) return selectionMask;
 269
 270   Double_t eta = track->Eta();
 271   if ( eta > -4. && eta < -2.5 ) selectionMask |= kMuEta;
 272
 273   Double_t thetaAbsEndDeg = AliAnalysisMuonUtility::GetThetaAbsDeg(track);
 274
 275   if ( thetaAbsEndDeg > 2. && thetaAbsEndDeg < 10. ) selectionMask |= kMuThetaAbs;
 276
 277   Int_t matchTrig = AliAnalysisMuonUtility::GetMatchTrigger(track);
 278   Int_t cutLevel[3] = {kMuMatchApt, kMuMatchLpt, kMuMatchHpt};
 279   Double_t pt = track->Pt();
 280   for ( Int_t ilevel=0; ilevel<3; ilevel++ ) {
 281     if ( matchTrig < ilevel+1 ) break;
 282     if ( fSharpPtCut && pt < fOADBParam.GetSharpPtCut(ilevel) ) break;
 283     selectionMask |= cutLevel[ilevel];
 284   }
 285
 286   if ( AliAnalysisMuonUtility::GetChi2perNDFtracker(track) < fOADBParam.GetChi2NormCut() ) selectionMask |= kMuTrackChiSquare;
 287
 288   TVector3 dcaAtVz = GetCorrectedDCA(track);
 289   Double_t pTotMean = GetAverageMomentum(track);
 290
 291   Double_t pDca = pTotMean * dcaAtVz.Mag();
 292
 293   Double_t pTot = track->P();
 294
 295   Double_t sigmaPdca = IsThetaAbs23(track) ? fOADBParam.GetSigmaPdca23() : fOADBParam.GetSigmaPdca310();
 296
 297   // Momentum resolution only
 298   // The cut depends on the momentum resolution. In particular:
 299   // Sigma_pDCA = Sqrt( Sigma_pDCA_measured^2 + Sigma_p/p * pDCA)
 300   // The relative momentum distribution Sigma_p/p is estimated from the track Delta_p,
 301   // and from the error on sagitta Delta_s:
 302   // Delta_p/p = p*Delta_s/(1+p*Delta_s)
 303   // A cut at N sigmas requres a cut in N Delta_s, so
 304   // Sigma_p/p(N) = p*N*Delta_s/(1+p*N*Delta_s)
 305   //Double_t pResolutionEffect = pDca * pTot * GetRelPResolution() / ( 1. + GetNSigmaPdca() * GetRelPResolution()*pTot );
 306
 307   //Double_t pResolutionEffect = 0.4 * pTot;  // Values used in 2010 data
 308   //Double_t pResolutionEffect = 0.32 * pTot; // Values in 2011
 309   //Double_t sigmaPdcaWithRes = TMath::Sqrt( sigmaPdca*sigmaPdca + pResolutionEffect*pResolutionEffect );
 310
 311   // Momentum resolution and slope resolution
 312   // Due to the momentum resolution, the measured momentum is biased
 313   // Since we want to keep as much signal as possible, we want to avoid
 314   // that a measured pxDCA is rejected since the momentum is overestimated
 315   // p_true = p_meas - Delta_p
 316   // p_true = p_meas - N*Delta_s*p_meas / (1+n*Delta_s*p_meas)
 317   // Hence:
 318   // p_true x DCA < N * Sigma_pDCA_meas =>
 319   // p_meas x DCA < N * Sigma_pDCA_meas / ( 1 - N*Delta_s*p_meas / (1+n*Delta_s*p_meas))
 320   // Finally the cut value has to be summed in quadrature with the error on DCA,
 321   // which is given by the slope resolution
 322   // 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)
 323   Double_t nrp = fOADBParam.GetNSigmaPdca() * fOADBParam.GetRelPResolution() * pTot;
 324   Double_t pResolutionEffect = sigmaPdca / ( 1. - nrp / ( 1. + nrp ) );
 325   Double_t slopeResolutionEffect = 535. * fOADBParam.GetSlopeResolution() * pTot;
 326
 327   Double_t sigmaPdcaWithRes = TMath::Sqrt( pResolutionEffect*pResolutionEffect + slopeResolutionEffect*slopeResolutionEffect );
 328
 329   if ( pDca < fOADBParam.GetNSigmaPdca() * sigmaPdcaWithRes ) selectionMask |= kMuPdca;
 330
 331   AliDebug(1, Form("Selection mask 0x%x\n", selectionMask));
 332
 333   return selectionMask;
 334 }
 335
 336
 337 //________________________________________________________________________
 338 Bool_t AliMuonTrackCuts::IsSelected( TList* /* list */)
 339 {
 340   /// Not implemented
 341   AliError("Function not implemented: Use IsSelected(TObject*)");
 342   return kFALSE;
 343 }
 344
 345
 346 //________________________________________________________________________
 347 Bool_t AliMuonTrackCuts::IsThetaAbs23 ( const AliVParticle* track ) const
 348 {
 349   /// Check if theta_abs is smaller than 3 degrees
 350   return ( AliAnalysisMuonUtility::GetThetaAbsDeg(track) < 3. );
 351 }
 352
 353
 354 //________________________________________________________________________
 355 TVector3 AliMuonTrackCuts::GetCorrectedDCA ( const AliVParticle* track ) const
 356 {
 357   /// Get corrected DCA
 358
 359   TVector3 vertex(AliAnalysisMuonUtility::GetXatVertex(track), AliAnalysisMuonUtility::GetYatVertex(track), AliAnalysisMuonUtility::GetZatVertex(track));
 360
 361   TVector3 dcaTrack(AliAnalysisMuonUtility::GetXatDCA(track), AliAnalysisMuonUtility::GetYatDCA(track), AliAnalysisMuonUtility::GetZatDCA(track));
 362
 363   TVector3 dcaAtVz = dcaTrack - vertex - fOADBParam.GetMeanDCA();
 364
 365   return dcaAtVz;
 366 }
 367
 368 //________________________________________________________________________
 369 Double_t AliMuonTrackCuts::GetAverageMomentum ( const AliVParticle* track ) const
 370 {
 371   /// Get average momentum before and after the absorber
 372
 373   Double_t pTotMean = 0.;
 374   Double_t pTot = track->P();
 375   //if ( isESD ) pTotMean = 0.5 * ( pTot + ((AliESDMuonTrack*)track)->PUncorrected() );
 376   if ( ! AliAnalysisMuonUtility::IsAODTrack(track) ) pTotMean = ((AliESDMuonTrack*)track)->PUncorrected(); // Increased stability if using uncorrected value
 377   else {
 378     Double_t meanPcorr = IsThetaAbs23(track) ? fOADBParam.GetMeanPCorr23() : fOADBParam.GetMeanPCorr310();
 379     pTotMean = pTot - meanPcorr;
 380   }
 381
 382   return pTotMean;
 383 }
 384
 385
 386 //________________________________________________________________________
 387 void AliMuonTrackCuts::SetDefaultFilterMask ()
 388 {
 389   /// Standard cuts for single muon
 390   SetFilterMask ( kMuEta | kMuThetaAbs | kMuPdca | kMuMatchApt );
 391 }
 392
 393 //________________________________________________________________________
 394 Bool_t AliMuonTrackCuts::TrackPtCutMatchTrigClass ( const AliVParticle* track, const TArrayI ptCutFromClass) const
 395 {
 396   /// Check if track passes the trigger pt cut level used in the trigger class
 397   Int_t matchTrig = AliAnalysisMuonUtility::GetMatchTrigger(track);
 398   Bool_t matchTrackerPt = kTRUE;
 399   if ( IsApplySharpPtCutInMatching() ) {
 400     matchTrackerPt = ( track->Pt() >= fOADBParam.GetSharpPtCut(ptCutFromClass[0]-1,kFALSE) );
 401   }
 402   Bool_t passCut = ( ( matchTrig >= ptCutFromClass[0] ) && matchTrackerPt );
 403   AliDebug(1,Form("Class matchTrig %i %i  trackMatchTrig %i  trackPt %g (required %i)  passCut %i", ptCutFromClass[0], ptCutFromClass[1], matchTrig, track->Pt(), IsApplySharpPtCutInMatching(),passCut));
 404   return passCut;
 405 }
 406
 407
 408 //________________________________________________________________________
 409 void AliMuonTrackCuts::Print(Option_t* option) const
 410 {
 411   //
 412   /// Print info
 413   //
 414   TString sopt(option);
 415   sopt.ToLower();
 416   if ( sopt.IsNull() || sopt.Contains("*") || sopt.Contains("all") ) sopt = "mask param";
 417   UInt_t filterMask = GetFilterMask();
 418   Int_t cutLevel[3] = {kMuMatchApt, kMuMatchLpt, kMuMatchHpt};
 419   TString cutLevelName[3] = {"Apt", "Lpt", "Hpt"};
 420   if ( sopt.Contains("mask") ) {
 421     printf(" *** Muon track filter mask: *** \n");
 422     printf("  0x%x\n", filterMask);
 423     if ( filterMask & kMuEta ) printf("  -4 < eta < -2.5\n");
 424     if ( filterMask & kMuThetaAbs ) printf("  2 < theta_abs < 10 deg\n");
 425     if ( filterMask & kMuPdca ) printf("  pxDCA cut\n");
 426     for ( Int_t ilevel=0; ilevel<3; ilevel++ ) {
 427       if ( filterMask & cutLevel[ilevel] ) {
 428         printf("  match %s", cutLevelName[ilevel].Data());
 429         if ( fSharpPtCut ) printf(" && sharp pt from tracker");
 430         printf("\n");
 431       }
 432     }
 433     if ( filterMask & kMuTrackChiSquare ) printf("  Chi2 cut on track\n");
 434     printf(" ******************** \n");
 435   }
 436   if ( sopt.Contains("param") ) fOADBParam.Print();
 437 }