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 /// If Monte Carlo information is present, some basics checks are performed.
25 /// \author Diego Stocco
26 //-----------------------------------------------------------------------------
28 #define AliAnalysisTaskSingleMu_cxx
30 #include "AliAnalysisTaskSingleMu.h"
40 #include "TObjString.h"
41 #include "TObjArray.h"
44 //#include "TMCProcess.h"
46 #include "TPaveStats.h"
47 #include "TFitResultPtr.h"
49 #include "THashList.h"
52 #include "AliAODEvent.h"
53 #include "AliAODTrack.h"
54 #include "AliAODMCParticle.h"
55 #include "AliMCEvent.h"
56 #include "AliMCParticle.h"
57 #include "AliESDEvent.h"
58 #include "AliESDMuonTrack.h"
59 #include "AliVHeader.h"
60 #include "AliAODMCHeader.h"
64 #include "AliAnalysisManager.h"
67 #include "AliCFContainer.h"
68 #include "AliCFGridSparse.h"
69 #include "AliCFEffGrid.h"
72 #include "AliVAnalysisMuon.h"
73 #include "AliMergeableCollection.h"
74 #include "AliCounterCollection.h"
75 #include "AliMuonEventCuts.h"
76 #include "AliMuonTrackCuts.h"
77 #include "AliAnalysisMuonUtility.h"
81 ClassImp(AliAnalysisTaskSingleMu) // Class implementation in ROOT context
85 //________________________________________________________________________
86 AliAnalysisTaskSingleMu::AliAnalysisTaskSingleMu() :
94 //________________________________________________________________________
95 AliAnalysisTaskSingleMu::AliAnalysisTaskSingleMu(const char *name, const AliMuonTrackCuts& cuts) :
96 AliVAnalysisMuon(name, cuts),
103 TString thetaAbsKeys = "ThetaAbs23 ThetaAbs310";
104 fThetaAbsKeys = thetaAbsKeys.Tokenize(" ");
108 //________________________________________________________________________
109 AliAnalysisTaskSingleMu::~AliAnalysisTaskSingleMu()
115 delete fThetaAbsKeys;
118 //___________________________________________________________________________
119 void AliAnalysisTaskSingleMu::MyUserCreateOutputObjects()
123 TString histoName = "", histoTitle = "";
126 Double_t vzMin = -20., vzMax = 20.;
127 TString vzName("Vz"), vzTitle("Vz"), vzUnits("cm");
129 histoName = "hIpVtx";
130 histo = new TH1F(histoName.Data(), histoName.Data(), nVzBins, vzMin, vzMax);
131 histo->SetXTitle("v_{z} (cm)");
132 AddObjectToCollection(histo, kIPVz);
136 Double_t ptMin = 0., ptMax = 80.;
137 TString ptName("Pt"), ptTitle("p_{T}"), ptUnits("GeV/c");
140 Double_t etaMin = -4.5, etaMax = -2.;
141 TString etaName("Eta"), etaTitle("#eta"), etaUnits("");
144 Double_t phiMin = 0.; Double_t phiMax = 2.*TMath::Pi();
145 TString phiName("Phi"), phiTitle("#phi"), phiUnits("rad");
147 Int_t nChargeBins = 2;
148 Double_t chargeMin = -2, chargeMax = 2.;
149 TString chargeName("Charge"), chargeTitle("charge"), chargeUnits("e");
151 Int_t nThetaAbsEndBins = 2;
152 Double_t thetaAbsEndMin = -0.5, thetaAbsEndMax = 1.5;
153 TString thetaAbsEndName("ThetaAbsEnd"), thetaAbsEndTitle("#theta_{abs}"), thetaAbsEndUnits("a.u.");
155 Int_t nMotherTypeBins = kNtrackSources;
156 Double_t motherTypeMin = -0.5, motherTypeMax = (Double_t)kNtrackSources - 0.5;
157 TString motherType("MotherType"), motherTypeTitle("motherType"), motherTypeUnits("");
159 Int_t nbins[kNvars] = {nPtBins, nEtaBins, nPhiBins, nVzBins, nChargeBins, nThetaAbsEndBins, nMotherTypeBins};
160 Double_t xmin[kNvars] = {ptMin, etaMin, phiMin, vzMin, chargeMin, thetaAbsEndMin, motherTypeMin};
161 Double_t xmax[kNvars] = {ptMax, etaMax, phiMax, vzMax, chargeMax, thetaAbsEndMax, motherTypeMax};
162 TString axisTitle[kNvars] = {ptTitle, etaTitle, phiTitle, vzTitle, chargeTitle, thetaAbsEndTitle, motherTypeTitle};
163 TString axisUnits[kNvars] = {ptUnits, etaUnits, phiUnits, vzUnits, chargeUnits, thetaAbsEndUnits, motherTypeUnits};
165 AliCFContainer* cfContainer = new AliCFContainer("SingleMuContainer","Container for tracks",kNsteps,kNvars,nbins);
167 for ( Int_t idim = 0; idim<kNvars; idim++){
168 histoTitle = Form("%s (%s)", axisTitle[idim].Data(), axisUnits[idim].Data());
169 histoTitle.ReplaceAll("()","");
171 cfContainer->SetVarTitle(idim, histoTitle.Data());
172 cfContainer->SetBinLimits(idim, xmin[idim], xmax[idim]);
175 TString stepTitle[kNsteps] = {"reconstructed", "generated"};
177 TAxis* currAxis = 0x0;
178 for (Int_t istep=0; istep<kNsteps; istep++){
179 cfContainer->SetStepTitle(istep, stepTitle[istep].Data());
180 AliCFGridSparse* gridSparse = cfContainer->GetGrid(istep);
182 currAxis = gridSparse->GetAxis(kHvarMotherType);
183 for ( Int_t ibin=0; ibin<fSrcKeys->GetEntries(); ibin++ ) {
184 currAxis->SetBinLabel(ibin+1, fSrcKeys->At(ibin)->GetName());
188 AddObjectToCollection(cfContainer, kTrackContainer);
190 histoName = "hRecoDimu";
191 TH2* histoDimu = new TH2F(histoName.Data(), histoName.Data(), 24, 0., 120., 30, 0., 150.);
192 histoDimu->SetXTitle("min. p_{T}^{#mu} (GeV/c)");
193 histoDimu->SetYTitle("M_{#mu#mu} (GeV/c^{2})");
194 AddObjectToCollection(histoDimu, kNobjectTypes);
196 histoName = "hGeneratedZ";
197 AddObjectToCollection(histoDimu->Clone(histoName.Data()), kNobjectTypes+1);
200 histoName = "hZmuEtaCorr";
201 histoDimu = new TH2F(histoName.Data(), histoName.Data(), 160, -8., 8., 160,-8., 8.);
202 histoDimu->SetXTitle("#eta_{#mu-}");
203 histoDimu->SetYTitle("#eta_{#mu+}");
204 AddObjectToCollection(histoDimu, kNobjectTypes+2);
206 fMuonTrackCuts->Print("mask");
208 AliInfo(Form("Apply cut on dimuon (60<M_mumu<120 GeV/c^2) to reject Z contribution: %i", fCutOnDimu));
211 //________________________________________________________________________
212 void AliAnalysisTaskSingleMu::ProcessEvent(TString physSel, const TObjArray& selectTrigClasses, TString centrality)
215 /// Fill output objects
218 Double_t ipVz = AliAnalysisMuonUtility::GetVertexSPD(InputEvent())->GetZ();
219 Double_t ipVzMC = MCEvent() ? AliAnalysisMuonUtility::GetMCVertexZ(InputEvent(),MCEvent()) : 0.;
221 for ( Int_t itrig=0; itrig<selectTrigClasses.GetEntries(); ++itrig ) {
222 TString trigClassName = ((TObjString*)selectTrigClasses.At(itrig))->GetString();
223 ((TH1*)GetMergeableObject(physSel, trigClassName, centrality, "hIpVtx"))->Fill(ipVz);
226 // 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
227 // if ( isPileupFromSPD ) return; // UNCOMMENT TO REJECT PILEUP
229 Double_t containerInput[kNvars];
230 AliVParticle* track = 0x0;
232 Int_t nSteps = MCEvent() ? 2 : 1;
233 for ( Int_t istep = 0; istep<nSteps; ++istep ) {
234 Int_t nTracks = ( istep == kStepReconstructed ) ? AliAnalysisMuonUtility::GetNTracks(InputEvent()) : MCEvent()->GetNumberOfTracks();
236 TObjArray selectedTracks(nTracks);
237 TArrayI trackSources(nTracks);
238 for (Int_t itrack = 0; itrack < nTracks; itrack++) {
239 track = ( istep == kStepReconstructed ) ? AliAnalysisMuonUtility::GetTrack(itrack,InputEvent()) : MCEvent()->GetTrack(itrack);
241 // In case of MC we usually ask that the particle is a muon
242 // However, in W or Z simulations, Pythia stores both the initial muon
243 // (before ISR, FSR and kt kick) and the final state one.
244 // The first muon is of course there only for information and should be rejected.
245 // The Pythia code for initial state particles is 21
246 Bool_t isSelected = ( istep == kStepReconstructed ) ? fMuonTrackCuts->IsSelected(track) : ( TMath::Abs(track->PdgCode()) == 13 && AliAnalysisMuonUtility::GetStatusCode(track) != 21 );
247 if ( ! isSelected ) continue;
249 selectedTracks.AddAt(track,itrack);
250 trackSources[itrack] = GetParticleType(track);
254 // Loop on selected tracks
255 TArrayI rejectTrack(nTracks);
256 for ( Int_t itrack=0; itrack<nTracks; itrack++) {
257 track = static_cast<AliVParticle*>(selectedTracks.At(itrack));
258 if ( ! track ) continue;
260 if ( istep == kStepGeneratedMC || fCutOnDimu ) {
262 for ( Int_t jtrack=itrack+1; jtrack<nTracks; jtrack++ ) {
263 AliVParticle* auxTrack = static_cast<AliVParticle*>(selectedTracks.At(jtrack));
264 if ( ! auxTrack ) continue;
265 if ( track->Charge() * auxTrack->Charge() >= 0 ) continue;
267 TLorentzVector dimuPair = AliAnalysisMuonUtility::GetTrackPair(track,auxTrack);
268 Double_t ptMin = TMath::Min(track->Pt(),auxTrack->Pt());
269 Double_t invMass = dimuPair.M();
270 if ( invMass > 60. && invMass < 120. ) {
271 rejectTrack[itrack] = 1;
272 rejectTrack[jtrack] = 1;
274 if ( istep == kStepReconstructed ) {
275 for ( Int_t itrig=0; itrig<selectTrigClasses.GetEntries(); ++itrig ) {
276 TString trigClassName = ((TObjString*)selectTrigClasses.At(itrig))->GetString();
277 if ( ! fMuonTrackCuts->TrackPtCutMatchTrigClass(track, fMuonEventCuts->GetTrigClassPtCutLevel(trigClassName)) || ! fMuonTrackCuts->TrackPtCutMatchTrigClass(auxTrack, fMuonEventCuts->GetTrigClassPtCutLevel(trigClassName)) ) continue;
278 ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, "hRecoDimu"))->Fill(ptMin,invMass);
279 } // loop on triggers
282 if ( trackSources[itrack] == kZbosonMu && trackSources[jtrack] == kZbosonMu ) {
283 Bool_t isAccepted = kTRUE;
284 AliVParticle* muDaughter[2] = {0x0, 0x0};
285 for ( Int_t imu=0; imu<2; imu++ ) {
286 AliVParticle* currPart = ( imu == 0 ) ? track : auxTrack;
287 if ( currPart->Charge() < 0. ) muDaughter[0] = currPart;
288 else muDaughter[1] = currPart;
289 if ( currPart->Eta() < -4.5 || currPart->Eta() > -2. ) {
292 } // loop on muons in the pair
294 Double_t pairRapidity = dimuPair.Rapidity();
295 if ( pairRapidity < -4. || pairRapidity > -2.5 ) isAccepted = kFALSE;
296 //printf("Rapidity Z %g pair %g\n",track->Y(), pairRapidity);
298 for ( Int_t itrig=0; itrig<selectTrigClasses.GetEntries(); ++itrig ) {
299 TString trigClassName = ((TObjString*)selectTrigClasses.At(itrig))->GetString();
300 if ( isAccepted ) ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, "hGeneratedZ"))->Fill(ptMin,invMass);
301 ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, "hZmuEtaCorr"))->Fill(muDaughter[0]->Eta(),muDaughter[1]->Eta());
302 } // loop on selected trig
303 } // both muons from Z
304 } // kStepGeneratedMC
305 } // loop on auxiliary tracks
306 } // apply cut on dimu
307 if ( rejectTrack[itrack] > 0 ) continue;
309 Double_t thetaAbsEndDeg = 0;
310 if ( istep == kStepReconstructed ) {
311 thetaAbsEndDeg = AliAnalysisMuonUtility::GetThetaAbsDeg(track);
314 thetaAbsEndDeg = ( TMath::Pi()-track->Theta() ) * TMath::RadToDeg();
316 Int_t thetaAbsBin = ( thetaAbsEndDeg < 3. ) ? kThetaAbs23 : kThetaAbs310;
318 containerInput[kHvarPt] = track->Pt();
319 containerInput[kHvarEta] = track->Eta();
320 containerInput[kHvarPhi] = track->Phi();
321 containerInput[kHvarVz] = ( istep == kStepReconstructed ) ? ipVz : ipVzMC;
322 containerInput[kHvarCharge] = track->Charge()/3.;
323 containerInput[kHvarThetaAbs] = (Double_t)thetaAbsBin;
324 containerInput[kHvarMotherType] = (Double_t)trackSources[itrack];
326 for ( Int_t itrig=0; itrig<selectTrigClasses.GetEntries(); ++itrig ) {
327 TString trigClassName = ((TObjString*)selectTrigClasses.At(itrig))->GetString();
328 if ( istep == kStepReconstructed && ! fMuonTrackCuts->TrackPtCutMatchTrigClass(track, fMuonEventCuts->GetTrigClassPtCutLevel(trigClassName)) ) continue;
329 ((AliCFContainer*)GetMergeableObject(physSel, trigClassName, centrality, "SingleMuContainer"))->Fill(containerInput,istep);
330 } // loop on selected trigger classes
332 } // loop on container steps
336 //________________________________________________________________________
337 void AliAnalysisTaskSingleMu::Terminate(Option_t *) {
339 /// Draw some histograms at the end.
342 AliVAnalysisMuon::Terminate("");
344 if ( ! fMergeableCollection ) return;
346 TString physSel = fTerminateOptions->At(0)->GetName();
347 TString trigClassName = fTerminateOptions->At(1)->GetName();
348 TString centralityRange = fTerminateOptions->At(2)->GetName();
349 TString furtherOpt = fTerminateOptions->At(3)->GetName();
351 TString minBiasTrig = "";
352 TObjArray* optArr = furtherOpt.Tokenize(" ");
353 TString currName = "";
354 for ( Int_t iopt=0; iopt<optArr->GetEntries(); iopt++ ) {
355 currName = optArr->At(iopt)->GetName();
356 if ( currName.Contains("-B-") || currName.Contains("ANY") ) minBiasTrig = currName;
360 furtherOpt.ToUpper();
361 Bool_t plotChargeAsymmetry = furtherOpt.Contains("ASYM");
363 AliCFContainer* cfContainer = static_cast<AliCFContainer*> ( GetSum(physSel,trigClassName,centralityRange,"SingleMuContainer") );
364 if ( ! cfContainer ) return;
366 AliCFEffGrid* effSparse = new AliCFEffGrid(Form("eff%s", cfContainer->GetName()),Form("Efficiency %s", cfContainer->GetTitle()),*cfContainer);
367 effSparse->CalculateEfficiency(kStepReconstructed, kStepGeneratedMC);
369 AliCFGridSparse* gridSparseArray[3] = {effSparse->GetNum(), effSparse->GetDen(), effSparse};
370 TString gridSparseName[3] = {cfContainer->GetStepTitle(kStepReconstructed), cfContainer->GetStepTitle(kStepGeneratedMC), "Efficiency"};
372 Int_t srcColors[kNtrackSources] = {kBlack, kRed, kSpring, kTeal, kBlue, kViolet, kMagenta, kOrange, kGray};
373 // TString allSrcNames = "";
374 // for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
375 // if ( ! allSrcNames.IsNull() ) allSrcNames.Append(" ");
376 // allSrcNames += fSrcKeys->At(isrc)->GetName();
385 Bool_t isMC = furtherOpt.Contains("MC");
387 TAxis* srcAxis = gridSparseArray[0]->GetAxis(kHvarMotherType);
388 Int_t unIdBin = srcAxis->GetNbins();
389 for ( Int_t ibin=1; ibin<=srcAxis->GetNbins(); ibin++ ) {
390 currName = srcAxis->GetBinLabel(ibin);
391 if ( currName.Contains("Unidentified") ) unIdBin = ibin;
394 Int_t firstSrcBin = ( isMC ) ? 1 : unIdBin;
395 Int_t lastSrcBin = ( isMC ) ? srcAxis->GetNbins() - 1 : unIdBin;
396 if ( ! isMC ) srcColors[unIdBin-1] = 1;
401 TString chargeNames[3] = {fChargeKeys->At(0)->GetName(), fChargeKeys->At(1)->GetName(), "Total"};
402 THashList histoList[3];
403 for ( Int_t icharge=0; icharge<3; icharge++ ) {
404 histoList[icharge].SetName(chargeNames[icharge].Data());
406 for ( Int_t isrc = firstSrcBin; isrc <= lastSrcBin; ++isrc ) {
407 for ( Int_t icharge=0; icharge<3; ++icharge ) {
408 Int_t icharge1 = ( icharge < 2 ) ? icharge : 0;
409 Int_t icharge2 = ( icharge < 2 ) ? icharge : 1;
410 for ( Int_t igrid=0; igrid<3; ++igrid ) {
411 if ( gridSparseArray[igrid]->GetEntries() == 0. ) break;
412 if ( gridSparseArray[igrid]->IsA() != AliCFEffGrid::Class() ) {
413 SetSparseRange(gridSparseArray[igrid], kHvarEta, "", -3.999, -2.501);
414 SetSparseRange(gridSparseArray[igrid], kHvarMotherType, "", isrc, isrc, "USEBIN");
415 SetSparseRange(gridSparseArray[igrid], kHvarCharge, "", icharge1+1, icharge2+1, "USEBIN");
417 TH1* histo = gridSparseArray[igrid]->Project(kHvarPt, kHvarEta);
418 histo->SetName(Form("hPtEta_%s_%s_%s", gridSparseName[igrid].Data(), srcAxis->GetBinLabel(isrc), chargeNames[icharge].Data()));
419 if ( histo->Integral() > 0 ) histoList[icharge].Add(histo);
420 for ( Int_t iproj=0; iproj<4; ++iproj ) {
421 histo = gridSparseArray[igrid]->Project(iproj);
422 histo->SetName(Form("proj%i_%s_%s_%s", iproj, gridSparseName[igrid].Data(), srcAxis->GetBinLabel(isrc), chargeNames[icharge].Data()));
423 if ( histo->Integral() > 0 ) histoList[icharge].Add(histo);
424 } // loop on projections
425 } // loop on grid sparse
427 } // loop on track sources
429 // Get charge asymmetry or mu+/mu-
430 THashList histoListRatio;
431 TString basePlotName = plotChargeAsymmetry ? "ChargeAsymmetry" : "ChargeRatio";
432 histoListRatio.SetName(basePlotName.Data());
433 Int_t baseCharge = 1;
434 Int_t auxCharge = 1-baseCharge;
435 for ( Int_t ihisto=0; ihisto<histoList[baseCharge].GetEntries(); ihisto++ ) {
436 TObject* obj = histoList[baseCharge].At(ihisto);
437 TString histoName = obj->GetName();
438 if ( histoName.Contains(gridSparseName[2].Data()) ) continue;
439 TString searchName = histoName;
440 searchName.ReplaceAll(fChargeKeys->At(baseCharge)->GetName(), fChargeKeys->At(auxCharge)->GetName());
441 TH1* auxHisto = static_cast<TH1*> (histoList[auxCharge].FindObject(searchName.Data()));
442 if ( ! auxHisto ) continue;
443 histoName.ReplaceAll(fChargeKeys->At(baseCharge)->GetName(),basePlotName.Data());
444 TH1* histo = static_cast<TH1*> (obj->Clone(histoName.Data()));
445 if ( plotChargeAsymmetry ) {
446 histo->Add(auxHisto, -1.);
447 // h2 + h1 = 2xh2 + (h1-h2)
448 auxHisto->Add(auxHisto, histo, 2.);
450 histo->Divide(auxHisto);
451 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());
452 axisTitle.ReplaceAll("MuPlus","#mu^{+}");
453 axisTitle.ReplaceAll("MuMinus","#mu^{-}");
454 histo->GetYaxis()->SetTitle(axisTitle.Data());
455 histo->SetStats(kFALSE);
456 histoListRatio.Add(histo);
460 TString histoName = "", drawOpt = "";
461 for ( Int_t itype=0; itype<3; itype++ ) {
462 THashList* currList = 0x0;
464 if ( itype == 1 ) currList = &histoListRatio;
465 else if ( itype == 2 ) currList = &histoList[2];
467 for ( Int_t igrid=0; igrid<3; ++igrid ) {
469 TCanvas* canKine = 0x0;
470 TLegend* legKine = 0x0;
471 for ( Int_t iproj=0; iproj<4; ++iproj ) {
472 for ( Int_t isrc = firstSrcBin; isrc <= lastSrcBin; ++isrc ) {
473 for ( Int_t icharge=0; icharge<nCharges; ++icharge ) {
474 if ( itype == 0 ) currList = &histoList[icharge];
475 histoName = Form("proj%i_%s_%s_%s", iproj, gridSparseName[igrid].Data(), srcAxis->GetBinLabel(isrc), currList->GetName());
476 TH1* histo = static_cast<TH1*>(currList->FindObject(histoName.Data()));
477 if ( ! histo ) continue;
481 currName = Form("%s_%s_%s", GetName(), currList->GetName(), gridSparseName[igrid].Data());
482 canKine = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
483 canKine->Divide(2,2);
484 legKine = new TLegend(0.6, 0.6, 0.8, 0.8);
486 canKine->cd(iproj+1);
488 if ( ( iproj == kHvarPt || iproj == kHvarVz ) && gridSparseArray[igrid]->IsA() != AliCFEffGrid::Class() ) gPad->SetLogy();
490 Bool_t isFirst = ( gPad->GetListOfPrimitives()->GetEntries() == 0 );
491 drawOpt = isFirst ? "e" : "esames";
492 histo->SetLineColor(srcColors[isrc-1]);
493 histo->SetMarkerColor(srcColors[isrc-1]);
494 histo->SetMarkerStyle(20+4*icharge);
495 histo->Draw(drawOpt.Data());
496 TPaveStats* paveStats = (TPaveStats*)histo->FindObject("stats");
497 if ( paveStats ) paveStats->SetTextColor(srcColors[isrc-1]);
499 TString legEntry = ( itype == 0 ) ? fChargeKeys->At(icharge)->GetName() : "";
500 if ( isMC ) legEntry += Form(" %s", srcAxis->GetBinLabel(isrc));
501 if ( ! legEntry.IsNull() ) legKine->AddEntry(histo,legEntry.Data(), "lp");
503 } // loop on mu charge
504 } // loop on track sources
505 } // loop on projections
506 if ( legKine && legKine->GetNRows() > 0 ) {
508 legKine->Draw("same");
510 } // loop on grid sparse
514 for ( Int_t igrid=0; igrid<3; igrid++ ) {
515 if ( gridSparseArray[igrid]->IsA() == AliCFEffGrid::Class() ) continue;
516 SetSparseRange(gridSparseArray[igrid], kHvarCharge, "", 1, gridSparseArray[igrid]->GetAxis(kHvarCharge)->GetNbins(), "USEBIN"); // Reset range
517 } // loop on container steps
519 //////////////////////
520 // Event statistics //
521 //////////////////////
522 printf("\nTotal analyzed events:\n");
523 TString evtSel = Form("trigger:%s", trigClassName.Data());
524 fEventCounters->PrintSum(evtSel.Data());
525 printf("Physics selected analyzed events:\n");
526 evtSel = Form("trigger:%s/selected:yes", trigClassName.Data());
527 fEventCounters->PrintSum(evtSel.Data());
529 TString countPhysSel = "any";
530 if ( physSel.Contains(fPhysSelKeys->At(kPhysSelPass)->GetName()) ) countPhysSel = "yes";
531 else if ( physSel.Contains(fPhysSelKeys->At(kPhysSelReject)->GetName()) ) countPhysSel="no";
532 countPhysSel.Prepend("selected:");
533 printf("Analyzed events vs. centrality:\n");
534 evtSel = Form("trigger:%s/%s", trigClassName.Data(), countPhysSel.Data());
535 fEventCounters->Print("centrality",evtSel.Data(),kTRUE);
538 TString outFilename = Form("/tmp/out%s.root", GetName());
539 TFile* outFile = new TFile(outFilename.Data(),"RECREATE");
540 for ( Int_t icharge=0; icharge<3; icharge++ ) {
541 histoList[icharge].Write();
543 histoListRatio.Write();
545 printf("\nCreating file %s\n", outFilename.Data(
551 if ( ! furtherOpt.Contains("VERTEX") ) return;
552 Int_t firstMother = ( isMC ) ? 0 : kUnidentified;
553 Int_t lastMother = ( isMC ) ? kNtrackSources - 1 : kUnidentified;
555 TH1* eventVertex = (TH1*)GetSum(physSel, minBiasTrig, centralityRange, "hIpVtx");
556 if ( ! eventVertex ) return;
557 Double_t minZ = -9.99, maxZ = 9.99;
558 Double_t meanZ = 0., sigmaZ = 4.;
559 Double_t nSigma = 2.;
560 TString fitOpt = "R0S";
561 Bool_t fixFitRange = kFALSE;
562 TString fitFormula = Form("[0]+[1]*(x+[2])");
565 if ( eventVertex->GetSumw2N() == 0 ) eventVertex->Sumw2();
566 Double_t eventVtxIntegral = eventVertex->Integral(0,eventVertex->GetNbinsX()+1); // Include under/overflow
567 printf("Event vertex integral %.0f\n\n", eventVtxIntegral);
568 if ( eventVtxIntegral <= 0. ) return;
569 eventVertex->Scale(1./eventVtxIntegral);
570 printf("\nFit MB vertex\n");
571 eventVertex->Fit("gaus",fitOpt.Data(),"",minZ,maxZ);
572 TF1* vtxFit = (TF1*)eventVertex->GetListOfFunctions()->FindObject("gaus");
573 currName = "vtxIntegrated";
574 can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
577 vtxFit->Draw("same");
580 enum {kRecoHF, kRecoBkg, kInputHF, kInputDecay, kRecoAll, kNrecoHistos};
581 TString baseRecoName[kNrecoHistos] = {"RecoHF", "RecoBkg", "InputHF", "InputDecay", "RecoAll"};
582 TArrayI sumMothers[kNrecoHistos];
583 sumMothers[kRecoHF].Set(0);
584 sumMothers[kRecoBkg].Set(0);
585 sumMothers[kInputHF].Set(3);
586 sumMothers[kInputHF][0] = kCharmMu;
587 sumMothers[kInputHF][1] = kBeautyMu;
588 sumMothers[kInputHF][2] = kQuarkoniumMu;
589 sumMothers[kInputDecay].Set(1);
590 sumMothers[kInputDecay][0] = kDecayMu;
591 sumMothers[kRecoAll].Set(srcAxis->GetNbins());
592 for ( Int_t isrc=1; isrc<srcAxis->GetNbins(); ++isrc ) {
593 sumMothers[kRecoAll][isrc-1] = isrc;
596 meanZ = vtxFit->GetParameter(1);
597 sigmaZ = vtxFit->GetParameter(2);
599 Double_t minZfit = ( fixFitRange ) ? minZ : meanZ - nSigma*sigmaZ;
600 Double_t maxZfit = ( fixFitRange ) ? maxZ : meanZ + nSigma*sigmaZ;
602 TF1* fitFunc = new TF1("fitFunc", fitFormula.Data(), minZ, maxZ);
603 fitFunc->SetLineColor(2);
604 fitFunc->SetParNames("Line norm", "Line slope", "Free path");
605 const Double_t kFreePath = 153.; // 150.; // 130.; // cm
606 //fitFunc->SetParameters(0.,1.);
607 fitFunc->FixParameter(2, kFreePath);
609 AliCFGridSparse* gridSparse = cfContainer->GetGrid(kStepReconstructed);
610 TAxis* ptAxis = gridSparse->GetAxis(kHvarPt);
613 Double_t limitNorm = 0., limitSlope = 0.;
614 Int_t firstPtBin = 0, lastPtBin = 0;
616 gStyle->SetOptFit(1111);
618 for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) {
620 SetSparseRange(gridSparse, kHvarThetaAbs, "", itheta+1, itheta+1, "USEBIN");
621 SetSparseRange(gridSparse, kHvarPt, "", 1, ptAxis->GetNbins(), "USEBIN");
622 TH1* recoHisto[kNrecoHistos];
623 for ( Int_t ireco=0; ireco<kNrecoHistos; ++ireco ) {
624 recoHisto[ireco] = gridSparse->Project(kHvarPt);
625 histoName = Form("%sMuon_%s", baseRecoName[ireco].Data(), fThetaAbsKeys->At(itheta)->GetName());
626 recoHisto[ireco]->SetName(histoName.Data());
627 recoHisto[ireco]->SetTitle(histoName.Data());
628 recoHisto[ireco]->Reset();
629 recoHisto[ireco]->Sumw2();
630 for ( Int_t isrc=0; isrc<sumMothers[ireco].GetSize(); ++isrc ) {
631 SetSparseRange(gridSparse, kHvarMotherType, "", sumMothers[ireco][isrc], sumMothers[ireco][isrc], "USEBIN");
632 TH1* auxHisto = gridSparse->Project(kHvarPt);
633 recoHisto[ireco]->Add(auxHisto);
637 SetSparseRange(gridSparse, kHvarMotherType, "", firstMother+1, lastMother+1, "USEBIN");
640 for ( Int_t ibinpt=0; ibinpt<=ptAxis->GetNbins(); ++ibinpt ) {
642 lastPtBin = ( ibinpt == 0 ) ? ptAxis->GetNbins() : ibinpt;
643 SetSparseRange(gridSparse, kHvarPt, "", firstPtBin, lastPtBin, "USEBIN");
644 TH1* histo = gridSparse->Project(kHvarVz);
645 histo->SetName(Form("hVtx_%s_%s_ptBin%i", cfContainer->GetStepTitle(kStepReconstructed), fThetaAbsKeys->At(itheta)->GetName(), ibinpt));
646 if ( histo->Integral() < 100. ) break;
647 printf("\nFit %.2f < pt < %.2f (entries %g)\n", ptAxis->GetBinLowEdge(firstPtBin), ptAxis->GetBinUpEdge(lastPtBin), histo->GetEntries());
648 histo->Divide(eventVertex);
649 Double_t norm = histo->GetBinContent(histo->FindBin(0.));
650 histo->GetYaxis()->SetTitle("#frac{dN_{#mu}}{dv_{z}} / #left(#frac{1}{N_{MB}}#frac{dN_{MB}}{dv_{z}}#right)");
651 slope = ( histo->GetBinContent(histo->FindBin(meanZ+sigmaZ)) -
652 histo->GetBinContent(histo->FindBin(meanZ-sigmaZ)) ) / ( 2. * sigmaZ );
654 if ( slope < 0. ) slope = norm/kFreePath;
656 // Try to fit twice: it fit fails the first time
657 // set some limits on parameters
658 for ( Int_t itry=0; itry<2; itry++ ) {
659 fitFunc->SetParameter(0, norm);
660 fitFunc->SetParameter(1, slope);
662 limitNorm = 2.*histo->Integral();
663 limitSlope = 2.*histo->Integral()/kFreePath;
664 //fitFunc->SetParLimits(0, 0., limitNorm); // REMEMBER TO CHECK
665 fitFunc->SetParLimits(1, 0., limitSlope); // REMEMBER TO CHECK
666 printf("Norm 0. < %f < %f slope 0. < %f < %f\n", norm, limitNorm, slope, limitSlope);
668 TFitResultPtr fitRes = histo->Fit(fitFunc, fitOpt.Data(), "", minZfit, maxZfit);
670 // if ( gMinuit->fCstatu.Contains("CONVERGED") &&
671 if ( ((Int_t)fitRes) == 0 &&
672 fitFunc->GetParameter(0) > 0. &&
673 fitFunc->GetParameter(1) > 0. )
675 else if ( furtherOpt.Contains("REFIT") ) printf("Re-fit with limits\n");
677 printf("Warning: fit problems !!!\n");
682 Double_t p0 = fitFunc->GetParameter(0);
683 Double_t p0err = fitFunc->GetParError(0);
684 Double_t p1 = fitFunc->GetParameter(1);
685 Double_t p1err = fitFunc->GetParError(1);
687 Double_t nullVz = ( p1 != 0. ) ? -p0/p1 : 0.;
688 Double_t nullVzErr = ( p0 != 0. && p1 != 0. ) ? TMath::Abs(nullVz) * TMath::Sqrt(p0err*p0err/(p0*p0) + p1err*p1err/(p1*p1) ) : 0.;
690 printf("Null value at %f +- %f\n", nullVz - kFreePath, nullVzErr);
692 recoHisto[kRecoHF]->SetBinContent(ibinpt, p0);
693 recoHisto[kRecoHF]->SetBinError(ibinpt, p0err);
694 recoHisto[kRecoBkg]->SetBinContent(ibinpt, ( kFreePath + meanZ ) * p1);
695 recoHisto[kRecoBkg]->SetBinError(ibinpt, ( kFreePath + meanZ ) * p1err);
696 if ( currDraw%4 == 0 ){
697 currName = Form("vtx_%s_PtBin%i",fThetaAbsKeys->At(itheta)->GetName(), ibinpt);
698 can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
701 can->cd( currDraw%4 + 1 );
704 fitFunc->DrawCopy("same");
707 SetSparseRange(gridSparse, kHvarMotherType, "", firstMother+1, lastMother+1, "USEBIN");
708 currName = Form("recoPt_%s",fThetaAbsKeys->At(itheta)->GetName());
709 can = new TCanvas(currName.Data(),currName.Data(),(igroup1+1)*xshift,igroup2*yshift,600,600);
710 TLegend* leg = new TLegend(0.6, 0.6, 0.8, 0.8);
712 for ( Int_t ireco=0; ireco<kNrecoHistos-1; ++ireco ) {
713 if ( recoHisto[ireco]->GetEntries() == 0. ) continue;
714 TH1* ratio = (TH1*)recoHisto[ireco]->Clone(Form("%s_ratio", recoHisto[ireco]->GetName()));
715 ratio->Divide(recoHisto[kRecoAll]);
716 ratio->SetLineColor(srcColors[ireco]);
717 ratio->SetMarkerColor(srcColors[ireco]);
718 ratio->SetMarkerStyle(20+ireco);
719 ratio->GetYaxis()->SetTitle("fraction of total");
720 ratio->Draw(drawOpt.Data());
721 leg->AddEntry(ratio,baseRecoName[ireco].Data(), "lp");
725 } // loop on theta abs