]>
Commit | Line | Data |
---|---|---|
aad6618e | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-2007, 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 | ||
27de2dfb | 16 | /* $Id$ */ |
17 | ||
662e37fe | 18 | //----------------------------------------------------------------------------- |
19 | /// \class AliAnalysisTaskSingleMu | |
20 | /// Analysis task for single muons in the spectrometer. | |
a07db2fb | 21 | /// The output is a list of histograms and CF containers. |
b201705a | 22 | /// The macro class can run on AODs or ESDs. |
9728bcfd | 23 | /// If Monte Carlo information is present, some basics checks are performed. |
662e37fe | 24 | /// |
25 | /// \author Diego Stocco | |
26 | //----------------------------------------------------------------------------- | |
27 | ||
aad6618e | 28 | #define AliAnalysisTaskSingleMu_cxx |
29 | ||
a07db2fb | 30 | #include "AliAnalysisTaskSingleMu.h" |
31 | ||
aad6618e | 32 | // ROOT includes |
aad6618e | 33 | #include "TROOT.h" |
9728bcfd | 34 | #include "TH1.h" |
35 | #include "TH2.h" | |
9728bcfd | 36 | #include "TAxis.h" |
662e37fe | 37 | #include "TCanvas.h" |
a07db2fb | 38 | #include "TLegend.h" |
9728bcfd | 39 | #include "TMath.h" |
589e3f71 | 40 | #include "TObjString.h" |
55065f3f | 41 | #include "TObjArray.h" |
a07db2fb | 42 | #include "TF1.h" |
43 | #include "TStyle.h" | |
44 | //#include "TMCProcess.h" | |
a07db2fb | 45 | #include "TArrayI.h" |
261f31ad | 46 | #include "TArrayD.h" |
c6e0141f | 47 | #include "TPaveStats.h" |
12e33589 | 48 | #include "TFitResultPtr.h" |
ac2e110b | 49 | #include "TFile.h" |
50 | #include "THashList.h" | |
aad6618e | 51 | |
52 | // STEER includes | |
aad6618e | 53 | #include "AliAODEvent.h" |
54 | #include "AliAODTrack.h" | |
a07db2fb | 55 | #include "AliAODMCParticle.h" |
75008b31 | 56 | #include "AliMCEvent.h" |
9728bcfd | 57 | #include "AliMCParticle.h" |
b201705a | 58 | #include "AliESDEvent.h" |
59 | #include "AliESDMuonTrack.h" | |
a07db2fb | 60 | #include "AliVHeader.h" |
61 | #include "AliAODMCHeader.h" | |
62 | #include "AliStack.h" | |
9728bcfd | 63 | |
589e3f71 | 64 | // ANALYSIS includes |
65 | #include "AliAnalysisManager.h" | |
589e3f71 | 66 | |
67 | // CORRFW includes | |
589e3f71 | 68 | #include "AliCFContainer.h" |
69 | #include "AliCFGridSparse.h" | |
c6e0141f | 70 | #include "AliCFEffGrid.h" |
589e3f71 | 71 | |
c6e0141f | 72 | // PWG includes |
a07db2fb | 73 | #include "AliVAnalysisMuon.h" |
74 | #include "AliMergeableCollection.h" | |
75 | #include "AliCounterCollection.h" | |
ebba9ef0 | 76 | #include "AliMuonEventCuts.h" |
c6e0141f | 77 | #include "AliMuonTrackCuts.h" |
ebba9ef0 | 78 | #include "AliAnalysisMuonUtility.h" |
a07db2fb | 79 | |
aad6618e | 80 | |
662e37fe | 81 | /// \cond CLASSIMP |
82 | ClassImp(AliAnalysisTaskSingleMu) // Class implementation in ROOT context | |
83 | /// \endcond | |
aad6618e | 84 | |
9728bcfd | 85 | |
aad6618e | 86 | //________________________________________________________________________ |
a07db2fb | 87 | AliAnalysisTaskSingleMu::AliAnalysisTaskSingleMu() : |
88 | AliVAnalysisMuon(), | |
47b31d8b CP |
89 | fThetaAbsKeys(0x0), |
90 | fCutOnDimu(kFALSE) | |
aad6618e | 91 | { |
a07db2fb | 92 | /// Default ctor. |
aad6618e | 93 | } |
94 | ||
9728bcfd | 95 | //________________________________________________________________________ |
a07db2fb | 96 | AliAnalysisTaskSingleMu::AliAnalysisTaskSingleMu(const char *name, const AliMuonTrackCuts& cuts) : |
97 | AliVAnalysisMuon(name, cuts), | |
47b31d8b CP |
98 | fThetaAbsKeys(0x0), |
99 | fCutOnDimu(kFALSE) | |
aad6618e | 100 | { |
101 | // | |
a07db2fb | 102 | /// Constructor. |
9728bcfd | 103 | // |
a07db2fb | 104 | TString thetaAbsKeys = "ThetaAbs23 ThetaAbs310"; |
105 | fThetaAbsKeys = thetaAbsKeys.Tokenize(" "); | |
589e3f71 | 106 | } |
107 | ||
108 | ||
a07db2fb | 109 | //________________________________________________________________________ |
110 | AliAnalysisTaskSingleMu::~AliAnalysisTaskSingleMu() | |
589e3f71 | 111 | { |
112 | // | |
a07db2fb | 113 | /// Destructor |
589e3f71 | 114 | // |
115 | ||
a07db2fb | 116 | delete fThetaAbsKeys; |
aad6618e | 117 | } |
118 | ||
119 | //___________________________________________________________________________ | |
a07db2fb | 120 | void AliAnalysisTaskSingleMu::MyUserCreateOutputObjects() |
aad6618e | 121 | { |
589e3f71 | 122 | |
a07db2fb | 123 | TH1* histo = 0x0; |
124 | TString histoName = "", histoTitle = ""; | |
9728bcfd | 125 | |
589e3f71 | 126 | Int_t nVzBins = 40; |
127 | Double_t vzMin = -20., vzMax = 20.; | |
a07db2fb | 128 | TString vzName("Vz"), vzTitle("Vz"), vzUnits("cm"); |
9728bcfd | 129 | |
a07db2fb | 130 | histoName = "hIpVtx"; |
131 | histo = new TH1F(histoName.Data(), histoName.Data(), nVzBins, vzMin, vzMax); | |
132 | histo->SetXTitle("v_{z} (cm)"); | |
133 | AddObjectToCollection(histo, kIPVz); | |
589e3f71 | 134 | |
a07db2fb | 135 | |
136 | Int_t nPtBins = 80; | |
137 | Double_t ptMin = 0., ptMax = 80.; | |
171652d3 | 138 | TString ptName("Pt"), ptTitle("p_{T}"), ptUnits("GeV/c"); |
a07db2fb | 139 | |
140 | Int_t nEtaBins = 25; | |
141 | Double_t etaMin = -4.5, etaMax = -2.; | |
142 | TString etaName("Eta"), etaTitle("#eta"), etaUnits(""); | |
143 | ||
144 | Int_t nPhiBins = 36; | |
12e33589 | 145 | Double_t phiMin = 0.; Double_t phiMax = 2.*TMath::Pi(); |
a07db2fb | 146 | TString phiName("Phi"), phiTitle("#phi"), phiUnits("rad"); |
147 | ||
589e3f71 | 148 | Int_t nChargeBins = 2; |
149 | Double_t chargeMin = -2, chargeMax = 2.; | |
150 | TString chargeName("Charge"), chargeTitle("charge"), chargeUnits("e"); | |
9728bcfd | 151 | |
a07db2fb | 152 | Int_t nThetaAbsEndBins = 2; |
153 | Double_t thetaAbsEndMin = -0.5, thetaAbsEndMax = 1.5; | |
154 | TString thetaAbsEndName("ThetaAbsEnd"), thetaAbsEndTitle("#theta_{abs}"), thetaAbsEndUnits("a.u."); | |
155 | ||
589e3f71 | 156 | Int_t nMotherTypeBins = kNtrackSources; |
157 | Double_t motherTypeMin = -0.5, motherTypeMax = (Double_t)kNtrackSources - 0.5; | |
158 | TString motherType("MotherType"), motherTypeTitle("motherType"), motherTypeUnits(""); | |
a07db2fb | 159 | |
160 | Int_t nbins[kNvars] = {nPtBins, nEtaBins, nPhiBins, nVzBins, nChargeBins, nThetaAbsEndBins, nMotherTypeBins}; | |
161 | Double_t xmin[kNvars] = {ptMin, etaMin, phiMin, vzMin, chargeMin, thetaAbsEndMin, motherTypeMin}; | |
162 | Double_t xmax[kNvars] = {ptMax, etaMax, phiMax, vzMax, chargeMax, thetaAbsEndMax, motherTypeMax}; | |
163 | TString axisTitle[kNvars] = {ptTitle, etaTitle, phiTitle, vzTitle, chargeTitle, thetaAbsEndTitle, motherTypeTitle}; | |
164 | TString axisUnits[kNvars] = {ptUnits, etaUnits, phiUnits, vzUnits, chargeUnits, thetaAbsEndUnits, motherTypeUnits}; | |
9728bcfd | 165 | |
a07db2fb | 166 | AliCFContainer* cfContainer = new AliCFContainer("SingleMuContainer","Container for tracks",kNsteps,kNvars,nbins); |
4ab8d5a6 | 167 | |
589e3f71 | 168 | for ( Int_t idim = 0; idim<kNvars; idim++){ |
169 | histoTitle = Form("%s (%s)", axisTitle[idim].Data(), axisUnits[idim].Data()); | |
170 | histoTitle.ReplaceAll("()",""); | |
a07db2fb | 171 | |
172 | cfContainer->SetVarTitle(idim, histoTitle.Data()); | |
173 | cfContainer->SetBinLimits(idim, xmin[idim], xmax[idim]); | |
9728bcfd | 174 | } |
a07db2fb | 175 | |
176 | TString stepTitle[kNsteps] = {"reconstructed", "generated"}; | |
589e3f71 | 177 | |
a07db2fb | 178 | TAxis* currAxis = 0x0; |
589e3f71 | 179 | for (Int_t istep=0; istep<kNsteps; istep++){ |
a07db2fb | 180 | cfContainer->SetStepTitle(istep, stepTitle[istep].Data()); |
181 | AliCFGridSparse* gridSparse = cfContainer->GetGrid(istep); | |
182 | ||
183 | currAxis = gridSparse->GetAxis(kHvarMotherType); | |
184 | for ( Int_t ibin=0; ibin<fSrcKeys->GetEntries(); ibin++ ) { | |
185 | currAxis->SetBinLabel(ibin+1, fSrcKeys->At(ibin)->GetName()); | |
9728bcfd | 186 | } |
187 | } | |
4ab8d5a6 | 188 | |
a07db2fb | 189 | AddObjectToCollection(cfContainer, kTrackContainer); |
190 | ||
47b31d8b CP |
191 | histoName = "hRecoDimu"; |
192 | TH2* histoDimu = new TH2F(histoName.Data(), histoName.Data(), 24, 0., 120., 30, 0., 150.); | |
193 | histoDimu->SetXTitle("min. p_{T}^{#mu} (GeV/c)"); | |
194 | histoDimu->SetYTitle("M_{#mu#mu} (GeV/c^{2})"); | |
195 | AddObjectToCollection(histoDimu, kNobjectTypes); | |
196 | ||
197 | histoName = "hGeneratedZ"; | |
3575043a | 198 | TH2* histoGenZ = static_cast<TH2*>(histoDimu->Clone(histoName.Data())); |
199 | histoGenZ->SetTitle(histoName.Data()); | |
200 | AddObjectToCollection(histoGenZ, kNobjectTypes+1); | |
47b31d8b CP |
201 | |
202 | ||
203 | histoName = "hZmuEtaCorr"; | |
3575043a | 204 | histoDimu = new TH2F(histoName.Data(), histoName.Data(), 160, -8., 8., 160, -8., 8.); |
47b31d8b CP |
205 | histoDimu->SetXTitle("#eta_{#mu-}"); |
206 | histoDimu->SetYTitle("#eta_{#mu+}"); | |
207 | AddObjectToCollection(histoDimu, kNobjectTypes+2); | |
208 | ||
c6e0141f | 209 | fMuonTrackCuts->Print("mask"); |
56e01f1b | 210 | |
211 | AliInfo(Form("Apply cut on dimuon (60<M_mumu<120 GeV/c^2) to reject Z contribution: %i", fCutOnDimu)); | |
aad6618e | 212 | } |
213 | ||
214 | //________________________________________________________________________ | |
a07db2fb | 215 | void AliAnalysisTaskSingleMu::ProcessEvent(TString physSel, const TObjArray& selectTrigClasses, TString centrality) |
aad6618e | 216 | { |
217 | // | |
a07db2fb | 218 | /// Fill output objects |
aad6618e | 219 | // |
220 | ||
ebba9ef0 | 221 | Double_t ipVz = AliAnalysisMuonUtility::GetVertexSPD(InputEvent())->GetZ(); |
ac2e110b | 222 | Double_t ipVzMC = MCEvent() ? AliAnalysisMuonUtility::GetMCVertexZ(InputEvent(),MCEvent()) : 0.; |
a07db2fb | 223 | |
224 | for ( Int_t itrig=0; itrig<selectTrigClasses.GetEntries(); ++itrig ) { | |
225 | TString trigClassName = ((TObjString*)selectTrigClasses.At(itrig))->GetString(); | |
226 | ((TH1*)GetMergeableObject(physSel, trigClassName, centrality, "hIpVtx"))->Fill(ipVz); | |
9728bcfd | 227 | } |
5936324e | 228 | |
d5a197f5 | 229 | // Bool_t isPileupFromSPD = ( fAODEvent && ! fAODEvent->GetTracklets() ) ? InputEvent()->IsPileupFromSPD(3, 0.8, 3., 2., 5.) : InputEvent()->IsPileupFromSPDInMultBins(); // Avoid break when reading Muon AODs (tracklet info is not present and IsPileupFromSPDInMultBins crashes // UNCOMMENT TO REJECT PILEUP |
230 | // if ( isPileupFromSPD ) return; // UNCOMMENT TO REJECT PILEUP | |
a07db2fb | 231 | |
232 | Double_t containerInput[kNvars]; | |
233 | AliVParticle* track = 0x0; | |
93d6def1 | 234 | |
ac2e110b | 235 | Int_t nSteps = MCEvent() ? 2 : 1; |
236 | for ( Int_t istep = 0; istep<nSteps; ++istep ) { | |
237 | Int_t nTracks = ( istep == kStepReconstructed ) ? AliAnalysisMuonUtility::GetNTracks(InputEvent()) : MCEvent()->GetNumberOfTracks(); | |
47b31d8b CP |
238 | |
239 | TObjArray selectedTracks(nTracks); | |
240 | TArrayI trackSources(nTracks); | |
261f31ad | 241 | TArrayD trackWgt(nTracks); |
5936324e | 242 | trackWgt.Reset(1.); |
a07db2fb | 243 | for (Int_t itrack = 0; itrack < nTracks; itrack++) { |
ac2e110b | 244 | track = ( istep == kStepReconstructed ) ? AliAnalysisMuonUtility::GetTrack(itrack,InputEvent()) : MCEvent()->GetTrack(itrack); |
a07db2fb | 245 | |
47b31d8b CP |
246 | // In case of MC we usually ask that the particle is a muon |
247 | // However, in W or Z simulations, Pythia stores both the initial muon | |
248 | // (before ISR, FSR and kt kick) and the final state one. | |
249 | // The first muon is of course there only for information and should be rejected. | |
250 | // The Pythia code for initial state particles is 21 | |
f7852207 | 251 | // When running with POWHEG, Pythia puts the hard process input of POWHEG in the stack |
252 | // with state 21, and then re-add it to stack before applying ISR, FSR and kt kick. | |
253 | // This muon produces the final state muon, and its status code is 11 | |
254 | // To avoid all problems, keep only final state muon (status code <10) | |
255 | // FIXME: is the convention valid for other generators as well? | |
256 | Bool_t isSelected = ( istep == kStepReconstructed ) ? fMuonTrackCuts->IsSelected(track) : ( TMath::Abs(track->PdgCode()) == 13 && AliAnalysisMuonUtility::GetStatusCode(track) < 10 ); | |
a07db2fb | 257 | if ( ! isSelected ) continue; |
258 | ||
47b31d8b CP |
259 | selectedTracks.AddAt(track,itrack); |
260 | trackSources[itrack] = GetParticleType(track); | |
261 | ||
261f31ad | 262 | TObject* wgtObj = GetWeight(fSrcKeys->At(trackSources[itrack])->GetName()); |
263 | ||
264 | if ( wgtObj ) { | |
265 | AliVParticle* mcTrack = ( istep == kStepReconstructed ) ? MCEvent()->GetTrack(track->GetLabel()) : track; | |
266 | if ( wgtObj->IsA() == TF1::Class() ) trackWgt[itrack] = static_cast<TF1*>(wgtObj)->Eval(mcTrack->Pt()); | |
267 | else if ( wgtObj->IsA()->InheritsFrom(TH1::Class()) ) { | |
268 | TH1* wgtHisto = static_cast<TH1*>(wgtObj); | |
269 | trackWgt[itrack] = wgtHisto->GetBinContent(wgtHisto->GetXaxis()->FindBin(mcTrack->Pt())); | |
270 | } | |
271 | AliDebug(3,Form("Apply weights %s: pt %g gen pt %g weight %g",wgtObj->GetName(),track->Pt(),mcTrack->Pt(),trackWgt[itrack])); | |
5936324e | 272 | // Int_t iancestor = fUtilityMuonAncestor->GetAncestor(track,MCEvent()); |
273 | // AliVParticle* motherPart = MCEvent()->GetTrack(iancestor); | |
274 | // trackWgt[itrack] = beautyMuWgt->GetBinContent(beautyMuWgt->GetXaxis()->FindBin(motherPart->Pt())); | |
275 | } | |
47b31d8b CP |
276 | } // loop on tracks |
277 | ||
278 | // Loop on selected tracks | |
279 | TArrayI rejectTrack(nTracks); | |
280 | for ( Int_t itrack=0; itrack<nTracks; itrack++) { | |
281 | track = static_cast<AliVParticle*>(selectedTracks.At(itrack)); | |
282 | if ( ! track ) continue; | |
283 | ||
3575043a | 284 | // Check dimuons |
285 | for ( Int_t jtrack=itrack+1; jtrack<nTracks; jtrack++ ) { | |
286 | AliVParticle* auxTrack = static_cast<AliVParticle*>(selectedTracks.At(jtrack)); | |
287 | if ( ! auxTrack ) continue; | |
288 | if ( track->Charge() * auxTrack->Charge() >= 0 ) continue; | |
47b31d8b | 289 | |
3575043a | 290 | TLorentzVector dimuPair = AliAnalysisMuonUtility::GetTrackPair(track,auxTrack); |
291 | Double_t ptMin = TMath::Min(track->Pt(),auxTrack->Pt()); | |
292 | Double_t invMass = dimuPair.M(); | |
293 | if ( fCutOnDimu && invMass > 60. && invMass < 120. ) { | |
294 | rejectTrack[itrack] = 1; | |
295 | rejectTrack[jtrack] = 1; | |
296 | } | |
297 | Double_t dimuWgt = trackWgt[itrack] * trackWgt[jtrack]; | |
298 | if ( istep == kStepReconstructed ) { | |
299 | for ( Int_t itrig=0; itrig<selectTrigClasses.GetEntries(); ++itrig ) { | |
300 | TString trigClassName = ((TObjString*)selectTrigClasses.At(itrig))->GetString(); | |
301 | if ( ! fMuonTrackCuts->TrackPtCutMatchTrigClass(track, fMuonEventCuts->GetTrigClassPtCutLevel(trigClassName)) || ! fMuonTrackCuts->TrackPtCutMatchTrigClass(auxTrack, fMuonEventCuts->GetTrigClassPtCutLevel(trigClassName)) ) continue; | |
302 | ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, "hRecoDimu"))->Fill(ptMin,invMass,dimuWgt); | |
303 | } // loop on triggers | |
304 | } | |
305 | else { | |
306 | if ( trackSources[itrack] == kZbosonMu && trackSources[jtrack] == kZbosonMu ) { | |
307 | Bool_t isAccepted = kTRUE; | |
308 | AliVParticle* muDaughter[2] = {0x0, 0x0}; | |
309 | for ( Int_t imu=0; imu<2; imu++ ) { | |
310 | AliVParticle* currPart = ( imu == 0 ) ? track : auxTrack; | |
311 | if ( currPart->Charge() < 0. ) muDaughter[0] = currPart; | |
312 | else muDaughter[1] = currPart; | |
313 | if ( currPart->Eta() < -4.5 || currPart->Eta() > -2. ) { | |
314 | isAccepted = kFALSE; | |
315 | } | |
316 | } // loop on muons in the pair | |
47b31d8b | 317 | |
3575043a | 318 | Double_t pairRapidity = dimuPair.Rapidity(); |
319 | if ( pairRapidity < -4. || pairRapidity > -2.5 ) isAccepted = kFALSE; | |
320 | //printf("Rapidity Z %g pair %g\n",track->Y(), pairRapidity); | |
47b31d8b | 321 | |
3575043a | 322 | for ( Int_t itrig=0; itrig<selectTrigClasses.GetEntries(); ++itrig ) { |
323 | TString trigClassName = ((TObjString*)selectTrigClasses.At(itrig))->GetString(); | |
324 | if ( isAccepted ) ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, "hGeneratedZ"))->Fill(ptMin,invMass,dimuWgt); | |
325 | ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, "hZmuEtaCorr"))->Fill(muDaughter[0]->Eta(),muDaughter[1]->Eta(),dimuWgt); | |
326 | } // loop on selected trig | |
327 | } // both muons from Z | |
328 | } // kStepGeneratedMC | |
329 | } // loop on auxiliary tracks | |
47b31d8b | 330 | if ( rejectTrack[itrack] > 0 ) continue; |
a07db2fb | 331 | |
332 | Double_t thetaAbsEndDeg = 0; | |
333 | if ( istep == kStepReconstructed ) { | |
ebba9ef0 | 334 | thetaAbsEndDeg = AliAnalysisMuonUtility::GetThetaAbsDeg(track); |
93d6def1 | 335 | } |
a07db2fb | 336 | else { |
337 | thetaAbsEndDeg = ( TMath::Pi()-track->Theta() ) * TMath::RadToDeg(); | |
93d6def1 | 338 | } |
a07db2fb | 339 | Int_t thetaAbsBin = ( thetaAbsEndDeg < 3. ) ? kThetaAbs23 : kThetaAbs310; |
340 | ||
341 | containerInput[kHvarPt] = track->Pt(); | |
342 | containerInput[kHvarEta] = track->Eta(); | |
343 | containerInput[kHvarPhi] = track->Phi(); | |
344 | containerInput[kHvarVz] = ( istep == kStepReconstructed ) ? ipVz : ipVzMC; | |
345 | containerInput[kHvarCharge] = track->Charge()/3.; | |
346 | containerInput[kHvarThetaAbs] = (Double_t)thetaAbsBin; | |
47b31d8b | 347 | containerInput[kHvarMotherType] = (Double_t)trackSources[itrack]; |
a07db2fb | 348 | |
349 | for ( Int_t itrig=0; itrig<selectTrigClasses.GetEntries(); ++itrig ) { | |
350 | TString trigClassName = ((TObjString*)selectTrigClasses.At(itrig))->GetString(); | |
ebba9ef0 | 351 | if ( istep == kStepReconstructed && ! fMuonTrackCuts->TrackPtCutMatchTrigClass(track, fMuonEventCuts->GetTrigClassPtCutLevel(trigClassName)) ) continue; |
5936324e | 352 | ((AliCFContainer*)GetMergeableObject(physSel, trigClassName, centrality, "SingleMuContainer"))->Fill(containerInput,istep,trackWgt[itrack]); |
a07db2fb | 353 | } // loop on selected trigger classes |
354 | } // loop on tracks | |
355 | } // loop on container steps | |
589e3f71 | 356 | } |
662e37fe | 357 | |
662e37fe | 358 | |
9728bcfd | 359 | //________________________________________________________________________ |
a07db2fb | 360 | void AliAnalysisTaskSingleMu::Terminate(Option_t *) { |
589e3f71 | 361 | // |
a07db2fb | 362 | /// Draw some histograms at the end. |
589e3f71 | 363 | // |
33898693 | 364 | |
a07db2fb | 365 | AliVAnalysisMuon::Terminate(""); |
33898693 | 366 | |
a07db2fb | 367 | if ( ! fMergeableCollection ) return; |
368 | ||
369 | TString physSel = fTerminateOptions->At(0)->GetName(); | |
370 | TString trigClassName = fTerminateOptions->At(1)->GetName(); | |
371 | TString centralityRange = fTerminateOptions->At(2)->GetName(); | |
372 | TString furtherOpt = fTerminateOptions->At(3)->GetName(); | |
12e33589 | 373 | |
374 | TString minBiasTrig = ""; | |
375 | TObjArray* optArr = furtherOpt.Tokenize(" "); | |
376 | TString currName = ""; | |
377 | for ( Int_t iopt=0; iopt<optArr->GetEntries(); iopt++ ) { | |
378 | currName = optArr->At(iopt)->GetName(); | |
df1e3f7a | 379 | if ( currName.Contains("-B-") || currName.Contains("ANY") ) minBiasTrig = currName; |
12e33589 | 380 | } |
381 | delete optArr; | |
33898693 | 382 | |
a07db2fb | 383 | furtherOpt.ToUpper(); |
df1e3f7a | 384 | Bool_t plotChargeAsymmetry = furtherOpt.Contains("ASYM"); |
a07db2fb | 385 | |
9bb3925a | 386 | AliCFContainer* inContainer = static_cast<AliCFContainer*> ( GetSum(physSel,trigClassName,centralityRange,"SingleMuContainer") ); |
387 | if ( ! inContainer ) return; | |
388 | ||
389 | AliCFContainer* cfContainer = inContainer; | |
390 | ||
391 | if ( ! furtherOpt.Contains("GENINTRIGCLASS") && trigClassName != "ANY" ) { | |
392 | // The output container contains both the reconstructed and (in MC) | |
393 | // the generated muons in a specific trigger class. | |
394 | // Since the trigger pt level of the track is matched to the trigger class, | |
395 | // analyzing the MUHigh trigger (for example) is equivalent of analyzing | |
396 | // Hpt tracks. | |
397 | // However, in this way, the generated muons distribution depend | |
398 | // on a condition on the reconstructed track. | |
399 | // If we calculate the Acc.xEff. as the ratio of reconstructed and | |
400 | // generated muons IN A TRIGGER CLASS, we are biasing the final result. | |
401 | // The correct value of Acc.xEff. is hence obtained as the distribution | |
402 | // of reconstructed tracks in a muon trigger class, divided by the | |
403 | // total number of generated class (which is in the "class" ANY). | |
404 | // The following code sets the generated muons as the one in the class ANY. | |
405 | // The feature is the default. If you want the generated muons in the same | |
406 | // trigger class as the generated tracks, use option "GENINTRIGCLASS" | |
407 | AliCFContainer* fullMCcontainer = static_cast<AliCFContainer*> ( GetSum(physSel,"ANY",centralityRange,"SingleMuContainer") ); | |
408 | if ( fullMCcontainer ) { | |
409 | cfContainer = static_cast<AliCFContainer*>(cfContainer->Clone("SingleMuContainerCombo")); | |
410 | cfContainer->SetGrid(kStepGeneratedMC,fullMCcontainer->GetGrid(kStepGeneratedMC)); | |
411 | } | |
412 | } | |
c6e0141f | 413 | |
414 | AliCFEffGrid* effSparse = new AliCFEffGrid(Form("eff%s", cfContainer->GetName()),Form("Efficiency %s", cfContainer->GetTitle()),*cfContainer); | |
415 | effSparse->CalculateEfficiency(kStepReconstructed, kStepGeneratedMC); | |
416 | ||
417 | AliCFGridSparse* gridSparseArray[3] = {effSparse->GetNum(), effSparse->GetDen(), effSparse}; | |
418 | TString gridSparseName[3] = {cfContainer->GetStepTitle(kStepReconstructed), cfContainer->GetStepTitle(kStepGeneratedMC), "Efficiency"}; | |
33898693 | 419 | |
420 | Int_t srcColors[kNtrackSources] = {kBlack, kRed, kSpring, kTeal, kBlue, kViolet, kMagenta, kOrange, kGray}; | |
421 | // TString allSrcNames = ""; | |
422 | // for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) { | |
423 | // if ( ! allSrcNames.IsNull() ) allSrcNames.Append(" "); | |
424 | // allSrcNames += fSrcKeys->At(isrc)->GetName(); | |
425 | // } | |
426 | ||
a07db2fb | 427 | TCanvas* can = 0x0; |
428 | Int_t xshift = 100; | |
429 | Int_t yshift = 100; | |
430 | Int_t igroup1 = -1; | |
431 | Int_t igroup2 = 0; | |
432 | ||
433 | Bool_t isMC = furtherOpt.Contains("MC"); | |
33898693 | 434 | |
435 | TAxis* srcAxis = gridSparseArray[0]->GetAxis(kHvarMotherType); | |
8ac7181e | 436 | Int_t unIdBin = srcAxis->FindBin(fSrcKeys->At(kUnidentified)->GetName()); |
437 | if ( unIdBin < 1 ) unIdBin = srcAxis->GetNbins(); | |
33898693 | 438 | |
439 | Int_t firstSrcBin = ( isMC ) ? 1 : unIdBin; | |
440 | Int_t lastSrcBin = ( isMC ) ? srcAxis->GetNbins() - 1 : unIdBin; | |
441 | if ( ! isMC ) srcColors[unIdBin-1] = 1; | |
ac2e110b | 442 | |
a07db2fb | 443 | //////////////// |
444 | // Kinematics // | |
445 | //////////////// | |
ac2e110b | 446 | TString chargeNames[3] = {fChargeKeys->At(0)->GetName(), fChargeKeys->At(1)->GetName(), "Total"}; |
447 | THashList histoList[3]; | |
448 | for ( Int_t icharge=0; icharge<3; icharge++ ) { | |
449 | histoList[icharge].SetName(chargeNames[icharge].Data()); | |
450 | } | |
33898693 | 451 | for ( Int_t isrc = firstSrcBin; isrc <= lastSrcBin; ++isrc ) { |
ac2e110b | 452 | for ( Int_t icharge=0; icharge<3; ++icharge ) { |
453 | Int_t icharge1 = ( icharge < 2 ) ? icharge : 0; | |
454 | Int_t icharge2 = ( icharge < 2 ) ? icharge : 1; | |
c6e0141f | 455 | for ( Int_t igrid=0; igrid<3; ++igrid ) { |
456 | if ( gridSparseArray[igrid]->GetEntries() == 0. ) break; | |
457 | if ( gridSparseArray[igrid]->IsA() != AliCFEffGrid::Class() ) { | |
458 | SetSparseRange(gridSparseArray[igrid], kHvarEta, "", -3.999, -2.501); | |
33898693 | 459 | SetSparseRange(gridSparseArray[igrid], kHvarMotherType, "", isrc, isrc, "USEBIN"); |
ac2e110b | 460 | SetSparseRange(gridSparseArray[igrid], kHvarCharge, "", icharge1+1, icharge2+1, "USEBIN"); |
c6e0141f | 461 | } |
ac2e110b | 462 | TH1* histo = gridSparseArray[igrid]->Project(kHvarPt, kHvarEta); |
33898693 | 463 | histo->SetName(Form("hPtEta_%s_%s_%s", gridSparseName[igrid].Data(), srcAxis->GetBinLabel(isrc), chargeNames[icharge].Data())); |
f7852207 | 464 | histo->SetDirectory(0); |
ac2e110b | 465 | if ( histo->Integral() > 0 ) histoList[icharge].Add(histo); |
c6e0141f | 466 | for ( Int_t iproj=0; iproj<4; ++iproj ) { |
ac2e110b | 467 | histo = gridSparseArray[igrid]->Project(iproj); |
33898693 | 468 | histo->SetName(Form("proj%i_%s_%s_%s", iproj, gridSparseName[igrid].Data(), srcAxis->GetBinLabel(isrc), chargeNames[icharge].Data())); |
f7852207 | 469 | histo->SetDirectory(0); |
ac2e110b | 470 | if ( histo->Integral() > 0 ) histoList[icharge].Add(histo); |
471 | } // loop on projections | |
472 | } // loop on grid sparse | |
473 | } // loop on charge | |
c6e0141f | 474 | } // loop on track sources |
475 | ||
ac2e110b | 476 | // Get charge asymmetry or mu+/mu- |
477 | THashList histoListRatio; | |
478 | TString basePlotName = plotChargeAsymmetry ? "ChargeAsymmetry" : "ChargeRatio"; | |
479 | histoListRatio.SetName(basePlotName.Data()); | |
480 | Int_t baseCharge = 1; | |
481 | Int_t auxCharge = 1-baseCharge; | |
482 | for ( Int_t ihisto=0; ihisto<histoList[baseCharge].GetEntries(); ihisto++ ) { | |
483 | TObject* obj = histoList[baseCharge].At(ihisto); | |
484 | TString histoName = obj->GetName(); | |
485 | if ( histoName.Contains(gridSparseName[2].Data()) ) continue; | |
486 | TString searchName = histoName; | |
487 | searchName.ReplaceAll(fChargeKeys->At(baseCharge)->GetName(), fChargeKeys->At(auxCharge)->GetName()); | |
488 | TH1* auxHisto = static_cast<TH1*> (histoList[auxCharge].FindObject(searchName.Data())); | |
489 | if ( ! auxHisto ) continue; | |
490 | histoName.ReplaceAll(fChargeKeys->At(baseCharge)->GetName(),basePlotName.Data()); | |
491 | TH1* histo = static_cast<TH1*> (obj->Clone(histoName.Data())); | |
33898693 | 492 | if ( plotChargeAsymmetry ) { |
ac2e110b | 493 | histo->Add(auxHisto, -1.); |
494 | // h2 + h1 = 2xh2 + (h1-h2) | |
495 | auxHisto->Add(auxHisto, histo, 2.); | |
496 | } | |
497 | histo->Divide(auxHisto); | |
498 | TString axisTitle = plotChargeAsymmetry ? Form("(%s - %s) / (%s + %s)", fChargeKeys->At(1)->GetName(), fChargeKeys->At(0)->GetName(), fChargeKeys->At(1)->GetName(), fChargeKeys->At(0)->GetName()) : Form("%s / %s", fChargeKeys->At(1)->GetName(), fChargeKeys->At(0)->GetName()); | |
499 | axisTitle.ReplaceAll("MuPlus","#mu^{+}"); | |
500 | axisTitle.ReplaceAll("MuMinus","#mu^{-}"); | |
501 | histo->GetYaxis()->SetTitle(axisTitle.Data()); | |
502 | histo->SetStats(kFALSE); | |
f7852207 | 503 | histo->SetDirectory(0); |
ac2e110b | 504 | histoListRatio.Add(histo); |
505 | } | |
506 | ||
507 | // Plot kinematics | |
508 | TString histoName = "", drawOpt = ""; | |
509 | for ( Int_t itype=0; itype<3; itype++ ) { | |
510 | THashList* currList = 0x0; | |
511 | Int_t nCharges = 1; | |
512 | if ( itype == 1 ) currList = &histoListRatio; | |
513 | else if ( itype == 2 ) currList = &histoList[2]; | |
514 | else nCharges = 2; | |
515 | for ( Int_t igrid=0; igrid<3; ++igrid ) { | |
516 | igroup1 = igrid; | |
517 | TCanvas* canKine = 0x0; | |
518 | TLegend* legKine = 0x0; | |
519 | for ( Int_t iproj=0; iproj<4; ++iproj ) { | |
33898693 | 520 | for ( Int_t isrc = firstSrcBin; isrc <= lastSrcBin; ++isrc ) { |
ac2e110b | 521 | for ( Int_t icharge=0; icharge<nCharges; ++icharge ) { |
522 | if ( itype == 0 ) currList = &histoList[icharge]; | |
33898693 | 523 | histoName = Form("proj%i_%s_%s_%s", iproj, gridSparseName[igrid].Data(), srcAxis->GetBinLabel(isrc), currList->GetName()); |
ac2e110b | 524 | TH1* histo = static_cast<TH1*>(currList->FindObject(histoName.Data())); |
525 | if ( ! histo ) continue; | |
526 | if ( ! canKine ) { | |
527 | igroup2 = itype; | |
528 | igroup1 = igrid; | |
529 | currName = Form("%s_%s_%s", GetName(), currList->GetName(), gridSparseName[igrid].Data()); | |
530 | canKine = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600); | |
531 | canKine->Divide(2,2); | |
532 | legKine = new TLegend(0.6, 0.6, 0.8, 0.8); | |
533 | } | |
534 | canKine->cd(iproj+1); | |
535 | if ( itype != 1 ) { | |
536 | if ( ( iproj == kHvarPt || iproj == kHvarVz ) && gridSparseArray[igrid]->IsA() != AliCFEffGrid::Class() ) gPad->SetLogy(); | |
537 | } | |
538 | Bool_t isFirst = ( gPad->GetListOfPrimitives()->GetEntries() == 0 ); | |
539 | drawOpt = isFirst ? "e" : "esames"; | |
33898693 | 540 | histo->SetLineColor(srcColors[isrc-1]); |
541 | histo->SetMarkerColor(srcColors[isrc-1]); | |
ac2e110b | 542 | histo->SetMarkerStyle(20+4*icharge); |
543 | histo->Draw(drawOpt.Data()); | |
544 | TPaveStats* paveStats = (TPaveStats*)histo->FindObject("stats"); | |
33898693 | 545 | if ( paveStats ) paveStats->SetTextColor(srcColors[isrc-1]); |
ac2e110b | 546 | if ( iproj == 0 ) { |
547 | TString legEntry = ( itype == 0 ) ? fChargeKeys->At(icharge)->GetName() : ""; | |
33898693 | 548 | if ( isMC ) legEntry += Form(" %s", srcAxis->GetBinLabel(isrc)); |
ac2e110b | 549 | if ( ! legEntry.IsNull() ) legKine->AddEntry(histo,legEntry.Data(), "lp"); |
550 | } | |
551 | } // loop on mu charge | |
552 | } // loop on track sources | |
553 | } // loop on projections | |
554 | if ( legKine && legKine->GetNRows() > 0 ) { | |
555 | canKine->cd(1); | |
556 | legKine->Draw("same"); | |
557 | } | |
558 | } // loop on grid sparse | |
559 | } // loop on types | |
ac2e110b | 560 | |
c6e0141f | 561 | |
562 | for ( Int_t igrid=0; igrid<3; igrid++ ) { | |
c6e0141f | 563 | if ( gridSparseArray[igrid]->IsA() == AliCFEffGrid::Class() ) continue; |
564 | SetSparseRange(gridSparseArray[igrid], kHvarCharge, "", 1, gridSparseArray[igrid]->GetAxis(kHvarCharge)->GetNbins(), "USEBIN"); // Reset range | |
a07db2fb | 565 | } // loop on container steps |
a07db2fb | 566 | |
567 | ////////////////////// | |
568 | // Event statistics // | |
569 | ////////////////////// | |
570 | printf("\nTotal analyzed events:\n"); | |
571 | TString evtSel = Form("trigger:%s", trigClassName.Data()); | |
572 | fEventCounters->PrintSum(evtSel.Data()); | |
573 | printf("Physics selected analyzed events:\n"); | |
574 | evtSel = Form("trigger:%s/selected:yes", trigClassName.Data()); | |
575 | fEventCounters->PrintSum(evtSel.Data()); | |
576 | ||
577 | TString countPhysSel = "any"; | |
578 | if ( physSel.Contains(fPhysSelKeys->At(kPhysSelPass)->GetName()) ) countPhysSel = "yes"; | |
579 | else if ( physSel.Contains(fPhysSelKeys->At(kPhysSelReject)->GetName()) ) countPhysSel="no"; | |
580 | countPhysSel.Prepend("selected:"); | |
581 | printf("Analyzed events vs. centrality:\n"); | |
582 | evtSel = Form("trigger:%s/%s", trigClassName.Data(), countPhysSel.Data()); | |
583 | fEventCounters->Print("centrality",evtSel.Data(),kTRUE); | |
584 | ||
ac2e110b | 585 | TString outFilename = Form("/tmp/out%s.root", GetName()); |
f7852207 | 586 | printf("\nCreating file %s\n", outFilename.Data()); |
ac2e110b | 587 | TFile* outFile = new TFile(outFilename.Data(),"RECREATE"); |
588 | for ( Int_t icharge=0; icharge<3; icharge++ ) { | |
589 | histoList[icharge].Write(); | |
590 | } | |
591 | histoListRatio.Write(); | |
592 | outFile->Close(); | |
ac2e110b | 593 | |
a07db2fb | 594 | /////////////////// |
595 | // Vertex method // | |
596 | /////////////////// | |
597 | if ( ! furtherOpt.Contains("VERTEX") ) return; | |
a07db2fb | 598 | igroup1++; |
12e33589 | 599 | TH1* eventVertex = (TH1*)GetSum(physSel, minBiasTrig, centralityRange, "hIpVtx"); |
a07db2fb | 600 | if ( ! eventVertex ) return; |
601 | Double_t minZ = -9.99, maxZ = 9.99; | |
602 | Double_t meanZ = 0., sigmaZ = 4.; | |
603 | Double_t nSigma = 2.; | |
12e33589 | 604 | TString fitOpt = "R0S"; |
a07db2fb | 605 | Bool_t fixFitRange = kFALSE; |
606 | TString fitFormula = Form("[0]+[1]*(x+[2])"); | |
33898693 | 607 | |
608 | // Get vertex shape | |
a07db2fb | 609 | if ( eventVertex->GetSumw2N() == 0 ) eventVertex->Sumw2(); |
610 | Double_t eventVtxIntegral = eventVertex->Integral(0,eventVertex->GetNbinsX()+1); // Include under/overflow | |
611 | printf("Event vertex integral %.0f\n\n", eventVtxIntegral); | |
612 | if ( eventVtxIntegral <= 0. ) return; | |
613 | eventVertex->Scale(1./eventVtxIntegral); | |
614 | printf("\nFit MB vertex\n"); | |
615 | eventVertex->Fit("gaus",fitOpt.Data(),"",minZ,maxZ); | |
616 | TF1* vtxFit = (TF1*)eventVertex->GetListOfFunctions()->FindObject("gaus"); | |
617 | currName = "vtxIntegrated"; | |
618 | can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600); | |
619 | can->SetLogy(); | |
620 | eventVertex->Draw(); | |
621 | vtxFit->Draw("same"); | |
33898693 | 622 | |
a07db2fb | 623 | |
624 | enum {kRecoHF, kRecoBkg, kInputHF, kInputDecay, kRecoAll, kNrecoHistos}; | |
625 | TString baseRecoName[kNrecoHistos] = {"RecoHF", "RecoBkg", "InputHF", "InputDecay", "RecoAll"}; | |
626 | TArrayI sumMothers[kNrecoHistos]; | |
627 | sumMothers[kRecoHF].Set(0); | |
628 | sumMothers[kRecoBkg].Set(0); | |
629 | sumMothers[kInputHF].Set(3); | |
8ac7181e | 630 | |
631 | sumMothers[kInputHF][0] = srcAxis->FindBin(fSrcKeys->At(kCharmMu)->GetName()); | |
632 | sumMothers[kInputHF][1] = srcAxis->FindBin(fSrcKeys->At(kBeautyMu)->GetName()); | |
633 | sumMothers[kInputHF][2] = srcAxis->FindBin(fSrcKeys->At(kQuarkoniumMu)->GetName()); | |
a07db2fb | 634 | sumMothers[kInputDecay].Set(1); |
8ac7181e | 635 | sumMothers[kInputDecay][0] = srcAxis->FindBin(fSrcKeys->At(kDecayMu)->GetName()); |
33898693 | 636 | sumMothers[kRecoAll].Set(srcAxis->GetNbins()); |
8ac7181e | 637 | for ( Int_t isrc=0; isrc<srcAxis->GetNbins(); ++isrc ) { |
638 | sumMothers[kRecoAll][isrc] = isrc+1; | |
b201705a | 639 | } |
a07db2fb | 640 | |
641 | meanZ = vtxFit->GetParameter(1); | |
642 | sigmaZ = vtxFit->GetParameter(2); | |
643 | ||
644 | Double_t minZfit = ( fixFitRange ) ? minZ : meanZ - nSigma*sigmaZ; | |
645 | Double_t maxZfit = ( fixFitRange ) ? maxZ : meanZ + nSigma*sigmaZ; | |
646 | ||
647 | TF1* fitFunc = new TF1("fitFunc", fitFormula.Data(), minZ, maxZ); | |
648 | fitFunc->SetLineColor(2); | |
649 | fitFunc->SetParNames("Line norm", "Line slope", "Free path"); | |
650 | const Double_t kFreePath = 153.; // 150.; // 130.; // cm | |
651 | //fitFunc->SetParameters(0.,1.); | |
652 | fitFunc->FixParameter(2, kFreePath); | |
33898693 | 653 | |
a07db2fb | 654 | AliCFGridSparse* gridSparse = cfContainer->GetGrid(kStepReconstructed); |
655 | TAxis* ptAxis = gridSparse->GetAxis(kHvarPt); | |
656 | ||
657 | Double_t slope = 0.; | |
658 | Double_t limitNorm = 0., limitSlope = 0.; | |
12e33589 | 659 | Int_t firstPtBin = 0, lastPtBin = 0; |
660 | ||
661 | gStyle->SetOptFit(1111); | |
33898693 | 662 | |
a07db2fb | 663 | for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) { |
664 | igroup2++; | |
665 | SetSparseRange(gridSparse, kHvarThetaAbs, "", itheta+1, itheta+1, "USEBIN"); | |
666 | SetSparseRange(gridSparse, kHvarPt, "", 1, ptAxis->GetNbins(), "USEBIN"); | |
667 | TH1* recoHisto[kNrecoHistos]; | |
668 | for ( Int_t ireco=0; ireco<kNrecoHistos; ++ireco ) { | |
669 | recoHisto[ireco] = gridSparse->Project(kHvarPt); | |
670 | histoName = Form("%sMuon_%s", baseRecoName[ireco].Data(), fThetaAbsKeys->At(itheta)->GetName()); | |
671 | recoHisto[ireco]->SetName(histoName.Data()); | |
672 | recoHisto[ireco]->SetTitle(histoName.Data()); | |
673 | recoHisto[ireco]->Reset(); | |
674 | recoHisto[ireco]->Sumw2(); | |
675 | for ( Int_t isrc=0; isrc<sumMothers[ireco].GetSize(); ++isrc ) { | |
33898693 | 676 | SetSparseRange(gridSparse, kHvarMotherType, "", sumMothers[ireco][isrc], sumMothers[ireco][isrc], "USEBIN"); |
a07db2fb | 677 | TH1* auxHisto = gridSparse->Project(kHvarPt); |
678 | recoHisto[ireco]->Add(auxHisto); | |
679 | delete auxHisto; | |
680 | } | |
b201705a | 681 | } |
be436bdc | 682 | SetSparseRange(gridSparse, kHvarMotherType, "", firstPtBin, lastSrcBin, "USEBIN"); |
a07db2fb | 683 | Int_t currDraw = 0; |
33898693 | 684 | |
a07db2fb | 685 | for ( Int_t ibinpt=0; ibinpt<=ptAxis->GetNbins(); ++ibinpt ) { |
686 | firstPtBin = ibinpt; | |
687 | lastPtBin = ( ibinpt == 0 ) ? ptAxis->GetNbins() : ibinpt; | |
688 | SetSparseRange(gridSparse, kHvarPt, "", firstPtBin, lastPtBin, "USEBIN"); | |
689 | TH1* histo = gridSparse->Project(kHvarVz); | |
690 | histo->SetName(Form("hVtx_%s_%s_ptBin%i", cfContainer->GetStepTitle(kStepReconstructed), fThetaAbsKeys->At(itheta)->GetName(), ibinpt)); | |
ac2e110b | 691 | if ( histo->Integral() < 100. ) break; |
a07db2fb | 692 | printf("\nFit %.2f < pt < %.2f (entries %g)\n", ptAxis->GetBinLowEdge(firstPtBin), ptAxis->GetBinUpEdge(lastPtBin), histo->GetEntries()); |
693 | histo->Divide(eventVertex); | |
694 | Double_t norm = histo->GetBinContent(histo->FindBin(0.)); | |
695 | histo->GetYaxis()->SetTitle("#frac{dN_{#mu}}{dv_{z}} / #left(#frac{1}{N_{MB}}#frac{dN_{MB}}{dv_{z}}#right)"); | |
b84f935c | 696 | histo->SetTitle(Form("%g < p_{T} (GeV/c) < %g",ptAxis->GetBinLowEdge(firstPtBin), ptAxis->GetBinUpEdge(lastPtBin))); |
33898693 | 697 | slope = ( histo->GetBinContent(histo->FindBin(meanZ+sigmaZ)) - |
a07db2fb | 698 | histo->GetBinContent(histo->FindBin(meanZ-sigmaZ)) ) / ( 2. * sigmaZ ); |
699 | ||
700 | if ( slope < 0. ) slope = norm/kFreePath; | |
701 | ||
702 | // Try to fit twice: it fit fails the first time | |
703 | // set some limits on parameters | |
704 | for ( Int_t itry=0; itry<2; itry++ ) { | |
705 | fitFunc->SetParameter(0, norm); | |
706 | fitFunc->SetParameter(1, slope); | |
707 | if ( itry > 0 ) { | |
708 | limitNorm = 2.*histo->Integral(); | |
709 | limitSlope = 2.*histo->Integral()/kFreePath; | |
710 | //fitFunc->SetParLimits(0, 0., limitNorm); // REMEMBER TO CHECK | |
711 | fitFunc->SetParLimits(1, 0., limitSlope); // REMEMBER TO CHECK | |
712 | printf("Norm 0. < %f < %f slope 0. < %f < %f\n", norm, limitNorm, slope, limitSlope); | |
713 | } | |
12e33589 | 714 | TFitResultPtr fitRes = histo->Fit(fitFunc, fitOpt.Data(), "", minZfit, maxZfit); |
a07db2fb | 715 | |
33898693 | 716 | // if ( gMinuit->fCstatu.Contains("CONVERGED") && |
12e33589 | 717 | if ( ((Int_t)fitRes) == 0 && |
33898693 | 718 | fitFunc->GetParameter(0) > 0. && |
a07db2fb | 719 | fitFunc->GetParameter(1) > 0. ) |
720 | break; | |
7e7619d1 | 721 | else if ( furtherOpt.Contains("REFIT") ) printf("Re-fit with limits\n"); |
a07db2fb | 722 | else { |
7e7619d1 | 723 | printf("Warning: fit problems !!!\n"); |
724 | break; | |
a07db2fb | 725 | } |
726 | } | |
727 | ||
728 | Double_t p0 = fitFunc->GetParameter(0); | |
729 | Double_t p0err = fitFunc->GetParError(0); | |
730 | Double_t p1 = fitFunc->GetParameter(1); | |
731 | Double_t p1err = fitFunc->GetParError(1); | |
732 | ||
733 | Double_t nullVz = ( p1 != 0. ) ? -p0/p1 : 0.; | |
734 | Double_t nullVzErr = ( p0 != 0. && p1 != 0. ) ? TMath::Abs(nullVz) * TMath::Sqrt(p0err*p0err/(p0*p0) + p1err*p1err/(p1*p1) ) : 0.; | |
735 | ||
736 | printf("Null value at %f +- %f\n", nullVz - kFreePath, nullVzErr); | |
737 | ||
738 | recoHisto[kRecoHF]->SetBinContent(ibinpt, p0); | |
739 | recoHisto[kRecoHF]->SetBinError(ibinpt, p0err); | |
740 | recoHisto[kRecoBkg]->SetBinContent(ibinpt, ( kFreePath + meanZ ) * p1); | |
741 | recoHisto[kRecoBkg]->SetBinError(ibinpt, ( kFreePath + meanZ ) * p1err); | |
742 | if ( currDraw%4 == 0 ){ | |
743 | currName = Form("vtx_%s_PtBin%i",fThetaAbsKeys->At(itheta)->GetName(), ibinpt); | |
744 | can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600); | |
745 | can->Divide(2,2); | |
746 | } | |
747 | can->cd( currDraw%4 + 1 ); | |
748 | can->SetLogy(); | |
749 | histo->Draw(); | |
750 | fitFunc->DrawCopy("same"); | |
751 | currDraw++; | |
752 | } // loop on pt bins | |
be436bdc | 753 | SetSparseRange(gridSparse, kHvarMotherType, "", firstSrcBin, lastSrcBin, "USEBIN"); |
a07db2fb | 754 | currName = Form("recoPt_%s",fThetaAbsKeys->At(itheta)->GetName()); |
755 | can = new TCanvas(currName.Data(),currName.Data(),(igroup1+1)*xshift,igroup2*yshift,600,600); | |
756 | TLegend* leg = new TLegend(0.6, 0.6, 0.8, 0.8); | |
757 | drawOpt = "e"; | |
758 | for ( Int_t ireco=0; ireco<kNrecoHistos-1; ++ireco ) { | |
759 | if ( recoHisto[ireco]->GetEntries() == 0. ) continue; | |
760 | TH1* ratio = (TH1*)recoHisto[ireco]->Clone(Form("%s_ratio", recoHisto[ireco]->GetName())); | |
761 | ratio->Divide(recoHisto[kRecoAll]); | |
762 | ratio->SetLineColor(srcColors[ireco]); | |
763 | ratio->SetMarkerColor(srcColors[ireco]); | |
764 | ratio->SetMarkerStyle(20+ireco); | |
765 | ratio->GetYaxis()->SetTitle("fraction of total"); | |
766 | ratio->Draw(drawOpt.Data()); | |
767 | leg->AddEntry(ratio,baseRecoName[ireco].Data(), "lp"); | |
768 | drawOpt = "esame"; | |
b201705a | 769 | } |
a07db2fb | 770 | leg->Draw("same"); |
771 | } // loop on theta abs | |
589e3f71 | 772 | } |