1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 #include "AliAnalysisMuMuFromESD.h"
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"
30 #include "TDirectory.h"
33 #include "TLorentzVector.h"
36 #include "TObjArray.h"
37 #include "TPaveText.h"
39 #include "AliPhysicsSelection.h"
40 #include "AliPhysicsSelectionTask.h"
41 #include "AliBackgroundSelection.h"
42 #include "AliCounterCollection.h"
43 #include "THashList.h"
48 /// AliAnalysisMuMuFromESD
50 /// cuts must be added through AddSingleCut and AddPaircut methods
54 ClassImp(AliAnalysisMuMuFromESD)
56 //_____________________________________________________________________________
57 AliAnalysisMuMuFromESD::AliAnalysisMuMuFromESD()
58 : AliAnalysisMuMu(), fVertex(0x0)
63 //_____________________________________________________________________________
64 AliAnalysisMuMuFromESD::AliAnalysisMuMuFromESD(Bool_t aa)
65 :AliAnalysisMuMu(kTRUE,aa), fVertex(0x0)
71 //_____________________________________________________________________________
72 AliAnalysisMuMuFromESD::AliAnalysisMuMuFromESD(TList* triggerClasses)
73 :AliAnalysisMuMu(kTRUE,triggerClasses), fVertex(0x0)
79 //_____________________________________________________________________________
80 void AliAnalysisMuMuFromESD::Ctor()
86 AddGlobalEventSelection("SD2");
90 AddGlobalEventSelection("ALL");
93 AddGlobalEventSelection("PS");
95 // fBranchNames = "ESD:AliESDRun.,AliESDHeader.,MuonTracks";
98 //_____________________________________________________________________________
99 AliAnalysisMuMuFromESD::~AliAnalysisMuMuFromESD()
105 //_____________________________________________________________________________
106 void AliAnalysisMuMuFromESD::MuUserExec(Option_t*)
109 // Called for each event
114 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
121 // consider only events with OSM2 fired
123 UInt_t sd2 = (1<<12);
124 UInt_t trigger = esd->GetHeader()->GetL0TriggerInputs();
126 Bool_t ok = ( ( trigger & sd2 ) == sd2);
133 TString runNumber(RunNumber(*esd));
135 TString triggers = esd->GetFiredTriggerClasses();
137 if ( IsDynamicTriggerClasses() ) AddTriggerClasses(triggers.Data());
139 TObjArray* a = triggers.Tokenize(" ");
143 AliInputEventHandler* inputEventHandler = static_cast<AliInputEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
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 );
151 TString centrality(DefaultCentralityName());
153 fVertex = new AliESDVertex(*esd->GetPrimaryVertexSPD());
157 // FIXME : should get it from centrality task...
159 // Double_t cent(-100);
160 // Double_t v0mult = vzero->GetMTotV0A()+vzero->GetMTotV0C();
162 // for ( std::vector<double>::size_type i = 0 ; i < fCentralityLimits.size()-1 && cent < -10 ; ++i )
164 // if ( v0mult > fCentralityAmplitude[i] )
166 // cent = fCentralityLimits[i];
172 // centrality = CentralityName(cent);
176 TString eventtype("all");
179 if ( fAA ) eventtype="sd2";
183 while ( ( tname = static_cast<TObjString*>(next()) ) )
185 if ( !fTriggerClasses->FindObject(tname->GetName()) )
187 // skip the trigger classes we're not supposed to analyse
192 fEventCounters->Count(Form("event:%s/trigger:%s/run:%d", eventtype.Data(),tname->GetName(), runNumber.Atoi()));
194 AssertHistogramCollection(et.Data(),tname->String().Data());
196 FillHistos(et.Data(),tname->String().Data(),centrality.Data(),*esd);
198 if ( isPhysicsSelected )
201 fEventCounters->Count(Form("event:ps/trigger:%s/run:%d", tname->GetName(), runNumber.Atoi()));
203 AssertHistogramCollection("PS",tname->String().Data());
205 FillHistos("PS",tname->String().Data(),centrality.Data(),*esd);
213 //_____________________________________________________________________________
215 AliAnalysisMuMuFromESD::RunNumber(const AliESDEvent& esd) const
217 /// Get the current runnumber
218 return Form("%09d",esd.GetRunNumber());
221 //_____________________________________________________________________________
222 Bool_t AliAnalysisMuMuFromESD::TrackMatchCut(const AliESDMuonTrack& track) const
224 /// whether the track match the trigger (any pt)
225 return track.GetMatchTrigger() > 0;
228 //_____________________________________________________________________________
229 Bool_t AliAnalysisMuMuFromESD::TrackMatchLowCut(const AliESDMuonTrack& track) const
231 /// whether the track match the trigger (low pt)
232 return track.GetMatchTrigger() > 1;
235 //_____________________________________________________________________________
236 Bool_t AliAnalysisMuMuFromESD::TrackMatchHighCut(const AliESDMuonTrack& track) const
238 /// whether the track match the trigger (high pt)
239 return track.GetMatchTrigger() > 2;
242 //_____________________________________________________________________________
243 Bool_t AliAnalysisMuMuFromESD::TrackRabsCut(const AliESDMuonTrack& track) const
245 /// Cut between 2 and 9 degrees
247 Double_t angle = TMath::ATan(track.GetRAtAbsorberEnd()/AbsZEnd());
249 return ( angle > Deg2() && angle < Deg10() );
252 //_____________________________________________________________________________
253 Bool_t AliAnalysisMuMuFromESD::TrackPtCut(const AliESDMuonTrack& track) const
255 /// whether the track pass the pt cut
256 return track.Pt() > 1.0;
259 //_____________________________________________________________________________
260 Bool_t AliAnalysisMuMuFromESD::TrackChi2(const AliESDMuonTrack& track) const
262 /// whether the track pass the chi2 cut
263 return track.GetNormalizedChi2() < 4.0;
266 //_____________________________________________________________________________
267 Double_t AliAnalysisMuMuFromESD::CorrectedDCA(const AliESDMuonTrack& track) const
269 /// Get corrected DCA of the track
271 TVector3 eventVertex(fVertex->GetX(),fVertex->GetY(),fVertex->GetZ());
273 TVector3 trackDcaAtVz(track.GetNonBendingCoorAtDCA(), track.GetBendingCoorAtDCA(), fVertex->GetZ());
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().
280 //_____________________________________________________________________________
281 Double_t AliAnalysisMuMuFromESD::PDCACutValue(const AliESDMuonTrack& track) const
283 /// Get P x DCA cut value of the track
285 TVector3 pTot(track.Px(),track.Py(),track.Pz());
287 Double_t theta = TMath::ATan(track.GetRAtAbsorberEnd()/AbsZEnd());
288 Double_t cutValue = (theta > Deg3()) ? 63. : 120.;
290 // AliInfo(Form("theta %e cutValue %e",theta*TMath::RadToDeg(),cutValue));
292 return 5*TMath::Sqrt(cutValue*cutValue + 0.4*0.4*pTot.Mag2());
295 //_____________________________________________________________________________
296 Bool_t AliAnalysisMuMuFromESD::TrackDCACut(const AliESDMuonTrack& track) const
298 /// whether the track pass the P x DCA cut
300 if (!fVertex) return kFALSE;
302 TVector3 pTot(track.Px(),track.Py(),track.Pz());
304 TVector3 pTotUncorr(track.PxUncorrected(), track.PyUncorrected(), track.PzUncorrected());
306 Double_t pTotMean = 0.5 * ( pTot.Mag() + pTotUncorr.Mag() );
308 Double_t correctedDca = CorrectedDCA(track);
310 Double_t cutVariable = pTotMean * correctedDca;
312 Double_t cutValue = PDCACutValue(track);
314 return ( cutVariable < cutValue );
317 //_____________________________________________________________________________
318 Bool_t AliAnalysisMuMuFromESD::TrackEtaCut(const AliESDMuonTrack& track) const
320 /// whether the track passes the eta cut
321 return track.Eta() < -2.5 && track.Eta() > -4;
324 //_____________________________________________________________________________
325 UInt_t AliAnalysisMuMuFromESD::GetTrackMask(const AliESDMuonTrack& track) const
327 /// Get the mask of all cuts this track pass
331 if ( TrackRabsCut(track) ) m |= kRabs;
333 if ( TrackPtCut(track) ) m |= kPt;
335 if ( TrackMatchCut(track) ) m |= kMatched;
337 if ( TrackMatchLowCut(track) ) m |= kMatchedLow;
339 if ( TrackMatchHighCut(track) ) m |= kMatchedHigh;
341 if ( TrackEtaCut(track) ) m |= kEta;
343 if ( TrackChi2(track) ) m |= kChi2;
345 if ( TrackDCACut(track) ) m |= kDCA;
351 //_____________________________________________________________________________
352 void AliAnalysisMuMuFromESD::FillHistosForTrack(const char* physics,
353 const char* triggerClassName,
354 const char* centrality,
355 const AliESDMuonTrack& track,
356 const char* /*runNumber*/)
358 /// Fill histograms for one track
360 if ( track.ContainTrackerData() )
366 TString charge("Plus");
368 if ( track.Charge() < 0 ) charge = "Minus";
370 UInt_t mask = GetTrackMask(track);
372 Double_t theta = TMath::ATan(track.GetRAtAbsorberEnd()/AbsZEnd());
374 Double_t xdca = track.GetNonBendingCoorAtDCA();
375 Double_t ydca = track.GetBendingCoorAtDCA();
376 Double_t dca = TMath::Sqrt(xdca*xdca+ydca*ydca);
378 TIter next(fSingleTrackCutNames);
381 while ( ( str = static_cast<TObjString*>(next()) ) )
383 Bool_t test = ( ( str->GetUniqueID() & mask ) == str->GetUniqueID() );
388 TVector3 pTot(track.Px(),track.Py(),track.Pz());
390 TVector3 pTotUncorr(track.PxUncorrected(), track.PyUncorrected(), track.PzUncorrected());
392 Double_t pTotMean = 0.5 * ( pTot.Mag() + pTotUncorr.Mag() );
394 Double_t pdcacorr = pTotMean*CorrectedDCA(track);
395 Double_t pdcacut = PDCACutValue(track);
397 Histo(physics,triggerClassName,centrality,str->String().Data(),Form("PtEtaMu%s",charge.Data()))->Fill(p.Eta(),p.Pt());
399 Histo(physics,triggerClassName,centrality,str->String().Data(),Form("PEtaMu%s",charge.Data()))->Fill(p.Eta(),p.P());
402 Histo(physics,triggerClassName,centrality,str->String().Data(),Form("XYdcaMu%s",charge.Data()))->Fill(ydca,xdca);
404 Histo(physics,triggerClassName,centrality,str->String().Data(),Form("Chi2Mu%s",charge.Data()))->Fill(track.GetNormalizedChi2());
406 if ( theta >= Deg2() && theta < Deg3() )
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);
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);
416 Histo(physics,triggerClassName,centrality,str->String().Data(),Form("dcaPwPtCut23Mu%s",charge.Data()))->Fill(p.P(),dca);
419 else if ( theta >= Deg3() && theta < Deg10() )
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);
427 Histo(physics,triggerClassName,centrality,str->String().Data(),Form("dcaPwPtCut310Mu%s",charge.Data()))->Fill(p.P(),dca);
435 //_____________________________________________________________________________
436 void AliAnalysisMuMuFromESD::FillHistos(const char* physics,
437 const char* triggerClassName,
438 const char* centrality,
439 const AliESDEvent& esd)
441 /// Fill histograms for /physics/triggerClassName/centrality
443 Histo(physics,triggerClassName,centrality,"BCX")->Fill(1.0*esd.GetBunchCrossNumber());
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());
452 AliESDVZERO* vzero = esd.GetVZEROData();
456 Histo(physics,triggerClassName,centrality,"V0Time")->Fill(vzero->GetV0ATime(),vzero->GetV0CTime());
458 Double_t v0mult = vzero->GetMTotV0A()+vzero->GetMTotV0C();
460 Histo(physics,triggerClassName,centrality,"V0Amplitude")->Fill(v0mult);
465 TString runNumber(RunNumber(esd));
467 for (Int_t i = 0; i < esd.GetNumberOfMuonTracks(); ++i)
469 AliESDMuonTrack* tracki = esd.GetMuonTrack(i);
471 if ( tracki->ContainTrackerData() )
473 UInt_t maski = GetTrackMask(*tracki);
475 FillHistosForTrack(physics,triggerClassName,centrality,*tracki,runNumber.Data());
479 tracki->LorentzP(pi);
481 for (Int_t j = i+1; j < esd.GetNumberOfMuonTracks(); ++j)
483 AliESDMuonTrack* trackj = esd.GetMuonTrack(j);
485 if ( trackj->ContainTrackerData() )
487 UInt_t maskj = GetTrackMask(*trackj);
490 trackj->LorentzP(pj);
494 TIter next(fPairTrackCutNames);
497 while ( ( str = static_cast<TObjString*>(next()) ) )
499 UInt_t singleTrackMask(0);
502 DecodePairCutMask(str->GetUniqueID(),singleTrackMask,pairMask);
504 Bool_t testi = ( ( maski & singleTrackMask ) == singleTrackMask ) ;
505 Bool_t testj = ( ( maskj & singleTrackMask ) == singleTrackMask ) ;
506 Bool_t testij(kTRUE);
510 testij = ( (maski & pairMask) == pairMask ) &&
511 ( (maskj & pairMask) == pairMask );
514 if ( ( testi || testj ) && testij )
516 if ( tracki->Charge() != trackj->Charge() )
518 Histo(physics,triggerClassName,centrality,str->String().Data(),"MinvUSPt")->Fill(pj.Pt(),pj.M());
520 else if ( tracki->Charge() > 0 && trackj->Charge() > 0 )
522 Histo(physics,triggerClassName,centrality,str->String().Data(),"MinvPPPt")->Fill(pj.Pt(),pj.M());
526 Histo(physics,triggerClassName,centrality,str->String().Data(),"MinvMMPt")->Fill(pj.Pt(),pj.M());