PWG/muon/AliAnalysisMuMuMinv.cxx

   1 #include "AliAnalysisMuMuMinv.h"
   2
   3 /**
   4  *
   5  * \ingroup pwg-muon-mumu
   6  *
   7  * \class AliAnalysisMuMuMinv
   8  *
   9  * Analysis which fills a bunch of histograms for invariant mass analysis of J/psi
  10  *
  11  * Can be used on real data and/or simulated (for instance to get Acc x Eff)
  12  *
  13  * Can optionally use as input an already computed Acc x Eff matrix that will be applied
  14  * when filling the invariant mass histograms.
  15  *
  16  */
  17
  18 #include "TH2F.h"
  19 #include "AliLog.h"
  20 #include "TObjArray.h"
  21 #include "AliAnalysisMuMuBinning.h"
  22 #include "TString.h"
  23 #include "TLorentzVector.h"
  24 #include "TString.h"
  25 #include "AliMCEvent.h"
  26 #include "AliMergeableCollection.h"
  27 #include "AliAnalysisMuonUtility.h"
  28 #include "TParameter.h"
  29 #include <cassert>
  30
  31 ClassImp(AliAnalysisMuMuMinv)
  32
  33 //_____________________________________________________________________________
  34 AliAnalysisMuMuMinv::AliAnalysisMuMuMinv()
  35 : AliAnalysisMuMuBase(),
  36 fAccEffHisto(0x0)
  37 {
  38   // FIXME : find the AccxEff histogram from HistogramCollection()->Histo("/EXCHANGE/JpsiAccEff")
  39
  40 //  if ( accEff )
  41 //  {
  42 //    fAccEffHisto = static_cast<TH2F*>(accEff->Clone());
  43 //    fAccEffHisto->SetDirectory(0);
  44 //  }
  45 }
  46
  47 //_____________________________________________________________________________
  48 AliAnalysisMuMuMinv::~AliAnalysisMuMuMinv()
  49 {
  50   /// dtor
  51   delete fAccEffHisto;
  52 }
  53
  54 //_____________________________________________________________________________
  55 void
  56 AliAnalysisMuMuMinv::DefineHistogramCollection(const char* eventSelection,
  57                                                const char* triggerClassName,
  58                                                const char* centrality)
  59 {
  60   /// Define the histograms this analysis will use
  61
  62   if ( Histo(eventSelection,triggerClassName,centrality,"AliAnalysisMuMuMinv") )
  63   {
  64     return;
  65   }
  66
  67   // dummy histogram to signal that we already defined all our histograms (see above)
  68   CreateEventHistos(kHistoForData | kHistoForMCInput,eventSelection,triggerClassName,centrality,"AliAnalysisMuMuMinv","Dummy semaphore",1,0,1);
  69
  70   /// Create invariant mass histograms
  71
  72   Double_t minvMin = 0;
  73   Double_t minvMax = 16;
  74   Int_t nMinvBins = GetNbins(minvMin,minvMax,0.025);
  75
  76   Int_t nMCMinvBins = GetNbins(minvMin,minvMax,0.1);
  77
  78   //  TObjArray* bins = fBinning->CreateBinObjArray("psi","y vs pt,integrated,pt,y,phi","");
  79   TObjArray* bins = Binning()->CreateBinObjArray("psi");
  80
  81   CreatePairHistos(eventSelection,triggerClassName,centrality,"Pt","#mu+#mu- Pt distribution",
  82                   200,0,20);
  83
  84   Int_t nbinsy = 6;
  85   Double_t ymin = -4.0;
  86   Double_t ymax = -2.5;
  87
  88   CreatePairHistos(eventSelection,triggerClassName,centrality,"y","#mu+#mu- y distribution",nbinsy,ymin,ymax);
  89
  90   //  CreatePairHistos(eventSelection,triggerClassName,centrality,"BinFlowPt","#mu+#mu- BinFlowPt distribution",
  91   //                  200,0,20);
  92
  93   CreatePairHistos(eventSelection,triggerClassName,centrality,"PtRecVsSim","#mu+#mu- Pt distribution rec vs sim",
  94                   200,0,20,200,0,20);
  95
  96   if (!Histo("INPUT","ALL","Pt"))
  97   {
  98     HistogramCollection()->Adopt("/INPUT/ALL",new TH1F("Pt","Pt",200,0,20));
  99     HistogramCollection()->Adopt("/INPUT/INYRANGE",new TH1F("Pt","Pt",200,0,20));
 100     HistogramCollection()->Adopt("/INPUT/ALL",new TH1F("Y","Y",nbinsy,ymin,ymax));
 101     HistogramCollection()->Adopt("/INPUT/INYRANGE",new TH1F("Y","Y",nbinsy,ymin,ymax));
 102   }
 103
 104   TIter next(bins);
 105   AliAnalysisMuMuBinning::Range* r;
 106   Int_t nb(0);
 107
 108   while ( ( r = static_cast<AliAnalysisMuMuBinning::Range*>(next()) ) )
 109   {
 110     TString minvName(GetMinvHistoName(*r,ShouldCorrectDimuonForAccEff()));
 111
 112     if ( IsHistogramDisabled(minvName.Data()) ) continue;
 113
 114     ++nb;
 115
 116     AliDebug(1,Form("bin %d %s histoname = %s",nb,r->AsString().Data(),minvName.Data()));
 117
 118     CreatePairHistos(eventSelection,triggerClassName,centrality,minvName.Data(),
 119                      Form("#mu+#mu- inv. mass %s %s;M_{#mu^+#mu^-} (GeV/c^2)",
 120                           ShouldCorrectDimuonForAccEff() ? "(AccxEff corrected)":"",
 121                           r->AsString().Data()),
 122                      nMinvBins,minvMin,minvMax);
 123
 124     TString hname(GetMinvHistoName(*r,kFALSE));
 125
 126     TH1* h = HistogramCollection()->Histo("/INPUT/ALL",hname.Data());
 127     if (!h)
 128     {
 129       h = new TH1F(hname.Data(),Form("MC #mu+#mu- inv. mass %s",r->AsString().Data()),
 130                    nMCMinvBins,minvMin,minvMax);
 131
 132       HistogramCollection()->Adopt("/INPUT/ALL",h);
 133
 134       HistogramCollection()->Adopt("/INPUT/INYRANGE",static_cast<TH1*>(h->Clone()));
 135     }
 136   }
 137
 138   delete bins;
 139 }
 140
 141 //_____________________________________________________________________________
 142 void AliAnalysisMuMuMinv::FillHistosForPair(const char* eventSelection, const char* triggerClassName,
 143                                             const char* centrality, const char* pairCutName,
 144                                             const AliVParticle& tracki,
 145                                             const AliVParticle& trackj)
 146 {
 147   /** Fill histograms for muon track pairs
 148    */
 149
 150   if (!AliAnalysisMuonUtility::IsMuonTrack(&tracki) ) return;
 151   if (!AliAnalysisMuonUtility::IsMuonTrack(&trackj) ) return;
 152
 153   TLorentzVector pi(tracki.Px(),tracki.Py(),tracki.Pz(),
 154                     TMath::Sqrt(AliAnalysisMuonUtility::MuonMass2()+tracki.P()*tracki.P()));
 155
 156
 157   TLorentzVector pair4Momentum(trackj.Px(),trackj.Py(),trackj.Pz(),
 158                                TMath::Sqrt(AliAnalysisMuonUtility::MuonMass2()+trackj.P()*trackj.P()));
 159
 160   pair4Momentum += pi;
 161
 162
 163   if (!IsHistogramDisabled("Chi12"))
 164   {
 165     Histo(eventSelection,triggerClassName,centrality,pairCutName,"Chi12")
 166     ->Fill(
 167            AliAnalysisMuonUtility::GetChi2perNDFtracker(&tracki),
 168            AliAnalysisMuonUtility::GetChi2perNDFtracker(&trackj));
 169   }
 170
 171   if (!IsHistogramDisabled("Rabs12"))
 172   {
 173     Histo(eventSelection,triggerClassName,centrality,pairCutName,"Rabs12")
 174     ->Fill(AliAnalysisMuonUtility::GetRabs(&tracki),
 175            AliAnalysisMuonUtility::GetRabs(&trackj));
 176   }
 177
 178   if ( ( tracki.Charge() != trackj.Charge() ) )
 179   {
 180     if ( !IsHistogramDisabled("Pt") )
 181     {
 182       Histo(eventSelection,triggerClassName,centrality,pairCutName,"Pt")->Fill(pair4Momentum.Pt());
 183     }
 184
 185     if ( !IsHistogramDisabled("y") )
 186     {
 187       Histo(eventSelection,triggerClassName,centrality,pairCutName,"y")->Fill(pair4Momentum.Y());
 188     }
 189
 190     if ( HasMC() )
 191     {
 192       Int_t labeli = tracki.GetLabel();
 193       Int_t labelj = trackj.GetLabel();
 194
 195       if ( labeli < 0 || labelj < 0 )
 196       {
 197         AliError("Got negative labels!");
 198       }
 199       else
 200       {
 201         AliVParticle* mcTracki = MCEvent()->GetTrack(labeli);
 202         AliVParticle* mcTrackj = MCEvent()->GetTrack(labelj);
 203
 204         TLorentzVector mcpi(mcTracki->Px(),mcTracki->Py(),mcTracki->Pz(),
 205                             TMath::Sqrt(AliAnalysisMuonUtility::MuonMass2()+mcTracki->P()*mcTracki->P()));
 206         TLorentzVector mcpj(mcTrackj->Px(),mcTrackj->Py(),mcTrackj->Pz(),
 207                             TMath::Sqrt(AliAnalysisMuonUtility::MuonMass2()+mcTrackj->P()*mcTrackj->P()));
 208
 209         mcpj += mcpi;
 210
 211         Histo(eventSelection,triggerClassName,centrality,pairCutName,"PtRecVsSim")->Fill(mcpj.Pt(),pair4Momentum.Pt());
 212
 213       }
 214     }
 215
 216
 217     TObjArray* bins = Binning()->CreateBinObjArray("psi","ptvsy,yvspt,pt,y,phi","");
 218     TIter nextBin(bins);
 219     AliAnalysisMuMuBinning::Range* r;
 220
 221     while ( ( r = static_cast<AliAnalysisMuMuBinning::Range*>(nextBin()) ) )
 222     {
 223       Bool_t ok(kFALSE);
 224
 225       if ( r->IsIntegrated() )
 226       {
 227         ok = kTRUE;
 228       }
 229       else if ( r->Is2D() )
 230       {
 231         if ( r->AsString().BeginsWith("PTVSY") )
 232         {
 233           ok = r->IsInRange(pair4Momentum.Rapidity(),pair4Momentum.Pt());
 234         }
 235         else if ( r->AsString().BeginsWith("YVSPT") )
 236         {
 237           ok = r->IsInRange(pair4Momentum.Pt(),pair4Momentum.Rapidity());
 238         }
 239         else
 240         {
 241           AliError(Form("Don't know how to deal with 2D bin %s",r->AsString().Data()));
 242         }
 243       }
 244       else
 245       {
 246         if ( r->Quantity() == "PT" )
 247         {
 248           ok = r->IsInRange(pair4Momentum.Pt());
 249         }
 250         else if ( r->Quantity() == "Y" )
 251         {
 252           ok = r->IsInRange(pair4Momentum.Rapidity());
 253         }
 254         else if ( r->Quantity() == "PHI" )
 255         {
 256           ok = r->IsInRange(pair4Momentum.Phi());
 257         }
 258         else if ( r->Quantity() == "NCH" )
 259         {
 260           TList* list = static_cast<TList*>(Event()->FindListObject("NCH"));
 261           if (list)
 262           {
 263             TIter next(list);
 264             TParameter<Double_t>* p;
 265
 266             while ( ( p = static_cast<TParameter<Double_t>*>(next())) )
 267             {
 268               if (TString(eventSelection).Contains(p->GetName()))
 269               {
 270                 ok = r->IsInRange(p->GetVal());
 271                 break;
 272               }
 273             }
 274           }
 275         }
 276       }
 277
 278       if ( ok )
 279       {
 280         TString minvName = GetMinvHistoName(*r,kFALSE);
 281
 282         if (!IsHistogramDisabled(minvName.Data()))
 283         {
 284           TH1* h = Histo(eventSelection,triggerClassName,centrality,pairCutName,minvName.Data());
 285
 286           if (!h)
 287           {
 288             AliError(Form("Could not get %s",minvName.Data()));
 289             continue;
 290           }
 291           h->Fill(pair4Momentum.M());
 292
 293           if ( fAccEffHisto )
 294           {
 295             // FIXME : fill Minv with weight = 1/AccxEff
 296           }
 297         }
 298       }
 299     }
 300
 301     delete bins;
 302   }
 303 }
 304
 305 //_____________________________________________________________________________
 306 void AliAnalysisMuMuMinv::FillHistosForMCEvent(const char* /*eventSelection*/,
 307                                                const char* /*triggerClassName*/,
 308                                                const char* /*centrality*/)
 309 {
 310   /** Fill histograms for MC event
 311    *
 312    * FIXME: this is here we should streamline a bit the code to carefully
 313    * compute the measured and truth values for Pt, Y, and (Pt,Y), in order
 314    * to be able to investigate unfolding techniques later on...
 315    */
 316
 317   Int_t nMCTracks = MCEvent()->GetNumberOfTracks();
 318
 319   for ( Int_t i = 0; i < nMCTracks; ++i )
 320   {
 321     AliVParticle* part = MCEvent()->GetTrack(i);
 322
 323     if  ( AliAnalysisMuonUtility::IsPrimary(part,MCEvent()) &&
 324           AliAnalysisMuonUtility::GetMotherIndex(part)==-1 )
 325     {
 326             Histo("INPUT","ALL","Pt")->Fill(part->Pt());
 327             Histo("INPUT","ALL","Y")->Fill(part->Y());
 328             if ( part->Y() < -2.5 && part->Y() > -4.0 )
 329             {
 330               Histo("INPUT","INYRANGE","Pt")->Fill(part->Pt());
 331               Histo("INPUT","INYRANGE","Y")->Fill(part->Y());
 332             }
 333
 334     }
 335
 336   }
 337
 338 //  Int_t nTracks = AliAnalysisMuonUtility::GetNTracks(Event());
 339 //
 340 //  Double_t measuredY(0.0);
 341 //
 342 //  for (Int_t i = 0; i < nTracks; ++i)
 343 //  {
 344 //    AliVParticle* tracki = AliAnalysisMuonUtility::GetTrack(i,Event());
 345 //
 346 //    if (!AliAnalysisMuonUtility::IsMuonTrack(tracki) ) continue;
 347 //
 348 //    for (Int_t j = i+1; j < nTracks; ++j )
 349 //    {
 350 //      AliVParticle* trackj = AliAnalysisMuonUtility::GetTrack(j,Event());
 351 //
 352 //      if (!AliAnalysisMuonUtility::IsMuonTrack(trackj) ) continue;
 353 //
 354 //      TLorentzVector pi(tracki->Px(),tracki->Py(),tracki->Pz(),
 355 //                        TMath::Sqrt(AliAnalysisMuonUtility::MuonMass2()+tracki->P()*tracki->P()));
 356 //
 357 //      TLorentzVector pair4Momentum(trackj->Px(),trackj->Py(),trackj->Pz(),
 358 //                                   TMath::Sqrt(AliAnalysisMuonUtility::MuonMass2()+trackj->P()*trackj->P()));
 359 //
 360 //      pair4Momentum += pi;
 361 //
 362 //      measuredY = pi.Y();
 363 //
 364 //      rur->Fill(measuredY,trueY); // FIXME : is this working if more than one pair is found in the reconstructed tracks ???
 365 //      // should we try to find the closest one in minv and assign the other(s)
 366 //      // to rur->Fake() ?
 367 //    }
 368 //
 369 //  }
 370 //
 371 //  if ( ! ( measuredY < 0.0 ) )
 372 //  {
 373 //    rur->Miss(trueY);
 374 //  }
 375 }
 376
 377 ////_____________________________________________________________________________
 378 //void AliAnalysisMuMuMinv::FillMC()
 379 //{
 380 //  // Fill the input Monte-Carlo histograms related to muons
 381 //
 382 //  // Specific things for MC
 383 //  if (!Histo("INPUT","ALL","Pt"))
 384 //  {
 385 //    HistogramCollection()->Adopt("/INPUT/ALL",new TH1F("Pt","Pt",200,0,20));
 386 //    HistogramCollection()->Adopt("/INPUT/INYRANGE",new TH1F("Pt","Pt",200,0,20));
 387 //
 388 //    Double_t rapidityMin = -5;
 389 //    Double_t rapidityMax = -2;
 390 //    Int_t nbinsRapidity = GetNbins(rapidityMin,rapidityMax,0.05);
 391 //
 392 //    HistogramCollection()->Adopt("/INPUT/ALL",new TH1F("Y","Y",nbinsRapidity,rapidityMin,rapidityMax));
 393 //    HistogramCollection()->Adopt("/INPUT/INYRANGE",new TH1F("Y","Y",nbinsRapidity,rapidityMin,rapidityMax));
 394 //
 395 //    Double_t etaMin = -5;
 396 //    Double_t etaMax = -2;
 397 //    Int_t nbinsEta = GetNbins(etaMin,etaMax,0.05);
 398 //
 399 //    HistogramCollection()->Adopt("/INPUT/ALL",new TH1F("Eta","Eta",nbinsEta,etaMin,etaMax));
 400 //    HistogramCollection()->Adopt("/INPUT/INYRANGE",new TH1F("Eta","Eta",nbinsEta,etaMin,etaMax));
 401 //  }
 402 //
 403 //  Int_t nMCTracks = MCEvent()->GetNumberOfTracks();
 404 //
 405 //  TObjArray* bins = Binning()->CreateBinObjArray("psi","ptvsy,yvspt,pt,y,phi","");
 406 //  TIter nextBin(bins);
 407 //  AliAnalysisMuMuBinning::Range* r;
 408 //
 409 //  for ( Int_t i = 0; i < nMCTracks; ++i )
 410 //  {
 411 //    AliVParticle* part = MCEvent()->GetTrack(i);
 412 //
 413 //    //    std::cout << "part " << i << " isprimary=" << AliAnalysisMuonUtility::IsPrimary(part,MCEvent()) << " motherindex=" << AliAnalysisMuonUtility::GetMotherIndex(part) << std::endl;
 414 //    //
 415 //    //    part->Print();
 416 //
 417 //    if  (AliAnalysisMuonUtility::IsPrimary(part,MCEvent()) &&
 418 //         AliAnalysisMuonUtility::GetMotherIndex(part)==-1)
 419 //    {
 420 //
 421 //      Histo("INPUT","ALL","Pt")->Fill(part->Pt());
 422 //      Histo("INPUT","ALL","Y")->Fill(part->Y());
 423 //      Histo("INPUT","ALL","Eta")->Fill(part->Eta());
 424 //
 425 //      if ( part->Y() < -2.5 && part->Y() > -4.0 )
 426 //      {
 427 //        Histo("INPUT","INYRANGE","Pt")->Fill(part->Pt());
 428 //        Histo("INPUT","INYRANGE","Y")->Fill(part->Y());
 429 //        Histo("INPUT","INYRANGE","Eta")->Fill(part->Eta());
 430 //      }
 431 //
 432 //      nextBin.Reset();
 433 //
 434 //      while ( ( r = static_cast<AliAnalysisMuMuBinning::Range*>(nextBin()) ) )
 435 //      {
 436 //        Bool_t ok(kFALSE);
 437 //
 438 //        if ( r->IsIntegrated() )
 439 //        {
 440 //          ok = kTRUE;
 441 //        }
 442 //        else if ( r->Is2D() )
 443 //        {
 444 //          if ( r->AsString().BeginsWith("PTVSY") )
 445 //          {
 446 //            ok = r->IsInRange(part->Y(),part->Pt());
 447 //          }
 448 //          else if ( r->AsString().BeginsWith("YVSPT") )
 449 //          {
 450 //            ok = r->IsInRange(part->Pt(),part->Y());
 451 //          }
 452 //          else
 453 //          {
 454 //            AliError(Form("Don't know how to deal with 2D bin %s",r->AsString().Data()));
 455 //          }
 456 //        }
 457 //        else
 458 //        {
 459 //          if ( r->Quantity() == "PT" )
 460 //          {
 461 //            ok = r->IsInRange(part->Pt());
 462 //          }
 463 //          else if ( r->Quantity() == "Y" )
 464 //          {
 465 //            ok = r->IsInRange(part->Y());
 466 //          }
 467 //          else if ( r->Quantity() == "PHI" )
 468 //          {
 469 //            ok = r->IsInRange(part->Phi());
 470 //          }
 471 //        }
 472 //
 473 //        if ( ok )
 474 //        {
 475 //          TString hname = GetMinvHistoName(*r,kFALSE);
 476 //
 477 //          if (!IsHistogramDisabled(hname.Data()))
 478 //          {
 479 //
 480 //            TH1* h = Histo("INPUT","ALL",hname.Data());
 481 //
 482 //            if (!h)
 483 //            {
 484 //              AliError(Form("Could not get ALL %s",hname.Data()));
 485 //              continue;
 486 //            }
 487 //
 488 //            h->Fill(part->M());
 489 //
 490 //            if ( part->Y() < -2.5 && part->Y() > -4.0 )
 491 //            {
 492 //              h = Histo("INPUT","INYRANGE",hname.Data());
 493 //              if (!h)
 494 //              {
 495 //                AliError(Form("Could not get INYRANGE %s",hname.Data()));
 496 //                continue;
 497 //              }
 498 //              h->Fill(part->M());
 499 //            }
 500 //
 501 //          }
 502 //
 503 //        }
 504 //      }
 505 //    }
 506 //  }
 507 //
 508 //  delete bins;
 509 //}
 510 //
 511
 512 //_____________________________________________________________________________
 513 TString AliAnalysisMuMuMinv::GetMinvHistoName(const AliAnalysisMuMuBinning::Range& r, Bool_t accEffCorrected) const
 514 {
 515   /// Get the invariant mass histogram name
 516   return TString::Format("MinvUS%s%s",r.AsString().Data(),
 517                          accEffCorrected ? "_AccEffCorr" : "");
 518 }
 519
 520 //_____________________________________________________________________________
 521 Bool_t AliAnalysisMuMuMinv::IsRapidityInRange(const AliVParticle& t1, const AliVParticle& t2, Double_t& ymin, Double_t& ymax) const
 522 {
 523   /// Whether the particle pair has its rapidity within [ymin,ymax[
 524
 525   TLorentzVector p1(t1.Px(),t1.Py(),t1.Pz(),TMath::Sqrt(AliAnalysisMuonUtility::MuonMass2()+t1.P()*t1.P()));
 526   TLorentzVector p2(t2.Px(),t2.Py(),t2.Pz(),TMath::Sqrt(AliAnalysisMuonUtility::MuonMass2()+t2.P()*t2.P()));
 527
 528   TLorentzVector total(p1+p2);
 529
 530   Double_t y = total.Rapidity();
 531
 532   return  ( y < ymax && y > ymin );
 533 }
 534
 535 //_____________________________________________________________________________
 536 void AliAnalysisMuMuMinv::NameOfIsRapidityInRange(TString& name, Double_t& ymin, Double_t& ymax) const
 537 {
 538   /** Get the name of the rapidity range (making the IsRapidityInRange method useable as an
 539    * AliAnalysisMuMuCutElement
 540    */
 541
 542   name.Form("PAIRYIN%2.1f-%2.1f",ymin,ymax);
 543 }
 544