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