PWG/muon/AliAnalysisMuMuFromAOD.cxx

   1 /**************************************************************************
   2 * Copyright(c) 1998-1999, 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$
  17
  18 #include "AliAnalysisMuMuFromAOD.h"
  19
  20 #include "AliAODEvent.h"
  21 #include "AliAODHeader.h"
  22 #include "AliAODMCHeader.h"
  23 #include "AliAODMCParticle.h"
  24 #include "AliAODTrack.h"
  25 #include "AliAnalysisManager.h"
  26 #include "AliCentrality.h"
  27 #include "AliCodeTimer.h"
  28 #include "AliCounterCollection.h"
  29 #include "AliHistogramCollection.h"
  30 #include "AliLog.h"
  31 #include "AliMCEvent.h"
  32 #include "AliVParticle.h"
  33 #include "Riostream.h"
  34 #include "TCanvas.h"
  35 #include "TDatabasePDG.h"
  36 #include "TDirectory.h"
  37 #include "TH1F.h"
  38 #include "TH2F.h"
  39 #include "THashList.h"
  40 #include "TLorentzVector.h"
  41 #include "TMap.h"
  42 #include "TMath.h"
  43 #include "TObjArray.h"
  44 #include "TPaveText.h"
  45 #include <cassert>
  46 #include <set>
  47
  48 ///
  49 /// AliAnalysisMuMuFromAOD
  50 ///
  51 /// The list of cut considered is entered using the AddSingleCut and
  52 /// AddPairCut methods.
  53 ///
  54 ///
  55
  56 ClassImp(AliAnalysisMuMuFromAOD)
  57
  58 //_____________________________________________________________________________
  59 AliAnalysisMuMuFromAOD::AliAnalysisMuMuFromAOD()
  60 : AliAnalysisMuMu(), fVertex(0), fGlobalEventSelectionName("ALL")
  61 {
  62   // default ctor
  63 }
  64
  65 //_____________________________________________________________________________
  66 AliAnalysisMuMuFromAOD::AliAnalysisMuMuFromAOD(Bool_t aa)
  67 : AliAnalysisMuMu(kFALSE,aa), fVertex(0), fGlobalEventSelectionName("ALL")
  68 {
  69   // Constructor
  70   // Must define here the single track and track pair cuts,
  71   // and also the global event selection "names"
  72   Ctor(fGlobalEventSelectionName);
  73 }
  74
  75 //_____________________________________________________________________________
  76 AliAnalysisMuMuFromAOD::AliAnalysisMuMuFromAOD(TList* triggerClasses)
  77 : AliAnalysisMuMu(kFALSE,triggerClasses), fVertex(0), fGlobalEventSelectionName("ALL")
  78 {
  79   // Constructor
  80   // Must define here the single track and track pair cuts,
  81   // and also the global event selection "names"
  82
  83   Ctor(fGlobalEventSelectionName);
  84 }
  85
  86 //_____________________________________________________________________________
  87 void AliAnalysisMuMuFromAOD::Ctor(const char* globaleventselectionname)
  88 {
  89   /// Common ctor
  90
  91   if ( fAA )
  92   {
  93     fGlobalEventSelectionName="SD2";
  94   }
  95   else
  96   {
  97     AddGlobalEventSelection(globaleventselectionname);
  98   }
  99 }
 100
 101 //_____________________________________________________________________________
 102 AliAnalysisMuMuFromAOD::~AliAnalysisMuMuFromAOD()
 103 {
 104   // dtor
 105 }
 106
 107 //_____________________________________________________________________________
 108 void AliAnalysisMuMuFromAOD::DumpMC(const AliAODEvent& aod)
 109 {
 110   // dump mc part (if present)
 111
 112   //  MC header
 113   AliAODMCHeader *mcHeader = static_cast<AliAODMCHeader*>(aod.FindListObject(AliAODMCHeader::StdBranchName()));
 114   if(mcHeader)
 115   {
 116     AliInfo(Form("================================ Generator %s x %e y %e z %e",
 117                  mcHeader->GetGeneratorName(),
 118                  mcHeader->GetVtxX(),
 119                  mcHeader->GetVtxY(),
 120                  mcHeader->GetVtxZ()));
 121   }
 122
 123   // MC particles
 124   TClonesArray* mcParticles = static_cast<TClonesArray*>(aod.FindListObject(AliAODMCParticle::StdBranchName()));
 125   if (mcParticles)
 126   {
 127     AliInfo(Form("%d mcparticles in that event",mcParticles->GetEntries()));
 128
 129     // loop on muon tracks are find back their decay chain
 130
 131     TIter next(aod.GetTracks());
 132     AliAODTrack* t;
 133 //    Int_t nmu(0);
 134 //
 135 //    AliInfo("Muon tracks");
 136 //
 137 //    while ( ( t = static_cast<AliAODTrack*>(next()) ) )
 138 //    {
 139 //      if ( t->IsMuonTrack() )
 140 //      {
 141 //        AliInfo(Form("mutrack #%2d Pt %e Eta %e Charge %d Label %d",
 142 //                     nmu,t->Pt(),t->Eta(),t->Charge(),t->GetLabel()));
 143 //        ++nmu;
 144 //      }
 145 //    }
 146 //
 147 //    AliInfo("MC particles...");
 148     TIter nextMC(mcParticles);
 149     AliAODMCParticle* p;
 150 //    Int_t nmc(0);
 151 //
 152 //    while ( ( p = static_cast<AliAODMCParticle*>(nextMC()) ) )
 153 //    {
 154 //      AliInfo(Form("mcpart #%4d",nmc));
 155 //      p->Print();
 156 //      ++nmc;
 157 //    }
 158 //    next.Reset();
 159
 160     AliInfo("Muon track ancestors...");
 161
 162     while ( ( t = static_cast<AliAODTrack*>(next()) ) )
 163     {
 164       if ( t->IsMuonTrack() )
 165       {
 166         AliInfo(Form("Track label %d",t->GetLabel()));
 167         Int_t label = t->GetLabel();
 168
 169         while ( label >= 0 )
 170         {
 171           p = static_cast<AliAODMCParticle*>(mcParticles->At(label));
 172           if (p)
 173           {
 174             p->Print();
 175             label = p->GetMother();
 176           }
 177           else
 178           {
 179             label = -1;
 180           }
 181         }
 182       }
 183     }
 184
 185   }
 186
 187 }
 188
 189 //_____________________________________________________________________________
 190 void AliAnalysisMuMuFromAOD::MuUserExec(Option_t*)
 191 {
 192   // Main loop
 193   // Called for each event
 194
 195   AliAODEvent* aod = dynamic_cast<AliAODEvent*>(InputEvent());
 196   if ( !aod ) return;
 197
 198   if ( fAA )
 199   {
 200     // consider only events with OSM2 fired
 201     UInt_t trigger = aod->GetHeader()->GetL0TriggerInputs();
 202     UInt_t sd2 = (1<<12);
 203     Bool_t ok = ( ( trigger & sd2 ) == sd2);
 204     if (!ok) return;
 205   }
 206
 207   if ( IsDynamicTriggerClasses() ) AddTriggerClasses(aod->GetFiredTriggerClasses());
 208
 209   TObjArray* a = aod->GetFiredTriggerClasses().Tokenize(" ");
 210
 211   TIter next(a);
 212   TObjString* tname;
 213
 214   TString centrality(DefaultCentralityName());
 215
 216   if ( fAA )
 217   {
 218     AliCentrality* acent = aod->GetCentrality();
 219
 220     Float_t fcent(-1);
 221
 222     if (acent) fcent = acent->GetCentralityPercentile("V0M");
 223
 224     Double_t cent(-100);
 225
 226     if ( fcent > 0 )
 227     {
 228       for ( std::vector<double>::size_type i = 0 ; i < fCentralityLimits.size() && cent < 0 ; ++i )
 229       {
 230         if ( fcent < fCentralityLimits[i] )
 231         {
 232           cent = fCentralityLimits[i];
 233         }
 234       }
 235     }
 236
 237     if ( cent > -1 )
 238     {
 239       centrality = CentralityName(cent);
 240     }
 241   }
 242
 243   fVertex = aod->GetPrimaryVertex();
 244
 245 //  Bool_t hasPileUp = aod->IsPileupFromSPD(3,0.8);
 246 //  Bool_t hasPileUp2 = aod->IsPileupFromSPD(5,0.8);
 247 //  Bool_t isIsolated = ( aod->GetBunchCrossNumber() > 1000 && aod->GetBunchCrossNumber() < 2900 );
 248
 249 //  Bool_t hasPileUp2 = aod->IsPileupFromSPDInMultBins();
 250
 251 //  TIter nextV(aod->GetVertices());
 252 //  AliAODVertex* v;
 253 //  while ( ( v = static_cast<AliAODVertex*>(nextV())) && !hasPileUp2 )
 254 //  {
 255 //    if ( v->GetType() == AliAODVertex::kPileupSPD ) hasPileUp2 = kTRUE;
 256 //  }
 257
 258   TString eventtype(fGlobalEventSelectionName);
 259   eventtype.ToLower();
 260
 261   while ( ( tname = static_cast<TObjString*>(next()) ) )
 262   {
 263     if ( !fTriggerClasses->FindObject(tname->GetName()) )
 264     {
 265       continue;
 266       // skip the trigger classes we're not supposed to analyse
 267     }
 268
 269
 270     fEventCounters->Count(Form("event:%s/trigger:%s/run:%d", eventtype.Data(), tname->GetName(), fCurrentRunNumber));
 271
 272     AssertHistogramCollection(fGlobalEventSelectionName.Data(),tname->String().Data());
 273
 274     FillHistos(fGlobalEventSelectionName.Data(),tname->String().Data(),centrality.Data(),*aod);
 275
 276 //    if ( hasPileUp )
 277 //    {
 278 //      fEventCounters->Count(Form("event:PILEUP/trigger:%s/run:%d", tname->GetName(), fCurrentRunNumber));
 279 //
 280 //      AssertHistogramCollection("PILEUP",tname->String().Data());
 281 //
 282 //      FillHistos("PILEUP",tname->String().Data(),centrality.Data(),*aod);
 283 //    }
 284 //
 285 //    if ( hasPileUp2 )
 286 //    {
 287 //      fEventCounters->Count(Form("event:PILEUP2/trigger:%s/run:%d", tname->GetName(), fCurrentRunNumber));
 288 //
 289 //      AssertHistogramCollection("PILEUP2",tname->String().Data());
 290 //
 291 //      FillHistos("PILEUP2",tname->String().Data(),centrality.Data(),*aod);
 292 //    }
 293 //
 294   }
 295
 296   delete a;
 297 }
 298
 299 //_____________________________________________________________________________
 300 Bool_t AliAnalysisMuMuFromAOD::TrackChi2(const AliAODTrack& track) const
 301 {
 302   /// Whether the track pass the chi2 cut
 303   return track.Chi2perNDF() < 3.5;
 304 }
 305
 306 //_____________________________________________________________________________
 307 Bool_t AliAnalysisMuMuFromAOD::TrackEtaCut(const AliAODTrack& track) const
 308 {
 309   /// Whether the track pass the eta cut
 310   return track.Eta() < -2.5 && track.Eta() > -4;
 311 }
 312
 313 //_____________________________________________________________________________
 314 Bool_t AliAnalysisMuMuFromAOD::TrackMatchCut(const AliAODTrack& track) const
 315 {
 316   /// Whether the track matched the trigger (any pt)
 317
 318   return track.GetMatchTrigger() > 0;
 319 }
 320
 321 //_____________________________________________________________________________
 322 Bool_t AliAnalysisMuMuFromAOD::TrackMatchLowCut(const AliAODTrack& track) const
 323 {
 324   /// Whether the track matched the trigger (low pt)
 325   return track.GetMatchTrigger() > 1;
 326 }
 327
 328 //_____________________________________________________________________________
 329 Bool_t AliAnalysisMuMuFromAOD::TrackMatchHighCut(const AliAODTrack& track) const
 330 {
 331   /// Whether the track matched the trigger (high pt)
 332
 333   return track.GetMatchTrigger() > 2;
 334 }
 335
 336 //_____________________________________________________________________________
 337 Bool_t AliAnalysisMuMuFromAOD::TrackRabsCut(const AliAODTrack& track) const
 338 {
 339   /// Cut between 2 and 10 degrees
 340
 341   Double_t angle = TMath::ATan(track.GetRAtAbsorberEnd()/AbsZEnd());
 342
 343   return ( angle > Deg2() && angle < Deg10() );
 344 }
 345
 346 //_____________________________________________________________________________
 347 Bool_t AliAnalysisMuMuFromAOD::TrackPtCut(const AliAODTrack& track) const
 348 {
 349   /// Whether the track is above pt cut
 350   return track.Pt() > 1.0;
 351 }
 352
 353 //_____________________________________________________________________________
 354 Bool_t AliAnalysisMuMuFromAOD::TrackBelowPtCut(const AliAODTrack& track) const
 355 {
 356   /// Whether the track is below pt cut
 357   return track.Pt() < 4.0;
 358 }
 359
 360 //_____________________________________________________________________________
 361 Bool_t AliAnalysisMuMuFromAOD::TrackDCACut(const AliAODTrack& track) const
 362 {
 363   /// Whether the track pass the P x DCA cut
 364
 365   if (!fVertex) return kFALSE;
 366
 367   TLorentzVector p(track.Px(),track.Py(),track.Pz(),
 368                    TMath::Sqrt(MuonMass2()+track.P()*track.P()));
 369
 370   Double_t theta = TMath::ATan(track.GetRAtAbsorberEnd()/AbsZEnd());
 371
 372   Double_t meanPcorr = ( theta < ( 180. - 3. ) * TMath::DegToRad() ) ? 1.2 : 1.49;
 373   Double_t pTotMean = p.P() - meanPcorr;
 374
 375   TVector3 eventVertex(fVertex->GetX(),fVertex->GetY(),fVertex->GetZ());
 376
 377   TVector3 trackDcaAtVz(track.XAtDCA(), track.YAtDCA(), fVertex->GetZ());
 378   TVector3 meanDca(-0.46, -0.92, 0.); // LHC10h1
 379   TVector3 dcaAtVz = trackDcaAtVz - eventVertex - meanDca;
 380   Double_t correctedDca = dcaAtVz.Mag(); // it should also be equal to dcaAtVz.Pt().
 381   Double_t cutVariable = pTotMean * correctedDca;
 382
 383   Double_t cutValue = (theta > Deg3()) ? 63. : 120.;
 384   cutValue = TMath::Sqrt(cutValue*cutValue + 0.4*0.4*p.P()*p.P());
 385
 386   return ( cutVariable < 5*cutValue );
 387 }
 388
 389 //_____________________________________________________________________________
 390 UInt_t AliAnalysisMuMuFromAOD::GetTrackMask(const AliAODTrack& track) const
 391 {
 392   /// Get the mask of all the cuts this track pass
 393
 394   UInt_t m(kAll);
 395
 396   if ( TrackRabsCut(track) ) m |= kRabs;
 397
 398   if ( TrackPtCut(track) ) m |= kPt;
 399
 400   if ( TrackMatchCut(track) ) m |= kMatched;
 401
 402   if ( TrackMatchLowCut(track) ) m |= kMatchedLow;
 403
 404   if ( TrackMatchHighCut(track) ) m |= kMatchedHigh;
 405
 406   if ( TrackEtaCut(track) ) m |= kEta;
 407
 408   if ( TrackChi2(track) ) m |= kChi2;
 409
 410   if ( TrackDCACut(track) ) m |= kDCA;
 411
 412   if ( TrackBelowPtCut(track) ) m |= kBelowPt;
 413
 414   return m;
 415 }
 416
 417 //_____________________________________________________________________________
 418 Bool_t AliAnalysisMuMuFromAOD::PairRapidityCut(const AliAODTrack& t1, const AliAODTrack& t2) const
 419 {
 420   /// Whether the pair passes the rapidity cut
 421
 422   TLorentzVector p1(t1.Px(),t1.Py(),t1.Pz(),TMath::Sqrt(MuonMass2()+t1.P()*t1.P()));
 423   TLorentzVector p2(t2.Px(),t2.Py(),t2.Pz(),TMath::Sqrt(MuonMass2()+t2.P()*t2.P()));
 424
 425   TLorentzVector total(p1+p2);
 426
 427   Double_t y = total.Rapidity();
 428
 429   Bool_t ok = ( y < -2.5 && y > -4.0 );
 430
 431   return ok;
 432 }
 433
 434 //_____________________________________________________________________________
 435 void AliAnalysisMuMuFromAOD::GetPairMask(const AliAODTrack& t1, const AliAODTrack& t2,
 436                                          UInt_t& mask1, UInt_t& mask2,
 437                                          UInt_t& mask12) const
 438 {
 439   /// Get the mask of the track pair
 440
 441   mask1 = GetTrackMask(t1);
 442   mask2 = GetTrackMask(t2);
 443
 444   mask12 = mask1 | mask2;
 445
 446   if ( PairRapidityCut(t1,t2) ) mask12 |= kPairRapidity;
 447 }
 448
 449 //_____________________________________________________________________________
 450 void AliAnalysisMuMuFromAOD::FillHistosForTrack(const char* physics,
 451                                                 const char* triggerClassName,
 452                                                 const char* centrality,
 453                                                 const AliAODTrack& track)
 454 {
 455   /// Fill histograms for one track
 456
 457   TLorentzVector p(track.Px(),track.Py(),track.Pz(),
 458                    TMath::Sqrt(MuonMass2()+track.P()*track.P()));
 459
 460
 461   TString charge("Plus");
 462
 463   if ( track.Charge() < 0 ) charge = "Minus";
 464
 465   UInt_t mask = GetTrackMask(track);
 466
 467   Double_t xdca = track.XAtDCA();
 468   Double_t ydca = track.YAtDCA();
 469   Double_t dca = TMath::Sqrt(xdca*xdca+ydca*ydca);
 470   Double_t theta = TMath::ATan(track.GetRAtAbsorberEnd()/AbsZEnd());
 471
 472   TIter next(fSingleTrackCutNames);
 473   TObjString* str;
 474
 475   while ( ( str = static_cast<TObjString*>(next()) ) )
 476   {
 477     Bool_t test = ( ( str->GetUniqueID() & mask ) == str->GetUniqueID() );
 478
 479     if ( test )
 480     {
 481       Histo(physics,triggerClassName,centrality,str->String().Data(),Form("PtEtaMu%s",charge.Data()))->Fill(p.Eta(),p.Pt());
 482
 483       Histo(physics,triggerClassName,centrality,str->String().Data(),Form("PEtaMu%s",charge.Data()))->Fill(p.Eta(),p.P());
 484
 485
 486       Histo(physics,triggerClassName,centrality,str->String().Data(),Form("XYdcaMu%s",charge.Data()))->Fill(ydca,xdca);
 487
 488       Histo(physics,triggerClassName,centrality,str->String().Data(),Form("Chi2Mu%s",charge.Data()))->Fill(track.Chi2perNDF());
 489
 490       if ( theta >= Deg2() && theta < Deg3() )
 491       {
 492         Histo(physics,triggerClassName,centrality,str->String().Data(),Form("dcaP23Mu%s",charge.Data()))->Fill(p.P(),dca);
 493         if ( p.Pt() > 2 )
 494         {
 495           Histo(physics,triggerClassName,centrality,str->String().Data(),Form("dcaPwPtCut23Mu%s",charge.Data()))->Fill(p.P(),dca);
 496         }
 497       }
 498       else if ( theta >= Deg3() && theta < Deg10() )
 499       {
 500         Histo(physics,triggerClassName,centrality,str->String().Data(),Form("dcaP310Mu%s",charge.Data()))->Fill(p.P(),dca);
 501         if ( p.Pt() > 2 )
 502         {
 503           Histo(physics,triggerClassName,centrality,str->String().Data(),Form("dcaPwPtCut310Mu%s",charge.Data()))->Fill(p.P(),dca);
 504         }
 505       }
 506     }
 507   }
 508
 509
 510 }
 511
 512 //_____________________________________________________________________________
 513 void AliAnalysisMuMuFromAOD::FillHistos(const char* physics, const char* triggerClassName,
 514                                         const char* centrality, const AliAODEvent& aod)
 515 {
 516   /// Fill histograms for /physics/triggerClassName/centrality
 517
 518   Histo(physics,triggerClassName,centrality,"BCX")->Fill(1.0*aod.GetBunchCrossNumber());
 519
 520   Histo(physics,triggerClassName,centrality,"Nevents")->Fill(1.0);
 521
 522   if ( fVertex )
 523   {
 524     Histo(physics,triggerClassName,centrality,"Xvertex")->Fill(fVertex->GetX());
 525     Histo(physics,triggerClassName,centrality,"Yvertex")->Fill(fVertex->GetY());
 526     Histo(physics,triggerClassName,centrality,"Zvertex")->Fill(fVertex->GetZ());
 527
 528     Histo(physics,triggerClassName,centrality,"YXvertex")->Fill(fVertex->GetX(),fVertex->GetY());
 529
 530 //    TIter nextV(aod.GetVertices());
 531 //    AliAODVertex* v;
 532 //    while ( ( v = static_cast<AliAODVertex*>(nextV())) )
 533 //    {
 534 //      if ( v->GetType() == AliAODVertex::kPileupSPD )
 535 //      {
 536 //        Histo(physics,triggerClassName,centrality,"PileUpXvertex")->Fill(v->GetX());
 537 //        Histo(physics,triggerClassName,centrality,"PileUpYvertex")->Fill(v->GetY());
 538 //
 539 //        Histo(physics,triggerClassName,centrality,"PileUpZvertex")->Fill(v->GetZ() - fVertex->GetZ());
 540 //
 541 //        Histo(physics,triggerClassName,centrality,"PileUpYXvertex")->Fill(v->GetX(),v->GetY());
 542 //
 543 //      }
 544 //    }
 545   }
 546
 547
 548   // Track loop
 549
 550     for (Int_t i = 0; i < aod.GetNumberOfTracks(); ++i)
 551   {
 552     AliAODTrack* tracki = aod.GetTrack(i);
 553
 554     if (!tracki->IsMuonTrack()) continue;
 555
 556     FillHistosForTrack(physics,triggerClassName,centrality,*tracki);
 557
 558     TLorentzVector pi(tracki->Px(),tracki->Py(),tracki->Pz(),
 559                       TMath::Sqrt(MuonMass2()+tracki->P()*tracki->P()));
 560
 561     for (Int_t j = i+1; j < aod.GetNumberOfTracks(); ++j)
 562     {
 563       AliAODTrack* trackj = aod.GetTrack(j);
 564
 565       if (!trackj->IsMuonTrack()) continue;
 566
 567       TLorentzVector pj(trackj->Px(),trackj->Py(),trackj->Pz(),
 568                        TMath::Sqrt(MuonMass2()+trackj->P()*trackj->P()));
 569
 570       pj += pi;
 571
 572       TIter next(fPairTrackCutNames);
 573       TObjString* str;
 574
 575       UInt_t maski(0),maskj(0),maskij(0);
 576
 577       GetPairMask(*trackj,*tracki,maski,maskj,maskij);
 578
 579       while ( ( str = static_cast<TObjString*>(next()) ) )
 580       {
 581         UInt_t singleTrackMask(0);
 582         UInt_t pairMask(0);
 583
 584         DecodePairCutMask(str->GetUniqueID(),singleTrackMask,pairMask);
 585
 586         Bool_t testi = ( ( maski & singleTrackMask ) == singleTrackMask ) ;
 587         Bool_t testj = ( ( maskj & singleTrackMask ) == singleTrackMask ) ;
 588         Bool_t testij(kTRUE);
 589
 590         if (pairMask>0) testij = ( ( maskij & pairMask ) == pairMask ) ;
 591
 592         if ( testi && testj && testij )
 593         {
 594           if ( tracki->Charge() != trackj->Charge() )
 595           {
 596             Histo(physics,triggerClassName,centrality,str->String().Data(),"MinvUSPt")->Fill(pj.Pt(),pj.M());
 597           }
 598           else if ( tracki->Charge() > 0 && trackj->Charge() > 0 )
 599           {
 600             Histo(physics,triggerClassName,centrality,str->String().Data(),"MinvPPPt")->Fill(pj.Pt(),pj.M());
 601           }
 602           else
 603           {
 604             Histo(physics,triggerClassName,centrality,str->String().Data(),"MinvMMPt")->Fill(pj.Pt(),pj.M());
 605           }
 606         }
 607       }
 608     }
 609   } //track loop
 610 }