1 /**************************************************************************
2 * Copyright(c) 1998-2007, 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 //-----------------------------------------------------------------------------
19 /// \class AliAnalysisTaskSingleMu
20 /// Analysis task for single muons in the spectrometer.
21 /// The output is a list of histograms and CF containers.
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.
24 /// If Monte Carlo information is present, some basics checks are performed.
26 /// \author Diego Stocco
27 //-----------------------------------------------------------------------------
29 #define AliAnalysisTaskSingleMu_cxx
31 #include "AliAnalysisTaskSingleMu.h"
41 #include "TObjString.h"
42 #include "TObjArray.h"
45 //#include "TMCProcess.h"
48 #include "TPaveStats.h"
51 #include "AliAODEvent.h"
52 #include "AliAODTrack.h"
53 #include "AliAODMCParticle.h"
54 #include "AliMCEvent.h"
55 #include "AliMCParticle.h"
56 #include "AliESDEvent.h"
57 #include "AliESDMuonTrack.h"
58 #include "AliVHeader.h"
59 #include "AliAODMCHeader.h"
63 #include "AliAnalysisManager.h"
66 #include "AliCFContainer.h"
67 #include "AliCFGridSparse.h"
68 #include "AliCFEffGrid.h"
71 #include "AliVAnalysisMuon.h"
72 #include "AliMergeableCollection.h"
73 #include "AliCounterCollection.h"
74 #include "AliMuonTrackCuts.h"
78 ClassImp(AliAnalysisTaskSingleMu) // Class implementation in ROOT context
82 //________________________________________________________________________
83 AliAnalysisTaskSingleMu::AliAnalysisTaskSingleMu() :
85 fMinNvtxContirbutors(0),
91 //________________________________________________________________________
92 AliAnalysisTaskSingleMu::AliAnalysisTaskSingleMu(const char *name, const AliMuonTrackCuts& cuts) :
93 AliVAnalysisMuon(name, cuts),
94 fMinNvtxContirbutors(1),
100 TString thetaAbsKeys = "ThetaAbs23 ThetaAbs310";
101 fThetaAbsKeys = thetaAbsKeys.Tokenize(" ");
105 //________________________________________________________________________
106 AliAnalysisTaskSingleMu::~AliAnalysisTaskSingleMu()
112 delete fThetaAbsKeys;
115 //___________________________________________________________________________
116 void AliAnalysisTaskSingleMu::MyUserCreateOutputObjects()
120 TString histoName = "", histoTitle = "";
123 Double_t vzMin = -20., vzMax = 20.;
124 TString vzName("Vz"), vzTitle("Vz"), vzUnits("cm");
126 histoName = "hIpVtx";
127 histo = new TH1F(histoName.Data(), histoName.Data(), nVzBins, vzMin, vzMax);
128 histo->SetXTitle("v_{z} (cm)");
129 AddObjectToCollection(histo, kIPVz);
133 Double_t ptMin = 0., ptMax = 80.;
134 TString ptName("Pt"), ptTitle("p_{t}"), ptUnits("GeV/c");
137 Double_t etaMin = -4.5, etaMax = -2.;
138 TString etaName("Eta"), etaTitle("#eta"), etaUnits("");
141 Double_t phiMin = 0.; Double_t phiMax = 2*TMath::Pi();
142 TString phiName("Phi"), phiTitle("#phi"), phiUnits("rad");
144 Int_t nChargeBins = 2;
145 Double_t chargeMin = -2, chargeMax = 2.;
146 TString chargeName("Charge"), chargeTitle("charge"), chargeUnits("e");
148 Int_t nThetaAbsEndBins = 2;
149 Double_t thetaAbsEndMin = -0.5, thetaAbsEndMax = 1.5;
150 TString thetaAbsEndName("ThetaAbsEnd"), thetaAbsEndTitle("#theta_{abs}"), thetaAbsEndUnits("a.u.");
152 Int_t nMotherTypeBins = kNtrackSources;
153 Double_t motherTypeMin = -0.5, motherTypeMax = (Double_t)kNtrackSources - 0.5;
154 TString motherType("MotherType"), motherTypeTitle("motherType"), motherTypeUnits("");
156 Int_t nbins[kNvars] = {nPtBins, nEtaBins, nPhiBins, nVzBins, nChargeBins, nThetaAbsEndBins, nMotherTypeBins};
157 Double_t xmin[kNvars] = {ptMin, etaMin, phiMin, vzMin, chargeMin, thetaAbsEndMin, motherTypeMin};
158 Double_t xmax[kNvars] = {ptMax, etaMax, phiMax, vzMax, chargeMax, thetaAbsEndMax, motherTypeMax};
159 TString axisTitle[kNvars] = {ptTitle, etaTitle, phiTitle, vzTitle, chargeTitle, thetaAbsEndTitle, motherTypeTitle};
160 TString axisUnits[kNvars] = {ptUnits, etaUnits, phiUnits, vzUnits, chargeUnits, thetaAbsEndUnits, motherTypeUnits};
162 AliCFContainer* cfContainer = new AliCFContainer("SingleMuContainer","Container for tracks",kNsteps,kNvars,nbins);
164 for ( Int_t idim = 0; idim<kNvars; idim++){
165 histoTitle = Form("%s (%s)", axisTitle[idim].Data(), axisUnits[idim].Data());
166 histoTitle.ReplaceAll("()","");
168 cfContainer->SetVarTitle(idim, histoTitle.Data());
169 cfContainer->SetBinLimits(idim, xmin[idim], xmax[idim]);
172 TString stepTitle[kNsteps] = {"reconstructed", "generated"};
174 TAxis* currAxis = 0x0;
175 for (Int_t istep=0; istep<kNsteps; istep++){
176 cfContainer->SetStepTitle(istep, stepTitle[istep].Data());
177 AliCFGridSparse* gridSparse = cfContainer->GetGrid(istep);
179 currAxis = gridSparse->GetAxis(kHvarMotherType);
180 for ( Int_t ibin=0; ibin<fSrcKeys->GetEntries(); ibin++ ) {
181 currAxis->SetBinLabel(ibin+1, fSrcKeys->At(ibin)->GetName());
185 AddObjectToCollection(cfContainer, kTrackContainer);
187 fMuonTrackCuts->Print("mask");
190 //________________________________________________________________________
191 void AliAnalysisTaskSingleMu::ProcessEvent(TString physSel, const TObjArray& selectTrigClasses, TString centrality)
194 /// Fill output objects
197 AliVVertex* primaryVertex = ( fAODEvent ) ? (AliVVertex*)fAODEvent->GetPrimaryVertexSPD() : (AliVVertex*)fESDEvent->GetPrimaryVertexSPD();
198 if ( primaryVertex->GetNContributors() < fMinNvtxContirbutors ) return;
200 Double_t ipVz = primaryVertex->GetZ();
203 if ( fMCEvent ) ipVzMC = fMCEvent->GetPrimaryVertex()->GetZ();
204 else if ( fAODEvent ) {
205 AliAODMCHeader* aodMCHeader = (AliAODMCHeader *)fAODEvent->FindListObject(AliAODMCHeader::StdBranchName());
206 if ( aodMCHeader ) ipVzMC = aodMCHeader->GetVtxZ();
210 for ( Int_t itrig=0; itrig<selectTrigClasses.GetEntries(); ++itrig ) {
211 TString trigClassName = ((TObjString*)selectTrigClasses.At(itrig))->GetString();
212 ((TH1*)GetMergeableObject(physSel, trigClassName, centrality, "hIpVtx"))->Fill(ipVz);
215 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
216 if ( isPileupFromSPD ) return;
218 Double_t containerInput[kNvars];
219 AliVParticle* track = 0x0;
221 for ( Int_t istep = 0; istep<2; ++istep ) {
222 Int_t nTracks = ( istep == kStepReconstructed ) ? GetNTracks() : GetNMCTracks();
223 for (Int_t itrack = 0; itrack < nTracks; itrack++) {
224 track = ( istep == kStepReconstructed ) ? GetTrack(itrack) : GetMCTrack(itrack);
226 Bool_t isSelected = ( istep == kStepReconstructed ) ? fMuonTrackCuts->IsSelected(track) : ( TMath::Abs(track->PdgCode()) == 13 );
227 if ( ! isSelected ) continue;
229 // In W simulations with Pythia, sometimes muon is stored twice.
230 // Remove muon in case it has another muon as daugther
231 if ( istep == kStepGeneratedMC ) {
232 Int_t firstDaughter = GetDaughterIndex(track, 0);
233 if ( firstDaughter >= 0 ) {
234 Bool_t hasMuonDaughter = kFALSE;
235 Int_t lastDaughter = GetDaughterIndex(track, 1);
236 for ( Int_t idaugh=firstDaughter; idaugh<=lastDaughter; idaugh++ ) {
237 AliVParticle* currTrack = GetMCTrack(idaugh);
238 if ( currTrack->PdgCode() == track->PdgCode() ) {
239 hasMuonDaughter = kTRUE;
243 if ( hasMuonDaughter ) {
244 AliDebug(1, Form("Current muon (%i) has muon daughter: rejecting it", itrack));
250 Int_t trackSrc = ( istep == kStepReconstructed ) ? GetParticleType(track) : RecoTrackMother(track);
252 Double_t thetaAbsEndDeg = 0;
253 if ( istep == kStepReconstructed ) {
254 Double_t rAbsEnd = ( fAODEvent ) ? ((AliAODTrack*)track)->GetRAtAbsorberEnd(): ((AliESDMuonTrack*)track)->GetRAtAbsorberEnd();
255 thetaAbsEndDeg = TMath::ATan( rAbsEnd / 505. ) * TMath::RadToDeg();
258 thetaAbsEndDeg = ( TMath::Pi()-track->Theta() ) * TMath::RadToDeg();
260 Int_t thetaAbsBin = ( thetaAbsEndDeg < 3. ) ? kThetaAbs23 : kThetaAbs310;
262 containerInput[kHvarPt] = track->Pt();
263 containerInput[kHvarEta] = track->Eta();
264 containerInput[kHvarPhi] = track->Phi();
265 containerInput[kHvarVz] = ( istep == kStepReconstructed ) ? ipVz : ipVzMC;
266 containerInput[kHvarCharge] = track->Charge()/3.;
267 containerInput[kHvarThetaAbs] = (Double_t)thetaAbsBin;
268 containerInput[kHvarMotherType] = (Double_t)trackSrc;
270 for ( Int_t itrig=0; itrig<selectTrigClasses.GetEntries(); ++itrig ) {
271 TString trigClassName = ((TObjString*)selectTrigClasses.At(itrig))->GetString();
272 if ( istep == kStepReconstructed && ! TrackPtCutMatchTrigClass(track, trigClassName) ) continue;
273 ((AliCFContainer*)GetMergeableObject(physSel, trigClassName, centrality, "SingleMuContainer"))->Fill(containerInput,istep);
274 } // loop on selected trigger classes
276 } // loop on container steps
280 //________________________________________________________________________
281 void AliAnalysisTaskSingleMu::Terminate(Option_t *) {
283 /// Draw some histograms at the end.
286 AliVAnalysisMuon::Terminate("");
288 if ( ! fMergeableCollection ) return;
290 TString physSel = fTerminateOptions->At(0)->GetName();
291 TString trigClassName = fTerminateOptions->At(1)->GetName();
292 TString centralityRange = fTerminateOptions->At(2)->GetName();
293 TString furtherOpt = fTerminateOptions->At(3)->GetName();
294 furtherOpt.ToUpper();
296 AliCFContainer* cfContainer = static_cast<AliCFContainer*> ( GetSum(physSel,trigClassName,centralityRange,"SingleMuContainer") );
297 if ( ! cfContainer ) return;
299 AliCFEffGrid* effSparse = new AliCFEffGrid(Form("eff%s", cfContainer->GetName()),Form("Efficiency %s", cfContainer->GetTitle()),*cfContainer);
300 effSparse->CalculateEfficiency(kStepReconstructed, kStepGeneratedMC);
302 AliCFGridSparse* gridSparseArray[3] = {effSparse->GetNum(), effSparse->GetDen(), effSparse};
303 TString gridSparseName[3] = {cfContainer->GetStepTitle(kStepReconstructed), cfContainer->GetStepTitle(kStepGeneratedMC), "Efficiency"};
305 Int_t srcColors[kNtrackSources] = {kBlack, kRed, kSpring, kTeal, kBlue, kViolet, kMagenta, kOrange};
306 // TString allSrcNames = "";
307 // for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
308 // if ( ! allSrcNames.IsNull() ) allSrcNames.Append(" ");
309 // allSrcNames += fSrcKeys->At(isrc)->GetName();
318 Bool_t isMC = furtherOpt.Contains("MC");
319 Int_t firstSrc = ( isMC ) ? 0 : kUnidentified;
320 Int_t lastSrc = ( isMC ) ? kNtrackSources - 1 : kUnidentified;
321 if ( ! isMC ) srcColors[kUnidentified] = 1;
323 TString histoName = "", currName = "", histoPattern = "", drawOpt = "";
327 TCanvas* canKine[3] = {0x0, 0x0, 0x0};
328 TLegend* legKine[3] = {0x0, 0x0, 0x0};
329 for ( Int_t isrc = firstSrc; isrc <= lastSrc; ++isrc ) {
330 for ( Int_t icharge=0; icharge<2; ++icharge ) {
331 for ( Int_t igrid=0; igrid<3; ++igrid ) {
332 if ( gridSparseArray[igrid]->GetEntries() == 0. ) break;
333 if ( gridSparseArray[igrid]->IsA() != AliCFEffGrid::Class() ) {
334 SetSparseRange(gridSparseArray[igrid], kHvarEta, "", -3.999, -2.501);
335 SetSparseRange(gridSparseArray[igrid], kHvarMotherType, "", isrc+1, isrc+1, "USEBIN");
336 SetSparseRange(gridSparseArray[igrid], kHvarCharge, "", icharge+1, icharge+1, "USEBIN");
338 if ( ! canKine[igrid] ) {
341 currName = Form("%s_proj_%s", GetName(), gridSparseName[igrid].Data());
342 canKine[igrid] = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
343 canKine[igrid]->Divide(2,2);
344 legKine[igrid] = new TLegend(0.6, 0.6, 0.8, 0.8);
347 for ( Int_t iproj=0; iproj<4; ++iproj ) {
348 canKine[igrid]->cd(iproj+1);
349 if ( ( iproj == kHvarPt || iproj == kHvarVz ) && gridSparseArray[igrid]->IsA() != AliCFEffGrid::Class() ) gPad->SetLogy();
350 TH1* projHisto = gridSparseArray[igrid]->Project(iproj);
351 projHisto->SetName(Form("proj%i_%s_src%i_charge%i", iproj, gridSparseName[igrid].Data(), isrc, icharge));
352 if ( projHisto->GetEntries() == 0 ) continue;
353 Bool_t isFirst = ( gPad->GetListOfPrimitives()->GetEntries() == 0 );
354 drawOpt = isFirst ? "e" : "esames";
355 //if ( isrc == kUnidentified && ! drawOpt.Contains("same") ) isMC = kFALSE;
356 //if ( ! isMC ) srcColors[kUnidentified] = 1;
357 projHisto->SetLineColor(srcColors[isrc]);
358 projHisto->SetMarkerColor(srcColors[isrc]);
359 projHisto->SetMarkerStyle(20+4*icharge);
360 projHisto->Draw(drawOpt.Data());
362 TPaveStats* paveStats = (TPaveStats*)projHisto->FindObject("stats");
363 if ( paveStats ) paveStats->SetTextColor(srcColors[isrc]);
365 TString legEntry = fChargeKeys->At(icharge)->GetName();
366 if ( isMC ) legEntry += Form(" %s", fSrcKeys->At(isrc)->GetName());
367 legKine[igrid]->AddEntry(projHisto,legEntry.Data(), "lp");
369 } // loop on grid sparse
370 } // loop on projections
371 } // loop on mu charge
372 } // loop on track sources
375 for ( Int_t igrid=0; igrid<3; igrid++ ) {
376 if ( ! canKine[igrid] ) continue;
377 canKine[igrid]->cd(1);
378 legKine[igrid]->Draw("same");
379 if ( gridSparseArray[igrid]->IsA() == AliCFEffGrid::Class() ) continue;
380 SetSparseRange(gridSparseArray[igrid], kHvarCharge, "", 1, gridSparseArray[igrid]->GetAxis(kHvarCharge)->GetNbins(), "USEBIN"); // Reset range
381 } // loop on container steps
384 //////////////////////
385 // Event statistics //
386 //////////////////////
387 printf("\nTotal analyzed events:\n");
388 TString evtSel = Form("trigger:%s", trigClassName.Data());
389 fEventCounters->PrintSum(evtSel.Data());
390 printf("Physics selected analyzed events:\n");
391 evtSel = Form("trigger:%s/selected:yes", trigClassName.Data());
392 fEventCounters->PrintSum(evtSel.Data());
394 TString countPhysSel = "any";
395 if ( physSel.Contains(fPhysSelKeys->At(kPhysSelPass)->GetName()) ) countPhysSel = "yes";
396 else if ( physSel.Contains(fPhysSelKeys->At(kPhysSelReject)->GetName()) ) countPhysSel="no";
397 countPhysSel.Prepend("selected:");
398 printf("Analyzed events vs. centrality:\n");
399 evtSel = Form("trigger:%s/%s", trigClassName.Data(), countPhysSel.Data());
400 fEventCounters->Print("centrality",evtSel.Data(),kTRUE);
406 if ( ! furtherOpt.Contains("VERTEX") ) return;
407 Int_t firstMother = kUnidentified, lastMother = kUnidentified;
409 TH1* eventVertex = (TH1*)GetSum(physSel, "CINT7-I-NOPF-ALLNOTRD", centralityRange, "hIpVtx");
410 if ( ! eventVertex ) return;
411 Double_t minZ = -9.99, maxZ = 9.99;
412 Double_t meanZ = 0., sigmaZ = 4.;
413 Double_t nSigma = 2.;
414 TString fitOpt = "R0";
415 Bool_t fixFitRange = kFALSE;
416 TString fitFormula = Form("[0]+[1]*(x+[2])");
419 if ( eventVertex->GetSumw2N() == 0 ) eventVertex->Sumw2();
420 Double_t eventVtxIntegral = eventVertex->Integral(0,eventVertex->GetNbinsX()+1); // Include under/overflow
421 printf("Event vertex integral %.0f\n\n", eventVtxIntegral);
422 if ( eventVtxIntegral <= 0. ) return;
423 eventVertex->Scale(1./eventVtxIntegral);
424 printf("\nFit MB vertex\n");
425 eventVertex->Fit("gaus",fitOpt.Data(),"",minZ,maxZ);
426 TF1* vtxFit = (TF1*)eventVertex->GetListOfFunctions()->FindObject("gaus");
427 currName = "vtxIntegrated";
428 can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
431 vtxFit->Draw("same");
434 enum {kRecoHF, kRecoBkg, kInputHF, kInputDecay, kRecoAll, kNrecoHistos};
435 TString baseRecoName[kNrecoHistos] = {"RecoHF", "RecoBkg", "InputHF", "InputDecay", "RecoAll"};
436 TArrayI sumMothers[kNrecoHistos];
437 sumMothers[kRecoHF].Set(0);
438 sumMothers[kRecoBkg].Set(0);
439 sumMothers[kInputHF].Set(3);
440 sumMothers[kInputHF][0] = kCharmMu;
441 sumMothers[kInputHF][1] = kBeautyMu;
442 sumMothers[kInputHF][2] = kQuarkoniumMu;
443 sumMothers[kInputDecay].Set(1);
444 sumMothers[kInputDecay][0] = kDecayMu;
445 sumMothers[kRecoAll].Set(kNtrackSources);
446 for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
447 sumMothers[kRecoAll][isrc] = isrc;
450 meanZ = vtxFit->GetParameter(1);
451 sigmaZ = vtxFit->GetParameter(2);
453 Double_t minZfit = ( fixFitRange ) ? minZ : meanZ - nSigma*sigmaZ;
454 Double_t maxZfit = ( fixFitRange ) ? maxZ : meanZ + nSigma*sigmaZ;
456 TF1* fitFunc = new TF1("fitFunc", fitFormula.Data(), minZ, maxZ);
457 fitFunc->SetLineColor(2);
458 fitFunc->SetParNames("Line norm", "Line slope", "Free path");
459 const Double_t kFreePath = 153.; // 150.; // 130.; // cm
460 //fitFunc->SetParameters(0.,1.);
461 fitFunc->FixParameter(2, kFreePath);
463 AliCFGridSparse* gridSparse = cfContainer->GetGrid(kStepReconstructed);
464 TAxis* ptAxis = gridSparse->GetAxis(kHvarPt);
467 Double_t limitNorm = 0., limitSlope = 0.;
468 Int_t firstPtBin = 0, lastPtBin = 0;
470 for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) {
472 SetSparseRange(gridSparse, kHvarThetaAbs, "", itheta+1, itheta+1, "USEBIN");
473 SetSparseRange(gridSparse, kHvarPt, "", 1, ptAxis->GetNbins(), "USEBIN");
474 TH1* recoHisto[kNrecoHistos];
475 for ( Int_t ireco=0; ireco<kNrecoHistos; ++ireco ) {
476 recoHisto[ireco] = gridSparse->Project(kHvarPt);
477 histoName = Form("%sMuon_%s", baseRecoName[ireco].Data(), fThetaAbsKeys->At(itheta)->GetName());
478 recoHisto[ireco]->SetName(histoName.Data());
479 recoHisto[ireco]->SetTitle(histoName.Data());
480 recoHisto[ireco]->Reset();
481 recoHisto[ireco]->Sumw2();
482 for ( Int_t isrc=0; isrc<sumMothers[ireco].GetSize(); ++isrc ) {
483 SetSparseRange(gridSparse, kHvarMotherType, "", sumMothers[ireco][isrc]+1, sumMothers[ireco][isrc]+1, "USEBIN");
484 TH1* auxHisto = gridSparse->Project(kHvarPt);
485 recoHisto[ireco]->Add(auxHisto);
489 SetSparseRange(gridSparse, kHvarMotherType, "", firstMother+1, lastMother+1, "USEBIN");
492 for ( Int_t ibinpt=0; ibinpt<=ptAxis->GetNbins(); ++ibinpt ) {
494 lastPtBin = ( ibinpt == 0 ) ? ptAxis->GetNbins() : ibinpt;
495 SetSparseRange(gridSparse, kHvarPt, "", firstPtBin, lastPtBin, "USEBIN");
496 TH1* histo = gridSparse->Project(kHvarVz);
497 histo->SetName(Form("hVtx_%s_%s_ptBin%i", cfContainer->GetStepTitle(kStepReconstructed), fThetaAbsKeys->At(itheta)->GetName(), ibinpt));
498 if ( histo->GetEntries() < 100. ) break;
499 printf("\nFit %.2f < pt < %.2f (entries %g)\n", ptAxis->GetBinLowEdge(firstPtBin), ptAxis->GetBinUpEdge(lastPtBin), histo->GetEntries());
500 histo->Divide(eventVertex);
501 Double_t norm = histo->GetBinContent(histo->FindBin(0.));
502 histo->GetYaxis()->SetTitle("#frac{dN_{#mu}}{dv_{z}} / #left(#frac{1}{N_{MB}}#frac{dN_{MB}}{dv_{z}}#right)");
503 slope = ( histo->GetBinContent(histo->FindBin(meanZ+sigmaZ)) -
504 histo->GetBinContent(histo->FindBin(meanZ-sigmaZ)) ) / ( 2. * sigmaZ );
506 if ( slope < 0. ) slope = norm/kFreePath;
508 // Try to fit twice: it fit fails the first time
509 // set some limits on parameters
510 for ( Int_t itry=0; itry<2; itry++ ) {
511 fitFunc->SetParameter(0, norm);
512 fitFunc->SetParameter(1, slope);
514 limitNorm = 2.*histo->Integral();
515 limitSlope = 2.*histo->Integral()/kFreePath;
516 //fitFunc->SetParLimits(0, 0., limitNorm); // REMEMBER TO CHECK
517 fitFunc->SetParLimits(1, 0., limitSlope); // REMEMBER TO CHECK
518 printf("Norm 0. < %f < %f slope 0. < %f < %f\n", norm, limitNorm, slope, limitSlope);
520 histo->Fit(fitFunc, fitOpt.Data(), "", minZfit, maxZfit);
522 if ( gMinuit->fCstatu.Contains("CONVERGED") &&
523 fitFunc->GetParameter(0) > 0. &&
524 fitFunc->GetParameter(1) > 0. )
526 else if ( furtherOpt.Contains("REFIT") ) printf("Re-fit with limits\n");
528 printf("Warning: fit problems !!!\n");
533 Double_t p0 = fitFunc->GetParameter(0);
534 Double_t p0err = fitFunc->GetParError(0);
535 Double_t p1 = fitFunc->GetParameter(1);
536 Double_t p1err = fitFunc->GetParError(1);
538 Double_t nullVz = ( p1 != 0. ) ? -p0/p1 : 0.;
539 Double_t nullVzErr = ( p0 != 0. && p1 != 0. ) ? TMath::Abs(nullVz) * TMath::Sqrt(p0err*p0err/(p0*p0) + p1err*p1err/(p1*p1) ) : 0.;
541 printf("Null value at %f +- %f\n", nullVz - kFreePath, nullVzErr);
543 recoHisto[kRecoHF]->SetBinContent(ibinpt, p0);
544 recoHisto[kRecoHF]->SetBinError(ibinpt, p0err);
545 recoHisto[kRecoBkg]->SetBinContent(ibinpt, ( kFreePath + meanZ ) * p1);
546 recoHisto[kRecoBkg]->SetBinError(ibinpt, ( kFreePath + meanZ ) * p1err);
547 if ( currDraw%4 == 0 ){
548 currName = Form("vtx_%s_PtBin%i",fThetaAbsKeys->At(itheta)->GetName(), ibinpt);
549 can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
552 can->cd( currDraw%4 + 1 );
555 fitFunc->DrawCopy("same");
558 SetSparseRange(gridSparse, kHvarMotherType, "", firstMother+1, lastMother+1, "USEBIN");
559 TH1* totalPtDistrib = gridSparse->Project(kHvarPt);
560 totalPtDistrib->SetName(Form("totalPtDistrib_%s", fThetaAbsKeys->GetName()));
561 currName = Form("recoPt_%s",fThetaAbsKeys->At(itheta)->GetName());
562 can = new TCanvas(currName.Data(),currName.Data(),(igroup1+1)*xshift,igroup2*yshift,600,600);
563 TLegend* leg = new TLegend(0.6, 0.6, 0.8, 0.8);
565 for ( Int_t ireco=0; ireco<kNrecoHistos-1; ++ireco ) {
566 if ( recoHisto[ireco]->GetEntries() == 0. ) continue;
567 TH1* ratio = (TH1*)recoHisto[ireco]->Clone(Form("%s_ratio", recoHisto[ireco]->GetName()));
568 ratio->Divide(recoHisto[kRecoAll]);
569 ratio->SetLineColor(srcColors[ireco]);
570 ratio->SetMarkerColor(srcColors[ireco]);
571 ratio->SetMarkerStyle(20+ireco);
572 ratio->GetYaxis()->SetTitle("fraction of total");
573 ratio->Draw(drawOpt.Data());
574 leg->AddEntry(ratio,baseRecoName[ireco].Data(), "lp");
578 } // loop on theta abs