Include standard cuts in analysis. Make it more suitable for proof analysis: eliminat...
[u/mrichter/AliRoot.git] / PWG / muon / AliAnalysisTaskSingleMu.cxx
CommitLineData
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.
23/// In the latter case a flag can be activated to produce a tree as output.
9728bcfd 24/// If Monte Carlo information is present, some basics checks are performed.
662e37fe 25///
26/// \author Diego Stocco
27//-----------------------------------------------------------------------------
28
aad6618e 29#define AliAnalysisTaskSingleMu_cxx
30
a07db2fb 31#include "AliAnalysisTaskSingleMu.h"
32
aad6618e 33// ROOT includes
aad6618e 34#include "TROOT.h"
9728bcfd 35#include "TH1.h"
36#include "TH2.h"
9728bcfd 37#include "TAxis.h"
662e37fe 38#include "TCanvas.h"
a07db2fb 39#include "TLegend.h"
9728bcfd 40#include "TMath.h"
589e3f71 41#include "TObjString.h"
55065f3f 42#include "TObjArray.h"
a07db2fb 43#include "TF1.h"
44#include "TStyle.h"
45//#include "TMCProcess.h"
46#include "TMinuit.h"
47#include "TArrayI.h"
aad6618e 48
49// STEER includes
aad6618e 50#include "AliAODEvent.h"
51#include "AliAODTrack.h"
a07db2fb 52#include "AliAODMCParticle.h"
75008b31 53#include "AliMCEvent.h"
9728bcfd 54#include "AliMCParticle.h"
b201705a 55#include "AliESDEvent.h"
56#include "AliESDMuonTrack.h"
a07db2fb 57#include "AliVHeader.h"
58#include "AliAODMCHeader.h"
59#include "AliStack.h"
9728bcfd 60
589e3f71 61// ANALYSIS includes
62#include "AliAnalysisManager.h"
589e3f71 63
64// CORRFW includes
589e3f71 65#include "AliCFContainer.h"
66#include "AliCFGridSparse.h"
589e3f71 67
a07db2fb 68// PWG3 includes
69#include "AliVAnalysisMuon.h"
70#include "AliMergeableCollection.h"
71#include "AliCounterCollection.h"
72
aad6618e 73
662e37fe 74/// \cond CLASSIMP
75ClassImp(AliAnalysisTaskSingleMu) // Class implementation in ROOT context
76/// \endcond
aad6618e 77
9728bcfd 78
aad6618e 79//________________________________________________________________________
a07db2fb 80AliAnalysisTaskSingleMu::AliAnalysisTaskSingleMu() :
81 AliVAnalysisMuon(),
82 fThetaAbsKeys(0x0)
aad6618e 83{
a07db2fb 84 /// Default ctor.
aad6618e 85}
86
9728bcfd 87//________________________________________________________________________
a07db2fb 88AliAnalysisTaskSingleMu::AliAnalysisTaskSingleMu(const char *name, const AliMuonTrackCuts& cuts) :
89 AliVAnalysisMuon(name, cuts),
90 fThetaAbsKeys(0x0)
aad6618e 91{
92 //
a07db2fb 93 /// Constructor.
9728bcfd 94 //
a07db2fb 95 TString thetaAbsKeys = "ThetaAbs23 ThetaAbs310";
96 fThetaAbsKeys = thetaAbsKeys.Tokenize(" ");
589e3f71 97}
98
99
a07db2fb 100//________________________________________________________________________
101AliAnalysisTaskSingleMu::~AliAnalysisTaskSingleMu()
589e3f71 102{
103 //
a07db2fb 104 /// Destructor
589e3f71 105 //
106
a07db2fb 107 delete fThetaAbsKeys;
aad6618e 108}
109
110//___________________________________________________________________________
a07db2fb 111void AliAnalysisTaskSingleMu::MyUserCreateOutputObjects()
aad6618e 112{
589e3f71 113
a07db2fb 114 TH1* histo = 0x0;
115 TString histoName = "", histoTitle = "";
9728bcfd 116
589e3f71 117 Int_t nVzBins = 40;
118 Double_t vzMin = -20., vzMax = 20.;
a07db2fb 119 TString vzName("Vz"), vzTitle("Vz"), vzUnits("cm");
9728bcfd 120
a07db2fb 121 histoName = "hIpVtx";
122 histo = new TH1F(histoName.Data(), histoName.Data(), nVzBins, vzMin, vzMax);
123 histo->SetXTitle("v_{z} (cm)");
124 AddObjectToCollection(histo, kIPVz);
589e3f71 125
a07db2fb 126
127 Int_t nPtBins = 80;
128 Double_t ptMin = 0., ptMax = 80.;
129 TString ptName("Pt"), ptTitle("p_{t}"), ptUnits("GeV/c");
130
131 Int_t nEtaBins = 25;
132 Double_t etaMin = -4.5, etaMax = -2.;
133 TString etaName("Eta"), etaTitle("#eta"), etaUnits("");
134
135 Int_t nPhiBins = 36;
136 Double_t phiMin = 0.; Double_t phiMax = 2*TMath::Pi();
137 TString phiName("Phi"), phiTitle("#phi"), phiUnits("rad");
138
589e3f71 139 Int_t nChargeBins = 2;
140 Double_t chargeMin = -2, chargeMax = 2.;
141 TString chargeName("Charge"), chargeTitle("charge"), chargeUnits("e");
9728bcfd 142
a07db2fb 143 Int_t nThetaAbsEndBins = 2;
144 Double_t thetaAbsEndMin = -0.5, thetaAbsEndMax = 1.5;
145 TString thetaAbsEndName("ThetaAbsEnd"), thetaAbsEndTitle("#theta_{abs}"), thetaAbsEndUnits("a.u.");
146
589e3f71 147 Int_t nMotherTypeBins = kNtrackSources;
148 Double_t motherTypeMin = -0.5, motherTypeMax = (Double_t)kNtrackSources - 0.5;
149 TString motherType("MotherType"), motherTypeTitle("motherType"), motherTypeUnits("");
a07db2fb 150
151 Int_t nbins[kNvars] = {nPtBins, nEtaBins, nPhiBins, nVzBins, nChargeBins, nThetaAbsEndBins, nMotherTypeBins};
152 Double_t xmin[kNvars] = {ptMin, etaMin, phiMin, vzMin, chargeMin, thetaAbsEndMin, motherTypeMin};
153 Double_t xmax[kNvars] = {ptMax, etaMax, phiMax, vzMax, chargeMax, thetaAbsEndMax, motherTypeMax};
154 TString axisTitle[kNvars] = {ptTitle, etaTitle, phiTitle, vzTitle, chargeTitle, thetaAbsEndTitle, motherTypeTitle};
155 TString axisUnits[kNvars] = {ptUnits, etaUnits, phiUnits, vzUnits, chargeUnits, thetaAbsEndUnits, motherTypeUnits};
9728bcfd 156
a07db2fb 157 AliCFContainer* cfContainer = new AliCFContainer("SingleMuContainer","Container for tracks",kNsteps,kNvars,nbins);
4ab8d5a6 158
589e3f71 159 for ( Int_t idim = 0; idim<kNvars; idim++){
160 histoTitle = Form("%s (%s)", axisTitle[idim].Data(), axisUnits[idim].Data());
161 histoTitle.ReplaceAll("()","");
a07db2fb 162
163 cfContainer->SetVarTitle(idim, histoTitle.Data());
164 cfContainer->SetBinLimits(idim, xmin[idim], xmax[idim]);
9728bcfd 165 }
a07db2fb 166
167 TString stepTitle[kNsteps] = {"reconstructed", "generated"};
589e3f71 168
a07db2fb 169 TAxis* currAxis = 0x0;
589e3f71 170 for (Int_t istep=0; istep<kNsteps; istep++){
a07db2fb 171 cfContainer->SetStepTitle(istep, stepTitle[istep].Data());
172 AliCFGridSparse* gridSparse = cfContainer->GetGrid(istep);
173
174 currAxis = gridSparse->GetAxis(kHvarMotherType);
175 for ( Int_t ibin=0; ibin<fSrcKeys->GetEntries(); ibin++ ) {
176 currAxis->SetBinLabel(ibin+1, fSrcKeys->At(ibin)->GetName());
9728bcfd 177 }
178 }
4ab8d5a6 179
a07db2fb 180 AddObjectToCollection(cfContainer, kTrackContainer);
181
182 fMuonTrackCuts.Print("mask");
aad6618e 183}
184
185//________________________________________________________________________
a07db2fb 186void AliAnalysisTaskSingleMu::ProcessEvent(TString physSel, const TObjArray& selectTrigClasses, TString centrality)
aad6618e 187{
188 //
a07db2fb 189 /// Fill output objects
aad6618e 190 //
191
a07db2fb 192 AliVVertex* primaryVertex = ( fAODEvent ) ? (AliVVertex*)fAODEvent->GetPrimaryVertexSPD() : (AliVVertex*)fESDEvent->GetPrimaryVertexSPD();
193 if ( primaryVertex->GetNContributors() < 1 ) return;
b201705a 194
a07db2fb 195 Double_t ipVz = primaryVertex->GetZ();
196 Double_t ipVzMC = 0;
197 if ( IsMC() ) {
198 if ( fMCEvent ) ipVzMC = fMCEvent->GetPrimaryVertex()->GetZ();
b201705a 199 else {
a07db2fb 200 AliAODMCHeader* aodMCHeader = (AliAODMCHeader *)fAODEvent->FindListObject(AliAODMCHeader::StdBranchName());
201 if ( aodMCHeader ) ipVzMC = aodMCHeader->GetVtxZ();
55065f3f 202 }
589e3f71 203 }
a07db2fb 204
205 for ( Int_t itrig=0; itrig<selectTrigClasses.GetEntries(); ++itrig ) {
206 TString trigClassName = ((TObjString*)selectTrigClasses.At(itrig))->GetString();
207 ((TH1*)GetMergeableObject(physSel, trigClassName, centrality, "hIpVtx"))->Fill(ipVz);
9728bcfd 208 }
aad6618e 209
a07db2fb 210 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
211 if ( isPileupFromSPD ) return;
212
213 Double_t containerInput[kNvars];
214 AliVParticle* track = 0x0;
93d6def1 215
a07db2fb 216 for ( Int_t istep = 0; istep<2; ++istep ) {
217 Int_t nTracks = ( istep == kStepReconstructed ) ? GetNTracks() : GetNMCTracks();
218 for (Int_t itrack = 0; itrack < nTracks; itrack++) {
219 track = ( istep == kStepReconstructed ) ? GetTrack(itrack) : GetMCTrack(itrack);
220
221 Bool_t isSelected = ( istep == kStepReconstructed ) ? fMuonTrackCuts.IsSelected(track) : ( TMath::Abs(track->PdgCode()) == 13 );
222 if ( ! isSelected ) continue;
223
224 Int_t trackSrc = ( istep == kStepReconstructed ) ? GetParticleType(track) : RecoTrackMother(track);
225
226 Double_t thetaAbsEndDeg = 0;
227 if ( istep == kStepReconstructed ) {
228 Double_t rAbsEnd = ( fAODEvent ) ? ((AliAODTrack*)track)->GetRAtAbsorberEnd(): ((AliESDMuonTrack*)track)->GetRAtAbsorberEnd();
229 thetaAbsEndDeg = TMath::ATan( rAbsEnd / 505. ) * TMath::RadToDeg();
93d6def1 230 }
a07db2fb 231 else {
232 thetaAbsEndDeg = ( TMath::Pi()-track->Theta() ) * TMath::RadToDeg();
93d6def1 233 }
a07db2fb 234 Int_t thetaAbsBin = ( thetaAbsEndDeg < 3. ) ? kThetaAbs23 : kThetaAbs310;
235
236 containerInput[kHvarPt] = track->Pt();
237 containerInput[kHvarEta] = track->Eta();
238 containerInput[kHvarPhi] = track->Phi();
239 containerInput[kHvarVz] = ( istep == kStepReconstructed ) ? ipVz : ipVzMC;
240 containerInput[kHvarCharge] = track->Charge()/3.;
241 containerInput[kHvarThetaAbs] = (Double_t)thetaAbsBin;
242 containerInput[kHvarMotherType] = (Double_t)trackSrc;
243
244 for ( Int_t itrig=0; itrig<selectTrigClasses.GetEntries(); ++itrig ) {
245 TString trigClassName = ((TObjString*)selectTrigClasses.At(itrig))->GetString();
246 if ( istep == kStepReconstructed && ! TrackPtCutMatchTrigClass(track, trigClassName) ) continue;
247 ((AliCFContainer*)GetMergeableObject(physSel, trigClassName, centrality, "SingleMuContainer"))->Fill(containerInput,istep);
248 } // loop on selected trigger classes
249 } // loop on tracks
250 } // loop on container steps
589e3f71 251}
662e37fe 252
662e37fe 253
9728bcfd 254//________________________________________________________________________
a07db2fb 255void AliAnalysisTaskSingleMu::Terminate(Option_t *) {
589e3f71 256 //
a07db2fb 257 /// Draw some histograms at the end.
589e3f71 258 //
589e3f71 259
a07db2fb 260 AliVAnalysisMuon::Terminate("");
b201705a 261
a07db2fb 262 if ( ! fMergeableCollection ) return;
263
264 TString physSel = fTerminateOptions->At(0)->GetName();
265 TString trigClassName = fTerminateOptions->At(1)->GetName();
266 TString centralityRange = fTerminateOptions->At(2)->GetName();
267 TString furtherOpt = fTerminateOptions->At(3)->GetName();
268 furtherOpt.ToUpper();
269
270 AliCFContainer* cfContainer = static_cast<AliCFContainer*> ( GetSum(physSel,trigClassName,centralityRange,"SingleMuContainer") );
271 if ( ! cfContainer ) return;
272
273 Int_t srcColors[kNtrackSources] = {kBlack, kRed, kGreen, kBlue, kViolet, 7, kOrange};
274// TString allSrcNames = "";
275// for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
276// if ( ! allSrcNames.IsNull() ) allSrcNames.Append(" ");
277// allSrcNames += fSrcKeys->At(isrc)->GetName();
278// }
279
280 TCanvas* can = 0x0;
281 Int_t xshift = 100;
282 Int_t yshift = 100;
283 Int_t igroup1 = -1;
284 Int_t igroup2 = 0;
285
286 Bool_t isMC = furtherOpt.Contains("MC");
287 Int_t firstSrc = ( isMC ) ? 0 : kUnidentified;
288 Int_t lastSrc = ( isMC ) ? kNtrackSources - 1 : kUnidentified;
289 if ( ! isMC ) srcColors[kUnidentified] = 1;
290
291 TString histoName = "", currName = "", histoPattern = "", drawOpt = "";
292 ////////////////
293 // Kinematics //
294 ////////////////
295 for ( Int_t istep=0; istep<kNsteps; ++istep ) {
296 igroup1++;
297 igroup2 = 0;
298 AliCFGridSparse* gridSparse = cfContainer->GetGrid(istep);
299 if ( gridSparse->GetEntries() == 0. ) continue;
300 SetSparseRange(gridSparse, kHvarEta, "", -3.999, -2.501);
301 currName = Form("%s_proj_%s", GetName(), cfContainer->GetStepTitle(istep));
302 can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
303 can->Divide(2,2);
304 TLegend* leg = new TLegend(0.6, 0.6, 0.8, 0.8);
305 igroup2++;
306 for ( Int_t iproj=0; iproj<4; ++iproj ) {
307 can->cd(iproj+1);
308 if ( iproj == kHvarPt || iproj == kHvarVz ) gPad->SetLogy();
309 for ( Int_t isrc = firstSrc; isrc <= lastSrc; ++isrc ) {
310 SetSparseRange(gridSparse, kHvarMotherType, "", isrc+1, isrc+1, "USEBIN");
311 for ( Int_t icharge=0; icharge<2; ++icharge ) {
312 SetSparseRange(gridSparse, kHvarCharge, "", icharge+1, icharge+1, "USEBIN");
313 TH1* projHisto = gridSparse->Project(iproj);
314 projHisto->SetName(Form("proj%i_%s_src%i_charge%i", iproj, cfContainer->GetStepTitle(istep), isrc, icharge));
315 if ( projHisto->GetEntries() == 0 ) continue;
316 Bool_t isFirst = ( gPad->GetListOfPrimitives()->GetEntries() == 0 );
317 drawOpt = isFirst ? "e" : "esame";
318 //if ( isrc == kUnidentified && ! drawOpt.Contains("same") ) isMC = kFALSE;
319 //if ( ! isMC ) srcColors[kUnidentified] = 1;
320 projHisto->SetLineColor(srcColors[isrc]);
321 projHisto->SetMarkerColor(srcColors[isrc]);
322 projHisto->SetMarkerStyle(20+4*icharge);
323 projHisto->Draw(drawOpt.Data());
324 if ( iproj == 0 ) {
325 TString legEntry = fChargeKeys->At(icharge)->GetName();
326 if ( isMC ) legEntry += Form(" %s", fSrcKeys->At(isrc)->GetName());
327 leg->AddEntry(projHisto,legEntry.Data(), "lp");
328 }
329 } // loop on mu charge
330 } // loop on track sources
331 if ( iproj == 0 ) leg->Draw("same");
332 } // loop on projections
333 SetSparseRange(gridSparse, kHvarCharge, "", 1, gridSparse->GetAxis(kHvarCharge)->GetNbins(), "USEBIN"); // Reset range
334 } // loop on container steps
335
336
337 //////////////////////
338 // Event statistics //
339 //////////////////////
340 printf("\nTotal analyzed events:\n");
341 TString evtSel = Form("trigger:%s", trigClassName.Data());
342 fEventCounters->PrintSum(evtSel.Data());
343 printf("Physics selected analyzed events:\n");
344 evtSel = Form("trigger:%s/selected:yes", trigClassName.Data());
345 fEventCounters->PrintSum(evtSel.Data());
346
347 TString countPhysSel = "any";
348 if ( physSel.Contains(fPhysSelKeys->At(kPhysSelPass)->GetName()) ) countPhysSel = "yes";
349 else if ( physSel.Contains(fPhysSelKeys->At(kPhysSelReject)->GetName()) ) countPhysSel="no";
350 countPhysSel.Prepend("selected:");
351 printf("Analyzed events vs. centrality:\n");
352 evtSel = Form("trigger:%s/%s", trigClassName.Data(), countPhysSel.Data());
353 fEventCounters->Print("centrality",evtSel.Data(),kTRUE);
354
355
356 ///////////////////
357 // Vertex method //
358 ///////////////////
359 if ( ! furtherOpt.Contains("VERTEX") ) return;
360 Int_t firstMother = kUnidentified, lastMother = kUnidentified;
361 igroup1++;
362 TH1* eventVertex = (TH1*)GetSum(physSel, "CINT7-I-NOPF-ALLNOTRD", centralityRange, "hIpVtx");
363 if ( ! eventVertex ) return;
364 Double_t minZ = -9.99, maxZ = 9.99;
365 Double_t meanZ = 0., sigmaZ = 4.;
366 Double_t nSigma = 2.;
367 TString fitOpt = "R0";
368 Bool_t fixFitRange = kFALSE;
369 TString fitFormula = Form("[0]+[1]*(x+[2])");
370
371 // Get vertex shape
372 if ( eventVertex->GetSumw2N() == 0 ) eventVertex->Sumw2();
373 Double_t eventVtxIntegral = eventVertex->Integral(0,eventVertex->GetNbinsX()+1); // Include under/overflow
374 printf("Event vertex integral %.0f\n\n", eventVtxIntegral);
375 if ( eventVtxIntegral <= 0. ) return;
376 eventVertex->Scale(1./eventVtxIntegral);
377 printf("\nFit MB vertex\n");
378 eventVertex->Fit("gaus",fitOpt.Data(),"",minZ,maxZ);
379 TF1* vtxFit = (TF1*)eventVertex->GetListOfFunctions()->FindObject("gaus");
380 currName = "vtxIntegrated";
381 can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
382 can->SetLogy();
383 eventVertex->Draw();
384 vtxFit->Draw("same");
b201705a 385
a07db2fb 386
387 enum {kRecoHF, kRecoBkg, kInputHF, kInputDecay, kRecoAll, kNrecoHistos};
388 TString baseRecoName[kNrecoHistos] = {"RecoHF", "RecoBkg", "InputHF", "InputDecay", "RecoAll"};
389 TArrayI sumMothers[kNrecoHistos];
390 sumMothers[kRecoHF].Set(0);
391 sumMothers[kRecoBkg].Set(0);
392 sumMothers[kInputHF].Set(3);
393 sumMothers[kInputHF][0] = kCharmMu;
394 sumMothers[kInputHF][1] = kBeautyMu;
395 sumMothers[kInputHF][2] = kQuarkoniumMu;
396 sumMothers[kInputDecay].Set(1);
397 sumMothers[kInputDecay][0] = kDecayMu;
398 sumMothers[kRecoAll].Set(kNtrackSources);
399 for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
400 sumMothers[kRecoAll][isrc] = isrc;
b201705a 401 }
a07db2fb 402
403 meanZ = vtxFit->GetParameter(1);
404 sigmaZ = vtxFit->GetParameter(2);
405
406 Double_t minZfit = ( fixFitRange ) ? minZ : meanZ - nSigma*sigmaZ;
407 Double_t maxZfit = ( fixFitRange ) ? maxZ : meanZ + nSigma*sigmaZ;
408
409 TF1* fitFunc = new TF1("fitFunc", fitFormula.Data(), minZ, maxZ);
410 fitFunc->SetLineColor(2);
411 fitFunc->SetParNames("Line norm", "Line slope", "Free path");
412 const Double_t kFreePath = 153.; // 150.; // 130.; // cm
413 //fitFunc->SetParameters(0.,1.);
414 fitFunc->FixParameter(2, kFreePath);
415
416 AliCFGridSparse* gridSparse = cfContainer->GetGrid(kStepReconstructed);
417 TAxis* ptAxis = gridSparse->GetAxis(kHvarPt);
418
419 Double_t slope = 0.;
420 Double_t limitNorm = 0., limitSlope = 0.;
421 Int_t firstPtBin = 0, lastPtBin = 0;
422
423 for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) {
424 igroup2++;
425 SetSparseRange(gridSparse, kHvarThetaAbs, "", itheta+1, itheta+1, "USEBIN");
426 SetSparseRange(gridSparse, kHvarPt, "", 1, ptAxis->GetNbins(), "USEBIN");
427 TH1* recoHisto[kNrecoHistos];
428 for ( Int_t ireco=0; ireco<kNrecoHistos; ++ireco ) {
429 recoHisto[ireco] = gridSparse->Project(kHvarPt);
430 histoName = Form("%sMuon_%s", baseRecoName[ireco].Data(), fThetaAbsKeys->At(itheta)->GetName());
431 recoHisto[ireco]->SetName(histoName.Data());
432 recoHisto[ireco]->SetTitle(histoName.Data());
433 recoHisto[ireco]->Reset();
434 recoHisto[ireco]->Sumw2();
435 for ( Int_t isrc=0; isrc<sumMothers[ireco].GetSize(); ++isrc ) {
436 SetSparseRange(gridSparse, kHvarMotherType, "", sumMothers[ireco][isrc]+1, sumMothers[ireco][isrc]+1, "USEBIN");
437 TH1* auxHisto = gridSparse->Project(kHvarPt);
438 recoHisto[ireco]->Add(auxHisto);
439 delete auxHisto;
440 }
b201705a 441 }
a07db2fb 442 SetSparseRange(gridSparse, kHvarMotherType, "", firstMother+1, lastMother+1, "USEBIN");
443 Int_t currDraw = 0;
444
445 for ( Int_t ibinpt=0; ibinpt<=ptAxis->GetNbins(); ++ibinpt ) {
446 firstPtBin = ibinpt;
447 lastPtBin = ( ibinpt == 0 ) ? ptAxis->GetNbins() : ibinpt;
448 SetSparseRange(gridSparse, kHvarPt, "", firstPtBin, lastPtBin, "USEBIN");
449 TH1* histo = gridSparse->Project(kHvarVz);
450 histo->SetName(Form("hVtx_%s_%s_ptBin%i", cfContainer->GetStepTitle(kStepReconstructed), fThetaAbsKeys->At(itheta)->GetName(), ibinpt));
451 if ( histo->GetEntries() < 100. ) break;
452 printf("\nFit %.2f < pt < %.2f (entries %g)\n", ptAxis->GetBinLowEdge(firstPtBin), ptAxis->GetBinUpEdge(lastPtBin), histo->GetEntries());
453 histo->Divide(eventVertex);
454 Double_t norm = histo->GetBinContent(histo->FindBin(0.));
455 histo->GetYaxis()->SetTitle("#frac{dN_{#mu}}{dv_{z}} / #left(#frac{1}{N_{MB}}#frac{dN_{MB}}{dv_{z}}#right)");
456 slope = ( histo->GetBinContent(histo->FindBin(meanZ+sigmaZ)) -
457 histo->GetBinContent(histo->FindBin(meanZ-sigmaZ)) ) / ( 2. * sigmaZ );
458
459 if ( slope < 0. ) slope = norm/kFreePath;
460
461 // Try to fit twice: it fit fails the first time
462 // set some limits on parameters
463 for ( Int_t itry=0; itry<2; itry++ ) {
464 fitFunc->SetParameter(0, norm);
465 fitFunc->SetParameter(1, slope);
466 if ( itry > 0 ) {
467 limitNorm = 2.*histo->Integral();
468 limitSlope = 2.*histo->Integral()/kFreePath;
469 //fitFunc->SetParLimits(0, 0., limitNorm); // REMEMBER TO CHECK
470 fitFunc->SetParLimits(1, 0., limitSlope); // REMEMBER TO CHECK
471 printf("Norm 0. < %f < %f slope 0. < %f < %f\n", norm, limitNorm, slope, limitSlope);
472 }
473 histo->Fit(fitFunc, fitOpt.Data(), "", minZfit, maxZfit);
474
475 if ( gMinuit->fCstatu.Contains("CONVERGED") &&
476 fitFunc->GetParameter(0) > 0. &&
477 fitFunc->GetParameter(1) > 0. )
478 break;
479 else {
480 printf("Warning: fit problems !!!\n"); // COMMENT TO REFIT
481 break; // COMMENT TO REFIT
482 //printf("Re-fit with limits\n"); //UNCOMMENT TO REFIT
483 }
484 }
485
486 Double_t p0 = fitFunc->GetParameter(0);
487 Double_t p0err = fitFunc->GetParError(0);
488 Double_t p1 = fitFunc->GetParameter(1);
489 Double_t p1err = fitFunc->GetParError(1);
490
491 Double_t nullVz = ( p1 != 0. ) ? -p0/p1 : 0.;
492 Double_t nullVzErr = ( p0 != 0. && p1 != 0. ) ? TMath::Abs(nullVz) * TMath::Sqrt(p0err*p0err/(p0*p0) + p1err*p1err/(p1*p1) ) : 0.;
493
494 printf("Null value at %f +- %f\n", nullVz - kFreePath, nullVzErr);
495
496 recoHisto[kRecoHF]->SetBinContent(ibinpt, p0);
497 recoHisto[kRecoHF]->SetBinError(ibinpt, p0err);
498 recoHisto[kRecoBkg]->SetBinContent(ibinpt, ( kFreePath + meanZ ) * p1);
499 recoHisto[kRecoBkg]->SetBinError(ibinpt, ( kFreePath + meanZ ) * p1err);
500 if ( currDraw%4 == 0 ){
501 currName = Form("vtx_%s_PtBin%i",fThetaAbsKeys->At(itheta)->GetName(), ibinpt);
502 can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
503 can->Divide(2,2);
504 }
505 can->cd( currDraw%4 + 1 );
506 can->SetLogy();
507 histo->Draw();
508 fitFunc->DrawCopy("same");
509 currDraw++;
510 } // loop on pt bins
511 SetSparseRange(gridSparse, kHvarMotherType, "", firstMother+1, lastMother+1, "USEBIN");
512 TH1* totalPtDistrib = gridSparse->Project(kHvarPt);
513 totalPtDistrib->SetName(Form("totalPtDistrib_%s", fThetaAbsKeys->GetName()));
514 currName = Form("recoPt_%s",fThetaAbsKeys->At(itheta)->GetName());
515 can = new TCanvas(currName.Data(),currName.Data(),(igroup1+1)*xshift,igroup2*yshift,600,600);
516 TLegend* leg = new TLegend(0.6, 0.6, 0.8, 0.8);
517 drawOpt = "e";
518 for ( Int_t ireco=0; ireco<kNrecoHistos-1; ++ireco ) {
519 if ( recoHisto[ireco]->GetEntries() == 0. ) continue;
520 TH1* ratio = (TH1*)recoHisto[ireco]->Clone(Form("%s_ratio", recoHisto[ireco]->GetName()));
521 ratio->Divide(recoHisto[kRecoAll]);
522 ratio->SetLineColor(srcColors[ireco]);
523 ratio->SetMarkerColor(srcColors[ireco]);
524 ratio->SetMarkerStyle(20+ireco);
525 ratio->GetYaxis()->SetTitle("fraction of total");
526 ratio->Draw(drawOpt.Data());
527 leg->AddEntry(ratio,baseRecoName[ireco].Data(), "lp");
528 drawOpt = "esame";
b201705a 529 }
a07db2fb 530 leg->Draw("same");
531 } // loop on theta abs
589e3f71 532}