PWG/muon/AliAnalysisMuMu.cxx

   1 #include "AliAnalysisMuMu.h"
   2
   3 #include "AliAnalysisManager.h"
   4 #include "AliCounterCollection.h"
   5 #include "AliHistogramCollection.h"
   6 #include "AliLog.h"
   7 #include "AliLog.h"
   8 #include "Riostream.h"
   9 #include "TCanvas.h"
  10 #include "TDatabasePDG.h"
  11 #include "TH1.h"
  12 #include "TH2.h"
  13 #include "THashList.h"
  14 #include "TList.h"
  15 #include "TMath.h"
  16 #include "TObjString.h"
  17 #include "TPaveText.h"
  18 #include "TROOT.h"
  19 #include <algorithm>
  20
  21 //
  22 // AliAnalysisMuMu : base class for mu pairs analysis
  23 //
  24 // Mainly invariant mass (for J/psi and Upsilon) but also
  25 // some single control histograms.
  26 //
  27 // This base class contains common things for ESD-based
  28 // and AOD-based analysis
  29 //
  30 // The output contains an AliHistogramCollection and
  31 // an AliCounterCollection
  32 //
  33 // author: L. Aphecetche (Subatech)
  34 //
  35
  36 ClassImp(AliAnalysisMuMu)
  37
  38 namespace
  39 {
  40   Int_t GetNbins(Double_t xmin, Double_t xmax, Double_t xstep)
  41   {
  42     if ( TMath::AreEqualRel(xstep,0.0,1E-9) ) return 1;
  43
  44     return TMath::Nint(TMath::Abs((xmax-xmin)/xstep));
  45   }
  46 }
  47
  48 //_____________________________________________________________________________
  49 AliAnalysisMuMu* AliAnalysisMuMu::Create(const char* inputDataType, TList* triggerClassesToConsider)
  50 {
  51   /// Create the right implementation depending on intputDataType
  52
  53   TString sinputDataType(inputDataType);
  54   sinputDataType.ToUpper();
  55
  56   if ( sinputDataType != "AOD" && sinputDataType != "ESD" )
  57   {
  58     AliErrorClass(Form("Unknown input data type : %s",inputDataType));
  59     return 0x0;
  60   }
  61
  62   return reinterpret_cast<AliAnalysisMuMu*>(gROOT->ProcessLineFast(Form("new AliAnalysisMuMuFrom%s((TList*)%p)",
  63                                                                         sinputDataType.Data(),
  64                                                                         triggerClassesToConsider)));
  65 }
  66
  67 //_____________________________________________________________________________
  68 AliAnalysisMuMu* AliAnalysisMuMu::Create(const char* inputDataType, Bool_t aa)
  69 {
  70   /// Create the right implementation depending on intputDataType
  71
  72   TString sinputDataType(inputDataType);
  73   sinputDataType.ToUpper();
  74
  75   if ( sinputDataType != "AOD" && sinputDataType != "ESD" )
  76   {
  77     AliErrorClass(Form("Unknown input data type : %s",inputDataType));
  78     return 0x0;
  79   }
  80
  81   return reinterpret_cast<AliAnalysisMuMu*>(gROOT->ProcessLineFast(Form("new AliAnalysisMuMuFrom%s(%d)",
  82                                                                         sinputDataType.Data(),aa)));
  83 }
  84
  85 //_____________________________________________________________________________
  86 AliAnalysisMuMu::AliAnalysisMuMu() : AliAnalysisTaskSE("AliAnalysisMuMu"),
  87 fHistogramCollection(0),
  88 fOutput(0),
  89 fTriggerClasses(0),
  90 fEventCounters(0),
  91 fIsFromESD(kFALSE),
  92 fSingleTrackCutNames(0x0),
  93 fPairTrackCutNames(0x0),
  94 fCentralityLimits(),
  95 fCentralityNames(0x0),
  96 fGlobalEventSelectionNames(0x0),
  97 fAA(kFALSE),
  98 fIsDynamicTriggerClasses(kFALSE)
  99 {
 100   /// default ctor
 101 }
 102
 103 //_____________________________________________________________________________
 104 AliAnalysisMuMu::AliAnalysisMuMu(Bool_t fromESD, Bool_t aa)
 105 : AliAnalysisTaskSE(Form("AliAnalysisMuMu-from%s",fromESD ? "ESD":"AOD")),
 106 fHistogramCollection(0),
 107 fOutput(0),
 108 fTriggerClasses(new THashList),
 109 fEventCounters(0),
 110 fIsFromESD(fromESD),
 111 fSingleTrackCutNames(0x0),
 112 fPairTrackCutNames(0x0),
 113 fCentralityLimits(),
 114 fCentralityNames(new TObjArray),
 115 fGlobalEventSelectionNames(0x0),
 116 fAA(aa),
 117 fIsDynamicTriggerClasses(kTRUE)
 118 {
 119   /// Constructor
 120   /// The list of triggers to be considered will be updated on the fly
 121   /// (see method AddTriggerClasses)
 122
 123   fTriggerClasses->SetOwner(kTRUE);
 124
 125   DefineOutput(1, TList::Class());
 126
 127   DefineCentralityClasses();
 128
 129 }
 130
 131 //_____________________________________________________________________________
 132 AliAnalysisMuMu::AliAnalysisMuMu(Bool_t fromESD, TList* triggerClasses)
 133 : AliAnalysisTaskSE(Form("AliAnalysisMuMu-from%s",fromESD ? "ESD":"AOD")),
 134 fHistogramCollection(0),
 135 fOutput(0),
 136 fTriggerClasses(new THashList),
 137 fEventCounters(0),
 138 fIsFromESD(fromESD),
 139 fSingleTrackCutNames(0x0),
 140 fPairTrackCutNames(0x0),
 141 fCentralityLimits(),
 142 fCentralityNames(new TObjArray),
 143 fGlobalEventSelectionNames(0x0),
 144 fAA(kFALSE),
 145 fIsDynamicTriggerClasses(kFALSE)
 146 {
 147   /// Constructor with a predefined list of triggers to consider
 148
 149   fTriggerClasses->SetOwner(kTRUE);
 150
 151   DefineOutput(1, TList::Class());
 152
 153   TObjString* tname;
 154   TIter next(triggerClasses);
 155
 156   while ( ( tname = static_cast<TObjString*>(next()) ) )
 157   {
 158     fTriggerClasses->Add(new TObjString(*tname));
 159     if ( tname->String().BeginsWith("CMB") )
 160     {
 161       fAA = kTRUE;
 162     }
 163   }
 164
 165   DefineCentralityClasses();
 166 }
 167
 168 //_____________________________________________________________________________
 169 AliAnalysisMuMu::~AliAnalysisMuMu()
 170 {
 171   /// dtor
 172   delete fTriggerClasses;
 173   delete fSingleTrackCutNames;
 174   delete fPairTrackCutNames;
 175   if(fOutput && ! AliAnalysisManager::GetAnalysisManager()->IsProofMode())
 176   {
 177     delete fOutput;
 178   }
 179   delete fCentralityNames;
 180   delete fGlobalEventSelectionNames;
 181 }
 182
 183 //_____________________________________________________________________________
 184 void AliAnalysisMuMu::AddGlobalEventSelection(const char* name)
 185 {
 186   /// Add a name to the list of global event event selection ones.
 187
 188   if (!fGlobalEventSelectionNames)
 189   {
 190     fGlobalEventSelectionNames = new TObjArray;
 191     fGlobalEventSelectionNames->SetOwner(kTRUE);
 192   }
 193   fGlobalEventSelectionNames->Add(new TObjString(name));
 194 }
 195
 196 //_____________________________________________________________________________
 197 UInt_t AliAnalysisMuMu::CodePairCutMask(UInt_t maskForOneOrBothTrack, UInt_t maskForTrackPair) const
 198 {
 199   /// Encode two masks (obviously each one should be 16 bits only) into one
 200
 201   return ( ( maskForOneOrBothTrack << 16 ) | maskForTrackPair );
 202 }
 203
 204 //_____________________________________________________________________________
 205 void AliAnalysisMuMu::DecodePairCutMask(UInt_t code, UInt_t& maskForOneOrBothTrack, UInt_t& maskForTrackPair) const
 206 {
 207   /// Decode the pair cut mask
 208   maskForOneOrBothTrack = ( ( code & 0xFFFF0000 ) >> 16 );
 209   maskForTrackPair = ( code & 0xFFFF );
 210 }
 211
 212 //_____________________________________________________________________________
 213 void AliAnalysisMuMu::DefineCentralityClasses()
 214 {
 215   /// Define the default centrality classes that will be used.
 216
 217   if ( fAA )
 218   {
 219     fCentralityLimits.push_back(10.0);
 220     fCentralityLimits.push_back(20.0);
 221     fCentralityLimits.push_back(30.0);
 222     fCentralityLimits.push_back(40.0);
 223     fCentralityLimits.push_back(50.0);
 224     fCentralityLimits.push_back(60.0);
 225     fCentralityLimits.push_back(70.0);
 226     fCentralityLimits.push_back(80.0);
 227     fCentralityLimits.push_back(90.0);
 228   }
 229 //  fCentralityLimits.push_back(100.0);
 230
 231   for ( std::vector<double>::size_type i = 0; i < fCentralityLimits.size(); ++i )
 232   {
 233     Double_t limit = fCentralityLimits[i];
 234     fCentralityNames->Add(new TObjString(CentralityName(limit)));
 235   }
 236
 237   fCentralityNames->Add(new TObjString(DefaultCentralityName()));
 238   fCentralityNames->SetOwner(kTRUE);
 239 }
 240
 241 //_____________________________________________________________________________
 242 void AliAnalysisMuMu::AddTriggerClasses(const char* triggerlist)
 243 {
 244   /// Given a list of trigger names (triggerlist) separated by spaces
 245   /// add those which we don't know yet (selecting only the few ones
 246   /// we're interested in, e.g CINT* CMU* CMB*
 247   ///
 248
 249   TString slist(triggerlist);
 250
 251   TObjArray* a = slist.Tokenize(" ");
 252   TObjString* s;
 253   TIter next(a);
 254
 255   while ( ( s = static_cast<TObjString*>(next()) ) )
 256   {
 257     TString trigger(s->String());
 258     Bool_t add(kFALSE);
 259
 260     if ( trigger.BeginsWith("CMB") && trigger.Contains("-B-") ) add = kTRUE;
 261     if ( trigger.BeginsWith("CINT") && trigger.Contains("-B-") && !trigger.Contains("WU") ) add = kTRUE;
 262     if ( trigger.BeginsWith("CMU") && trigger.Contains("-B-") ) add = kTRUE;
 263
 264     if ( trigger.BeginsWith("CMUP") ) add = kFALSE;
 265
 266     if ( add && !fTriggerClasses->FindObject(trigger.Data()) )
 267     {
 268       AliInfo(Form("Adding %s to considered trigger classes",trigger.Data()));
 269       fTriggerClasses->Add(new TObjString(trigger));
 270     }
 271
 272   }
 273
 274   delete a;
 275 }
 276
 277 //_____________________________________________________________________________
 278 void AliAnalysisMuMu::AddPairCut(const char* cutName, UInt_t maskForOneOrBothTrack, UInt_t maskForTrackPair)
 279 {
 280   /// Add a cut for a pair.
 281   /// maskForOneOrBothTrack is the mask of cuts that at least one track must satisfy
 282   /// maskForTrackPair is the mask of cuts that *both* tracks must satisfy.
 283   /// if maskForTrackPair is 0, then no (extra) condition is applied to the pair
 284
 285   if ( !fPairTrackCutNames )
 286   {
 287     fPairTrackCutNames = new TObjArray;
 288     fPairTrackCutNames->SetOwner(kTRUE);
 289   }
 290   TObjString* oname = new TObjString(Form("p%s",cutName));
 291   if ( !maskForTrackPair ) maskForTrackPair = maskForOneOrBothTrack;
 292   oname->SetUniqueID(CodePairCutMask(maskForOneOrBothTrack,maskForTrackPair));
 293   fPairTrackCutNames->Add(oname);
 294 }
 295
 296 //_____________________________________________________________________________
 297 void AliAnalysisMuMu::AddSingleCut(const char* name, UInt_t mask)
 298 {
 299   /// Add a cut for single tracks
 300   if ( !fSingleTrackCutNames )
 301   {
 302     fSingleTrackCutNames = new TObjArray;
 303     fSingleTrackCutNames->SetOwner(kTRUE);
 304   }
 305   TObjString* oname = new TObjString(Form("s%s",name));
 306   oname->SetUniqueID(mask);
 307   fSingleTrackCutNames->Add(oname);
 308 }
 309
 310 //_____________________________________________________________________________
 311 void AliAnalysisMuMu::BeautifyHistos()
 312 {
 313   /// Put the titles, marker sizes, color, etc...
 314
 315   TIter next(fHistogramCollection->CreateIterator());
 316   TH1* h;
 317
 318   while ( ( h = static_cast<TH1*>(next()) ) )
 319   {
 320     TString name(h->GetName());
 321
 322     if ( name.Contains("Plus") )
 323     {
 324       h->SetMarkerStyle(kFullCircle);
 325       h->SetMarkerColor(kBlue);
 326       h->SetLineColor(kBlue);
 327     }
 328     else
 329     {
 330       h->SetMarkerStyle(kOpenCircle);
 331       h->SetMarkerColor(kRed);
 332       h->SetLineColor(kRed);
 333     }
 334   }
 335 }
 336
 337 //_____________________________________________________________________________
 338 void
 339 AliAnalysisMuMu::CreateSingleHisto(const char* physics,
 340                                    const char* triggerClassName,
 341                                    const char* hname, const char* htitle,
 342                                    Int_t nbinsx, Double_t xmin, Double_t xmax,
 343                                    Int_t nbinsy, Double_t ymin, Double_t ymax) const
 344 {
 345   /// Append histograms for single track to our histogram collection
 346   CreateHisto(fSingleTrackCutNames,physics,triggerClassName,hname,htitle,
 347               nbinsx,xmin,xmax,nbinsy,ymin,ymax);
 348 }
 349
 350 //_____________________________________________________________________________
 351 void
 352 AliAnalysisMuMu::CreatePairHisto(const char* physics,
 353                                  const char* triggerClassName,
 354                                  const char* hname, const char* htitle,
 355                                  Int_t nbinsx, Double_t xmin, Double_t xmax,
 356                                  Int_t nbinsy, Double_t ymin, Double_t ymax) const
 357 {
 358   /// Append histograms for track pairs to our histogram collection
 359
 360   CreateHisto(fPairTrackCutNames,physics,triggerClassName,hname,htitle,
 361               nbinsx,xmin,xmax,nbinsy,ymin,ymax);
 362 }
 363
 364 //_____________________________________________________________________________
 365 void
 366 AliAnalysisMuMu::CreateEventHisto(const char* physics,
 367                                   const char* triggerClassName,
 368                                   const char* hname, const char* htitle,
 369                                   Int_t nbinsx, Double_t xmin, Double_t xmax,
 370                                   Int_t nbinsy, Double_t ymin, Double_t ymax) const
 371 {
 372   /// Append histograms at the event level
 373
 374   TIter next(fCentralityNames);
 375   TObjString* cent;
 376
 377   while ( ( cent = static_cast<TObjString*>(next()) ) )
 378   {
 379     TH1* h(0x0);
 380
 381     if ( nbinsy > 0 )
 382     {
 383       h = new TH2F(hname,htitle,nbinsx,xmin,xmax,nbinsy,ymin,ymax);
 384     }
 385     else
 386     {
 387       h = new TH1F(hname,htitle,nbinsx,xmin,xmax);
 388     }
 389
 390     fHistogramCollection->Adopt(physics,triggerClassName,cent->String().Data(),h);
 391   }
 392 }
 393
 394 //_____________________________________________________________________________
 395 void
 396 AliAnalysisMuMu::CreateHisto(TObjArray* array,
 397                              const char* physics,
 398                              const char* triggerClassName,
 399                              const char* hname, const char* htitle,
 400                              Int_t nbinsx, Double_t xmin, Double_t xmax,
 401                              Int_t nbinsy, Double_t ymin, Double_t ymax) const
 402 {
 403   /// Create a bunch of histograms for all centralities
 404   /// FIXME: have a way to specify the histo precision (i.e. F vs I vs S ...)
 405
 406   TIter next(array);
 407   TObjString* tcn;
 408   while ( ( tcn = static_cast<TObjString*>(next()) ) )
 409   {
 410     TIter nextCent(fCentralityNames);
 411     TObjString* cent;
 412
 413     while ( ( cent = static_cast<TObjString*>(nextCent()) ) )
 414     {
 415       TH1* h(0x0);
 416
 417       if ( nbinsy > 0 )
 418       {
 419         h = new TH2F(hname,htitle,nbinsx,xmin,xmax,nbinsy,ymin,ymax);
 420       }
 421       else
 422       {
 423         h = new TH1F(hname,htitle,nbinsx,xmin,xmax);
 424       }
 425
 426       fHistogramCollection->Adopt(physics,triggerClassName,cent->String().Data(),tcn->String().Data(),h);
 427     }
 428   }
 429 }
 430
 431 //_____________________________________________________________________________
 432 const char*
 433 AliAnalysisMuMu::DefaultCentralityName() const
 434 {
 435   /// Get default centrality name
 436   if ( fAA ) return "CENTX";
 437   else return "PP";
 438 }
 439
 440 //_____________________________________________________________________________
 441 const char*
 442 AliAnalysisMuMu::CentralityName(Double_t centrality) const
 443 {
 444   /// Get centrality name corresponding to the floating ^point value
 445
 446   if ( centrality > 0 && centrality <= 100.0 )
 447   {
 448     return Form("CENT%02d",TMath::Nint(centrality));
 449   }
 450   else
 451   {
 452     return DefaultCentralityName();
 453   }
 454 }
 455
 456 //_____________________________________________________________________________
 457 void AliAnalysisMuMu::AssertHistogramCollection(const char* physics, const char* triggerClassName)
 458 {
 459   // insure that a given set of histogram is created
 460   TH1* test = fHistogramCollection->Histo(physics,triggerClassName,DefaultCentralityName(),"Zvertex");
 461   if (!test) FillHistogramCollection(physics,triggerClassName);
 462 }
 463
 464 //_____________________________________________________________________________
 465 void AliAnalysisMuMu::FillHistogramCollection(const char* physics, const char* triggerClassName)
 466 {
 467   /// Actually create the histograms for phyics/triggerClassName
 468
 469   AliDebug(1,Form("(%s,%s)",physics,triggerClassName));
 470
 471   Double_t ptMin = 0;
 472   Double_t ptMax = 50;
 473   Int_t nbinsPt = GetNbins(ptMin,ptMax,0.25);
 474   Double_t pMin = 0;
 475   Double_t pMax = 300;
 476   Int_t nbinsP = GetNbins(pMin,pMax,1.0);
 477   Double_t etaMin = -6;
 478   Double_t etaMax = 1;
 479   Int_t nbinsEta = GetNbins(etaMin,etaMax,0.1);
 480
 481   CreateSingleHisto(physics,triggerClassName,"PtEtaMuPlus", "#mu+ P_{T} distribution vs Eta", nbinsEta,etaMin,etaMax, nbinsPt,ptMin,ptMax);
 482   CreateSingleHisto(physics,triggerClassName,"PtEtaMuMinus", "#mu- P_{T} distribution vs Eta",nbinsEta,etaMin,etaMax, nbinsPt,ptMin,ptMax);
 483
 484   CreateSingleHisto(physics,triggerClassName,"PEtaMuPlus", "#mu+ P distribution",nbinsEta,etaMin,etaMax,nbinsP,pMin,pMax);
 485   CreateSingleHisto(physics,triggerClassName,"PEtaMuMinus", "#mu- P distribution",nbinsEta,etaMin,etaMax,nbinsP,pMin,pMax);
 486
 487   Double_t chi2min = 0;
 488   Double_t chi2max = 20;
 489   Int_t nbinchi2 = GetNbins(chi2min,chi2max,0.1);
 490
 491   CreateSingleHisto(physics, triggerClassName, "Chi2MuPlus", "chisquare per NDF #mu+", nbinchi2, chi2min, chi2max);
 492   CreateSingleHisto(physics, triggerClassName, "Chi2MuMinus", "chisquare per NDF #mu-", nbinchi2, chi2min, chi2max);
 493
 494   Double_t minvMin = 0;
 495   Double_t minvMax = 16;
 496   Int_t nMinvBins = GetNbins(minvMin,minvMax,0.05);
 497
 498   CreatePairHisto(physics,triggerClassName,"MinvUSPt", "#mu+#mu- inv. mass vs Pt",nbinsPt,ptMin,ptMax,nMinvBins,minvMin,minvMax);
 499   CreatePairHisto(physics,triggerClassName,"MinvPPPt", "#mu+#mu+ inv. mass vs Pt",nbinsPt,ptMin,ptMax,nMinvBins,minvMin,minvMax);
 500   CreatePairHisto(physics,triggerClassName,"MinvMMPt", "#mu-#mu- inv. mass vs Pt",nbinsPt,ptMin,ptMax,nMinvBins,minvMin,minvMax);
 501
 502
 503   Double_t xmin = -40;
 504   Double_t xmax = +40;
 505   Int_t nbins = GetNbins(xmin,xmax,0.5);
 506
 507   CreateEventHisto(physics,triggerClassName,"Zvertex","z vertex",nbins,xmin,xmax);
 508
 509   xmin = -5;
 510   xmax = 5;
 511   nbins = GetNbins(xmin,xmax,0.01);
 512
 513   CreateEventHisto(physics,triggerClassName,"Xvertex","x vertex",nbins,xmin,xmax);
 514   CreateEventHisto(physics,triggerClassName,"Yvertex","y vertex",nbins,xmin,xmax);
 515   CreateEventHisto(physics,triggerClassName,"YXvertex","y vs x vertex",nbins,xmin,xmax,
 516                    nbins,xmin,xmax);
 517
 518   if (!fIsFromESD)
 519   {
 520
 521     CreateEventHisto(physics,triggerClassName,"PileUpZvertex","pileup z vertex",nbins,xmin,xmax);
 522
 523     CreateEventHisto(physics,triggerClassName,"PileUpXvertex","pileup x vertex",nbins,xmin,xmax);
 524     CreateEventHisto(physics,triggerClassName,"PileUpYvertex","pileup y vertex",nbins,xmin,xmax);
 525
 526     CreateEventHisto(physics,triggerClassName,"PileUpYXvertex","pileup y vs x vertex",nbins,xmin,xmax,
 527                      nbins,xmin,xmax);
 528
 529   }
 530   CreateEventHisto(physics,triggerClassName,"Nevents","number of events",2,-0.5,1.5);
 531
 532   xmin = 0;
 533   xmax = 3564;
 534   nbins = GetNbins(xmin,xmax,1.0);
 535
 536   CreateEventHisto(physics,triggerClassName,"BCX","bunch-crossing ids",nbins,xmin-0.5,xmax-0.5);
 537
 538   xmin = -200;
 539   xmax = +200;
 540   nbins = GetNbins(xmin,xmax,1.0);
 541
 542   CreateSingleHisto(physics,triggerClassName,"XYdcaMuPlus","#mu+ DCA non bending vs bending;dcaY;dcaX",nbins,xmin,xmax,nbins,xmin,xmax);
 543   CreateSingleHisto(physics,triggerClassName,"XYdcaMuMinus","#mu- DCA non bending vs bending;dcaY;dcaX",nbins,xmin,xmax,nbins,xmin,xmax);
 544
 545   xmin = 0;
 546   xmax = 300;
 547   nbins = GetNbins(xmin,xmax,1.0);
 548
 549   CreateSingleHisto(physics,triggerClassName,"dcaP23MuPlus","#mu+ DCA vs P for 2-3 degrees;P (GeV);DCA (cm)",nbinsP,pMin,pMax,nbins,xmin,xmax);
 550   CreateSingleHisto(physics,triggerClassName,"dcaP23MuMinus","#mu- DCA vs P for 2-3 degrees;P (GeV);DCA (cm)",nbinsP,pMin,pMax,nbins,xmin,xmax);
 551   CreateSingleHisto(physics,triggerClassName,"dcaP310MuPlus","#mu+ DCA vs P for 3-10 degrees;P (GeV);DCA (cm)",nbinsP,pMin,pMax,nbins,xmin,xmax);
 552   CreateSingleHisto(physics,triggerClassName,"dcaP310MuMinus","#mu- DCA vs P for 3-10 degrees;P (GeV);DCA (cm)",nbinsP,pMin,pMax,nbins,xmin,xmax);
 553
 554   CreateSingleHisto(physics,triggerClassName,"dcaPwPtCut23MuPlus","#mu+ DCA vs P for 2-3 degrees with Pt Cut;P (GeV);DCA (cm)",nbinsP,pMin,pMax,nbins,xmin,xmax);
 555   CreateSingleHisto(physics,triggerClassName,"dcaPwPtCut23MuMinus","#mu- DCA vs P for 2-3 degrees with Pt Cut;P (GeV);DCA (cm)",nbinsP,pMin,pMax,nbins,xmin,xmax);
 556   CreateSingleHisto(physics,triggerClassName,"dcaPwPtCut310MuPlus","#mu+ DCA vs P for 3-10 degrees with Pt Cut;P (GeV);DCA (cm)",nbinsP,pMin,pMax,nbins,xmin,xmax);
 557   CreateSingleHisto(physics,triggerClassName,"dcaPwPtCut310MuMinus","#mu- DCA vs P for 3-10 degrees with Pt Cut;P (GeV);DCA (cm)",nbinsP,pMin,pMax,nbins,xmin,xmax);
 558
 559   xmin = 0;
 560   xmax = 20000;
 561   nbins = GetNbins(xmin,xmax,100.0);
 562
 563   CreateSingleHisto(physics,triggerClassName,"PDCAcorrP23MuPlus","#mu+ PxDCAcorr vs P for 2-3 degrees;P (GeV/c);P.DCA (GeV.cm)",nbinsP,pMin,pMax,nbins,xmin,xmax);
 564   CreateSingleHisto(physics,triggerClassName,"PDCAcorrP23MuMinus","#mu- PxDCAcorr  vs P for 2-3 degrees;P (GeV/c);P.DCA (GeV.cm)",nbinsP,pMin,pMax,nbins,xmin,xmax);
 565   CreateSingleHisto(physics,triggerClassName,"PDCAcorrP310MuPlus","#mu+ PxDCAcorr  vs P for 3-10 degrees;P (GeV/c);P.DCA (GeV.cm)",nbinsP,pMin,pMax,nbins,xmin,xmax);
 566   CreateSingleHisto(physics,triggerClassName,"PDCAcorrP310MuMinus","#mu- PxDCAcorr  vs P for 3-10 degrees;P (GeV/c);P.DCA (GeV.cm)",nbinsP,pMin,pMax,nbins,xmin,xmax);
 567
 568   CreateSingleHisto(physics,triggerClassName,"PDCAcutP23MuPlus","#mu+ PxDCAcut vs P for 2-3 degrees;P (GeV/c);P.DCA (GeV.cm)",nbinsP,pMin,pMax,nbins,xmin,xmax);
 569   CreateSingleHisto(physics,triggerClassName,"PDCAcutP23MuMinus","#mu- PxDCAcut  vs P for 2-3 degrees;P (GeV/c);P.DCA (GeV.cm)",nbinsP,pMin,pMax,nbins,xmin,xmax);
 570   CreateSingleHisto(physics,triggerClassName,"PDCAcutP310MuPlus","#mu+ PxDCAcut  vs P for 3-10 degrees;P (GeV/c);P.DCA (GeV.cm)",nbinsP,pMin,pMax,nbins,xmin,xmax);
 571   CreateSingleHisto(physics,triggerClassName,"PDCAcutP310MuMinus","#mu- PxDCAcut  vs P for 3-10 degrees;P (GeV/c);P.DCA (GeV.cm)",nbinsP,pMin,pMax,nbins,xmin,xmax);
 572
 573
 574   if ( fIsFromESD )
 575   {
 576     xmin = -80;
 577     xmax = 80;
 578     nbins = GetNbins(xmin,xmax,0.1);
 579     CreateEventHisto(physics,triggerClassName,"V0Time","Mean Time V0C versus V0A;Time V0A (ns);Time V0C (ns)",nbins,xmin,xmax,nbins,xmin,xmax);
 580
 581     xmin = 0;
 582     xmax = 20000;
 583     nbins = GetNbins(xmin,xmax,10);
 584
 585     CreateEventHisto(physics,triggerClassName,"V0Amplitude","V0Amult+V0Cmult",nbins,xmin,xmax);
 586   }
 587
 588   if ( fAA )
 589   {
 590     Double_t* vbins = new Double_t[fCentralityLimits.size()+1];
 591
 592     vbins[0] = 0.0;
 593
 594     for ( std::vector<double>::size_type i = 0; i < fCentralityLimits.size(); ++i )
 595     {
 596       vbins[i+1] = fCentralityLimits[i];
 597     }
 598
 599     TH1* h = new TH1F("Centrality","Centrality",fCentralityLimits.size(),vbins);
 600
 601     delete[] vbins;
 602
 603     fHistogramCollection->Adopt(physics,triggerClassName,h);
 604
 605   }
 606
 607   xmin = 0;
 608   xmax = 5000;
 609   nbins = GetNbins(xmin,xmax,10);
 610
 611   CreateEventHisto(physics,triggerClassName,"Tracklets","Number of tracklets",nbins,xmin,xmax);
 612
 613 }
 614
 615
 616 //_____________________________________________________________________________
 617 Double_t AliAnalysisMuMu::MuonMass2() const
 618 {
 619   /// A usefull constant
 620   static Double_t m2 = 1.11636129640000012e-02; // using a constant here as the line below is a problem for CINT...
 621 //  static Double_t m2 = TDatabasePDG::Instance()->GetParticle("mu-")->Mass()*TDatabasePDG::Instance()->GetParticle("mu-")->Mass();
 622   return m2;
 623 }
 624
 625 //_____________________________________________________________________________
 626 void AliAnalysisMuMu::FinishTaskOutput()
 627 {
 628   /// prune empty histograms BEFORE mergin, in order to save some bytes...
 629
 630   if ( fHistogramCollection )
 631   {
 632     UInt_t size = fHistogramCollection->EstimateSize(kFALSE);
 633
 634     AliInfo(Form("size before prune histograms = %5.1f MB",size/1024.0/1024.0));
 635
 636     fHistogramCollection->PruneEmptyHistograms();
 637
 638     AliInfo(Form("size after prune histograms = %5.1f MB",size/1024.0/1024.0));
 639
 640   }
 641 }
 642
 643 //_____________________________________________________________________________
 644 TH1* AliAnalysisMuMu::Histo(const char* physics, const char* triggerClassName, const char* histoname)
 645 {
 646   /// Get one histo back
 647   return fHistogramCollection ? fHistogramCollection->Histo(physics,triggerClassName,histoname) : 0x0;
 648 }
 649
 650 //_____________________________________________________________________________
 651 TH1* AliAnalysisMuMu::Histo(const char* physics, const char* histoname)
 652 {
 653   /// Get one histo back
 654   return fHistogramCollection ? fHistogramCollection->Histo(physics,histoname) : 0x0;
 655 }
 656
 657 //_____________________________________________________________________________
 658 TH1* AliAnalysisMuMu::Histo(const char* physics,
 659                             const char* triggerClassName,
 660                             const char* what,
 661                             const char* histoname)
 662 {
 663   /// Get one histo back
 664   return fHistogramCollection ? fHistogramCollection->Histo(physics,triggerClassName,what,histoname) : 0x0;
 665 }
 666
 667 //_____________________________________________________________________________
 668 TH1* AliAnalysisMuMu::Histo(const char* physics,
 669                             const char* triggerClassName,
 670                             const char* cent,
 671                             const char* what,
 672                             const char* histoname)
 673 {
 674   /// Get one histo back
 675
 676   return fHistogramCollection ? fHistogramCollection->Histo(physics,triggerClassName,cent,what,histoname) : 0x0;
 677 }
 678
 679 //_____________________________________________________________________________
 680 void AliAnalysisMuMu::NotifyRun()
 681 {
 682   /// Called at each change of run ?
 683   AliInfo(Form("Run %09d File %s",fCurrentRunNumber,CurrentFileName()));
 684 }
 685
 686 //_____________________________________________________________________________
 687 void
 688 AliAnalysisMuMu::Print(Option_t* /*opt*/) const
 689 {
 690   /// Print the definition of this analysis
 691
 692   cout << ClassName() << " - " << GetName() << " - " << ( fAA ? "PbPb mode" : "pp mode" ) << endl;
 693
 694   if  ( !fSingleTrackCutNames || !fSingleTrackCutNames )
 695   {
 696     cout << "No single track cut defined yet" << endl;
 697   }
 698   else
 699   {
 700     TIter next(fSingleTrackCutNames);
 701     TObjString* str;
 702
 703     while ( ( str = static_cast<TObjString*>(next()) ) )
 704     {
 705       cout << Form("SINGLE CUT %20s MASK %x",str->String().Data(),str->GetUniqueID()) << endl;
 706     }
 707   }
 708
 709   if  ( !fPairTrackCutNames || !fPairTrackCutNames )
 710   {
 711     cout << "No track pair cut defined yet" << endl;
 712   }
 713   else
 714   {
 715     TIter next2(fPairTrackCutNames);
 716     TObjString* str;
 717
 718     while ( ( str = static_cast<TObjString*>(next2()) ) )
 719     {
 720       UInt_t single(0);
 721       UInt_t pair(0);
 722       DecodePairCutMask(str->GetUniqueID(),single,pair);
 723
 724       cout << Form("PAIR   CUT %20s UID %x SINGLE MASK %x PAIR MASK %x",str->String().Data(),str->GetUniqueID(),single,pair) << endl;
 725     }
 726   }
 727
 728   if ( !fTriggerClasses || !fTriggerClasses->First() )
 729   {
 730     cout << "No trigger classes defined yet" << endl;
 731   }
 732   else
 733   {
 734     cout << "Trigger classes that will be considered:" << endl;
 735     TIter next(fTriggerClasses);
 736     TObjString* s;
 737
 738     while ( ( s = static_cast<TObjString*>(next()) ) )
 739     {
 740       cout << s->String().Data() << endl;
 741     }
 742   }
 743 }
 744
 745 //_____________________________________________________________________________
 746 void
 747 AliAnalysisMuMu::Terminate(Option_t *)
 748 {
 749   /// Called once at the end of the query
 750   /// Just a simple printout of the stat we analyse and how many histograms
 751   /// we got
 752
 753   fOutput = dynamic_cast<TList*>(GetOutputData(1));
 754
 755   if (!fOutput) return;
 756
 757   fHistogramCollection = dynamic_cast<AliHistogramCollection*>(fOutput->FindObject("mumu"));
 758
 759   if (!fHistogramCollection)
 760   {
 761     AliError("Could not find back histogram collection in output...");
 762     return;
 763   }
 764
 765   fHistogramCollection->PruneEmptyHistograms();
 766
 767   UInt_t size2 = fHistogramCollection->EstimateSize();
 768
 769   AliInfo(Form("size after prune histograms = %5.1f MB",size2/1024.0/1024.0));
 770
 771   BeautifyHistos();
 772
 773   fHistogramCollection->Print("*Minv*");
 774
 775   fEventCounters = dynamic_cast<AliCounterCollection*>(fOutput->FindObject("eventCounters"));
 776
 777   if (!fEventCounters)
 778   {
 779     AliError("Could not find back counters in output...");
 780     return;
 781   }
 782
 783   fEventCounters->Print();
 784 }
 785
 786
 787 //_____________________________________________________________________________
 788 void AliAnalysisMuMu::UserExec(Option_t* opt)
 789 {
 790   /// Executed at each event
 791
 792   MuUserExec(opt);
 793
 794   // Post output data.
 795   PostData(1, fOutput);
 796 }
 797
 798 //_____________________________________________________________________________
 799 void AliAnalysisMuMu::UserCreateOutputObjects()
 800 {
 801   // Create histograms
 802   // Called once
 803
 804   fTriggerClasses->Print();
 805
 806   fOutput = new TList;
 807   fOutput->SetOwner(kTRUE);
 808
 809   fHistogramCollection = new AliHistogramCollection("mumu");
 810
 811   fOutput->Add(fHistogramCollection);
 812
 813   // initialize event counters
 814   fEventCounters = new AliCounterCollection("eventCounters");
 815
 816   TIter nextGSN(fGlobalEventSelectionNames);
 817   TObjString* str;
 818   TString eventRubric;
 819   while ( ( str = static_cast<TObjString*>(nextGSN()) ) )
 820   {
 821     if ( eventRubric.Length() > 0 ) eventRubric += "/";
 822     eventRubric += str->String();
 823   }
 824
 825   fEventCounters->AddRubric("event", eventRubric.Data());
 826
 827   fEventCounters->AddRubric("trigger", 100);
 828
 829   fEventCounters->AddRubric("run", 1000000);
 830
 831   fEventCounters->Init();
 832
 833   fOutput->Add(fEventCounters);
 834
 835   // Post output data.
 836   PostData(1, fOutput);
 837 }