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