Improve the utility to apply a different trigger pt cut to the track depending on...
[u/mrichter/AliRoot.git] / PWG / muon / AliVAnalysisMuon.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: AliVAnalysisMuon.cxx 47782 2011-02-24 18:37:31Z martinez $ */
17
18 //-----------------------------------------------------------------------------
19 /// \class AliVAnalysisMuon
20 /// Base class with utilities for muon analysis
21 ///
22 /// \author Diego Stocco
23 //-----------------------------------------------------------------------------
24
25 #include "AliVAnalysisMuon.h"
26
27 // ROOT includes
28 #include "TROOT.h"
29 #include "TH1.h"
30 #include "TH2.h"
31 #include "TAxis.h"
32 #include "TCanvas.h"
33 #include "TLegend.h"
34 #include "TMath.h"
35 #include "TObjString.h"
36 #include "TObjArray.h"
37 #include "THashList.h"
38 #include "TStyle.h"
39 //#include "TMCProcess.h"
40 #include "TLorentzVector.h"
41
42 // STEER includes
43 #include "AliInputEventHandler.h"
44 #include "AliCentrality.h"
45
46 #include "AliAODEvent.h"
47 #include "AliAODTrack.h"
48 #include "AliAODMCParticle.h"
49 #include "AliMCEvent.h"
50 #include "AliMCParticle.h"
51 //#include "AliStack.h"
52 #include "AliESDEvent.h"
53 #include "AliESDMuonTrack.h"
54 #include "AliCounterCollection.h"
55 #include "AliVVertex.h"
56
57 // ANALYSIS includes
58 #include "AliAnalysisManager.h"
59 #include "AliAnalysisTaskSE.h"
60 #include "AliAnalysisDataSlot.h"
61 #include "AliAnalysisDataContainer.h"
62
63 // CORRFW includes
64 #include "AliCFGridSparse.h"
65
66 // PWG3 includes
67 #include "AliMergeableCollection.h"
68 #include "AliMuonTrackCuts.h"
69 #include "AliMuonPairCuts.h"
70
71 /// \cond CLASSIMP
72 ClassImp(AliVAnalysisMuon) // Class implementation in ROOT context
73 /// \endcond
74
75
76 //________________________________________________________________________
77 AliVAnalysisMuon::AliVAnalysisMuon() :
78   AliAnalysisTaskSE(),
79   fMuonTrackCuts(0x0),
80   fMuonPairCuts(0x0),
81   fESDEvent(0x0),
82   fAODEvent(0x0),
83   fTerminateOptions(0x0),
84   fChargeKeys(0x0),
85   fSrcKeys(0x0),
86   fPhysSelKeys(0x0),
87   fTriggerClasses(0x0),
88   fCentralityClasses(0x0),
89   fEventCounters(0x0),
90   fMergeableCollection(0x0),
91   fOutputList(0x0),
92   fMinNvtxContirbutors(0),
93   fSelectedTrigPattern(0x0),
94   fRejectedTrigPattern(0x0),
95   fSelectedTrigLevel(0x0),
96   fOutputPrototypeList(0x0)
97 {
98   /// Default ctor.
99 }
100
101 //________________________________________________________________________
102 AliVAnalysisMuon::AliVAnalysisMuon(const char *name, const AliMuonTrackCuts& trackCuts, const AliMuonPairCuts& pairCuts) :
103   AliAnalysisTaskSE(name),
104   fMuonTrackCuts(new AliMuonTrackCuts(trackCuts)),
105   fMuonPairCuts(new AliMuonPairCuts(pairCuts)),
106   fESDEvent(0x0),
107   fAODEvent(0x0),
108   fTerminateOptions(0x0),
109   fChargeKeys(0x0),
110   fSrcKeys(0x0),
111   fPhysSelKeys(0x0),
112   fTriggerClasses(new THashList()),
113   fCentralityClasses(0x0),
114   fEventCounters(0x0),
115   fMergeableCollection(0x0),
116   fOutputList(0x0),
117   fMinNvtxContirbutors(1),
118   fSelectedTrigPattern(new TObjArray()),
119   fRejectedTrigPattern(new TObjArray()),
120   fSelectedTrigLevel(new TObjArray()),
121   fOutputPrototypeList(0x0)
122 {
123   //
124   /// Constructor.
125   //
126   
127   fTriggerClasses->SetOwner();
128   InitKeys();
129   SetTrigClassPatterns();
130   SetTrigClassLevels();
131   SetCentralityClasses();
132
133   DefineOutput(1, TObjArray::Class());
134 }
135
136 //________________________________________________________________________
137 AliVAnalysisMuon::AliVAnalysisMuon(const char *name, const AliMuonTrackCuts& trackCuts) :
138   AliAnalysisTaskSE(name),
139   fMuonTrackCuts(new AliMuonTrackCuts(trackCuts)),
140   fMuonPairCuts(0x0),
141   fESDEvent(0x0),
142   fAODEvent(0x0),
143   fTerminateOptions(0x0),
144   fChargeKeys(0x0),
145   fSrcKeys(0x0),
146   fPhysSelKeys(0x0),
147   fTriggerClasses(new THashList()),
148   fCentralityClasses(0x0),
149   fEventCounters(0x0),
150   fMergeableCollection(0x0),
151   fOutputList(0x0),
152   fMinNvtxContirbutors(1),
153   fSelectedTrigPattern(new TObjArray()),
154   fRejectedTrigPattern(new TObjArray()),
155   fSelectedTrigLevel(new TObjArray()),
156   fOutputPrototypeList(0x0)
157 {
158   //
159   /// Constructor.
160   //
161   
162   fTriggerClasses->SetOwner();
163   InitKeys();
164   SetTrigClassPatterns();
165   SetTrigClassLevels();
166   SetCentralityClasses();
167   
168   DefineOutput(1, TObjArray::Class());
169 }
170
171
172 //________________________________________________________________________
173 AliVAnalysisMuon::AliVAnalysisMuon(const char *name, const AliMuonPairCuts& pairCuts) :
174   AliAnalysisTaskSE(name),
175   fMuonTrackCuts(0x0),
176   fMuonPairCuts(new AliMuonPairCuts(pairCuts)),
177   fESDEvent(0x0),
178   fAODEvent(0x0),
179   fTerminateOptions(0x0),
180   fChargeKeys(0x0),
181   fSrcKeys(0x0),
182   fPhysSelKeys(0x0),
183   fTriggerClasses(new THashList()),
184   fCentralityClasses(0x0),
185   fEventCounters(0x0),
186   fMergeableCollection(0x0),
187   fOutputList(0x0),
188   fMinNvtxContirbutors(1),
189   fSelectedTrigPattern(new TObjArray()),
190   fRejectedTrigPattern(new TObjArray()),
191   fSelectedTrigLevel(new TObjArray()),
192   fOutputPrototypeList(0x0)
193 {
194   //
195   /// Constructor.
196   //
197   fTriggerClasses->SetOwner();
198   InitKeys();
199   SetTrigClassPatterns();
200   SetTrigClassLevels();
201   SetCentralityClasses();
202     
203   DefineOutput(1, TObjArray::Class());
204 }
205
206
207 //________________________________________________________________________
208 AliVAnalysisMuon::~AliVAnalysisMuon()
209 {
210   //
211   /// Destructor
212   //
213
214   delete fMuonTrackCuts;
215   delete fMuonPairCuts;
216   delete fTerminateOptions;
217   delete fChargeKeys;
218   delete fSrcKeys;
219   delete fPhysSelKeys;
220   delete fTriggerClasses;
221   delete fCentralityClasses;
222   delete fSelectedTrigPattern;
223   delete fRejectedTrigPattern;
224   delete fSelectedTrigLevel;
225   delete fOutputPrototypeList;
226
227
228   // For proof: do not delete output containers
229   if ( ! AliAnalysisManager::GetAnalysisManager() || ! AliAnalysisManager::GetAnalysisManager()->IsProofMode() ) {
230     delete fOutputList;
231   }
232 }
233
234 //___________________________________________________________________________
235 void AliVAnalysisMuon::FinishTaskOutput()
236 {
237   //
238   /// Remove empty histograms to reduce the number of histos to be merged
239   //
240
241
242   fMergeableCollection->PruneEmptyObjects();
243
244   TString objectName = "";
245   
246   // Add stat. info from physics selection
247   // (usefull when running on AODs)
248   if ( fInputHandler ) {
249     for ( Int_t istat=0; istat<2; istat++ ) {
250       TString statType = ( istat == 0 ) ? "ALL" : "BIN0";
251       TH2* hStat = dynamic_cast<TH2*>(fInputHandler->GetStatistics(statType.Data()));
252       if ( hStat ) {
253         objectName = Form("%s_%s", hStat->GetName(), GetName());
254         TH2* cloneStat = static_cast<TH2*>(hStat->Clone(objectName.Data()));
255         cloneStat->SetDirectory(0);
256         fOutputList->Add(cloneStat);
257       }
258       else {
259         AliWarning("Stat histogram not available");
260         break;
261       }
262     } // loop on stat type
263   }
264 }
265
266
267 //___________________________________________________________________________
268 void AliVAnalysisMuon::NotifyRun()
269 {
270   /// Set run number for cuts
271   if ( fMuonTrackCuts ) fMuonTrackCuts->SetRun(fCurrentRunNumber);
272   if ( fMuonPairCuts ) fMuonPairCuts->SetRun(fCurrentRunNumber);
273 }
274
275 //___________________________________________________________________________
276 void AliVAnalysisMuon::UserCreateOutputObjects() 
277 {
278   //
279   /// Create output objects
280   //
281   AliInfo(Form("   CreateOutputObjects of task %s\n", GetName()));
282   
283   fOutputList = new TObjArray();
284   fOutputList->SetOwner();
285
286   fEventCounters = new AliCounterCollection("eventCounters");
287
288   if ( ! fCentralityClasses ) SetCentralityClasses();
289   TString centralityClasses = "";
290   for ( Int_t icent=1; icent<=fCentralityClasses->GetNbins(); ++icent ) {
291     if ( ! centralityClasses.IsNull() ) centralityClasses += "/";
292     centralityClasses += fCentralityClasses->GetBinLabel(icent);
293   }
294   fEventCounters->AddRubric("selected", "yes/no");
295   fEventCounters->AddRubric("trigger", 100);
296   fEventCounters->AddRubric("centrality", centralityClasses);
297   fEventCounters->Init();
298   fOutputList->Add(fEventCounters);
299  
300   fMergeableCollection = new AliMergeableCollection("outputObjects");
301   fOutputList->Add(fMergeableCollection);
302
303   PostData(1, fOutputList);
304   
305   MyUserCreateOutputObjects();
306 }
307
308
309 //________________________________________________________________________
310 void AliVAnalysisMuon::UserExec(Option_t * /*option*/) 
311 {
312   //
313   /// Main loop
314   /// Called for each event
315   //
316
317   fAODEvent = dynamic_cast<AliAODEvent*> (InputEvent());
318   if ( ! fAODEvent ) 
319     fESDEvent = dynamic_cast<AliESDEvent*> (InputEvent());
320
321   if ( ! fAODEvent && ! fESDEvent ) {
322     AliError ("AOD or ESD event not found. Nothing done!");
323     return;
324   }
325
326   Int_t physSel = ( fInputHandler->IsEventSelected() & AliVEvent::kAny ) ? kPhysSelPass : kPhysSelReject;
327
328   //
329   // Global event info
330   //
331
332   TString firedTrigClasses = ( fAODEvent ) ? fAODEvent->GetFiredTriggerClasses() : fESDEvent->GetFiredTriggerClasses();
333   firedTrigClasses.Prepend("ANY ");
334   AliDebug(2, Form("Fired classes %s", firedTrigClasses.Data()));
335   TObjArray* selectTrigClasses = BuildTriggerClasses(firedTrigClasses);
336   if ( selectTrigClasses->GetEntries() == 0 ) {
337     delete selectTrigClasses;
338     return;
339   }
340
341   Double_t centrality = InputEvent()->GetCentrality()->GetCentralityPercentile("V0M");
342   Int_t centralityBin = fCentralityClasses->FindBin(centrality);
343   TString centralityBinLabel = fCentralityClasses->GetBinLabel(centralityBin);
344
345   TString selKey = ( physSel == kPhysSelPass ) ? "yes" : "no";
346   for ( Int_t itrig=0; itrig<selectTrigClasses->GetEntries(); ++itrig ) {
347     TString trigName = selectTrigClasses->At(itrig)->GetName();
348     fEventCounters->Count(Form("trigger:%s/selected:%s/centrality:%s", trigName.Data(), selKey.Data(), centralityBinLabel.Data()));
349   }
350
351   ProcessEvent(fPhysSelKeys->At(physSel)->GetName(), *selectTrigClasses, fCentralityClasses->GetBinLabel(centralityBin));
352
353   delete selectTrigClasses;
354
355   // Post final data. It will be written to a file with option "RECREATE"
356   PostData(1, fOutputList);
357 }
358
359 //________________________________________________________________________
360 void AliVAnalysisMuon::Terminate(Option_t *)
361 {
362   //
363   /// Draw some histogram at the end.
364   //
365   
366   if ( ! fTerminateOptions ) SetTerminateOptions();
367   
368   if ( gROOT->IsBatch() ) return;
369     
370   fOutputList = dynamic_cast<TObjArray*>(GetOutputData(1));
371   if ( ! fOutputList ) return;
372   fEventCounters = static_cast<AliCounterCollection*>(fOutputList->FindObject("eventCounters"));
373   fMergeableCollection = static_cast<AliMergeableCollection*>(fOutputList->FindObject("outputObjects"));
374   
375   if ( ! fMergeableCollection ) return;
376   AliInfo(Form("Mergeable collection size %g MB", fMergeableCollection->EstimateSize()/1024.0/1024.0));
377   if ( fTerminateOptions->At(3) ) {
378     TString sopt = fTerminateOptions->At(3)->GetName();
379     if ( sopt.Contains("verbose") ) fMergeableCollection->Print("*"); 
380   }
381 }
382
383
384
385 //________________________________________________________________________
386 Int_t AliVAnalysisMuon::GetNTracks()
387 {
388   //
389   /// Return the number of tracks in event
390   //
391   return ( fAODEvent ) ? fAODEvent->GetNTracks() : fESDEvent->GetNumberOfMuonTracks();
392 }
393
394
395 //________________________________________________________________________
396 AliVParticle* AliVAnalysisMuon::GetTrack(Int_t itrack)
397 {
398   //
399   /// Get the current track
400   //
401   AliVParticle* track = 0x0;
402   if ( fAODEvent ) track = fAODEvent->GetTrack(itrack);
403   else track = fESDEvent->GetMuonTrack(itrack);
404   return track;
405 }
406
407 //________________________________________________________________________
408 Double_t AliVAnalysisMuon::MuonMass2() const
409 {
410   /// A usefull constant
411   static Double_t m2 = 1.11636129640000012e-02;
412   return m2;
413 }
414
415 //________________________________________________________________________
416 TLorentzVector AliVAnalysisMuon::GetTrackPair(AliVParticle* track1, AliVParticle* track2) const
417 {
418   //
419   /// Get track pair
420   //
421   
422   AliVParticle* tracks[2] = {track1, track2};
423   
424   TLorentzVector vec[2];
425   for ( Int_t itrack=0; itrack<2; ++itrack ) {
426     Double_t trackP = tracks[itrack]->P();
427     Double_t energy = TMath::Sqrt(trackP*trackP + MuonMass2());
428     vec[itrack].SetPxPyPzE(tracks[itrack]->Px(), tracks[itrack]->Py(), tracks[itrack]->Pz(), energy);
429   }
430   
431   TLorentzVector vecPair = vec[0] + vec[1];
432   return vecPair;
433 }
434
435
436 //________________________________________________________________________
437 Int_t AliVAnalysisMuon::GetNMCTracks()
438 {
439   //
440   /// Return the number of MC tracks in event
441   //
442   Int_t nMCtracks = 0;
443   if ( fMCEvent ) nMCtracks = fMCEvent->GetNumberOfTracks();
444   else if ( fAODEvent ) {
445     TClonesArray* mcArray = (TClonesArray*)fAODEvent->GetList()->FindObject(AliAODMCParticle::StdBranchName());
446     if ( mcArray ) nMCtracks = mcArray->GetEntries();
447   }
448   return nMCtracks;
449 }
450
451 //________________________________________________________________________
452 AliVParticle* AliVAnalysisMuon::GetMCTrack(Int_t trackLabel)
453 {
454   //
455   /// MC information can be provided by the MC input handler
456   /// (mostly when analyising ESDs) or can be found inside AODs
457   /// This method returns the correct one
458   //
459   AliVParticle* mcTrack = 0x0;
460   if ( fMCEvent ) mcTrack = fMCEvent->GetTrack(trackLabel);
461   else if ( fAODEvent ) {
462     TClonesArray* mcArray = (TClonesArray*)fAODEvent->FindListObject(AliAODMCParticle::StdBranchName());
463     if ( mcArray ) mcTrack =  (AliVParticle*)mcArray->At(trackLabel);
464   }
465   if ( ! mcTrack ) AliWarning(Form("No track with label %i!", trackLabel));
466   return mcTrack;
467 }
468
469 //________________________________________________________________________
470 Int_t AliVAnalysisMuon::GetMotherIndex(AliVParticle* mcParticle)
471 {
472   //
473   /// Return the mother index
474   //
475   Int_t imother = ( fMCEvent ) ? ((AliMCParticle*)mcParticle)->GetMother() : ((AliAODMCParticle*)mcParticle)->GetMother();
476   return imother;
477 }
478
479 //________________________________________________________________________
480 Int_t AliVAnalysisMuon::GetDaughterIndex(AliVParticle* mcParticle, Int_t idaughter)
481 {
482   //
483   /// Return the daughter index
484   /// idaughter can be:
485   /// 0 -> first daughter
486   /// 1 -> last daughter
487   //
488   if ( idaughter < 0 || idaughter > 1 ) {
489     AliError(Form("Requested daughter %i Daughter index can be either 0 (first) or 1 (last)", idaughter));
490     return -1;
491   }
492   
493   if ( fMCEvent ) {
494     if ( idaughter == 0 ) return ((AliMCParticle*)mcParticle)->GetFirstDaughter();
495     else return ((AliMCParticle*)mcParticle)->GetLastDaughter();
496   }
497   
498   return ((AliAODMCParticle*)mcParticle)->GetDaughter(idaughter);
499 }
500
501
502
503 //________________________________________________________________________
504 Bool_t AliVAnalysisMuon::IsMC()
505 {
506   //
507   /// Contains MC info
508   //
509   return ( fMCEvent || ( fAODEvent && fAODEvent->FindListObject(AliAODMCParticle::StdBranchName()) ) );
510 }
511
512
513 //________________________________________________________________________
514 Int_t AliVAnalysisMuon::GetParticleType(AliVParticle* track)
515 {
516   //
517   /// Get particle type from mathced MC track
518   //
519   
520   Int_t trackSrc = kUnidentified;
521   Int_t trackLabel = track->GetLabel();
522   if ( trackLabel >= 0 ) {
523     AliVParticle* matchedMCTrack = GetMCTrack(trackLabel);
524     if ( matchedMCTrack ) trackSrc = RecoTrackMother(matchedMCTrack);
525   } // track has MC label
526   return trackSrc;
527 }
528
529
530 //________________________________________________________________________
531 Int_t AliVAnalysisMuon::RecoTrackMother(AliVParticle* mcParticle)
532 {
533   //
534   /// Find track mother from kinematics
535   //
536   
537   Int_t recoPdg = mcParticle->PdgCode();
538   
539   // Track is not a muon
540   if ( TMath::Abs(recoPdg) != 13 ) return kRecoHadron;
541   
542   Int_t imother = GetMotherIndex(mcParticle);
543   
544   Int_t den[3] = {100, 1000, 1};
545   
546   Int_t motherType = kDecayMu;
547   while ( imother >= 0 ) {
548     AliVParticle* part = GetMCTrack(imother);
549     //if ( ! part ) return motherType;
550     
551     Int_t absPdg = TMath::Abs(part->PdgCode());
552     
553     Bool_t isPrimary = ( fMCEvent ) ? ( imother < fMCEvent->GetNumberOfPrimaries() ) : ((AliAODMCParticle*)part)->IsPrimary();
554     
555     if ( isPrimary ) {
556       if ( absPdg == 24 ) return kWbosonMu;
557       
558       for ( Int_t idec=0; idec<3; idec++ ) {
559         Int_t flv = (absPdg%100000)/den[idec];
560         if ( flv > 0 && flv < 4 ) return kDecayMu;
561         else if ( flv == 0 || flv > 5 ) continue;
562         else {
563           if ( den[idec] == 100 ) motherType = kQuarkoniumMu;
564           else if ( flv == 4 ) motherType = kCharmMu;
565           else motherType = kBeautyMu;
566           break; // break loop on pdg code
567           // but continue the global loop to find higher mass HF
568         }
569       } // loop on pdg code
570       if ( absPdg < 10 ) break; // particle loop
571     } // is primary
572     else {
573       if ( part->Zv() < -90. ) {
574         // If hadronic process => secondary
575         //if ( part->GetUniqueID() == kPHadronic ) {
576         return kSecondaryMu;
577         //}
578       }
579     } // is secondary
580     
581     imother = GetMotherIndex(part);
582     
583   } // loop on mothers
584   
585   return motherType;
586 }
587
588
589 //________________________________________________________________________
590 AliVVertex* AliVAnalysisMuon::GetVertexSPD() const
591 {
592   //
593   /// Get vertex SPD
594   //
595   
596   AliVVertex* primaryVertex = ( fAODEvent ) ? (AliVVertex*)fAODEvent->GetPrimaryVertexSPD() : (AliVVertex*)fESDEvent->GetPrimaryVertexSPD();
597   return primaryVertex;
598 }
599
600
601 //________________________________________________________________________
602 Bool_t AliVAnalysisMuon::AddObjectToCollection(TObject* object, Int_t index)
603 {
604   //
605   /// Add object to collection
606   //
607   
608   if ( ! fOutputPrototypeList ) {
609     fOutputPrototypeList = new TObjArray();
610     fOutputPrototypeList->SetOwner();
611   }
612   if ( fOutputPrototypeList->FindObject(object->GetName() ) ) {
613     AliWarning(Form("Object with name %s already in the list", object->GetName()));
614     return kFALSE;
615   }
616   if ( index < 0 ) fOutputPrototypeList->Add(object);
617   else fOutputPrototypeList->AddAtAndExpand(object, index);
618   
619   return kTRUE;
620 }
621
622 //________________________________________________________________________
623 TObject* AliVAnalysisMuon::GetMergeableObject(TString physSel, TString trigClassName, TString centrality, TString objectName)
624 {
625   //
626   /// Get mergeable object
627   /// (create collection if necessary)
628   //
629   
630   TString identifier = Form("/%s/%s/%s/", physSel.Data(), trigClassName.Data(), centrality.Data());
631   
632   TObject* obj = fMergeableCollection->GetObject(identifier.Data(), objectName.Data());
633   if ( ! obj ) {
634     CreateMergeableObjects(physSel, trigClassName, centrality);
635     obj = fMergeableCollection->GetObject(identifier.Data(), objectName.Data());
636     AliInfo(Form("Mergeable object collection size %g MB", fMergeableCollection->EstimateSize()/1024.0/1024.0));
637   }
638   return obj;
639 }
640
641 //________________________________________________________________________
642 TObject* AliVAnalysisMuon::GetSum(TString physSel, TString trigClassNames, TString centrality, TString objectPattern)
643 {
644   //
645   /// Sum objects
646   /// - physSel, trigClassNames must be in the form: key1,key2
647   /// - centrality must be in the form minValue_maxValue
648   /// - objectPattern must be in the form match1&match2&match3,match4
649   ///   meaning that the object name must contain match1 and match2 and either one of match3 and match4
650   
651   if ( ! fMergeableCollection ) return 0x0;
652   
653   // Get centrality range
654   Int_t firstCentrality = 1;
655   Int_t lastCentrality = fCentralityClasses->GetNbins();
656   
657   TObjArray* centralityRange = centrality.Tokenize("_");
658   Float_t range[2] = {0., 100.};
659   if ( centralityRange->GetEntries() >= 2 ) {
660     for ( Int_t irange=0; irange<2; ++irange ) {
661       range[irange] = ((TObjString*)centralityRange->At(irange))->GetString().Atof();
662     }
663     firstCentrality = fCentralityClasses->FindBin(range[0]+0.0001);
664     lastCentrality = fCentralityClasses->FindBin(range[1]-0.0001);
665   }
666   delete centralityRange;
667   
668   TString sumCentralityString = "";
669   for ( Int_t icent=firstCentrality; icent<=lastCentrality; ++icent ) {
670     if ( ! sumCentralityString.IsNull() ) sumCentralityString += ",";
671     sumCentralityString += fCentralityClasses->GetBinLabel(icent);
672   }
673   
674   objectPattern.ReplaceAll(" ","");
675   TObjArray* objPatternList = objectPattern.Tokenize("&");
676   TObjArray objPatternMatrix(objPatternList->GetEntries());
677   objPatternMatrix.SetOwner();
678   for ( Int_t ikey=0; ikey<objPatternList->GetEntries(); ikey++ ) {
679     TObjArray* subKeyList = ((TObjString*)objPatternList->At(ikey))->GetString().Tokenize(",");
680     objPatternMatrix.AddAt(subKeyList, ikey);
681   }
682   delete objPatternList;
683   
684
685   TObjArray objectNameInCollection;
686   objectNameInCollection.SetOwner();
687   TObjArray* physSelList = physSel.Tokenize(",");
688   TObjArray* trigClassList = trigClassNames.Tokenize(",");
689   TObjArray* centralityList = sumCentralityString.Tokenize(",");
690   for ( Int_t isel=0; isel<physSelList->GetEntries(); isel++ ) {
691     for ( Int_t itrig = 0; itrig<trigClassList->GetEntries(); itrig++ ) {
692       for ( Int_t icent=0; icent<centralityList->GetEntries(); icent++ ) {
693         TString currId = Form("/%s/%s/%s/", physSelList->At(isel)->GetName(), trigClassList->At(itrig)->GetName(),centralityList->At(icent)->GetName());
694         TList* objNameList = fMergeableCollection->CreateListOfObjectNames(currId.Data());
695         for ( Int_t iobj=0; iobj<objNameList->GetEntries(); iobj++ ) {
696           TString objName = objNameList->At(iobj)->GetName();
697           if ( ! objectNameInCollection.FindObject(objName.Data()) ) objectNameInCollection.Add(new TObjString(objName.Data()));
698         }
699         delete objNameList;
700       }
701     }
702   }
703   delete physSelList;
704   delete trigClassList;
705   delete centralityList;
706
707   TString matchingObjectNames = "";
708   for ( Int_t iobj=0; iobj<objectNameInCollection.GetEntries(); iobj++ ) {
709     TString objName = objectNameInCollection.At(iobj)->GetName();
710     Bool_t matchAnd = kTRUE;
711     for ( Int_t ikey=0; ikey<objPatternMatrix.GetEntries(); ikey++ ) {
712       TObjArray*  subKeyList = (TObjArray*)objPatternMatrix.At(ikey);
713       Bool_t matchOr = kFALSE;
714       for ( Int_t isub=0; isub<subKeyList->GetEntries(); isub++ ) {
715         TString subKeyString = ((TObjString*)subKeyList->At(isub))->GetString();
716         if ( objName.Contains(subKeyString.Data()) ) {
717           matchOr = kTRUE;
718           break;
719         }
720       }
721       if ( ! matchOr ) {
722         matchAnd = kFALSE;
723         break;
724       }
725     }
726     if ( ! matchAnd ) continue;
727     if ( ! matchingObjectNames.IsNull() ) matchingObjectNames.Append(",");
728     matchingObjectNames += objName;
729   }
730
731   TString idPattern = Form("/%s/%s/%s/%s", physSel.Data(), trigClassNames.Data(), sumCentralityString.Data(), matchingObjectNames.Data());
732   idPattern.ReplaceAll(" ","");
733   
734   AliDebug(1,Form("Sum pattern %s", idPattern.Data()));
735   
736   return fMergeableCollection->GetSum(idPattern.Data());
737 }
738
739 //___________________________________________________________________________
740 void AliVAnalysisMuon::CreateMergeableObjects(TString physSel, TString trigClassName, TString centrality)
741 {
742   TObject* obj = 0x0;
743   TString objectName = "";
744   TString identifier = Form("/%s/%s/%s/", physSel.Data(), trigClassName.Data(), centrality.Data());
745   for ( Int_t iobj=0; iobj<fOutputPrototypeList->GetEntries(); ++iobj ) {
746     objectName = fOutputPrototypeList->At(iobj)->GetName();
747     obj = fOutputPrototypeList->At(iobj)->Clone(objectName.Data());
748     fMergeableCollection->Adopt(identifier, obj);
749   } // loop on histos
750 }
751
752
753 //_______________________________________________________________________
754 Bool_t AliVAnalysisMuon::SetSparseRange(AliCFGridSparse* gridSparse,
755                                         Int_t ivar, TString labelName,
756                                         Double_t varMin, Double_t varMax,
757                                         TString option)
758 {
759   //
760   /// Set range in a smart way.
761   /// Allows to select a bin from the label.
762   /// Check the bin limits.
763   //
764   
765   option.ToUpper();
766   Int_t minVarBin = -1, maxVarBin = -1;
767   TAxis* axis = gridSparse->GetAxis(ivar);
768   
769   if ( ! axis ) {
770     printf("Warning: Axis %i not found in %s", ivar, gridSparse->GetName());
771     return kFALSE;
772   }
773   
774   if ( ! labelName.IsNull() ) {
775     minVarBin = axis->FindBin(labelName.Data());
776     maxVarBin = minVarBin;
777     if ( minVarBin < 1 ) {
778       printf("Warning: %s: label %s not found. Nothing done", gridSparse->GetName(), labelName.Data());
779       return kFALSE;
780     }
781   }
782   else if ( option.Contains( "USEBIN" ) ) {
783     minVarBin = (Int_t)varMin;
784     maxVarBin = (Int_t)varMax;
785   }
786   else {
787     minVarBin = axis->FindBin(varMin);
788     maxVarBin = axis->FindBin(varMax);
789   }
790   
791   if ( axis->GetFirst() == minVarBin && axis->GetLast() == maxVarBin ) return kFALSE;
792   
793   gridSparse->SetRangeUser(ivar, axis->GetBinCenter(minVarBin), axis->GetBinCenter(maxVarBin));
794   
795   return kTRUE;
796 }
797
798 //________________________________________________________________________
799 void AliVAnalysisMuon::SetTrigClassPatterns(TString pattern)
800 {
801   /// Set trigger classes
802   ///
803   /// Classes are filled dynamically according to the pattern
804   /// - if name contains ! (without spaces): reject it
805   /// - otherwise, keep it
806   /// example:
807   /// SetTrigClassPatterns("CMBAC !ALLNOTRD")
808   /// keeps classes containing CMBAC, and not containing ALLNOTRD.
809   ///
810   /// CAVEAT: if you use an fCFContainer and you want an axis to contain the trigger classes,
811   ///         please be sure that each pattern matches only 1 trigger class, or triggers will be mixed up
812   ///         when merging different chuncks.
813
814   fSelectedTrigPattern->SetOwner();
815   if ( fSelectedTrigPattern->GetEntries() > 0 ) fSelectedTrigPattern->Delete();
816   fRejectedTrigPattern->SetOwner();
817   if ( fRejectedTrigPattern->GetEntries() > 0 ) fRejectedTrigPattern->Delete();
818
819   pattern.ReplaceAll("  "," ");
820   pattern.ReplaceAll("! ","!");
821
822   TObjArray* fullList = pattern.Tokenize(" ");
823
824   for ( Int_t ipat=0; ipat<fullList->GetEntries(); ++ipat ) {
825     TString currPattern = fullList->At(ipat)->GetName();
826     if ( currPattern.Contains("!") ) {
827       currPattern.ReplaceAll("!","");
828       fRejectedTrigPattern->AddLast(new TObjString(currPattern));
829     }
830     else fSelectedTrigPattern->AddLast(new TObjString(currPattern));
831   }
832   
833   delete fullList;
834 }
835
836 //________________________________________________________________________
837 void AliVAnalysisMuon::SetTrigClassLevels(TString pattern)
838 {
839   /// Set trigger cut level associated to the trigger class
840   ///
841   /// example:
842   /// SetTrigClassLevels("MSL:Lpt MSH:Hpt")
843   ///
844   /// For the trigger classes defined in SetTrigClassPatterns
845   /// it check if they contains the keywords MSL or MSH
846   /// Hence, in the analysis, the function
847   /// TrackPtCutMatchTrigClass(track, "CPBIMSL") returns true if track match Lpt
848   /// TrackPtCutMatchTrigClass(track, "CPBIMSH") returns true if track match Hpt
849   /// TrackPtCutMatchTrigClass(track, "CMBAC") always returns true
850   
851   fSelectedTrigLevel->SetOwner();
852   if ( fSelectedTrigLevel->GetEntries() > 0 ) fSelectedTrigLevel->Delete();
853   
854   pattern.ReplaceAll("  "," ");
855   pattern.ReplaceAll(" :",":");
856   
857   TObjArray* fullList = pattern.Tokenize(" ");
858   
859   for ( Int_t ipat=0; ipat<fullList->GetEntries(); ++ipat ) {
860     TString currPattern = fullList->At(ipat)->GetName();
861     TObjArray* arr = currPattern.Tokenize(":");
862     TObjString* trigClassPattern = new TObjString(arr->At(0)->GetName());
863     TString selTrigLevel = arr->At(1)->GetName();
864     selTrigLevel.ToUpper();
865     UInt_t trigLevel = 0;
866     if ( selTrigLevel.Contains("APT") ) trigLevel = 1;
867     else if ( selTrigLevel.Contains("LPT") ) trigLevel = 2;
868     else if ( selTrigLevel.Contains("HPT") ) trigLevel = 3;
869     trigClassPattern->SetUniqueID(trigLevel);
870     fSelectedTrigLevel->AddLast(trigClassPattern);
871     delete arr;
872   }
873   
874   delete fullList;
875 }
876
877
878 //________________________________________________________________________
879 void AliVAnalysisMuon::SetCentralityClasses(Int_t nCentralityBins, Double_t* centralityBins)
880 {
881   //
882   /// Set centrality classes
883   //
884   Double_t* bins = centralityBins;
885   Int_t nbins = nCentralityBins;
886   
887   Double_t defaultCentralityBins[] = {-5., 0., 5., 10., 20., 30., 40., 50., 60., 70., 80., 100., 105.};
888   if ( ! centralityBins ) {
889     bins = defaultCentralityBins;
890     nbins = sizeof(defaultCentralityBins)/sizeof(defaultCentralityBins[0])-1;
891   }
892
893   if ( fCentralityClasses ) delete fCentralityClasses;
894   fCentralityClasses = new TAxis(nbins, bins);
895   TString currClass = "";
896   for ( Int_t ibin=1; ibin<=fCentralityClasses->GetNbins(); ++ibin ){
897     currClass = Form("%.0f_%.0f",fCentralityClasses->GetBinLowEdge(ibin),fCentralityClasses->GetBinUpEdge(ibin));
898     fCentralityClasses->SetBinLabel(ibin, currClass.Data());
899   }
900 }
901
902 //________________________________________________________________________
903 void AliVAnalysisMuon::SetTerminateOptions(TString physSel, TString trigClass, TString centralityRange, TString furtherOpts)
904 {
905   //
906   /// Set terminate options
907   //
908   if ( ! fTerminateOptions ) {
909     fTerminateOptions = new TObjArray(4);
910     fTerminateOptions->SetOwner();
911   }
912   fTerminateOptions->AddAt(new TObjString(physSel), 0);
913   fTerminateOptions->AddAt(new TObjString(trigClass), 1);
914   fTerminateOptions->AddAt(new TObjString(centralityRange),2);
915   fTerminateOptions->AddLast(new TObjString(furtherOpts));
916 }
917
918 //________________________________________________________________________
919 void AliVAnalysisMuon::InitKeys()
920 {
921   TString chargeKeys = "MuMinus MuPlus";
922   fChargeKeys = chargeKeys.Tokenize(" ");
923   
924   TString srcKeys = "CharmMu BeautyMu QuarkoniumMu WbosonMu DecayMu SecondaryMu Hadron Unidentified";
925   fSrcKeys = srcKeys.Tokenize(" ");
926   
927   TString physSelKeys = "PhysSelPass PhysSelReject";
928   fPhysSelKeys = physSelKeys.Tokenize(" ");
929 }
930
931 //________________________________________________________________________
932 TObjArray* AliVAnalysisMuon::BuildTriggerClasses(TString firedTrigClasses)
933 {
934   //
935   /// Return the list of trigger classes to be considered
936   /// for current event. Update the global list if necessary
937   //
938
939   TObjArray* selectedTrigClasses = new TObjArray(0);
940   selectedTrigClasses->SetOwner();
941   
942   TObjArray* firedTrigClassesList = firedTrigClasses.Tokenize(" ");
943   
944   UInt_t trigLevel = 0;
945   for ( Int_t itrig=0; itrig<firedTrigClassesList->GetEntries(); ++itrig ) {
946     TString trigName = ((TObjString*)firedTrigClassesList->At(itrig))->GetString();
947     
948     TObject* foundTrig = fTriggerClasses->FindObject(trigName.Data());
949     if ( foundTrig ) trigLevel = foundTrig->GetUniqueID();
950     else {
951       Bool_t rejectTrig = kFALSE;
952       for ( Int_t ipat=0; ipat<fRejectedTrigPattern->GetEntries(); ++ipat ) {
953         if ( trigName.Contains(fRejectedTrigPattern->At(ipat)->GetName() ) ) {
954           rejectTrig = kTRUE;
955           break;
956         }
957       } // loop on reject pattern
958       if ( rejectTrig ) continue;
959
960       rejectTrig = kTRUE;
961       for ( Int_t ipat=0; ipat<fSelectedTrigPattern->GetEntries(); ++ipat ) {
962         if ( trigName.Contains(fSelectedTrigPattern->At(ipat)->GetName() ) ) {
963           rejectTrig = kFALSE;
964           break;
965         }
966       } // loop on keep pattern
967       if ( rejectTrig ) continue;
968     
969       trigLevel = 0;
970       for ( Int_t ipat=0; ipat<fSelectedTrigLevel->GetEntries(); ++ipat ) {
971         if ( trigName.Contains(fSelectedTrigLevel->At(ipat)->GetName() ) ) {
972           trigLevel = fSelectedTrigLevel->At(ipat)->GetUniqueID();
973           break;
974         }
975       } // loop on trig level patterns      
976     }
977     TObjString* currTrig = new TObjString(trigName);
978     currTrig->SetUniqueID(trigLevel);
979     selectedTrigClasses->AddLast(currTrig);
980     
981     if ( foundTrig ) continue;
982     AliInfo(Form("Adding %s (trig level %u) to considered trigger classes",trigName.Data(),trigLevel));
983     TObjString* addTrig = new TObjString(trigName);
984     addTrig->SetUniqueID(trigLevel);
985     fTriggerClasses->Add(addTrig);
986   } // loop on trigger classes
987
988   delete firedTrigClassesList;
989
990   return selectedTrigClasses;
991 }
992
993
994 //________________________________________________________________________
995 Bool_t AliVAnalysisMuon::TrackPtCutMatchTrigClass(AliVParticle* track, TString trigClassName)
996 {
997   /// Check if track passes the trigger pt cut level used in the trigger class
998   Int_t matchTrig = ( fAODEvent ) ? ((AliAODTrack*)track)->GetMatchTrigger() : ((AliESDMuonTrack*)track)->GetMatchTrigger();
999   Int_t classMatchLevel = GetTrigClassPtCutLevel(trigClassName);
1000   Bool_t matchTrackerPt = kTRUE;
1001   if ( fMuonTrackCuts && fMuonTrackCuts->IsApplySharpPtCutInMatching() ) {
1002     matchTrackerPt = ( track->Pt() >= fMuonTrackCuts->GetSharpPtCut(classMatchLevel-1,kFALSE) );
1003   }
1004   Bool_t passCut = ( ( matchTrig >= classMatchLevel ) && matchTrackerPt );
1005   AliDebug(1,Form("Class %s  matchTrig %i  trackMatchTrig %i  trackPt %g (required %i)  passCut %i", trigClassName.Data(), classMatchLevel, matchTrig, track->Pt(), fMuonTrackCuts->IsApplySharpPtCutInMatching(),passCut));
1006   return passCut;
1007 }
1008
1009
1010 //________________________________________________________________________
1011 Int_t AliVAnalysisMuon::GetTrigClassPtCutLevel(TString trigClassName)
1012 {
1013   /// Get trigger class pt cut level for tracking/trigger matching
1014   TObject* obj = fTriggerClasses->FindObject(trigClassName.Data());
1015   if ( ! obj ) {
1016     AliWarning(Form("Class %s not in the list!", trigClassName.Data()));
1017     return -1;
1018   }
1019   
1020   return obj->GetUniqueID();
1021 }