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