Update of Single Muon analysis to be included in the official ESD->AOD analysis train...
[u/mrichter/AliRoot.git] / PWG3 / 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 //-----------------------------------------------------------------------------
17 /// \class AliAnalysisTaskSingleMu
18 /// Analysis task for single muons in the spectrometer.
19 /// The output is a list of histograms.
20 /// The macro class can run on AOD or in the train with the ESD filter.
21 /// If Monte Carlo information is present, some basics checks are performed.
22 ///
23 /// \author Diego Stocco
24 //-----------------------------------------------------------------------------
25
26 //----------------------------------------------------------------------------
27 //    Implementation of the single muon analysis class
28 // An example of usage can be found in the macro RunSingleMuAnalysisFromAOD.C.
29 //----------------------------------------------------------------------------
30
31 #define AliAnalysisTaskSingleMu_cxx
32
33 // ROOT includes
34 #include "TROOT.h"
35 #include "TH1.h"
36 #include "TH2.h"
37 #include "TH3.h"
38 #include "TAxis.h"
39 #include "TList.h"
40 #include "TCanvas.h"
41 #include "TMath.h"
42
43 // STEER includes
44 #include "AliLog.h"
45
46 #include "AliAODEvent.h"
47 #include "AliAODTrack.h"
48 #include "AliAODVertex.h"
49
50 #include "AliMCEvent.h"
51 #include "AliMCParticle.h"
52 #include "AliStack.h"
53
54
55 #include "AliAnalysisTaskSE.h"
56 //#include "AliAnalysisDataSlot.h"
57 //#include "AliAnalysisManager.h"
58 #include "AliAnalysisTaskSingleMu.h"
59
60 /// \cond CLASSIMP
61 ClassImp(AliAnalysisTaskSingleMu) // Class implementation in ROOT context
62 /// \endcond
63
64
65 //________________________________________________________________________
66 AliAnalysisTaskSingleMu::AliAnalysisTaskSingleMu(const char *name) :
67   AliAnalysisTaskSE(name),
68   fUseMC(0),
69   fHistoList(0),
70   fHistoListMC(0)
71 {
72   //
73   /// Constructor.
74   //
75   DefineOutput(1, TList::Class());
76   DefineOutput(2, TList::Class());
77 }
78
79
80 //________________________________________________________________________
81 AliAnalysisTaskSingleMu::~AliAnalysisTaskSingleMu()
82 {
83   //
84   /// Destructor
85   //
86   if ( fHistoList ) delete fHistoList;
87   if ( fHistoListMC ) delete fHistoListMC;
88 }
89
90
91 //___________________________________________________________________________
92 void AliAnalysisTaskSingleMu::SetFlagsFromHandler()
93 {
94   //
95   /// Use the event handler information to correctly fill the analysis flags:
96   /// - check if Monte Carlo information is present
97   //  
98   
99   if ( MCEvent() ) {
100     fUseMC = kTRUE;
101     AliInfo("Monte Carlo information is present");
102   }
103   else {
104     AliInfo("No Monte Carlo information in run");
105   }
106
107 }
108
109
110 //___________________________________________________________________________
111 void AliAnalysisTaskSingleMu::UserCreateOutputObjects() 
112 {
113   //
114   /// Create output histograms
115   //
116   AliInfo(Form("   CreateOutputObjects of task %s\n", GetName()));
117
118   // initialize histogram lists
119   if(!fHistoList) fHistoList = new TList();
120   if(!fHistoListMC) fHistoListMC = new TList();
121   
122   Int_t nPtBins = 60;
123   Float_t ptMin = 0., ptMax = 30.;
124   TString ptName("Pt"), ptTitle("p_{t}"), ptUnits("GeV/c");
125   
126   Int_t nDcaBins = 100;
127   Float_t dcaMin = 0., dcaMax = 200.;
128   TString dcaName("DCA"), dcaTitle("DCA"), dcaUnits("cm");
129
130   Int_t nVzBins = 60;
131   Float_t vzMin = -30, vzMax = 30;
132   TString vzName("Vz"), vzTitle("Vz"), vzUnits("cm");
133   
134   Int_t nEtaBins = 25;
135   Float_t etaMin = -4.5, etaMax = -2.;
136   TString etaName("Eta"), etaTitle("#eta"), etaUnits("a.u.");
137   
138   Int_t nRapidityBins = 25;
139   Float_t rapidityMin = -4.5, rapidityMax = -2.;
140   TString rapidityName("Rapidity"), rapidityTitle("rapidity"), rapidityUnits("a.u.");
141
142   TString trigName[kNtrigCuts];
143   trigName[kNoMatchTrig] = "NoMatch";
144   trigName[kAllPtTrig]   = "AllPt";
145   trigName[kLowPtTrig]   = "LowPt";
146   trigName[kHighPtTrig]  = "HighPt";
147
148   TString srcName[kNtrackSources];
149   srcName[kCharmMu]     = "Charm";
150   srcName[kBeautyMu]    = "Beauty";
151   srcName[kPrimaryMu]   = "Decay";
152   srcName[kSecondaryMu] = "Secondary";
153   srcName[kRecoHadron]  = "Hadrons";
154   srcName[kUnknownPart] = "Unknown";
155
156   TH1F* histo1D = 0x0;
157   TH2F* histo2D = 0x0;
158   TString histoName, histoTitle;
159   Int_t histoIndex = 0;
160
161   // 1D histos
162   for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
163     histoName = Form("%sDistrib%sTrig", ptName.Data(), trigName[itrig].Data());
164     histoTitle = Form("%s distribution. Trigger: %s", ptTitle.Data(), trigName[itrig].Data());
165     histo1D = new TH1F(histoName.Data(), histoTitle.Data(), nPtBins, ptMin, ptMax);
166     histo1D->GetXaxis()->SetTitle(Form("%s (%s)", ptTitle.Data(), ptUnits.Data()));
167     histoIndex = GetHistoIndex(kHistoPt, itrig);
168     fHistoList->AddAt(histo1D, histoIndex);
169   }
170   
171   for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
172     histoName = Form("%sDistrib%sTrig", dcaName.Data(), trigName[itrig].Data());
173     histoTitle = Form("%s distribution. Trigger: %s", dcaTitle.Data(), trigName[itrig].Data());
174     histo1D = new TH1F(histoName.Data(), histoTitle.Data(), nDcaBins, dcaMin, dcaMax);
175     histo1D->GetXaxis()->SetTitle(Form("%s (%s)", dcaTitle.Data(), dcaUnits.Data()));
176     histoIndex = GetHistoIndex(kHistoDCA, itrig);
177     fHistoList->AddAt(histo1D, histoIndex);
178   }
179   
180   for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
181     histoName = Form("%sDistrib%sTrig", vzName.Data(), trigName[itrig].Data());
182     histoTitle = Form("%s distribution. Trigger: %s", vzTitle.Data(), trigName[itrig].Data());
183     histo1D = new TH1F(histoName.Data(), histoTitle.Data(), nVzBins, vzMin, vzMax);
184     histo1D->GetXaxis()->SetTitle(Form("%s (%s)", vzTitle.Data(), vzUnits.Data()));
185     histoIndex = GetHistoIndex(kHistoVz, itrig);
186     fHistoList->AddAt(histo1D, histoIndex);
187   }
188   
189   for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
190     histoName = Form("%sDistrib%sTrig", etaName.Data(), trigName[itrig].Data());
191     histoTitle = Form("%s distribution. Trigger: %s", etaTitle.Data(), trigName[itrig].Data());
192     histo1D = new TH1F(histoName.Data(), histoTitle.Data(), nEtaBins, etaMin, etaMax);
193     histo1D->GetXaxis()->SetTitle(Form("%s (%s)", etaTitle.Data(), etaUnits.Data()));
194     histoIndex = GetHistoIndex(kHistoEta, itrig);
195     fHistoList->AddAt(histo1D, histoIndex);
196   }
197   
198   for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
199     histoName = Form("%sDistrib%sTrig", rapidityName.Data(), trigName[itrig].Data());
200     histoTitle = Form("%s distribution. Trigger: %s", rapidityTitle.Data(), trigName[itrig].Data());
201     histo1D = new TH1F(histoName.Data(), histoTitle.Data(), nRapidityBins, rapidityMin, rapidityMax);
202     histo1D->GetXaxis()->SetTitle(Form("%s (%s)", rapidityTitle.Data(), rapidityUnits.Data()));
203     histoIndex = GetHistoIndex(kHistoRapidity, itrig);
204     fHistoList->AddAt(histo1D, histoIndex);
205   }
206
207   // 2D histos
208   for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
209     histoName = Form("%s%sDistrib%sTrig", dcaName.Data(), ptName.Data(), trigName[itrig].Data());
210     histoTitle = Form("%s vs. %s distribution. Trigger: %s", dcaTitle.Data(), ptTitle.Data(), trigName[itrig].Data());
211     histo2D = new TH2F(histoName.Data(), histoTitle.Data(),
212                        nPtBins, ptMin, ptMax,
213                        nDcaBins, dcaMin, dcaMax);
214     histo2D->GetXaxis()->SetTitle(Form("%s (%s)", ptTitle.Data(), ptUnits.Data()));
215     histo2D->GetYaxis()->SetTitle(Form("%s (%s)", dcaTitle.Data(), dcaUnits.Data()));
216     histoIndex = GetHistoIndex(kHistoPtDCA, itrig);
217     fHistoList->AddAt(histo2D, histoIndex);
218   }
219   
220   for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
221     histoName = Form("%s%sDistrib%sTrig", vzName.Data(), ptName.Data(), trigName[itrig].Data());
222     histoTitle = Form("%s vs. %s distribution. Trigger: %s", vzTitle.Data(), ptTitle.Data(), trigName[itrig].Data());
223     histo2D = new TH2F(histoName.Data(), histoTitle.Data(),
224                        nPtBins, ptMin, ptMax,
225                        nVzBins, vzMin, vzMax);
226     histo2D->GetXaxis()->SetTitle(Form("%s (%s)", ptTitle.Data(), ptUnits.Data()));
227     histo2D->GetYaxis()->SetTitle(Form("%s (%s)", vzTitle.Data(), vzUnits.Data()));
228     histoIndex = GetHistoIndex(kHistoPtVz, itrig);
229     fHistoList->AddAt(histo2D, histoIndex);
230   }
231   
232   for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
233     histoName = Form("%s%sDistrib%sTrig", rapidityName.Data(), ptName.Data(), trigName[itrig].Data());
234     histoTitle = Form("%s vs. %s distribution. Trigger: %s", rapidityTitle.Data(), ptTitle.Data(), trigName[itrig].Data());
235     histo2D = new TH2F(histoName.Data(), histoTitle.Data(),
236                      nPtBins, ptMin, ptMax,
237                      nRapidityBins, rapidityMin, rapidityMax);
238     histo2D->GetXaxis()->SetTitle(Form("%s (%s)", ptTitle.Data(), ptUnits.Data()));
239     histo2D->GetYaxis()->SetTitle(Form("%s (%s)", rapidityTitle.Data(), rapidityUnits.Data()));
240     histoIndex = GetHistoIndex(kHistoPtRapidity, itrig);
241     fHistoList->AddAt(histo2D, histoIndex);
242   }
243
244   // MC histos
245   for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
246     for (Int_t isrc = 0; isrc < kNtrackSources; isrc++) {
247       histoName = Form("%sResolution%sTrig%s", ptName.Data(), trigName[itrig].Data(), srcName[isrc].Data());
248       histoTitle = Form("%s resolution. Trigger: %s (%s)", ptTitle.Data(), trigName[itrig].Data(), srcName[isrc].Data());
249       histo1D = new TH1F(histoName.Data(), histoTitle.Data(), 100, -5., 5.);
250       histo1D->GetXaxis()->SetTitle(Form("%s^{reco} - %s^{MC} (%s)", ptTitle.Data(), ptTitle.Data(), ptUnits.Data()));
251       histoIndex = GetHistoIndex(kHistoPtResolution, itrig, isrc);
252       fHistoListMC->AddAt(histo1D, histoIndex);
253     }
254   }
255
256   for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
257     for (Int_t isrc = 0; isrc < kNtrackSources; isrc++) {
258       histoName = Form("%s%sDistrib%sTrig%s", dcaName.Data(), ptName.Data(),
259                        trigName[itrig].Data(), srcName[isrc].Data());
260       histoTitle = Form("%s vs. %s distribution. Trigger: %s (%s)", dcaTitle.Data(), ptTitle.Data(),
261                         trigName[itrig].Data(), srcName[isrc].Data());
262       histo2D = new TH2F(histoName.Data(), histoTitle.Data(),
263                          nPtBins, ptMin, ptMax,
264                          nDcaBins, dcaMin, dcaMax);
265       histo2D->GetXaxis()->SetTitle(Form("%s (%s)", ptTitle.Data(), ptUnits.Data()));
266       histo2D->GetYaxis()->SetTitle(Form("%s (%s)", dcaTitle.Data(), dcaUnits.Data()));
267       histoIndex = GetHistoIndex(kHistoPtDCAType, itrig, isrc);
268       fHistoListMC->AddAt(histo2D, histoIndex);
269     }
270   }
271   
272   for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
273     for (Int_t isrc = 0; isrc < kNtrackSources; isrc++) {
274       histoName = Form("%s%sDistrib%sTrig%s", vzName.Data(), ptName.Data(), 
275                        trigName[itrig].Data(), srcName[isrc].Data());
276       histoTitle = Form("%s vs. %s distribution. Trigger: %s (%s)", vzTitle.Data(), ptTitle.Data(),
277                         trigName[itrig].Data(), srcName[isrc].Data());
278       histo2D = new TH2F(histoName.Data(), histoTitle.Data(),
279                          nPtBins, ptMin, ptMax,
280                          nVzBins, vzMin, vzMax);
281       histo2D->GetXaxis()->SetTitle(Form("%s (%s)", ptTitle.Data(), ptUnits.Data()));
282       histo2D->GetYaxis()->SetTitle(Form("%s (%s)", vzTitle.Data(), vzUnits.Data()));
283       histoIndex = GetHistoIndex(kHistoPtVzType, itrig, isrc);
284       fHistoListMC->AddAt(histo2D, histoIndex);
285     }
286   }
287 }
288
289 //________________________________________________________________________
290 void AliAnalysisTaskSingleMu::UserExec(Option_t * /*option*/) 
291 {
292   //
293   /// Main loop
294   /// Called for each event
295   //
296
297   AliAODEvent* aodEvent = 0x0;
298   AliMCEvent*  mcEvent  = 0x0;
299
300   if ( Entry()==0 ) 
301     SetFlagsFromHandler();
302
303   aodEvent = dynamic_cast<AliAODEvent*> (InputEvent());
304   if ( ! aodEvent ) aodEvent = AODEvent();
305   if ( ! aodEvent ) {
306     AliError ("AOD event not found. Nothing done!");
307     return;
308   }
309
310   Int_t nTracks = aodEvent->GetNTracks();
311
312   if ( fUseMC ) mcEvent = MCEvent();
313
314   // Object declaration
315   AliAODTrack *aodTrack = 0x0;
316   AliMCParticle* mcTrack = 0x0;
317
318   const Float_t kDefaultFloat = -999.;
319
320   Float_t pt       = kDefaultFloat;
321   Float_t dca      = kDefaultFloat;
322   Float_t xAtDca   = kDefaultFloat;
323   Float_t yAtDca   = kDefaultFloat;
324   Float_t eta      = kDefaultFloat;
325   Float_t rapidity = kDefaultFloat;
326   Float_t vz       = kDefaultFloat;
327   Int_t trackLabel = -1;
328   Int_t matchTrig = -1;
329
330   for (Int_t itrack = 0; itrack < nTracks; itrack++) {
331
332     // Get variables
333     aodTrack = aodEvent->GetTrack(itrack);
334     if ( aodTrack->GetMostProbablePID() != AliAODTrack::kMuon ) continue;
335     pt = aodTrack->Pt();
336     xAtDca = aodTrack->XAtDCA();
337     yAtDca = aodTrack->YAtDCA();
338     dca = TMath::Sqrt( xAtDca * xAtDca + yAtDca * yAtDca );
339     eta = aodTrack->Eta();
340     rapidity = aodTrack->Y();
341     vz = aodTrack->Zv();
342     trackLabel = aodTrack->GetLabel();
343     matchTrig = aodTrack->GetMatchTrigger();
344
345     // Fill histograms
346     FillTriggerHistos(kHistoPt,       matchTrig, -1, pt);
347     FillTriggerHistos(kHistoDCA,      matchTrig, -1, dca);
348     FillTriggerHistos(kHistoVz,       matchTrig, -1, vz);
349     FillTriggerHistos(kHistoEta,      matchTrig, -1, eta);
350     FillTriggerHistos(kHistoRapidity, matchTrig, -1, rapidity);
351
352     FillTriggerHistos(kHistoPtDCA,      matchTrig, -1, pt, dca);
353     FillTriggerHistos(kHistoPtVz,       matchTrig, -1, pt, vz);
354     FillTriggerHistos(kHistoPtRapidity, matchTrig, -1, pt, rapidity);
355
356     // Monte Carlo part
357     if (! fUseMC ) continue;
358
359     Int_t motherType = kUnknownPart;
360
361     AliMCParticle* matchedMCTrack = 0x0;
362
363     Int_t nMCtracks = mcEvent->GetNumberOfTracks();
364     for(Int_t imctrack=0; imctrack<nMCtracks; imctrack++){
365       mcTrack = (AliMCParticle*)mcEvent->GetTrack(imctrack);
366       if ( trackLabel == mcTrack->GetLabel() ) {
367         matchedMCTrack = mcTrack;
368         break;
369       }
370     } // loop on MC tracks
371
372     if ( matchedMCTrack )
373       motherType = RecoTrackMother(matchedMCTrack, mcEvent);
374
375     AliDebug(1, Form("Found mother %i", motherType));
376
377     if ( motherType != kUnknownPart) {
378       FillTriggerHistos(kHistoPtResolution, matchTrig, motherType, pt - mcTrack->Pt());
379     }
380     FillTriggerHistos(kHistoPtDCAType, matchTrig, motherType, pt, dca);
381     FillTriggerHistos(kHistoPtVzType,  matchTrig, motherType, pt, vz);
382
383   } // loop on tracks
384   
385   // Post final data. It will be written to a file with option "RECREATE"
386   PostData(1, fHistoList);
387   if ( fUseMC ) PostData(2, fHistoListMC);
388 }
389
390 //________________________________________________________________________
391 void AliAnalysisTaskSingleMu::Terminate(Option_t *) {
392   //
393   /// Draw some histogram at the end.
394   //
395   if (!gROOT->IsBatch()) {
396     TCanvas *c1 = new TCanvas("c1","Vz vs Pt",10,10,310,310);
397     c1->SetFillColor(10); c1->SetHighLightColor(10);
398     c1->SetLeftMargin(0.15); c1->SetBottomMargin(0.15);
399     Int_t histoIndex = GetHistoIndex(kHistoPtVz, kAllPtTrig);
400     ((TH2F*)fHistoList->At(histoIndex))->DrawCopy("COLZ");
401   }
402 }
403
404 //________________________________________________________________________
405  void AliAnalysisTaskSingleMu::FillTriggerHistos(Int_t histoIndex, Int_t matchTrig, Int_t motherType,
406                                                 Float_t var1, Float_t var2, Float_t var3)
407 {
408   //
409   /// Fill all histograms passing the trigger cuts
410   //
411
412   Int_t nTrigs = TMath::Max(1, matchTrig);
413   TArrayI trigToFill(nTrigs);
414   trigToFill[0] = matchTrig;
415   for(Int_t itrig = 1; itrig < matchTrig; itrig++){
416     trigToFill[itrig] = itrig;
417   }
418
419   TList* histoList = (motherType < 0 ) ? fHistoList : fHistoListMC;
420
421   TString className;
422   for(Int_t itrig = 0; itrig < nTrigs; itrig++){
423     Int_t currIndex = GetHistoIndex(histoIndex, trigToFill[itrig], motherType);
424     className = histoList->At(currIndex)->ClassName();
425     AliDebug(3, Form("matchTrig %i  Fill %s", matchTrig, histoList->At(currIndex)->GetName()));
426     if (className.Contains("1"))
427       ((TH1F*)histoList->At(currIndex))->Fill(var1);
428     else if (className.Contains("2"))
429       ((TH2F*)histoList->At(currIndex))->Fill(var1, var2);
430     else if (className.Contains("3"))
431       ((TH3F*)histoList->At(currIndex))->Fill(var1, var2, var3);
432     else
433       AliWarning("Histogram not filled");
434   }
435 }
436
437 //________________________________________________________________________
438 Int_t AliAnalysisTaskSingleMu::RecoTrackMother(AliMCParticle* mcTrack, AliMCEvent* mcEvent)
439 {
440   //
441   /// Find track mother from kinematics
442   //
443
444   Int_t recoPdg = mcTrack->PdgCode();
445
446   Bool_t isMuon = (TMath::Abs(recoPdg) == 13) ? kTRUE : kFALSE;
447
448   // Track is not a muon
449   if ( ! isMuon ) return kRecoHadron;
450
451   Int_t imother = mcTrack->GetMother();
452
453   if ( imother<0 ) 
454     return kPrimaryMu; // Drell-Yan Muon
455
456   Int_t igrandma = imother;
457
458   AliMCParticle* motherPart = (AliMCParticle*)mcEvent->GetTrack(imother);
459   Int_t motherPdg = motherPart->PdgCode();
460
461   // Track is an heavy flavor muon
462   Int_t absPdg = TMath::Abs(motherPdg);
463   if(absPdg/100==5 || absPdg/1000==5) {
464     return kBeautyMu;
465   }
466   if(absPdg/100==4 || absPdg/1000==4){
467     Int_t newMother = -1;
468     igrandma = imother;
469     AliInfo("\nFound candidate B->C->mu. History:\n");
470     mcTrack->Print();
471     printf("- %2i   ", igrandma);
472     motherPart->Print();
473     Int_t absGrandMotherPdg = TMath::Abs(motherPart->PdgCode());
474     while ( absGrandMotherPdg > 10 ) {
475       igrandma = ((AliMCParticle*)mcEvent->GetTrack(igrandma))->GetMother();
476       if( igrandma < 0 ) break;
477       printf("- %2i   ", igrandma);
478       mcEvent->GetTrack(igrandma)->Print();
479       absGrandMotherPdg = TMath::Abs(((AliMCParticle*)mcEvent->GetTrack(igrandma))->PdgCode());
480     }
481
482     if (absGrandMotherPdg==5) newMother = kBeautyMu; // Charm from beauty
483     else if (absGrandMotherPdg==4) newMother = kCharmMu;
484
485     if(newMother<0) {
486       AliWarning("Mother not correctly found! Set to charm!\n");
487       newMother = kCharmMu;
488     }
489
490     return newMother;
491   }
492
493   Int_t nPrimaries = mcEvent->Stack()->GetNprimary();
494
495   // Track is a bkg. muon
496   if (imother<nPrimaries) {
497     return kPrimaryMu;
498   }
499   else {
500     return kSecondaryMu;
501   }
502 }
503
504 //________________________________________________________________________
505 Int_t AliAnalysisTaskSingleMu::GetHistoIndex(Int_t histoTypeIndex, Int_t trigIndex, Int_t srcIndex)
506  {
507    if ( srcIndex < 0 ) return kNtrigCuts * histoTypeIndex + trigIndex;
508
509    return 
510      kNtrackSources * kNtrigCuts * histoTypeIndex  + 
511      kNtrackSources * trigIndex  + 
512      srcIndex;
513 }