PWG/muon/AliAnalysisTaskMuMu.cxx

   1 #include "AliAnalysisTaskMuMu.h"
   2
   3 #include "AliAnalysisManager.h"
   4 #include "AliAnalysisMuMuBinning.h"
   5 #include "AliAnalysisMuonUtility.h"
   6 #include "AliAnalysisUtils.h"
   7 #include "AliAODEvent.h"
   8 #include "AliAODMCParticle.h"
   9 #include "AliAODTrack.h"
  10 #include "AliAODTZERO.h"
  11 #include "AliCentrality.h"
  12 #include "AliCodeTimer.h"
  13 #include "AliCounterCollection.h"
  14 #include "AliESDEvent.h"
  15 #include "AliESDTZERO.h"
  16 #include "AliInputEventHandler.h"
  17 #include "AliLog.h"
  18 #include "AliMCEvent.h"
  19 #include "AliMCEventHandler.h"
  20 #include "AliMergeableCollection.h"
  21 #include "AliMuonEventCuts.h"
  22 #include "AliMuonTrackCuts.h"
  23 #include "Riostream.h"
  24 #include "TCanvas.h"
  25 #include "TDatabasePDG.h"
  26 #include "TFormula.h"
  27 #include "TH1.h"
  28 #include "TH2.h"
  29 #include "THashList.h"
  30 #include "TList.h"
  31 #include "TMath.h"
  32 #include "TObjString.h"
  33 #include "TPaveText.h"
  34 #include "TProfile.h"
  35 #include "TRegexp.h"
  36 #include "TROOT.h"
  37 #include <algorithm>
  38 #include <cassert>
  39 #include "AliAnalysisMuMuBase.h"
  40 #include "AliInputEventHandler.h"
  41 #include "AliAnalysisMuMuCutRegistry.h"
  42 #include "AliAnalysisMuMuCutElement.h"
  43 #include "AliAnalysisMuMuCutCombination.h"
  44 #include <set>
  45
  46 /**
  47  * \class AliAnalysisTaskMuMu
  48  *
  49  * This class steers the work of one or more sub-analysis deriving from AliAnalysisMuMuBase
  50  * The output contains an AliHistogramCollection, an AliCounterCollection
  51  * and an AliAnalysisMuMuBinning
  52  * This task must be configured a bit before being used. For instance
  53  * you can select various event cuts, single muon track cuts and
  54  *  muon pairs cut, as well as defining various bins (for minv and mean pt
  55  * histograms) in pt,y,phi etc...
  56  *
  57  * Note that it's also possible to disable some (or all) histograms
  58  * (to save speed/memory), using DisableHistograms() method.
  59  *
  60  * For an example of such configuration, \see AddTaskMuMu.C
  61  */
  62
  63 using std::cout;
  64 using std::endl;
  65
  66 ClassImp(AliAnalysisTaskMuMu)
  67
  68 //_____________________________________________________________________________
  69 AliAnalysisTaskMuMu::AliAnalysisTaskMuMu()
  70 : AliAnalysisTaskSE("AliAnalysisTaskMuMu"),
  71 fHistogramCollection(0),
  72 fEventCounters(0),
  73 fBinning(0x0),
  74 fCutRegistry(0x0),
  75 fBeamYear(""),
  76 fHistogramToDisable(0x0),
  77 fSubAnalysisVector(0x0)
  78 {
  79   /// Constructor with a predefined list of triggers to consider
  80   /// Note that we take ownership of cutRegister
  81   ///
  82
  83 //  fBranchNames = "AOD:header,tracks,vertices,tracklets,AliAODTZERO,AliAODVZERO";
  84
  85   DefineOutput(1,AliMergeableCollection::Class());
  86   DefineOutput(2,AliCounterCollection::Class());
  87   DefineOutput(3,AliAnalysisMuMuBinning::Class());
  88 }
  89
  90 //_____________________________________________________________________________
  91 AliAnalysisTaskMuMu::~AliAnalysisTaskMuMu()
  92 {
  93   /// dtor
  94
  95   if (fHistogramCollection && ! AliAnalysisManager::GetAnalysisManager()->IsProofMode())
  96   {
  97     delete fHistogramCollection;
  98   }
  99
 100   if (fEventCounters && ! AliAnalysisManager::GetAnalysisManager()->IsProofMode())
 101   {
 102     delete fEventCounters;
 103   }
 104
 105   if (fBinning && ! AliAnalysisManager::GetAnalysisManager()->IsProofMode())
 106   {
 107     delete fBinning;
 108   }
 109
 110   delete fHistogramToDisable;
 111
 112   delete fCutRegistry;
 113
 114   delete fSubAnalysisVector;
 115 }
 116
 117 //_____________________________________________________________________________
 118 void AliAnalysisTaskMuMu::AdoptSubAnalysis(AliAnalysisMuMuBase* analysis)
 119 {
 120   if (!fSubAnalysisVector)
 121   {
 122     fSubAnalysisVector = new TObjArray;
 123     fSubAnalysisVector->SetOwner(kTRUE);
 124   }
 125   if ( !fSubAnalysisVector->FindObject(analysis) )
 126   {
 127     fSubAnalysisVector->Add(analysis);
 128   }
 129 }
 130
 131 //_____________________________________________________________________________
 132 AliAnalysisMuMuCutRegistry* AliAnalysisTaskMuMu::CutRegistry() const
 133 {
 134     /// Return (and create if not yet there) our cut registry
 135   if (!fCutRegistry)
 136   {
 137     fCutRegistry = new AliAnalysisMuMuCutRegistry;
 138   }
 139   return fCutRegistry;
 140 }
 141
 142 //_____________________________________________________________________________
 143 const char*
 144 AliAnalysisTaskMuMu::DefaultCentralityName() const
 145 {
 146   /// Get default centrality name
 147   if ( !fBeamYear.Contains("pp") ) return "CENTX";
 148   else return "PP";
 149 }
 150
 151 //_____________________________________________________________________________
 152 void AliAnalysisTaskMuMu::DisableHistograms(const char* pattern)
 153 {
 154   /// Disable the histogramming of all the histograms matching the pattern
 155
 156   TIter next(fSubAnalysisVector);
 157   AliAnalysisMuMuBase* a;
 158
 159   while ( ( a = static_cast<AliAnalysisMuMuBase*>(next()) ) )
 160   {
 161     a->DisableHistograms(pattern);
 162   }
 163 }
 164
 165 //_____________________________________________________________________________
 166 AliVEvent*
 167 AliAnalysisTaskMuMu::Event() const
 168 {
 169   // some const-dirty-dancing
 170   return const_cast<AliAnalysisTaskMuMu*>(this)->InputEvent();
 171 }
 172
 173 //_____________________________________________________________________________
 174 void AliAnalysisTaskMuMu::Fill(const char* eventSelection, const char* triggerClassName)
 175 {
 176   // Fill one set of histograms (only called for events which pass the eventSelection cut)
 177
 178   TString seventSelection(eventSelection);
 179   seventSelection.ToLower();
 180
 181   fEventCounters->Count(Form("event:%s/trigger:%s/centrality:%s/run:%d", seventSelection.Data(), triggerClassName, "ALL", fCurrentRunNumber));
 182
 183   if ( !IsHistogrammingDisabled() )
 184   {
 185     TObjArray* centralities = fBinning->CreateBinObjArray("centrality");
 186
 187     TIter next(centralities);
 188     AliAnalysisMuMuBinning::Range* r;
 189
 190     while ( ( r = static_cast<AliAnalysisMuMuBinning::Range*>(next()) ) )
 191     {
 192       TString estimator = r->Quantity();
 193
 194       Float_t fcent = Event()->GetCentrality()->GetCentralityPercentile(estimator.Data());
 195       if ( fcent < 0.) FillHistos(eventSelection,triggerClassName,"MV0");
 196       if ( fcent == 0.) FillHistos(eventSelection,triggerClassName,"0V0");
 197       if ( r->IsInRange(fcent) )
 198       {
 199         FillHistos(eventSelection,triggerClassName,r->AsString());
 200       }
 201     }
 202     delete centralities;
 203   }
 204 }
 205
 206 //_____________________________________________________________________________
 207 void AliAnalysisTaskMuMu::FillHistos(const char* eventSelection,
 208                                      const char* triggerClassName,
 209                                      const char* centrality)
 210 {
 211   /// Fill histograms for /physics/triggerClassName/centrality
 212
 213   AliCodeTimerAuto("",0);
 214
 215   TIter nextAnalysis(fSubAnalysisVector);
 216   AliAnalysisMuMuBase* analysis;
 217
 218   Int_t nTracks = AliAnalysisMuonUtility::GetNTracks(Event());
 219
 220   fEventCounters->Count(Form("event:%s/trigger:%s/centrality:%s/run:%d", eventSelection, triggerClassName, centrality, fCurrentRunNumber));
 221
 222   TIter nextTrackCut(fCutRegistry->GetCutCombinations(AliAnalysisMuMuCutElement::kTrack));
 223   TIter nextPairCut(fCutRegistry->GetCutCombinations(AliAnalysisMuMuCutElement::kTrackPair));
 224
 225   // loop on single tracks (whatever the type of tracks
 226   while ( ( analysis = static_cast<AliAnalysisMuMuBase*>(nextAnalysis()) ) )
 227   {
 228     analysis->DefineHistogramCollection(eventSelection,triggerClassName,centrality);
 229
 230     AliCodeTimerAuto(Form("%s (FillHistosForEvent)",analysis->ClassName()),1);
 231
 232     if ( MCEvent() != 0x0 )
 233     {
 234       analysis->FillHistosForMCEvent(eventSelection,triggerClassName,centrality);
 235     }
 236
 237     analysis->FillHistosForEvent(eventSelection,triggerClassName,centrality);
 238
 239     for (Int_t i = 0; i < nTracks; ++i)
 240     {
 241       AliVParticle* tracki = AliAnalysisMuonUtility::GetTrack(i,Event());
 242
 243       nextTrackCut.Reset();
 244       AliAnalysisMuMuCutCombination* trackCut;
 245
 246       while ( ( trackCut = static_cast<AliAnalysisMuMuCutCombination*>(nextTrackCut()) ) )
 247       {
 248         if ( trackCut->Pass(*tracki) )
 249         {
 250           analysis->FillHistosForTrack(eventSelection,triggerClassName,centrality,trackCut->GetName(),*tracki);
 251         }
 252       }
 253
 254       if (!AliAnalysisMuonUtility::IsMuonTrack(tracki) ) continue;
 255
 256       // loop on track pairs (here we only consider muon pairs)
 257
 258       for (Int_t j = i+1; j < nTracks; ++j)
 259       {
 260         AliVParticle* trackj = AliAnalysisMuonUtility::GetTrack(j,Event());
 261
 262         if (!AliAnalysisMuonUtility::IsMuonTrack(trackj) ) continue;
 263
 264         nextPairCut.Reset();
 265         AliAnalysisMuMuCutCombination* pairCut;
 266
 267         while ( ( pairCut = static_cast<AliAnalysisMuMuCutCombination*>(nextPairCut()) ) )
 268         {
 269           Bool_t testi = (pairCut->IsTrackCutter()) ? pairCut->Pass(*tracki) : kTRUE;
 270           Bool_t testj = (pairCut->IsTrackCutter()) ? pairCut->Pass(*trackj) : kTRUE;
 271           Bool_t testij = pairCut->Pass(*tracki,*trackj);
 272
 273           if ( ( testi || testj ) && testij )
 274           {
 275             analysis->FillHistosForPair(eventSelection,triggerClassName,centrality,pairCut->GetName(),*tracki,*trackj);
 276           }
 277         }
 278       }
 279     }
 280   }
 281 }
 282
 283 //_____________________________________________________________________________
 284 void AliAnalysisTaskMuMu::FinishTaskOutput()
 285 {
 286   /// prune empty histograms BEFORE mergin, in order to save some bytes...
 287
 288   if ( fHistogramCollection )
 289   {
 290     fHistogramCollection->PruneEmptyObjects();
 291   }
 292 }
 293
 294 //_____________________________________________________________________________
 295 void AliAnalysisTaskMuMu::GetSelectedTrigClassesInEvent(const AliVEvent* event, TObjArray& array)
 296 {
 297   /// Fills the array with a list of TObjString of the trigger classes that the various
 298   /// cuts accept for this event
 299
 300   array.Clear();
 301
 302   if (!event)
 303   {
 304     AliError("Will get a hard time selecting trigger classes with an empty event...");
 305     return;
 306   }
 307
 308   TString firedTriggerClasses = event->GetFiredTriggerClasses();
 309   UInt_t l0 = AliAnalysisMuonUtility::GetL0TriggerInputs(event);
 310   UInt_t l1 = AliAnalysisMuonUtility::GetL1TriggerInputs(event);
 311   UInt_t l2 = AliAnalysisMuonUtility::GetL2TriggerInputs(event);
 312
 313   std::set<std::string> tmpArray;
 314
 315   TIter nextCutCombination(CutRegistry()->GetCutCombinations(AliAnalysisMuMuCutElement::kTriggerClass));
 316   AliAnalysisMuMuCutCombination* cutCombination;
 317
 318   while ( ( cutCombination = static_cast<AliAnalysisMuMuCutCombination*>(nextCutCombination()) ) )
 319   {
 320     TString acceptedTriggerClasses;
 321
 322     if ( cutCombination->Pass(firedTriggerClasses,acceptedTriggerClasses,l0,l1,l2) )
 323     {
 324       TObjArray* split = acceptedTriggerClasses.Tokenize(" ");
 325       TIter next(split);
 326       TObjString* str;
 327       while ( ( str = static_cast<TObjString*>(next()) ) )
 328       {
 329         tmpArray.insert(str->String().Data());
 330       }
 331       delete split;
 332     }
 333   }
 334
 335   std::set<std::string>::const_iterator it;
 336
 337   for ( it = tmpArray.begin(); it != tmpArray.end(); ++it )
 338   {
 339     array.Add(new TObjString(it->c_str()));
 340   }
 341 }
 342
 343
 344 //_____________________________________________________________________________
 345 Bool_t AliAnalysisTaskMuMu::IsHistogramDisabled(const char* hname) const
 346 {
 347   /// Whether or not a given histogram (identified by its name)
 348   /// is disabled or not
 349
 350   TIter next(fSubAnalysisVector);
 351   AliAnalysisMuMuBase* analysis;
 352
 353   while ( ( analysis = static_cast<AliAnalysisMuMuBase*>(next()) ) )
 354   {
 355     if ( analysis->IsHistogramDisabled(hname) )
 356     {
 357       return kTRUE;
 358     }
 359   }
 360
 361   return kFALSE;
 362 }
 363
 364 //_____________________________________________________________________________
 365 Bool_t AliAnalysisTaskMuMu::IsHistogrammingDisabled() const
 366 {
 367   /// Whether or not *all* histograms are disabled
 368
 369   Bool_t disabled(kTRUE);
 370
 371   TIter next(fSubAnalysisVector);
 372   AliAnalysisMuMuBase* analysis;
 373
 374   while ( ( analysis = static_cast<AliAnalysisMuMuBase*>(next()) ) )
 375   {
 376     disabled = disabled && analysis->IsHistogrammingDisabled();
 377   }
 378
 379   return disabled;
 380 }
 381
 382 //_____________________________________________________________________________
 383 Bool_t AliAnalysisTaskMuMu::IsPP() const
 384 {
 385   // whether we're dealing with proton proton collisions
 386   return fBeamYear.Contains("pp");
 387 }
 388
 389 //_____________________________________________________________________________
 390 void AliAnalysisTaskMuMu::NotifyRun()
 391 {
 392   /// Called at each change of run
 393
 394   AliDebug(1,Form("Run %09d File %s",fCurrentRunNumber,CurrentFileName()));
 395
 396   TIter next(fSubAnalysisVector);
 397   AliAnalysisMuMuBase* analysis;
 398
 399   while ( ( analysis = static_cast<AliAnalysisMuMuBase*>(next()) ) )
 400   {
 401     analysis->SetRun(fInputHandler);
 402   }
 403 }
 404
 405 //_____________________________________________________________________________
 406 void
 407 AliAnalysisTaskMuMu::Print(Option_t* opt) const
 408 {
 409   /// Print the definition of this analysis
 410
 411   cout << ClassName() << " - " << GetName() << " - " << fBeamYear.Data() << endl;
 412
 413   TIter next(fSubAnalysisVector);
 414   AliAnalysisMuMuBase* analysis;
 415
 416   while ( ( analysis = static_cast<AliAnalysisMuMuBase*>(next()) ) )
 417   {
 418     analysis->Print(opt);
 419   }
 420
 421   fCutRegistry->Print("ALL");
 422
 423   if ( fBinning )
 424   {
 425     cout << "Binning" << endl;
 426     fBinning->Print();
 427   }
 428 }
 429
 430 //_____________________________________________________________________________
 431 void
 432 AliAnalysisTaskMuMu::Terminate(Option_t* opt)
 433 {
 434   /// Called once at the end of the query
 435   /// Just a simple printout of the stat we analyse and how many histograms
 436   /// we got
 437
 438   TIter next(fSubAnalysisVector);
 439   AliAnalysisMuMuBase* analysis;
 440
 441   while ( ( analysis = static_cast<AliAnalysisMuMuBase*>(next()) ) )
 442   {
 443     analysis->Terminate(opt);
 444   }
 445
 446   fHistogramCollection = dynamic_cast<AliMergeableCollection*>(GetOutputData(1));
 447
 448   if (!fHistogramCollection)
 449   {
 450     AliError("Could not find back histogram collection in output...");
 451   }
 452   else
 453   {
 454     // Removes empty objects and also the event histos of the Nch task
 455     fHistogramCollection->PruneEmptyObjects();
 456
 457     UInt_t size2 = fHistogramCollection->EstimateSize();
 458
 459     TIter nextHistogram(fHistogramCollection->CreateIterator());
 460     TObject* object;
 461
 462     while ( ( object = nextHistogram() ) )
 463     {
 464       if ( object->IsA()->InheritsFrom(TH1::Class()) )
 465       {
 466         TH1* h = static_cast<TH1*>(object);
 467         if ( h->GetXaxis()->GetLabels() )
 468         {
 469           h->LabelsDeflate("X");
 470         }
 471       }
 472     }
 473
 474     AliInfo(Form("size after prune histograms = %5.1f MB",size2/1024.0/1024.0));
 475
 476     fHistogramCollection->Print("-");
 477   }
 478
 479   fEventCounters = dynamic_cast<AliCounterCollection*>(GetOutputData(2));
 480
 481   if (!fEventCounters)
 482   {
 483     AliError("Could not find back counters in output...");
 484   }
 485   else
 486   {
 487     fEventCounters->Print("trigger/event");
 488   }
 489
 490   // post param container(s)
 491   PostData(3,fBinning);
 492 }
 493
 494 //_____________________________________________________________________________
 495 AliAnalysisMuMuBinning* AliAnalysisTaskMuMu::Binning() const
 496 {
 497   // Return our binning (making a default one if not already created
 498   if ( fBinning ) return fBinning;
 499
 500   fBinning = new AliAnalysisMuMuBinning("BIN");
 501
 502   return fBinning;
 503 }
 504
 505 //_____________________________________________________________________________
 506 void AliAnalysisTaskMuMu::UserExec(Option_t* /*opt*/)
 507 {
 508   /// Executed at each event
 509
 510 //  static Int_t n(0);
 511 //  AliInfo(Form("EVENT %10d Event()=%p MCEvent()=%p",n,Event(),MCEvent()));
 512 //  ++n;
 513 //
 514   AliCodeTimerAuto("",0);
 515
 516   Binning(); // insure we have a binning...
 517
 518   //  if ( MCEvent() )
 519   //  {
 520   TIter nextAnalysis(fSubAnalysisVector);
 521   AliAnalysisMuMuBase* analysis;
 522   while ( ( analysis = static_cast<AliAnalysisMuMuBase*>(nextAnalysis()) ) )
 523   {
 524     if ( MCEvent() ) // Set the MC flag for all analysis (prior to call anything from them
 525       // (e.g. any trigger class selection that might behave differently for
 526       // MC and real trigger classes)
 527     {
 528       analysis->SetMC();
 529     }
 530     analysis->SetEvent(Event(),MCEvent()); // Set the new event properties derived in the analysis
 531   }
 532   //  }
 533
 534
 535   TString firedTriggerClasses(AliAnalysisMuonUtility::GetFiredTriggerClasses(Event()));
 536
 537   // first loop to count things not associated to a specific trigger
 538   TIter nextEventCutCombination(CutRegistry()->GetCutCombinations(AliAnalysisMuMuCutElement::kEvent));
 539   AliAnalysisMuMuCutCombination* cutCombination;
 540
 541   while ( ( cutCombination = static_cast<AliAnalysisMuMuCutCombination*>(nextEventCutCombination())))
 542   {
 543     if ( cutCombination->Pass(*fInputHandler) )
 544     {
 545       fEventCounters->Count(Form("event:%s/trigger:%s/centrality:%s/run:%d", cutCombination->GetName(), "EVERYTHING",  "ALL", fCurrentRunNumber));
 546
 547       if ( firedTriggerClasses == "" )
 548       {
 549         fEventCounters->Count(Form("event:%s/trigger:%s/centrality:%s/run:%d", cutCombination->GetName(), "EMPTY", "ALL", fCurrentRunNumber));
 550       }
 551     }
 552   }
 553
 554   // second loop to count only the triggers we're interested in
 555   TObjArray selectedTriggerClasses;
 556
 557   GetSelectedTrigClassesInEvent(Event(),selectedTriggerClasses);
 558
 559   TIter next(&selectedTriggerClasses);
 560   TObjString* tname;
 561 //  Bool_t hasSetEventBeenCalled(kFALSE);
 562
 563   while ( ( tname = static_cast<TObjString*>(next()) ) )
 564   {
 565     nextEventCutCombination.Reset();
 566
 567     while ( ( cutCombination = static_cast<AliAnalysisMuMuCutCombination*>(nextEventCutCombination())) )
 568     {
 569       if ( cutCombination->Pass(*fInputHandler) )
 570       {
 571 //        if (!hasSetEventBeenCalled)
 572 //        {
 573 //          TIter nextAnalysis(fSubAnalysisVector);
 574 //          AliAnalysisMuMuBase* analysis;
 575 //
 576 //          while ( ( analysis = static_cast<AliAnalysisMuMuBase*>(nextAnalysis()) ) )
 577 //          {
 578 //            analysis->SetEvent(Event(),MCEvent());
 579 //          }
 580 //          hasSetEventBeenCalled = kTRUE;
 581 //        }
 582         Fill(cutCombination->GetName(),tname->String().Data());
 583       }
 584     }
 585   }
 586
 587   // Post output data.
 588   PostData(1, fHistogramCollection);
 589   PostData(2, fEventCounters);
 590   PostData(3, fBinning);
 591 }
 592
 593 //_____________________________________________________________________________
 594 void AliAnalysisTaskMuMu::UserCreateOutputObjects()
 595 {
 596   /// Create histograms
 597   /// Called once
 598
 599   OpenFile(1);
 600
 601   AliInfo(Form("fCutRegistry=%p",fCutRegistry));
 602
 603   if ( fCutRegistry )
 604   {
 605     fCutRegistry->Print();
 606   }
 607
 608   fHistogramCollection = new AliMergeableCollection("OC");
 609
 610   fEventCounters = new AliCounterCollection("CC");
 611
 612   // initialize event counters
 613
 614   TString eventRubric;
 615   TIter next(CutRegistry()->GetCutCombinations(AliAnalysisMuMuCutElement::kEvent));
 616   AliAnalysisMuMuCutCombination* cutCombination;
 617
 618   while ( ( cutCombination = static_cast<AliAnalysisMuMuCutCombination*>(next())) )
 619   {
 620     TString cutName = cutCombination->GetName();
 621     if ( eventRubric.Length() > 0 ) eventRubric += "/";
 622     eventRubric += cutName;
 623   }
 624
 625   fEventCounters->AddRubric("event", eventRubric.Data());
 626
 627   fEventCounters->AddRubric("trigger", 100);
 628
 629   fEventCounters->AddRubric("centrality", 100);
 630
 631   fEventCounters->AddRubric("run", 1000000);
 632
 633   // Initialize our subtasks, if any...
 634
 635   TIter nextAnalysis(fSubAnalysisVector);
 636   AliAnalysisMuMuBase* analysis;
 637
 638   while ( ( analysis = static_cast<AliAnalysisMuMuBase*>(nextAnalysis()) ) )
 639   {
 640     analysis->Init(*fEventCounters,*fHistogramCollection,*fBinning,*fCutRegistry);
 641   }
 642
 643   // finally end the counters initialization
 644   fEventCounters->Init();
 645
 646   // Post output data.
 647   PostData(1,fHistogramCollection);
 648   PostData(2,fEventCounters);
 649   PostData(3,fBinning);
 650 }