ba87c06ae81f06c1ee9606413be41b7c6fd75f74
[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()->GetCentralityPercentileUnchecked("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 ( ! fTerminateOptions ) SetTerminateOptions();
364   
365   if ( gROOT->IsBatch() ) return;
366     
367   fOutputList = dynamic_cast<TObjArray*>(GetOutputData(1));
368   if ( ! fOutputList ) return;
369   fEventCounters = static_cast<AliCounterCollection*>(fOutputList->FindObject("eventCounters"));
370   fMergeableCollection = static_cast<AliMergeableCollection*>(fOutputList->FindObject("outputObjects"));
371   
372   if ( ! fMergeableCollection ) return;
373   AliInfo(Form("Mergeable collection size %g MB", fMergeableCollection->EstimateSize()/1024.0/1024.0));
374   if ( fTerminateOptions->At(3) ) {
375     TString sopt = fTerminateOptions->At(3)->GetName();
376     if ( sopt.Contains("verbose") ) fMergeableCollection->Print("*"); 
377   }
378 }
379
380
381
382 //________________________________________________________________________
383 Int_t AliVAnalysisMuon::GetNTracks()
384 {
385   //
386   /// Return the number of tracks in event
387   //
388   return ( fAODEvent ) ? fAODEvent->GetNTracks() : fESDEvent->GetNumberOfMuonTracks();
389 }
390
391
392 //________________________________________________________________________
393 AliVParticle* AliVAnalysisMuon::GetTrack(Int_t itrack)
394 {
395   //
396   /// Get the current track
397   //
398   AliVParticle* track = 0x0;
399   if ( fAODEvent ) track = fAODEvent->GetTrack(itrack);
400   else track = fESDEvent->GetMuonTrack(itrack);
401   return track;
402 }
403
404 //________________________________________________________________________
405 Double_t AliVAnalysisMuon::MuonMass2() const
406 {
407   /// A usefull constant
408   static Double_t m2 = 1.11636129640000012e-02;
409   return m2;
410 }
411
412 //________________________________________________________________________
413 TLorentzVector AliVAnalysisMuon::GetTrackPair(AliVParticle* track1, AliVParticle* track2) const
414 {
415   //
416   /// Get track pair
417   //
418   
419   AliVParticle* tracks[2] = {track1, track2};
420   
421   TLorentzVector vec[2];
422   for ( Int_t itrack=0; itrack<2; ++itrack ) {
423     Double_t trackP = tracks[itrack]->P();
424     Double_t energy = TMath::Sqrt(trackP*trackP + MuonMass2());
425     vec[itrack].SetPxPyPzE(tracks[itrack]->Px(), tracks[itrack]->Py(), tracks[itrack]->Pz(), energy);
426   }
427   
428   TLorentzVector vecPair = vec[0] + vec[1];
429   return vecPair;
430 }
431
432
433 //________________________________________________________________________
434 Int_t AliVAnalysisMuon::GetNMCTracks()
435 {
436   //
437   /// Return the number of MC tracks in event
438   //
439   Int_t nMCtracks = 0;
440   if ( fMCEvent ) nMCtracks = fMCEvent->GetNumberOfTracks();
441   else if ( fAODEvent ) {
442     TClonesArray* mcArray = (TClonesArray*)fAODEvent->GetList()->FindObject(AliAODMCParticle::StdBranchName());
443     if ( mcArray ) nMCtracks = mcArray->GetEntries();
444   }
445   return nMCtracks;
446 }
447
448 //________________________________________________________________________
449 AliVParticle* AliVAnalysisMuon::GetMCTrack(Int_t trackLabel)
450 {
451   //
452   /// MC information can be provided by the MC input handler
453   /// (mostly when analyising ESDs) or can be found inside AODs
454   /// This method returns the correct one
455   //
456   AliVParticle* mcTrack = 0x0;
457   if ( fMCEvent ) mcTrack = fMCEvent->GetTrack(trackLabel);
458   else if ( fAODEvent ) {
459     TClonesArray* mcArray = (TClonesArray*)fAODEvent->FindListObject(AliAODMCParticle::StdBranchName());
460     if ( mcArray ) mcTrack =  (AliVParticle*)mcArray->At(trackLabel);
461   }
462   if ( ! mcTrack ) AliWarning(Form("No track with label %i!", trackLabel));
463   return mcTrack;
464 }
465
466 //________________________________________________________________________
467 Int_t AliVAnalysisMuon::GetMotherIndex(AliVParticle* mcParticle)
468 {
469   //
470   /// Return the mother index
471   //
472   Int_t imother = ( fMCEvent ) ? ((AliMCParticle*)mcParticle)->GetMother() : ((AliAODMCParticle*)mcParticle)->GetMother();
473   return imother;
474 }
475
476 //________________________________________________________________________
477 Int_t AliVAnalysisMuon::GetDaughterIndex(AliVParticle* mcParticle, Int_t idaughter)
478 {
479   //
480   /// Return the daughter index
481   /// idaughter can be:
482   /// 0 -> first daughter
483   /// 1 -> last daughter
484   //
485   if ( idaughter < 0 || idaughter > 1 ) {
486     AliError(Form("Requested daughter %i Daughter index can be either 0 (first) or 1 (last)", idaughter));
487     return -1;
488   }
489   
490   if ( fMCEvent ) {
491     if ( idaughter == 0 ) return ((AliMCParticle*)mcParticle)->GetFirstDaughter();
492     else return ((AliMCParticle*)mcParticle)->GetLastDaughter();
493   }
494   
495   return ((AliAODMCParticle*)mcParticle)->GetDaughter(idaughter);
496 }
497
498
499
500 //________________________________________________________________________
501 Bool_t AliVAnalysisMuon::IsMC()
502 {
503   //
504   /// Contains MC info
505   //
506   return ( fMCEvent || ( fAODEvent && fAODEvent->FindListObject(AliAODMCParticle::StdBranchName()) ) );
507 }
508
509
510 //________________________________________________________________________
511 Int_t AliVAnalysisMuon::GetParticleType(AliVParticle* track)
512 {
513   //
514   /// Get particle type from mathced MC track
515   //
516   
517   Int_t trackSrc = kUnidentified;
518   Int_t trackLabel = track->GetLabel();
519   if ( trackLabel >= 0 ) {
520     AliVParticle* matchedMCTrack = GetMCTrack(trackLabel);
521     if ( matchedMCTrack ) trackSrc = RecoTrackMother(matchedMCTrack);
522   } // track has MC label
523   return trackSrc;
524 }
525
526
527 //________________________________________________________________________
528 Int_t AliVAnalysisMuon::RecoTrackMother(AliVParticle* mcParticle)
529 {
530   //
531   /// Find track mother from kinematics
532   //
533   
534   Int_t recoPdg = mcParticle->PdgCode();
535   
536   // Track is not a muon
537   if ( TMath::Abs(recoPdg) != 13 ) return kRecoHadron;
538   
539   Int_t imother = GetMotherIndex(mcParticle);
540   
541   Int_t den[3] = {100, 1000, 1};
542   
543   Int_t motherType = kDecayMu;
544   while ( imother >= 0 ) {
545     AliVParticle* part = GetMCTrack(imother);
546     //if ( ! part ) return motherType;
547     
548     Int_t absPdg = TMath::Abs(part->PdgCode());
549     
550     Bool_t isPrimary = ( fMCEvent ) ? ( imother < fMCEvent->GetNumberOfPrimaries() ) : ((AliAODMCParticle*)part)->IsPrimary();
551     
552     if ( isPrimary ) {
553       if ( absPdg == 24 ) return kWbosonMu;
554       
555       for ( Int_t idec=0; idec<3; idec++ ) {
556         Int_t flv = (absPdg%100000)/den[idec];
557         if ( flv > 0 && flv < 4 ) return kDecayMu;
558         else if ( flv == 0 || flv > 5 ) continue;
559         else {
560           if ( den[idec] == 100 ) motherType = kQuarkoniumMu;
561           else if ( flv == 4 ) motherType = kCharmMu;
562           else motherType = kBeautyMu;
563           break; // break loop on pdg code
564           // but continue the global loop to find higher mass HF
565         }
566       } // loop on pdg code
567       if ( absPdg < 10 ) break; // particle loop
568     } // is primary
569     else {
570       if ( part->Zv() < -90. ) {
571         // If hadronic process => secondary
572         //if ( part->GetUniqueID() == kPHadronic ) {
573         return kSecondaryMu;
574         //}
575       }
576     } // is secondary
577     
578     imother = GetMotherIndex(part);
579     
580   } // loop on mothers
581   
582   return motherType;
583 }
584
585
586 //________________________________________________________________________
587 AliVVertex* AliVAnalysisMuon::GetVertexSPD() const
588 {
589   //
590   /// Get vertex SPD
591   //
592   
593   AliVVertex* primaryVertex = ( fAODEvent ) ? (AliVVertex*)fAODEvent->GetPrimaryVertexSPD() : (AliVVertex*)fESDEvent->GetPrimaryVertexSPD();
594   return primaryVertex;
595 }
596
597
598 //________________________________________________________________________
599 Bool_t AliVAnalysisMuon::AddObjectToCollection(TObject* object, Int_t index)
600 {
601   //
602   /// Add object to collection
603   //
604   
605   if ( ! fOutputPrototypeList ) {
606     fOutputPrototypeList = new TObjArray();
607     fOutputPrototypeList->SetOwner();
608   }
609   if ( fOutputPrototypeList->FindObject(object->GetName() ) ) {
610     AliWarning(Form("Object with name %s already in the list", object->GetName()));
611     return kFALSE;
612   }
613   if ( index < 0 ) fOutputPrototypeList->Add(object);
614   else fOutputPrototypeList->AddAtAndExpand(object, index);
615   
616   return kTRUE;
617 }
618
619 //________________________________________________________________________
620 TObject* AliVAnalysisMuon::GetMergeableObject(TString physSel, TString trigClassName, TString centrality, TString objectName)
621 {
622   //
623   /// Get mergeable object
624   /// (create collection if necessary)
625   //
626   
627   TString identifier = Form("/%s/%s/%s/", physSel.Data(), trigClassName.Data(), centrality.Data());
628   
629   TObject* obj = fMergeableCollection->GetObject(identifier.Data(), objectName.Data());
630   if ( ! obj ) {
631     CreateMergeableObjects(physSel, trigClassName, centrality);
632     obj = fMergeableCollection->GetObject(identifier.Data(), objectName.Data());
633     AliInfo(Form("Mergeable object collection size %g MB", fMergeableCollection->EstimateSize()/1024.0/1024.0));
634   }
635   return obj;
636 }
637
638 //________________________________________________________________________
639 TObject* AliVAnalysisMuon::GetSum(TString physSel, TString trigClassNames, TString centrality, TString objectPattern)
640 {
641   //
642   /// Sum objects
643   /// - physSel, trigClassNames must be in the form: key1,key2
644   /// - centrality must be in the form minValue_maxValue
645   /// - objectPattern must be in the form match1&match2&match3,match4
646   ///   meaning that the object name must contain match1 and match2 and either one of match3 and match4
647   
648   if ( ! fMergeableCollection ) return 0x0;
649   
650   // Get centrality range
651   Int_t firstCentrality = 1;
652   Int_t lastCentrality = fCentralityClasses->GetNbins();
653   
654   TObjArray* centralityRange = centrality.Tokenize("_");
655   Float_t range[2] = {0., 100.};
656   if ( centralityRange->GetEntries() >= 2 ) {
657     for ( Int_t irange=0; irange<2; ++irange ) {
658       range[irange] = ((TObjString*)centralityRange->At(irange))->GetString().Atof();
659     }
660     firstCentrality = fCentralityClasses->FindBin(range[0]+0.0001);
661     lastCentrality = fCentralityClasses->FindBin(range[1]-0.0001);
662   }
663   delete centralityRange;
664   
665   TString sumCentralityString = "";
666   for ( Int_t icent=firstCentrality; icent<=lastCentrality; ++icent ) {
667     if ( ! sumCentralityString.IsNull() ) sumCentralityString += ",";
668     sumCentralityString += fCentralityClasses->GetBinLabel(icent);
669   }
670   
671   objectPattern.ReplaceAll(" ","");
672   TObjArray* objPatternList = objectPattern.Tokenize("&");
673   TObjArray objPatternMatrix(objPatternList->GetEntries());
674   objPatternMatrix.SetOwner();
675   for ( Int_t ikey=0; ikey<objPatternList->GetEntries(); ikey++ ) {
676     TObjArray* subKeyList = ((TObjString*)objPatternList->At(ikey))->GetString().Tokenize(",");
677     objPatternMatrix.AddAt(subKeyList, ikey);
678   }
679   delete objPatternList;
680   
681
682   TObjArray objectNameInCollection;
683   objectNameInCollection.SetOwner();
684   TObjArray* physSelList = physSel.Tokenize(",");
685   TObjArray* trigClassList = trigClassNames.Tokenize(",");
686   TObjArray* centralityList = sumCentralityString.Tokenize(",");
687   for ( Int_t isel=0; isel<physSelList->GetEntries(); isel++ ) {
688     for ( Int_t itrig = 0; itrig<trigClassList->GetEntries(); itrig++ ) {
689       for ( Int_t icent=0; icent<centralityList->GetEntries(); icent++ ) {
690         TString currId = Form("/%s/%s/%s/", physSelList->At(isel)->GetName(), trigClassList->At(itrig)->GetName(),centralityList->At(icent)->GetName());
691         TList* objNameList = fMergeableCollection->CreateListOfObjectNames(currId.Data());
692         for ( Int_t iobj=0; iobj<objNameList->GetEntries(); iobj++ ) {
693           TString objName = objNameList->At(iobj)->GetName();
694           if ( ! objectNameInCollection.FindObject(objName.Data()) ) objectNameInCollection.Add(new TObjString(objName.Data()));
695         }
696         delete objNameList;
697       }
698     }
699   }
700   delete physSelList;
701   delete trigClassList;
702   delete centralityList;
703
704   TString matchingObjectNames = "";
705   for ( Int_t iobj=0; iobj<objectNameInCollection.GetEntries(); iobj++ ) {
706     TString objName = objectNameInCollection.At(iobj)->GetName();
707     Bool_t matchAnd = kTRUE;
708     for ( Int_t ikey=0; ikey<objPatternMatrix.GetEntries(); ikey++ ) {
709       TObjArray*  subKeyList = (TObjArray*)objPatternMatrix.At(ikey);
710       Bool_t matchOr = kFALSE;
711       for ( Int_t isub=0; isub<subKeyList->GetEntries(); isub++ ) {
712         TString subKeyString = ((TObjString*)subKeyList->At(isub))->GetString();
713         if ( objName.Contains(subKeyString.Data()) ) {
714           matchOr = kTRUE;
715           break;
716         }
717       }
718       if ( ! matchOr ) {
719         matchAnd = kFALSE;
720         break;
721       }
722     }
723     if ( ! matchAnd ) continue;
724     if ( ! matchingObjectNames.IsNull() ) matchingObjectNames.Append(",");
725     matchingObjectNames += objName;
726   }
727
728   TString idPattern = Form("/%s/%s/%s/%s", physSel.Data(), trigClassNames.Data(), sumCentralityString.Data(), matchingObjectNames.Data());
729   idPattern.ReplaceAll(" ","");
730   
731   AliDebug(1,Form("Sum pattern %s", idPattern.Data()));
732   
733   return fMergeableCollection->GetSum(idPattern.Data());
734 }
735
736 //___________________________________________________________________________
737 void AliVAnalysisMuon::CreateMergeableObjects(TString physSel, TString trigClassName, TString centrality)
738 {
739   TObject* obj = 0x0;
740   TString objectName = "";
741   TString identifier = Form("/%s/%s/%s/", physSel.Data(), trigClassName.Data(), centrality.Data());
742   for ( Int_t iobj=0; iobj<fOutputPrototypeList->GetEntries(); ++iobj ) {
743     objectName = fOutputPrototypeList->At(iobj)->GetName();
744     obj = fOutputPrototypeList->At(iobj)->Clone(objectName.Data());
745     fMergeableCollection->Adopt(identifier, obj);
746   } // loop on histos
747 }
748
749
750 //_______________________________________________________________________
751 Bool_t AliVAnalysisMuon::SetSparseRange(AliCFGridSparse* gridSparse,
752                                         Int_t ivar, TString labelName,
753                                         Double_t varMin, Double_t varMax,
754                                         TString option)
755 {
756   //
757   /// Set range in a smart way.
758   /// Allows to select a bin from the label.
759   /// Check the bin limits.
760   //
761   
762   option.ToUpper();
763   Int_t minVarBin = -1, maxVarBin = -1;
764   TAxis* axis = gridSparse->GetAxis(ivar);
765   
766   if ( ! axis ) {
767     printf("Warning: Axis %i not found in %s", ivar, gridSparse->GetName());
768     return kFALSE;
769   }
770   
771   if ( ! labelName.IsNull() ) {
772     minVarBin = axis->FindBin(labelName.Data());
773     maxVarBin = minVarBin;
774     if ( minVarBin < 1 ) {
775       printf("Warning: %s: label %s not found. Nothing done", gridSparse->GetName(), labelName.Data());
776       return kFALSE;
777     }
778   }
779   else if ( option.Contains( "USEBIN" ) ) {
780     minVarBin = (Int_t)varMin;
781     maxVarBin = (Int_t)varMax;
782   }
783   else {
784     minVarBin = axis->FindBin(varMin);
785     maxVarBin = axis->FindBin(varMax);
786   }
787   
788   if ( axis->GetFirst() == minVarBin && axis->GetLast() == maxVarBin ) return kFALSE;
789   
790   gridSparse->SetRangeUser(ivar, axis->GetBinCenter(minVarBin), axis->GetBinCenter(maxVarBin));
791   
792   return kTRUE;
793 }
794
795 //________________________________________________________________________
796 void AliVAnalysisMuon::SetTrigClassPatterns(TString pattern)
797 {
798   /// Set trigger classes
799   ///
800   /// Classes are filled dynamically according to the pattern
801   /// - if name contains ! (without spaces): reject it
802   /// - in the matching pattern it is also possible to specify the
803   ///   pt cut level associated to the trigger
804   /// example:
805   /// SetTrigClassPatterns("CMBAC CPBI1MSL:Lpt CPBI1MSH:Hpt !ALLNOTRD")
806   /// keeps classes containing CMBAC, CPBI1MSL and CPBI1MSH and not containing ALLNOTRD.
807   /// In addition, it knows that the class matching CPBI1MSL requires an Lpt trigger
808   /// and the one with CPBI1MSH requires a Hpt trigger.
809   /// Hence, in the analysis, the function
810   /// TrackPtCutMatchTrigClass(track, "CPBIMSL") returns true if track match Lpt
811   /// TrackPtCutMatchTrigClass(track, "CPBIMSL") returns true if track match Hpt
812   /// TrackPtCutMatchTrigClass(track, "CMBAC") always returns true
813   ///
814   /// CAVEAT: if you use an fCFContainer and you want an axis to contain the trigger classes,
815   ///         please be sure that each pattern matches only 1 trigger class, or triggers will be mixed up
816   ///         when merging different chuncks.
817
818   fSelectedTrigPattern->SetOwner();
819   if ( fSelectedTrigPattern->GetEntries() > 0 ) fSelectedTrigPattern->Delete();
820   fRejectedTrigPattern->SetOwner();
821   if ( fRejectedTrigPattern->GetEntries() > 0 ) fRejectedTrigPattern->Delete();
822   fSelectedTrigLevel->SetOwner();
823   if ( fSelectedTrigLevel->GetEntries() > 0 ) fSelectedTrigLevel->Delete();
824
825   pattern.ReplaceAll("  "," ");
826   pattern.ReplaceAll("! ","!");
827   pattern.ReplaceAll(" :",":");
828
829   TObjArray* fullList = pattern.Tokenize(" ");
830
831   for ( Int_t ipat=0; ipat<fullList->GetEntries(); ++ipat ) {
832     TString currPattern = fullList->At(ipat)->GetName();
833     if ( currPattern.Contains("!") ) {
834       currPattern.ReplaceAll("!","");
835       fRejectedTrigPattern->AddLast(new TObjString(currPattern));
836     }
837     else {
838       TObjArray* arr = currPattern.Tokenize(":");
839       fSelectedTrigPattern->AddLast(new TObjString(arr->At(0)->GetName()));
840       TString selTrigLevel = ( arr->At(1) ) ? arr->At(1)->GetName() : "none";
841       selTrigLevel.ToUpper();
842       fSelectedTrigLevel->AddLast(new TObjString(selTrigLevel));
843       delete arr;
844     }
845   }
846   
847   delete fullList;
848 }
849
850 //________________________________________________________________________
851 void AliVAnalysisMuon::SetCentralityClasses(Int_t nCentralityBins, Double_t* centralityBins)
852 {
853   //
854   /// Set centrality classes
855   //
856   Double_t* bins = centralityBins;
857   Int_t nbins = nCentralityBins;
858   
859   Double_t defaultCentralityBins[] = {-5., 0., 5., 10., 20., 30., 40., 50., 60., 70., 80., 100., 105.};
860   if ( ! centralityBins ) {
861     bins = defaultCentralityBins;
862     nbins = sizeof(defaultCentralityBins)/sizeof(defaultCentralityBins[0])-1;
863   }
864
865   if ( fCentralityClasses ) delete fCentralityClasses;
866   fCentralityClasses = new TAxis(nbins, bins);
867   TString currClass = "";
868   for ( Int_t ibin=1; ibin<=fCentralityClasses->GetNbins(); ++ibin ){
869     currClass = Form("%.0f_%.0f",fCentralityClasses->GetBinLowEdge(ibin),fCentralityClasses->GetBinUpEdge(ibin));
870     fCentralityClasses->SetBinLabel(ibin, currClass.Data());
871   }
872 }
873
874 //________________________________________________________________________
875 void AliVAnalysisMuon::SetTerminateOptions(TString physSel, TString trigClass, TString centralityRange, TString furtherOpts)
876 {
877   //
878   /// Set terminate options
879   //
880   if ( ! fTerminateOptions ) {
881     fTerminateOptions = new TObjArray(4);
882     fTerminateOptions->SetOwner();
883   }
884   fTerminateOptions->AddAt(new TObjString(physSel), 0);
885   fTerminateOptions->AddAt(new TObjString(trigClass), 1);
886   fTerminateOptions->AddAt(new TObjString(centralityRange),2);
887   fTerminateOptions->AddLast(new TObjString(furtherOpts));
888 }
889
890 //________________________________________________________________________
891 void AliVAnalysisMuon::InitKeys()
892 {
893   TString chargeKeys = "MuMinus MuPlus";
894   fChargeKeys = chargeKeys.Tokenize(" ");
895   
896   TString srcKeys = "CharmMu BeautyMu QuarkoniumMu WbosonMu DecayMu SecondaryMu Hadron Unidentified";
897   fSrcKeys = srcKeys.Tokenize(" ");
898   
899   TString physSelKeys = "PhysSelPass PhysSelReject";
900   fPhysSelKeys = physSelKeys.Tokenize(" ");
901 }
902
903 //________________________________________________________________________
904 TObjArray* AliVAnalysisMuon::BuildTriggerClasses(TString firedTrigClasses)
905 {
906   //
907   /// Return the list of trigger classes to be considered
908   /// for current event. Update the global list if necessary
909   //
910
911   TObjArray* selectedTrigClasses = new TObjArray(0);
912   selectedTrigClasses->SetOwner();
913   
914   TObjArray* firedTrigClassesList = firedTrigClasses.Tokenize(" ");
915
916   for ( Int_t itrig=0; itrig<firedTrigClassesList->GetEntries(); ++itrig ) {
917     TString trigName = ((TObjString*)firedTrigClassesList->At(itrig))->GetString();
918     Bool_t rejectTrig = kFALSE;
919     for ( Int_t ipat=0; ipat<fRejectedTrigPattern->GetEntries(); ++ipat ) {
920       if ( trigName.Contains(fRejectedTrigPattern->At(ipat)->GetName() ) ) {
921         rejectTrig = kTRUE;
922         break;
923       }
924     } // loop on reject pattern
925     if ( rejectTrig ) continue;
926
927     Int_t matchPatternIndex = -1;
928     for ( Int_t ipat=0; ipat<fSelectedTrigPattern->GetEntries(); ++ipat ) {
929       if ( trigName.Contains(fSelectedTrigPattern->At(ipat)->GetName() ) ) {
930         matchPatternIndex = ipat;
931         break;
932       }
933     } // loop on keep pattern
934     if ( matchPatternIndex < 0 ) continue;
935
936     selectedTrigClasses->AddLast(new TObjString(trigName));
937     if ( fTriggerClasses->FindObject(trigName.Data()) ) continue;
938     Int_t trigLevel = 0;
939     TString trigLevelString = fSelectedTrigLevel->At(matchPatternIndex)->GetName();
940     if ( trigLevelString.Contains("APT") ) trigLevel = 1;
941     else if ( trigLevelString.Contains("LPT") ) trigLevel = 2;
942     else if ( trigLevelString.Contains("HPT") ) trigLevel = 3;
943     AliInfo(Form("Adding %s to considered trigger classes",trigName.Data()));
944     TObjString* addTrig = new TObjString(trigName);
945     UInt_t uniqueId = trigLevel;
946     addTrig->SetUniqueID(uniqueId);
947     fTriggerClasses->Add(addTrig);
948   } // loop on trigger classes
949
950   delete firedTrigClassesList;
951
952   return selectedTrigClasses;
953 }
954
955
956 //________________________________________________________________________
957 Bool_t AliVAnalysisMuon::TrackPtCutMatchTrigClass(AliVParticle* track, TString trigClassName)
958 {
959   /// Check if track passes the trigger pt cut level used in the trigger class
960   Int_t matchTrig = ( fAODEvent ) ? ((AliAODTrack*)track)->GetMatchTrigger() : ((AliESDMuonTrack*)track)->GetMatchTrigger();
961   Int_t classMatchLevel = GetTrigClassPtCutLevel(trigClassName);
962   return matchTrig >= classMatchLevel;
963 }
964
965
966 //________________________________________________________________________
967 Int_t AliVAnalysisMuon::GetTrigClassPtCutLevel(TString trigClassName)
968 {
969   /// Get trigger class pt cut level for tracking/trigger matching
970   TObject* obj = fTriggerClasses->FindObject(trigClassName.Data());
971   if ( ! obj ) {
972     AliWarning(Form("Class %s not in the list!", trigClassName.Data()));
973     return -1;
974   }
975   
976   return obj->GetUniqueID();
977 }