]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG/muon/AliAnalysisMuMuFromESD.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Transition PWG3 --> PWGHF
[u/mrichter/AliRoot.git] / PWG / muon / AliAnalysisMuMuFromESD.cxx
CommitLineData
5eabea87 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// $Id$
17
18#include "AliAnalysisMuMuFromESD.h"
19
20#include "AliAnalysisManager.h"
21#include "AliInputEventHandler.h"
22#include "AliHistogramCollection.h"
23#include "AliAODEvent.h"
24#include "AliESDEvent.h"
25#include "AliESDMuonTrack.h"
26#include "AliESDVertex.h"
27#include "AliVParticle.h"
28#include "Riostream.h"
29#include "TCanvas.h"
30#include "TDirectory.h"
31#include "TH1F.h"
32#include "TH2F.h"
33#include "TLorentzVector.h"
34#include "TMap.h"
35#include "TMath.h"
36#include "TObjArray.h"
37#include "TPaveText.h"
38#include <cassert>
39#include "AliPhysicsSelection.h"
40#include "AliPhysicsSelectionTask.h"
41#include "AliBackgroundSelection.h"
42#include "AliCounterCollection.h"
43#include "THashList.h"
44
45#include "TGrid.h"
46
47///
48/// AliAnalysisMuMuFromESD
49///
50/// cuts must be added through AddSingleCut and AddPaircut methods
51///
52///
53
54ClassImp(AliAnalysisMuMuFromESD)
55
56//_____________________________________________________________________________
57AliAnalysisMuMuFromESD::AliAnalysisMuMuFromESD()
58: AliAnalysisMuMu(), fVertex(0x0)
59{
60 // default ctor
61}
62
63//_____________________________________________________________________________
64AliAnalysisMuMuFromESD::AliAnalysisMuMuFromESD(Bool_t aa)
65:AliAnalysisMuMu(kTRUE,aa), fVertex(0x0)
66{
67 // Constructor
68 Ctor();
69}
70
71//_____________________________________________________________________________
72AliAnalysisMuMuFromESD::AliAnalysisMuMuFromESD(TList* triggerClasses)
73:AliAnalysisMuMu(kTRUE,triggerClasses), fVertex(0x0)
74{
75 // Constructor
76 Ctor();
77}
78
79//_____________________________________________________________________________
80void AliAnalysisMuMuFromESD::Ctor()
81{
82 // Common ctor
83
84 if ( fAA )
85 {
86 AddGlobalEventSelection("SD2");
87 }
88 else
89 {
90 AddGlobalEventSelection("ALL");
91 }
92
93 AddGlobalEventSelection("PS");
94
95// fBranchNames = "ESD:AliESDRun.,AliESDHeader.,MuonTracks";
96}
97
98//_____________________________________________________________________________
99AliAnalysisMuMuFromESD::~AliAnalysisMuMuFromESD()
100{
101 /// dtor
102 delete fVertex;
103}
104
105//_____________________________________________________________________________
106void AliAnalysisMuMuFromESD::MuUserExec(Option_t*)
107{
108 // Main loop
109 // Called for each event
110
111 delete fVertex;
112 fVertex = 0x0;
113
114 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
115
116 if (esd)
117 {
118
119 if ( fAA )
120 {
121 // consider only events with OSM2 fired
122
123 UInt_t sd2 = (1<<12);
124 UInt_t trigger = esd->GetHeader()->GetL0TriggerInputs();
125
126 Bool_t ok = ( ( trigger & sd2 ) == sd2);
127
128
129 if (!ok) return;
130
131 }
132
133 TString runNumber(RunNumber(*esd));
134
135 TString triggers = esd->GetFiredTriggerClasses();
136
137 if ( IsDynamicTriggerClasses() ) AddTriggerClasses(triggers.Data());
138
139 TObjArray* a = triggers.Tokenize(" ");
140 TIter next(a);
141 TObjString* tname;
142
143 AliInputEventHandler* inputEventHandler = static_cast<AliInputEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
144
145 Bool_t isPhysicsSelected = ( inputEventHandler->IsEventSelected() & AliVEvent::kINT7 )
146 || ( inputEventHandler->IsEventSelected() & AliVEvent::kMUSH7 )
147 || ( inputEventHandler->IsEventSelected() & AliVEvent::kMUL7 )
148 || ( inputEventHandler->IsEventSelected() & AliVEvent::kMUS7 )
149 || ( inputEventHandler->IsEventSelected() & AliVEvent::kMUU7 );
150
151 TString centrality(DefaultCentralityName());
152
153 fVertex = new AliESDVertex(*esd->GetPrimaryVertexSPD());
154
155 if ( fAA )
156 {
157 // FIXME : should get it from centrality task...
158
159 // Double_t cent(-100);
160 // Double_t v0mult = vzero->GetMTotV0A()+vzero->GetMTotV0C();
161
162 // for ( std::vector<double>::size_type i = 0 ; i < fCentralityLimits.size()-1 && cent < -10 ; ++i )
163 // {
164 // if ( v0mult > fCentralityAmplitude[i] )
165 // {
166 // cent = fCentralityLimits[i];
167 // }
168 // }
169 //
170 // if ( cent > -1 )
171 // {
172 // centrality = CentralityName(cent);
173 // }
174 }
175
176 TString eventtype("all");
177 TString et;
178
179 if ( fAA ) eventtype="sd2";
180 et = eventtype;
181 et.ToUpper();
182
183 while ( ( tname = static_cast<TObjString*>(next()) ) )
184 {
185 if ( !fTriggerClasses->FindObject(tname->GetName()) )
186 {
187 // skip the trigger classes we're not supposed to analyse
188
189 continue;
190 }
191
192 fEventCounters->Count(Form("event:%s/trigger:%s/run:%d", eventtype.Data(),tname->GetName(), runNumber.Atoi()));
193
194 AssertHistogramCollection(et.Data(),tname->String().Data());
195
196 FillHistos(et.Data(),tname->String().Data(),centrality.Data(),*esd);
197
198 if ( isPhysicsSelected )
199 {
200
201 fEventCounters->Count(Form("event:ps/trigger:%s/run:%d", tname->GetName(), runNumber.Atoi()));
202
203 AssertHistogramCollection("PS",tname->String().Data());
204
205 FillHistos("PS",tname->String().Data(),centrality.Data(),*esd);
206 }
207 }
208
209 delete a;
210 }
211}
212
213//_____________________________________________________________________________
214const char*
215AliAnalysisMuMuFromESD::RunNumber(const AliESDEvent& esd) const
216{
217 /// Get the current runnumber
218 return Form("%09d",esd.GetRunNumber());
219}
220
221//_____________________________________________________________________________
222Bool_t AliAnalysisMuMuFromESD::TrackMatchCut(const AliESDMuonTrack& track) const
223{
224 /// whether the track match the trigger (any pt)
225 return track.GetMatchTrigger() > 0;
226}
227
228//_____________________________________________________________________________
229Bool_t AliAnalysisMuMuFromESD::TrackMatchLowCut(const AliESDMuonTrack& track) const
230{
231 /// whether the track match the trigger (low pt)
232 return track.GetMatchTrigger() > 1;
233}
234
235//_____________________________________________________________________________
236Bool_t AliAnalysisMuMuFromESD::TrackMatchHighCut(const AliESDMuonTrack& track) const
237{
238 /// whether the track match the trigger (high pt)
239 return track.GetMatchTrigger() > 2;
240}
241
242//_____________________________________________________________________________
243Bool_t AliAnalysisMuMuFromESD::TrackRabsCut(const AliESDMuonTrack& track) const
244{
245 /// Cut between 2 and 9 degrees
246
247 Double_t angle = TMath::ATan(track.GetRAtAbsorberEnd()/AbsZEnd());
248
249 return ( angle > Deg2() && angle < Deg10() );
250}
251
252//_____________________________________________________________________________
253Bool_t AliAnalysisMuMuFromESD::TrackPtCut(const AliESDMuonTrack& track) const
254{
255 /// whether the track pass the pt cut
256 return track.Pt() > 1.0;
257}
258
259//_____________________________________________________________________________
260Bool_t AliAnalysisMuMuFromESD::TrackChi2(const AliESDMuonTrack& track) const
261{
262 /// whether the track pass the chi2 cut
263 return track.GetNormalizedChi2() < 4.0;
264}
265
266//_____________________________________________________________________________
267Double_t AliAnalysisMuMuFromESD::CorrectedDCA(const AliESDMuonTrack& track) const
268{
269 /// Get corrected DCA of the track
270
271 TVector3 eventVertex(fVertex->GetX(),fVertex->GetY(),fVertex->GetZ());
272
273 TVector3 trackDcaAtVz(track.GetNonBendingCoorAtDCA(), track.GetBendingCoorAtDCA(), fVertex->GetZ());
274
275 TVector3 meanDca(-0.46, -0.92, 0.); // LHC10h1
276 TVector3 dcaAtVz = trackDcaAtVz - eventVertex - meanDca;
277 return dcaAtVz.Mag(); // it should also be equal to dcaAtVz.Pt().
278}
279
280//_____________________________________________________________________________
281Double_t AliAnalysisMuMuFromESD::PDCACutValue(const AliESDMuonTrack& track) const
282{
283 /// Get P x DCA cut value of the track
284
285 TVector3 pTot(track.Px(),track.Py(),track.Pz());
286
287 Double_t theta = TMath::ATan(track.GetRAtAbsorberEnd()/AbsZEnd());
288 Double_t cutValue = (theta > Deg3()) ? 63. : 120.;
289
290// AliInfo(Form("theta %e cutValue %e",theta*TMath::RadToDeg(),cutValue));
291
292 return 5*TMath::Sqrt(cutValue*cutValue + 0.4*0.4*pTot.Mag2());
293}
294
295//_____________________________________________________________________________
296Bool_t AliAnalysisMuMuFromESD::TrackDCACut(const AliESDMuonTrack& track) const
297{
298 /// whether the track pass the P x DCA cut
299
300 if (!fVertex) return kFALSE;
301
302 TVector3 pTot(track.Px(),track.Py(),track.Pz());
303
304 TVector3 pTotUncorr(track.PxUncorrected(), track.PyUncorrected(), track.PzUncorrected());
305
306 Double_t pTotMean = 0.5 * ( pTot.Mag() + pTotUncorr.Mag() );
307
308 Double_t correctedDca = CorrectedDCA(track);
309
310 Double_t cutVariable = pTotMean * correctedDca;
311
312 Double_t cutValue = PDCACutValue(track);
313
314 return ( cutVariable < cutValue );
315}
316
317//_____________________________________________________________________________
318Bool_t AliAnalysisMuMuFromESD::TrackEtaCut(const AliESDMuonTrack& track) const
319{
320 /// whether the track passes the eta cut
321 return track.Eta() < -2.5 && track.Eta() > -4;
322}
323
324//_____________________________________________________________________________
325UInt_t AliAnalysisMuMuFromESD::GetTrackMask(const AliESDMuonTrack& track) const
326{
327 /// Get the mask of all cuts this track pass
328
329 UInt_t m(kAll);
330
331 if ( TrackRabsCut(track) ) m |= kRabs;
332
333 if ( TrackPtCut(track) ) m |= kPt;
334
335 if ( TrackMatchCut(track) ) m |= kMatched;
336
337 if ( TrackMatchLowCut(track) ) m |= kMatchedLow;
338
339 if ( TrackMatchHighCut(track) ) m |= kMatchedHigh;
340
341 if ( TrackEtaCut(track) ) m |= kEta;
342
343 if ( TrackChi2(track) ) m |= kChi2;
344
345 if ( TrackDCACut(track) ) m |= kDCA;
346
347 return m;
348
349}
350
351//_____________________________________________________________________________
352void AliAnalysisMuMuFromESD::FillHistosForTrack(const char* physics,
353 const char* triggerClassName,
354 const char* centrality,
355 const AliESDMuonTrack& track,
356 const char* /*runNumber*/)
357{
358 /// Fill histograms for one track
359
360 if ( track.ContainTrackerData() )
361 {
362 TLorentzVector p;
363
364 track.LorentzP(p);
365
366 TString charge("Plus");
367
368 if ( track.Charge() < 0 ) charge = "Minus";
369
370 UInt_t mask = GetTrackMask(track);
371
372 Double_t theta = TMath::ATan(track.GetRAtAbsorberEnd()/AbsZEnd());
373
374 Double_t xdca = track.GetNonBendingCoorAtDCA();
375 Double_t ydca = track.GetBendingCoorAtDCA();
376 Double_t dca = TMath::Sqrt(xdca*xdca+ydca*ydca);
377
378 TIter next(fSingleTrackCutNames);
379 TObjString* str;
380
381 while ( ( str = static_cast<TObjString*>(next()) ) )
382 {
383 Bool_t test = ( ( str->GetUniqueID() & mask ) == str->GetUniqueID() );
384
385 if ( test )
386 {
387
388 TVector3 pTot(track.Px(),track.Py(),track.Pz());
389
390 TVector3 pTotUncorr(track.PxUncorrected(), track.PyUncorrected(), track.PzUncorrected());
391
392 Double_t pTotMean = 0.5 * ( pTot.Mag() + pTotUncorr.Mag() );
393
394 Double_t pdcacorr = pTotMean*CorrectedDCA(track);
395 Double_t pdcacut = PDCACutValue(track);
396
397 Histo(physics,triggerClassName,centrality,str->String().Data(),Form("PtEtaMu%s",charge.Data()))->Fill(p.Eta(),p.Pt());
398
399 Histo(physics,triggerClassName,centrality,str->String().Data(),Form("PEtaMu%s",charge.Data()))->Fill(p.Eta(),p.P());
400
401
402 Histo(physics,triggerClassName,centrality,str->String().Data(),Form("XYdcaMu%s",charge.Data()))->Fill(ydca,xdca);
403
404 Histo(physics,triggerClassName,centrality,str->String().Data(),Form("Chi2Mu%s",charge.Data()))->Fill(track.GetNormalizedChi2());
405
406 if ( theta >= Deg2() && theta < Deg3() )
407 {
408 // Histo(physics,triggerClassName,centrality,str->String().Data(),Form("PDCAP23Mu%s",charge.Data()))->Fill(p.P(),p.P()*dca);
409 Histo(physics,triggerClassName,centrality,str->String().Data(),Form("dcaP23Mu%s",charge.Data()))->Fill(p.P(),dca);
410
411 Histo(physics,triggerClassName,centrality,str->String().Data(),Form("PDCAcorrP23Mu%s",charge.Data()))->Fill(pTotMean,pdcacorr);
412 Histo(physics,triggerClassName,centrality,str->String().Data(),Form("PDCAcutP23Mu%s",charge.Data()))->Fill(p.P(),pdcacut);
413
414 if ( p.Pt() > 2 )
415 {
416 Histo(physics,triggerClassName,centrality,str->String().Data(),Form("dcaPwPtCut23Mu%s",charge.Data()))->Fill(p.P(),dca);
417 }
418 }
419 else if ( theta >= Deg3() && theta < Deg10() )
420 {
421 // Histo(physics,triggerClassName,centrality,str->String().Data(),Form("PDCAP310Mu%s",charge.Data()))->Fill(p.P(),p.P()*dca);
422 Histo(physics,triggerClassName,centrality,str->String().Data(),Form("dcaP310Mu%s",charge.Data()))->Fill(p.P(),dca);
423 Histo(physics,triggerClassName,centrality,str->String().Data(),Form("PDCAcorrP310Mu%s",charge.Data()))->Fill(pTotMean,pdcacorr);
424 Histo(physics,triggerClassName,centrality,str->String().Data(),Form("PDCAcutP310Mu%s",charge.Data()))->Fill(p.P(),pdcacut);
425 if ( p.Pt() > 2 )
426 {
427 Histo(physics,triggerClassName,centrality,str->String().Data(),Form("dcaPwPtCut310Mu%s",charge.Data()))->Fill(p.P(),dca);
428 }
429 }
430 }
431 }
432 }
433}
434
435//_____________________________________________________________________________
436void AliAnalysisMuMuFromESD::FillHistos(const char* physics,
437 const char* triggerClassName,
438 const char* centrality,
439 const AliESDEvent& esd)
440{
441 /// Fill histograms for /physics/triggerClassName/centrality
442
443 Histo(physics,triggerClassName,centrality,"BCX")->Fill(1.0*esd.GetBunchCrossNumber());
444
445 if ( fVertex )
446 {
447 Histo(physics,triggerClassName,centrality,"Xvertex")->Fill(fVertex->GetX());
448 Histo(physics,triggerClassName,centrality,"Yvertex")->Fill(fVertex->GetY());
449 Histo(physics,triggerClassName,centrality,"Zvertex")->Fill(fVertex->GetZ());
450 }
451
452 AliESDVZERO* vzero = esd.GetVZEROData();
453
454 if ( vzero )
455 {
456 Histo(physics,triggerClassName,centrality,"V0Time")->Fill(vzero->GetV0ATime(),vzero->GetV0CTime());
457
458 Double_t v0mult = vzero->GetMTotV0A()+vzero->GetMTotV0C();
459
460 Histo(physics,triggerClassName,centrality,"V0Amplitude")->Fill(v0mult);
461 }
462
463 // Track loop
464
465 TString runNumber(RunNumber(esd));
466
467 for (Int_t i = 0; i < esd.GetNumberOfMuonTracks(); ++i)
468 {
469 AliESDMuonTrack* tracki = esd.GetMuonTrack(i);
470
471 if ( tracki->ContainTrackerData() )
472 {
473 UInt_t maski = GetTrackMask(*tracki);
474
475 FillHistosForTrack(physics,triggerClassName,centrality,*tracki,runNumber.Data());
476
477 TLorentzVector pi;
478
479 tracki->LorentzP(pi);
480
481 for (Int_t j = i+1; j < esd.GetNumberOfMuonTracks(); ++j)
482 {
483 AliESDMuonTrack* trackj = esd.GetMuonTrack(j);
484
485 if ( trackj->ContainTrackerData() )
486 {
487 UInt_t maskj = GetTrackMask(*trackj);
488
489 TLorentzVector pj;
490 trackj->LorentzP(pj);
491
492 pj += pi;
493
494 TIter next(fPairTrackCutNames);
495 TObjString* str;
496
497 while ( ( str = static_cast<TObjString*>(next()) ) )
498 {
499 UInt_t singleTrackMask(0);
500 UInt_t pairMask(0);
501
502 DecodePairCutMask(str->GetUniqueID(),singleTrackMask,pairMask);
503
504 Bool_t testi = ( ( maski & singleTrackMask ) == singleTrackMask ) ;
505 Bool_t testj = ( ( maskj & singleTrackMask ) == singleTrackMask ) ;
506 Bool_t testij(kTRUE);
507
508 if ( pairMask > 0 )
509 {
510 testij = ( (maski & pairMask) == pairMask ) &&
511 ( (maskj & pairMask) == pairMask );
512 }
513
514 if ( ( testi || testj ) && testij )
515 {
516 if ( tracki->Charge() != trackj->Charge() )
517 {
518 Histo(physics,triggerClassName,centrality,str->String().Data(),"MinvUSPt")->Fill(pj.Pt(),pj.M());
519 }
520 else if ( tracki->Charge() > 0 && trackj->Charge() > 0 )
521 {
522 Histo(physics,triggerClassName,centrality,str->String().Data(),"MinvPPPt")->Fill(pj.Pt(),pj.M());
523 }
524 else
525 {
526 Histo(physics,triggerClassName,centrality,str->String().Data(),"MinvMMPt")->Fill(pj.Pt(),pj.M());
527 }
528 }
529 }
530 }
531 }
532 }
533 } //track loop
534}
535