In Terminate: get track source based on output (Diego)
[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 #include "TFile.h"
49 #include "THashList.h"
50
51 // STEER includes
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"
61 #include "AliStack.h"
62
63 // ANALYSIS includes
64 #include "AliAnalysisManager.h"
65
66 // CORRFW includes
67 #include "AliCFContainer.h"
68 #include "AliCFGridSparse.h"
69 #include "AliCFEffGrid.h"
70
71 // PWG includes
72 #include "AliVAnalysisMuon.h"
73 #include "AliMergeableCollection.h"
74 #include "AliCounterCollection.h"
75 #include "AliMuonEventCuts.h"
76 #include "AliMuonTrackCuts.h"
77 #include "AliAnalysisMuonUtility.h"
78
79
80 /// \cond CLASSIMP
81 ClassImp(AliAnalysisTaskSingleMu) // Class implementation in ROOT context
82 /// \endcond
83
84
85 //________________________________________________________________________
86 AliAnalysisTaskSingleMu::AliAnalysisTaskSingleMu() :
87   AliVAnalysisMuon(),
88   fThetaAbsKeys(0x0),
89   fCutOnDimu(kFALSE)
90 {
91   /// Default ctor.
92 }
93
94 //________________________________________________________________________
95 AliAnalysisTaskSingleMu::AliAnalysisTaskSingleMu(const char *name, const AliMuonTrackCuts& cuts) :
96   AliVAnalysisMuon(name, cuts),
97   fThetaAbsKeys(0x0),
98   fCutOnDimu(kFALSE)
99 {
100   //
101   /// Constructor.
102   //
103   TString thetaAbsKeys = "ThetaAbs23 ThetaAbs310";
104   fThetaAbsKeys = thetaAbsKeys.Tokenize(" ");
105 }
106
107
108 //________________________________________________________________________
109 AliAnalysisTaskSingleMu::~AliAnalysisTaskSingleMu()
110 {
111   //
112   /// Destructor
113   //
114
115   delete fThetaAbsKeys;
116 }
117
118 //___________________________________________________________________________
119 void AliAnalysisTaskSingleMu::MyUserCreateOutputObjects()
120 {
121
122   TH1* histo = 0x0;
123   TString histoName = "", histoTitle = "";
124   
125   Int_t nVzBins = 40;
126   Double_t vzMin = -20., vzMax = 20.;
127   TString vzName("Vz"), vzTitle("Vz"), vzUnits("cm");  
128   
129   histoName = "hIpVtx";
130   histo = new TH1F(histoName.Data(), histoName.Data(), nVzBins, vzMin, vzMax);
131   histo->SetXTitle("v_{z} (cm)");
132   AddObjectToCollection(histo, kIPVz);
133
134   
135   Int_t nPtBins = 80;
136   Double_t ptMin = 0., ptMax = 80.;
137   TString ptName("Pt"), ptTitle("p_{T}"), ptUnits("GeV/c");
138   
139   Int_t nEtaBins = 25;
140   Double_t etaMin = -4.5, etaMax = -2.;
141   TString etaName("Eta"), etaTitle("#eta"), etaUnits("");
142   
143   Int_t nPhiBins = 36;
144   Double_t phiMin = 0.; Double_t phiMax = 2.*TMath::Pi();
145   TString phiName("Phi"), phiTitle("#phi"), phiUnits("rad");
146     
147   Int_t nChargeBins = 2;
148   Double_t chargeMin = -2, chargeMax = 2.;
149   TString chargeName("Charge"), chargeTitle("charge"), chargeUnits("e");
150   
151   Int_t nThetaAbsEndBins = 2;
152   Double_t thetaAbsEndMin = -0.5, thetaAbsEndMax = 1.5;
153   TString thetaAbsEndName("ThetaAbsEnd"), thetaAbsEndTitle("#theta_{abs}"), thetaAbsEndUnits("a.u.");    
154   
155   Int_t nMotherTypeBins = kNtrackSources;
156   Double_t motherTypeMin = -0.5, motherTypeMax = (Double_t)kNtrackSources - 0.5;
157   TString motherType("MotherType"), motherTypeTitle("motherType"), motherTypeUnits("");
158     
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};
164
165   AliCFContainer* cfContainer = new AliCFContainer("SingleMuContainer","Container for tracks",kNsteps,kNvars,nbins);
166   
167   for ( Int_t idim = 0; idim<kNvars; idim++){
168     histoTitle = Form("%s (%s)", axisTitle[idim].Data(), axisUnits[idim].Data());
169     histoTitle.ReplaceAll("()","");
170     
171     cfContainer->SetVarTitle(idim, histoTitle.Data());
172     cfContainer->SetBinLimits(idim, xmin[idim], xmax[idim]);
173   }
174   
175   TString stepTitle[kNsteps] = {"reconstructed", "generated"};
176
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);
181         
182     currAxis = gridSparse->GetAxis(kHvarMotherType);
183     for ( Int_t ibin=0; ibin<fSrcKeys->GetEntries(); ibin++ ) {
184       currAxis->SetBinLabel(ibin+1, fSrcKeys->At(ibin)->GetName());
185     }
186   }
187   
188   AddObjectToCollection(cfContainer, kTrackContainer);
189   
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);
195   
196   histoName = "hGeneratedZ";
197   AddObjectToCollection(histoDimu->Clone(histoName.Data()), kNobjectTypes+1);
198   
199   
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);
205   
206   fMuonTrackCuts->Print("mask");
207   
208   AliInfo(Form("Apply cut on dimuon (60<M_mumu<120 GeV/c^2) to reject Z contribution: %i", fCutOnDimu));
209 }
210
211 //________________________________________________________________________
212 void AliAnalysisTaskSingleMu::ProcessEvent(TString physSel, const TObjArray& selectTrigClasses, TString centrality)
213 {
214   //
215   /// Fill output objects
216   //
217
218   Double_t ipVz = AliAnalysisMuonUtility::GetVertexSPD(InputEvent())->GetZ();
219   Double_t ipVzMC = MCEvent() ? AliAnalysisMuonUtility::GetMCVertexZ(InputEvent(),MCEvent()) : 0.;
220   
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);
224   }
225
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
228   
229   Double_t containerInput[kNvars];
230   AliVParticle* track = 0x0;
231
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();
235     
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);
240       
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;
248       
249       selectedTracks.AddAt(track,itrack);
250       trackSources[itrack] = GetParticleType(track);
251       
252     } // loop on tracks
253     
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;
259       
260       if ( istep == kStepGeneratedMC || fCutOnDimu ) {
261         // Check dimuons
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;
266           
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;
273           }
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
280           }
281           else {
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. ) {
290                   isAccepted = kFALSE;
291                 }
292               } // loop on muons in the pair
293               
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);
297               
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;
308       
309       Double_t thetaAbsEndDeg = 0;
310       if ( istep == kStepReconstructed ) {
311         thetaAbsEndDeg = AliAnalysisMuonUtility::GetThetaAbsDeg(track);
312       }
313       else {
314         thetaAbsEndDeg = ( TMath::Pi()-track->Theta() ) * TMath::RadToDeg();
315       }
316       Int_t thetaAbsBin = ( thetaAbsEndDeg < 3. ) ? kThetaAbs23 : kThetaAbs310;
317
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];
325       
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
331     } // loop on tracks
332   } // loop on container steps
333 }
334
335
336 //________________________________________________________________________
337 void AliAnalysisTaskSingleMu::Terminate(Option_t *) {
338   //
339   /// Draw some histograms at the end.
340   //
341   
342   AliVAnalysisMuon::Terminate("");
343   
344   if ( ! fMergeableCollection ) return;
345   
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();
350   
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;
357   }
358   delete optArr;
359   
360   furtherOpt.ToUpper();
361   Bool_t plotChargeAsymmetry = furtherOpt.Contains("ASYM");
362   
363   AliCFContainer* inContainer = static_cast<AliCFContainer*> ( GetSum(physSel,trigClassName,centralityRange,"SingleMuContainer") );
364   if ( ! inContainer ) return;
365   
366   AliCFContainer* cfContainer = inContainer;
367   
368   if ( ! furtherOpt.Contains("GENINTRIGCLASS") && trigClassName != "ANY" ) {
369     // The output container contains both the reconstructed and (in MC)
370     // the generated muons in a specific trigger class.
371     // Since the trigger pt level of the track is matched to the trigger class,
372     // analyzing the MUHigh trigger (for example) is equivalent of analyzing
373     // Hpt tracks.
374     // However, in this way, the generated muons distribution depend
375     // on a condition on the reconstructed track.
376     // If we calculate the Acc.xEff. as the ratio of reconstructed and
377     // generated muons IN A TRIGGER CLASS, we are biasing the final result.
378     // The correct value of Acc.xEff. is hence obtained as the distribution
379     // of reconstructed tracks in a muon trigger class, divided by the
380     // total number of generated class (which is in the "class" ANY).
381     // The following code sets the generated muons as the one in the class ANY.
382     // The feature is the default. If you want the generated muons in the same
383     // trigger class as the generated tracks, use option "GENINTRIGCLASS"
384     AliCFContainer* fullMCcontainer = static_cast<AliCFContainer*> ( GetSum(physSel,"ANY",centralityRange,"SingleMuContainer") );
385     if ( fullMCcontainer ) {
386       cfContainer = static_cast<AliCFContainer*>(cfContainer->Clone("SingleMuContainerCombo"));
387       cfContainer->SetGrid(kStepGeneratedMC,fullMCcontainer->GetGrid(kStepGeneratedMC));
388     }
389   }
390   
391   AliCFEffGrid* effSparse = new AliCFEffGrid(Form("eff%s", cfContainer->GetName()),Form("Efficiency %s", cfContainer->GetTitle()),*cfContainer);
392   effSparse->CalculateEfficiency(kStepReconstructed, kStepGeneratedMC);
393   
394   AliCFGridSparse* gridSparseArray[3] = {effSparse->GetNum(), effSparse->GetDen(), effSparse};
395   TString gridSparseName[3] = {cfContainer->GetStepTitle(kStepReconstructed), cfContainer->GetStepTitle(kStepGeneratedMC), "Efficiency"};
396   
397   Int_t srcColors[kNtrackSources] = {kBlack, kRed, kSpring, kTeal, kBlue, kViolet, kMagenta, kOrange, kGray};
398   //  TString allSrcNames = "";
399   //  for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
400   //    if ( ! allSrcNames.IsNull() ) allSrcNames.Append(" ");
401   //    allSrcNames += fSrcKeys->At(isrc)->GetName();
402   //  }
403   
404   TCanvas* can = 0x0;
405   Int_t xshift = 100;
406   Int_t yshift = 100;
407   Int_t igroup1 = -1;
408   Int_t igroup2 = 0;
409   
410   Bool_t isMC = furtherOpt.Contains("MC");
411   
412   TAxis* srcAxis = gridSparseArray[0]->GetAxis(kHvarMotherType);
413   Int_t unIdBin = srcAxis->GetNbins();
414   for ( Int_t ibin=1; ibin<=srcAxis->GetNbins(); ibin++ ) {
415     currName = srcAxis->GetBinLabel(ibin);
416     if ( currName.Contains("Unidentified") ) unIdBin = ibin;
417   }
418   
419   Int_t firstSrcBin = ( isMC ) ? 1 : unIdBin;
420   Int_t lastSrcBin  = ( isMC ) ? srcAxis->GetNbins() - 1 : unIdBin;
421   if ( ! isMC ) srcColors[unIdBin-1] = 1;
422   
423   ////////////////
424   // Kinematics //
425   ////////////////
426   TString chargeNames[3] = {fChargeKeys->At(0)->GetName(), fChargeKeys->At(1)->GetName(), "Total"};
427   THashList histoList[3];
428   for ( Int_t icharge=0; icharge<3; icharge++ ) {
429     histoList[icharge].SetName(chargeNames[icharge].Data());
430   }
431   for ( Int_t isrc = firstSrcBin; isrc <= lastSrcBin; ++isrc ) {
432     for ( Int_t icharge=0; icharge<3; ++icharge ) {
433       Int_t icharge1 = ( icharge < 2 ) ? icharge : 0;
434       Int_t icharge2 = ( icharge < 2 ) ? icharge : 1;
435       for ( Int_t igrid=0; igrid<3; ++igrid ) {
436         if ( gridSparseArray[igrid]->GetEntries() == 0. ) break;
437         if ( gridSparseArray[igrid]->IsA() != AliCFEffGrid::Class() ) {
438           SetSparseRange(gridSparseArray[igrid], kHvarEta, "", -3.999, -2.501);
439           SetSparseRange(gridSparseArray[igrid], kHvarMotherType, "", isrc, isrc, "USEBIN");
440           SetSparseRange(gridSparseArray[igrid], kHvarCharge, "", icharge1+1, icharge2+1, "USEBIN");
441         }
442         TH1* histo = gridSparseArray[igrid]->Project(kHvarPt, kHvarEta);
443         histo->SetName(Form("hPtEta_%s_%s_%s", gridSparseName[igrid].Data(), srcAxis->GetBinLabel(isrc), chargeNames[icharge].Data()));
444         if ( histo->Integral() > 0 ) histoList[icharge].Add(histo);
445         for ( Int_t iproj=0; iproj<4; ++iproj ) {
446           histo = gridSparseArray[igrid]->Project(iproj);
447           histo->SetName(Form("proj%i_%s_%s_%s", iproj, gridSparseName[igrid].Data(), srcAxis->GetBinLabel(isrc), chargeNames[icharge].Data()));
448           if ( histo->Integral() > 0 ) histoList[icharge].Add(histo);
449         } // loop on projections
450       } // loop on grid sparse
451     } // loop on charge
452   } // loop on track sources
453   
454   // Get charge asymmetry or mu+/mu-
455   THashList histoListRatio;
456   TString basePlotName = plotChargeAsymmetry ? "ChargeAsymmetry" : "ChargeRatio";
457   histoListRatio.SetName(basePlotName.Data());
458   Int_t baseCharge = 1;
459   Int_t auxCharge = 1-baseCharge;
460   for ( Int_t ihisto=0; ihisto<histoList[baseCharge].GetEntries(); ihisto++ ) {
461     TObject* obj = histoList[baseCharge].At(ihisto);
462     TString histoName = obj->GetName();
463     if ( histoName.Contains(gridSparseName[2].Data()) ) continue;
464     TString searchName = histoName;
465     searchName.ReplaceAll(fChargeKeys->At(baseCharge)->GetName(), fChargeKeys->At(auxCharge)->GetName());
466     TH1* auxHisto = static_cast<TH1*> (histoList[auxCharge].FindObject(searchName.Data()));
467     if ( ! auxHisto ) continue;
468     histoName.ReplaceAll(fChargeKeys->At(baseCharge)->GetName(),basePlotName.Data());
469     TH1* histo = static_cast<TH1*> (obj->Clone(histoName.Data()));
470     if ( plotChargeAsymmetry ) {
471       histo->Add(auxHisto, -1.);
472       // h2 + h1 = 2xh2 + (h1-h2)
473       auxHisto->Add(auxHisto, histo, 2.);
474     }
475     histo->Divide(auxHisto);
476     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());
477     axisTitle.ReplaceAll("MuPlus","#mu^{+}");
478     axisTitle.ReplaceAll("MuMinus","#mu^{-}");
479     histo->GetYaxis()->SetTitle(axisTitle.Data());
480     histo->SetStats(kFALSE);
481     histoListRatio.Add(histo);
482   }
483   
484   // Plot kinematics
485   TString histoName = "", drawOpt = "";
486   for ( Int_t itype=0; itype<3; itype++ ) {
487     THashList* currList = 0x0;
488     Int_t nCharges = 1;
489     if ( itype == 1 ) currList = &histoListRatio;
490     else if ( itype == 2 ) currList = &histoList[2];
491     else nCharges = 2;
492     for ( Int_t igrid=0; igrid<3; ++igrid ) {
493       igroup1 = igrid;
494       TCanvas* canKine = 0x0;
495       TLegend* legKine = 0x0;
496       for ( Int_t iproj=0; iproj<4; ++iproj ) {
497         for ( Int_t isrc = firstSrcBin; isrc <= lastSrcBin; ++isrc ) {
498           for ( Int_t icharge=0; icharge<nCharges; ++icharge ) {
499             if ( itype == 0 ) currList = &histoList[icharge];
500             histoName = Form("proj%i_%s_%s_%s", iproj, gridSparseName[igrid].Data(), srcAxis->GetBinLabel(isrc), currList->GetName());
501             TH1* histo = static_cast<TH1*>(currList->FindObject(histoName.Data()));
502             if ( ! histo ) continue;
503             if ( ! canKine ) {
504               igroup2 = itype;
505               igroup1 = igrid;
506               currName = Form("%s_%s_%s", GetName(), currList->GetName(), gridSparseName[igrid].Data());
507               canKine = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
508               canKine->Divide(2,2);
509               legKine = new TLegend(0.6, 0.6, 0.8, 0.8);
510             }
511             canKine->cd(iproj+1);
512             if ( itype != 1 ) {
513               if ( ( iproj == kHvarPt || iproj == kHvarVz ) && gridSparseArray[igrid]->IsA() != AliCFEffGrid::Class() ) gPad->SetLogy();
514             }
515             Bool_t isFirst = ( gPad->GetListOfPrimitives()->GetEntries() == 0 );
516             drawOpt = isFirst ? "e" : "esames";
517             histo->SetLineColor(srcColors[isrc-1]);
518             histo->SetMarkerColor(srcColors[isrc-1]);
519             histo->SetMarkerStyle(20+4*icharge);
520             histo->Draw(drawOpt.Data());
521             TPaveStats* paveStats = (TPaveStats*)histo->FindObject("stats");
522             if ( paveStats ) paveStats->SetTextColor(srcColors[isrc-1]);
523             if ( iproj == 0 ) {
524               TString legEntry = ( itype == 0 ) ? fChargeKeys->At(icharge)->GetName() : "";
525               if ( isMC ) legEntry += Form(" %s", srcAxis->GetBinLabel(isrc));
526               if ( ! legEntry.IsNull() ) legKine->AddEntry(histo,legEntry.Data(), "lp");
527             }
528           } // loop on mu charge
529         } // loop on track sources
530       } // loop on projections
531       if ( legKine && legKine->GetNRows() > 0 ) {
532         canKine->cd(1);
533         legKine->Draw("same");
534       }
535     } // loop on grid sparse
536   } // loop on types
537   
538   
539   for ( Int_t igrid=0; igrid<3; igrid++ ) {
540     if ( gridSparseArray[igrid]->IsA() == AliCFEffGrid::Class() ) continue;
541     SetSparseRange(gridSparseArray[igrid], kHvarCharge, "", 1, gridSparseArray[igrid]->GetAxis(kHvarCharge)->GetNbins(), "USEBIN"); // Reset range
542   } // loop on container steps
543   
544   //////////////////////
545   // Event statistics //
546   //////////////////////
547   printf("\nTotal analyzed events:\n");
548   TString evtSel = Form("trigger:%s", trigClassName.Data());
549   fEventCounters->PrintSum(evtSel.Data());
550   printf("Physics selected analyzed events:\n");
551   evtSel = Form("trigger:%s/selected:yes", trigClassName.Data());
552   fEventCounters->PrintSum(evtSel.Data());
553   
554   TString countPhysSel = "any";
555   if ( physSel.Contains(fPhysSelKeys->At(kPhysSelPass)->GetName()) ) countPhysSel = "yes";
556   else if ( physSel.Contains(fPhysSelKeys->At(kPhysSelReject)->GetName()) ) countPhysSel="no";
557   countPhysSel.Prepend("selected:");
558   printf("Analyzed events vs. centrality:\n");
559   evtSel = Form("trigger:%s/%s", trigClassName.Data(), countPhysSel.Data());
560   fEventCounters->Print("centrality",evtSel.Data(),kTRUE);
561   
562   
563   TString outFilename = Form("/tmp/out%s.root", GetName());
564   TFile* outFile = new TFile(outFilename.Data(),"RECREATE");
565   for ( Int_t icharge=0; icharge<3; icharge++ ) {
566     histoList[icharge].Write();
567   }
568   histoListRatio.Write();
569   outFile->Close();
570   printf("\nCreating file %s\n", outFilename.Data(
571          ));
572   
573   ///////////////////
574   // Vertex method //
575   ///////////////////
576   if ( ! furtherOpt.Contains("VERTEX") ) return;
577   igroup1++;
578   TH1* eventVertex = (TH1*)GetSum(physSel, minBiasTrig, centralityRange, "hIpVtx");
579   if ( ! eventVertex ) return;
580   Double_t minZ = -9.99, maxZ = 9.99;
581   Double_t meanZ = 0., sigmaZ = 4.;
582   Double_t nSigma = 2.;
583   TString fitOpt = "R0S";
584   Bool_t fixFitRange = kFALSE;
585   TString fitFormula = Form("[0]+[1]*(x+[2])");
586   
587   // Get vertex shape
588   if ( eventVertex->GetSumw2N() == 0 ) eventVertex->Sumw2();
589   Double_t eventVtxIntegral = eventVertex->Integral(0,eventVertex->GetNbinsX()+1); // Include under/overflow
590   printf("Event vertex integral %.0f\n\n", eventVtxIntegral);
591   if ( eventVtxIntegral <= 0. ) return;
592   eventVertex->Scale(1./eventVtxIntegral);
593   printf("\nFit MB vertex\n");
594   eventVertex->Fit("gaus",fitOpt.Data(),"",minZ,maxZ);
595   TF1* vtxFit = (TF1*)eventVertex->GetListOfFunctions()->FindObject("gaus");
596   currName = "vtxIntegrated";
597   can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
598   can->SetLogy();
599   eventVertex->Draw();
600   vtxFit->Draw("same");
601   
602   
603   enum {kRecoHF, kRecoBkg, kInputHF, kInputDecay, kRecoAll, kNrecoHistos};
604   TString baseRecoName[kNrecoHistos] = {"RecoHF", "RecoBkg", "InputHF", "InputDecay", "RecoAll"};
605   TArrayI sumMothers[kNrecoHistos];
606   sumMothers[kRecoHF].Set(0);
607   sumMothers[kRecoBkg].Set(0);
608   sumMothers[kInputHF].Set(3);
609   sumMothers[kInputHF][0] = kCharmMu;
610   sumMothers[kInputHF][1] = kBeautyMu;
611   sumMothers[kInputHF][2] = kQuarkoniumMu;
612   sumMothers[kInputDecay].Set(1);
613   sumMothers[kInputDecay][0] = kDecayMu;
614   sumMothers[kRecoAll].Set(srcAxis->GetNbins());
615   for ( Int_t isrc=1; isrc<srcAxis->GetNbins(); ++isrc ) {
616     sumMothers[kRecoAll][isrc-1] = isrc;
617   }
618   
619   meanZ = vtxFit->GetParameter(1);
620   sigmaZ = vtxFit->GetParameter(2);
621   
622   Double_t minZfit = ( fixFitRange ) ? minZ : meanZ - nSigma*sigmaZ;
623   Double_t maxZfit = ( fixFitRange ) ? maxZ : meanZ + nSigma*sigmaZ;
624   
625   TF1* fitFunc = new TF1("fitFunc", fitFormula.Data(), minZ, maxZ);
626   fitFunc->SetLineColor(2);
627   fitFunc->SetParNames("Line norm", "Line slope", "Free path");
628   const Double_t kFreePath = 153.; // 150.; // 130.; // cm
629   //fitFunc->SetParameters(0.,1.);
630   fitFunc->FixParameter(2, kFreePath);
631   
632   AliCFGridSparse* gridSparse = cfContainer->GetGrid(kStepReconstructed);
633   TAxis* ptAxis = gridSparse->GetAxis(kHvarPt);
634   
635   Double_t slope = 0.;
636   Double_t limitNorm = 0., limitSlope = 0.;
637   Int_t firstPtBin = 0, lastPtBin = 0;
638   
639   gStyle->SetOptFit(1111);
640   
641   for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) {
642     igroup2++;
643     SetSparseRange(gridSparse, kHvarThetaAbs, "", itheta+1, itheta+1, "USEBIN");
644     SetSparseRange(gridSparse, kHvarPt, "", 1, ptAxis->GetNbins(), "USEBIN");
645     TH1* recoHisto[kNrecoHistos];
646     for ( Int_t ireco=0; ireco<kNrecoHistos; ++ireco ) {
647       recoHisto[ireco] = gridSparse->Project(kHvarPt);
648       histoName = Form("%sMuon_%s", baseRecoName[ireco].Data(), fThetaAbsKeys->At(itheta)->GetName());
649       recoHisto[ireco]->SetName(histoName.Data());
650       recoHisto[ireco]->SetTitle(histoName.Data());
651       recoHisto[ireco]->Reset();
652       recoHisto[ireco]->Sumw2();
653       for ( Int_t isrc=0; isrc<sumMothers[ireco].GetSize(); ++isrc ) {
654         SetSparseRange(gridSparse, kHvarMotherType, "", sumMothers[ireco][isrc], sumMothers[ireco][isrc], "USEBIN");
655         TH1* auxHisto = gridSparse->Project(kHvarPt);
656         recoHisto[ireco]->Add(auxHisto);
657         delete auxHisto;
658       }
659     }
660     SetSparseRange(gridSparse, kHvarMotherType, "", firstPtBin, lastSrcBin, "USEBIN");
661     Int_t currDraw = 0;
662     
663     for ( Int_t ibinpt=0; ibinpt<=ptAxis->GetNbins(); ++ibinpt ) {
664       firstPtBin = ibinpt;
665       lastPtBin = ( ibinpt == 0 ) ? ptAxis->GetNbins() : ibinpt;
666       SetSparseRange(gridSparse, kHvarPt, "", firstPtBin, lastPtBin, "USEBIN");
667       TH1* histo = gridSparse->Project(kHvarVz);
668       histo->SetName(Form("hVtx_%s_%s_ptBin%i", cfContainer->GetStepTitle(kStepReconstructed), fThetaAbsKeys->At(itheta)->GetName(), ibinpt));
669       if ( histo->Integral() < 100. ) break;
670       printf("\nFit %.2f < pt < %.2f (entries %g)\n", ptAxis->GetBinLowEdge(firstPtBin), ptAxis->GetBinUpEdge(lastPtBin), histo->GetEntries());
671       histo->Divide(eventVertex);
672       Double_t norm = histo->GetBinContent(histo->FindBin(0.));
673       histo->GetYaxis()->SetTitle("#frac{dN_{#mu}}{dv_{z}} / #left(#frac{1}{N_{MB}}#frac{dN_{MB}}{dv_{z}}#right)");
674       slope = ( histo->GetBinContent(histo->FindBin(meanZ+sigmaZ)) -
675                histo->GetBinContent(histo->FindBin(meanZ-sigmaZ)) ) / ( 2. * sigmaZ );
676       
677       if ( slope < 0. ) slope = norm/kFreePath;
678       
679       // Try to fit twice: it fit fails the first time
680       // set some limits on parameters
681       for ( Int_t itry=0; itry<2; itry++ ) {
682         fitFunc->SetParameter(0, norm);
683         fitFunc->SetParameter(1, slope);
684         if ( itry > 0 ) {
685           limitNorm = 2.*histo->Integral();
686           limitSlope = 2.*histo->Integral()/kFreePath;
687           //fitFunc->SetParLimits(0, 0., limitNorm); // REMEMBER TO CHECK
688           fitFunc->SetParLimits(1, 0., limitSlope); // REMEMBER TO CHECK
689           printf("Norm 0. < %f < %f  slope  0. < %f < %f\n", norm, limitNorm, slope, limitSlope);
690         }
691         TFitResultPtr fitRes = histo->Fit(fitFunc, fitOpt.Data(), "", minZfit, maxZfit);
692         
693         //      if ( gMinuit->fCstatu.Contains("CONVERGED") &&
694         if ( ((Int_t)fitRes) == 0 &&
695             fitFunc->GetParameter(0) > 0. &&
696             fitFunc->GetParameter(1) > 0. )
697           break;
698         else if ( furtherOpt.Contains("REFIT") ) printf("Re-fit with limits\n");
699         else {
700           printf("Warning: fit problems !!!\n");
701           break;
702         }
703       }
704       
705       Double_t p0 = fitFunc->GetParameter(0);
706       Double_t p0err = fitFunc->GetParError(0);
707       Double_t p1 = fitFunc->GetParameter(1);
708       Double_t p1err = fitFunc->GetParError(1);
709       
710       Double_t nullVz = ( p1 != 0. ) ? -p0/p1 : 0.;
711       Double_t nullVzErr = ( p0 != 0. && p1 != 0. ) ? TMath::Abs(nullVz) * TMath::Sqrt(p0err*p0err/(p0*p0) + p1err*p1err/(p1*p1) ) : 0.;
712       
713       printf("Null value at %f +- %f\n", nullVz - kFreePath, nullVzErr);
714       
715       recoHisto[kRecoHF]->SetBinContent(ibinpt, p0);
716       recoHisto[kRecoHF]->SetBinError(ibinpt, p0err);
717       recoHisto[kRecoBkg]->SetBinContent(ibinpt, ( kFreePath + meanZ ) * p1);
718       recoHisto[kRecoBkg]->SetBinError(ibinpt, ( kFreePath + meanZ ) * p1err);
719       if ( currDraw%4 == 0 ){
720         currName = Form("vtx_%s_PtBin%i",fThetaAbsKeys->At(itheta)->GetName(), ibinpt);
721         can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
722         can->Divide(2,2);
723       }
724       can->cd( currDraw%4 + 1 );
725       can->SetLogy();
726       histo->Draw();
727       fitFunc->DrawCopy("same");
728       currDraw++;
729     } // loop on pt bins
730     SetSparseRange(gridSparse, kHvarMotherType, "", firstSrcBin, lastSrcBin, "USEBIN");
731     currName = Form("recoPt_%s",fThetaAbsKeys->At(itheta)->GetName());
732     can = new TCanvas(currName.Data(),currName.Data(),(igroup1+1)*xshift,igroup2*yshift,600,600);
733     TLegend* leg = new TLegend(0.6, 0.6, 0.8, 0.8);
734     drawOpt = "e";
735     for ( Int_t ireco=0; ireco<kNrecoHistos-1; ++ireco ) {
736       if ( recoHisto[ireco]->GetEntries() == 0. ) continue;
737       TH1* ratio = (TH1*)recoHisto[ireco]->Clone(Form("%s_ratio", recoHisto[ireco]->GetName()));
738       ratio->Divide(recoHisto[kRecoAll]);
739       ratio->SetLineColor(srcColors[ireco]);
740       ratio->SetMarkerColor(srcColors[ireco]);
741       ratio->SetMarkerStyle(20+ireco);
742       ratio->GetYaxis()->SetTitle("fraction of total");
743       ratio->Draw(drawOpt.Data());
744       leg->AddEntry(ratio,baseRecoName[ireco].Data(), "lp");
745       drawOpt = "esame";
746     }
747     leg->Draw("same");
748   } // loop on theta abs
749 }