PWG/muon/AliAnalysisMuMuFromESD.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 "AliAnalysisMuMuFromESD.h"
  19
  20 #include "AliAnalysisManager.h"
  21 #include "AliInputEventHandler.h"
  22 #include "AliHistogramCollection.h"
  23 #include "AliAODEvent.h"
  24 #include "AliESDEvent.h"
  25 #include "AliESDMuonTrack.h"
  26 #include "AliESDVertex.h"
  27 #include "AliVParticle.h"
  28 #include "Riostream.h"
  29 #include "TCanvas.h"
  30 #include "TDirectory.h"
  31 #include "TH1F.h"
  32 #include "TH2F.h"
  33 #include "TLorentzVector.h"
  34 #include "TMap.h"
  35 #include "TMath.h"
  36 #include "TObjArray.h"
  37 #include "TPaveText.h"
  38 #include <cassert>
  39 #include "AliPhysicsSelection.h"
  40 #include "AliPhysicsSelectionTask.h"
  41 #include "AliBackgroundSelection.h"
  42 #include "AliCounterCollection.h"
  43 #include "THashList.h"
  44
  45 #include "TGrid.h"
  46
  47 ///
  48 /// AliAnalysisMuMuFromESD
  49 ///
  50 /// cuts must be added through AddSingleCut and AddPaircut methods
  51 ///
  52 ///
  53
  54 ClassImp(AliAnalysisMuMuFromESD)
  55
  56 //_____________________________________________________________________________
  57 AliAnalysisMuMuFromESD::AliAnalysisMuMuFromESD()
  58 : AliAnalysisMuMu(), fVertex(0x0)
  59 {
  60   // default ctor
  61 }
  62
  63 //_____________________________________________________________________________
  64 AliAnalysisMuMuFromESD::AliAnalysisMuMuFromESD(Bool_t aa)
  65 :AliAnalysisMuMu(kTRUE,aa), fVertex(0x0)
  66 {
  67   // Constructor
  68   Ctor();
  69 }
  70
  71 //_____________________________________________________________________________
  72 AliAnalysisMuMuFromESD::AliAnalysisMuMuFromESD(TList* triggerClasses)
  73 :AliAnalysisMuMu(kTRUE,triggerClasses), fVertex(0x0)
  74 {
  75   // Constructor
  76   Ctor();
  77 }
  78
  79 //_____________________________________________________________________________
  80 void AliAnalysisMuMuFromESD::Ctor()
  81 {
  82   // Common ctor
  83
  84   if ( fAA )
  85   {
  86     AddGlobalEventSelection("SD2");
  87   }
  88   else
  89   {
  90     AddGlobalEventSelection("ALL");
  91   }
  92
  93   AddGlobalEventSelection("PS");
  94
  95 //  fBranchNames = "ESD:AliESDRun.,AliESDHeader.,MuonTracks";
  96 }
  97
  98 //_____________________________________________________________________________
  99 AliAnalysisMuMuFromESD::~AliAnalysisMuMuFromESD()
 100 {
 101   /// dtor
 102   delete fVertex;
 103 }
 104
 105 //_____________________________________________________________________________
 106 void AliAnalysisMuMuFromESD::MuUserExec(Option_t*)
 107 {
 108   // Main loop
 109   // Called for each event
 110
 111   delete fVertex;
 112   fVertex = 0x0;
 113
 114   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
 115
 116   if (esd)
 117   {
 118
 119     if ( fAA )
 120     {
 121     // consider only events with OSM2 fired
 122
 123       UInt_t sd2 = (1<<12);
 124       UInt_t trigger = esd->GetHeader()->GetL0TriggerInputs();
 125
 126       Bool_t ok = ( ( trigger & sd2 ) == sd2);
 127
 128
 129       if (!ok) return;
 130
 131     }
 132
 133     TString runNumber(RunNumber(*esd));
 134
 135     TString triggers = esd->GetFiredTriggerClasses();
 136
 137     if ( IsDynamicTriggerClasses() ) AddTriggerClasses(triggers.Data());
 138
 139     TObjArray* a = triggers.Tokenize(" ");
 140     TIter next(a);
 141     TObjString* tname;
 142
 143     AliInputEventHandler* inputEventHandler = static_cast<AliInputEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 144
 145     Bool_t isPhysicsSelected =  ( inputEventHandler->IsEventSelected() & AliVEvent::kINT7 )
 146     || ( inputEventHandler->IsEventSelected() & AliVEvent::kMUSH7 )
 147     || ( inputEventHandler->IsEventSelected() & AliVEvent::kMUL7 )
 148     || ( inputEventHandler->IsEventSelected() & AliVEvent::kMUS7 )
 149     || ( inputEventHandler->IsEventSelected() & AliVEvent::kMUU7 );
 150
 151     TString centrality(DefaultCentralityName());
 152
 153     fVertex = new AliESDVertex(*esd->GetPrimaryVertexSPD());
 154
 155     if ( fAA )
 156     {
 157       // FIXME : should get it from centrality task...
 158
 159       // Double_t cent(-100);
 160       //    Double_t v0mult = vzero->GetMTotV0A()+vzero->GetMTotV0C();
 161
 162       //    for ( std::vector<double>::size_type i = 0 ; i < fCentralityLimits.size()-1 && cent < -10 ; ++i )
 163       //    {
 164       //      if ( v0mult > fCentralityAmplitude[i] )
 165       //      {
 166       //        cent = fCentralityLimits[i];
 167       //      }
 168       //    }
 169       //
 170       // if ( cent > -1 )
 171       // {
 172       //   centrality = CentralityName(cent);
 173       // }
 174     }
 175
 176     TString eventtype("all");
 177     TString et;
 178
 179     if ( fAA ) eventtype="sd2";
 180     et = eventtype;
 181     et.ToUpper();
 182
 183     while ( ( tname = static_cast<TObjString*>(next()) ) )
 184     {
 185       if ( !fTriggerClasses->FindObject(tname->GetName()) )
 186       {
 187         // skip the trigger classes we're not supposed to analyse
 188
 189         continue;
 190       }
 191
 192       fEventCounters->Count(Form("event:%s/trigger:%s/run:%d", eventtype.Data(),tname->GetName(), runNumber.Atoi()));
 193
 194       AssertHistogramCollection(et.Data(),tname->String().Data());
 195
 196       FillHistos(et.Data(),tname->String().Data(),centrality.Data(),*esd);
 197
 198       if ( isPhysicsSelected )
 199       {
 200
 201         fEventCounters->Count(Form("event:ps/trigger:%s/run:%d", tname->GetName(), runNumber.Atoi()));
 202
 203         AssertHistogramCollection("PS",tname->String().Data());
 204
 205         FillHistos("PS",tname->String().Data(),centrality.Data(),*esd);
 206       }
 207     }
 208
 209     delete a;
 210   }
 211 }
 212
 213 //_____________________________________________________________________________
 214 const char*
 215 AliAnalysisMuMuFromESD::RunNumber(const AliESDEvent& esd) const
 216 {
 217   /// Get the current runnumber
 218   return Form("%09d",esd.GetRunNumber());
 219 }
 220
 221 //_____________________________________________________________________________
 222 Bool_t AliAnalysisMuMuFromESD::TrackMatchCut(const AliESDMuonTrack& track) const
 223 {
 224   /// whether the track match the trigger (any pt)
 225   return track.GetMatchTrigger() > 0;
 226 }
 227
 228 //_____________________________________________________________________________
 229 Bool_t AliAnalysisMuMuFromESD::TrackMatchLowCut(const AliESDMuonTrack& track) const
 230 {
 231   /// whether the track match the trigger (low pt)
 232   return track.GetMatchTrigger() > 1;
 233 }
 234
 235 //_____________________________________________________________________________
 236 Bool_t AliAnalysisMuMuFromESD::TrackMatchHighCut(const AliESDMuonTrack& track) const
 237 {
 238   /// whether the track match the trigger (high pt)
 239   return track.GetMatchTrigger() > 2;
 240 }
 241
 242 //_____________________________________________________________________________
 243 Bool_t AliAnalysisMuMuFromESD::TrackRabsCut(const AliESDMuonTrack& track) const
 244 {
 245   /// Cut between 2 and 9 degrees
 246
 247   Double_t angle = TMath::ATan(track.GetRAtAbsorberEnd()/AbsZEnd());
 248
 249   return ( angle > Deg2() && angle < Deg10() );
 250 }
 251
 252 //_____________________________________________________________________________
 253 Bool_t AliAnalysisMuMuFromESD::TrackPtCut(const AliESDMuonTrack& track) const
 254 {
 255   /// whether the track pass the pt cut
 256   return track.Pt() > 1.0;
 257 }
 258
 259 //_____________________________________________________________________________
 260 Bool_t AliAnalysisMuMuFromESD::TrackChi2(const AliESDMuonTrack& track) const
 261 {
 262   /// whether the track pass the chi2 cut
 263   return track.GetNormalizedChi2() < 4.0;
 264 }
 265
 266 //_____________________________________________________________________________
 267 Double_t AliAnalysisMuMuFromESD::CorrectedDCA(const AliESDMuonTrack& track) const
 268 {
 269   /// Get corrected DCA of the track
 270
 271   if (!fVertex) return -999.;
 272
 273   TVector3 eventVertex(fVertex->GetX(),fVertex->GetY(),fVertex->GetZ());
 274
 275   TVector3 trackDcaAtVz(track.GetNonBendingCoorAtDCA(), track.GetBendingCoorAtDCA(), fVertex->GetZ());
 276
 277   TVector3 meanDca(-0.46, -0.92, 0.); // LHC10h1
 278   TVector3 dcaAtVz = trackDcaAtVz - eventVertex - meanDca;
 279   return dcaAtVz.Mag(); // it should also be equal to dcaAtVz.Pt().
 280 }
 281
 282 //_____________________________________________________________________________
 283 Double_t AliAnalysisMuMuFromESD::PDCACutValue(const AliESDMuonTrack& track) const
 284 {
 285   /// Get P x DCA cut value of the track
 286
 287   TVector3 pTot(track.Px(),track.Py(),track.Pz());
 288
 289   Double_t theta = TMath::ATan(track.GetRAtAbsorberEnd()/AbsZEnd());
 290   Double_t cutValue = (theta > Deg3()) ? 63. : 120.;
 291
 292 //  AliInfo(Form("theta %e cutValue %e",theta*TMath::RadToDeg(),cutValue));
 293
 294   return 5*TMath::Sqrt(cutValue*cutValue + 0.4*0.4*pTot.Mag2());
 295 }
 296
 297 //_____________________________________________________________________________
 298 Bool_t AliAnalysisMuMuFromESD::TrackDCACut(const AliESDMuonTrack& track) const
 299 {
 300   /// whether the track pass the P x DCA cut
 301
 302   if (!fVertex) return kFALSE;
 303
 304   TVector3 pTot(track.Px(),track.Py(),track.Pz());
 305
 306   TVector3 pTotUncorr(track.PxUncorrected(), track.PyUncorrected(), track.PzUncorrected());
 307
 308   Double_t pTotMean = 0.5 * ( pTot.Mag() + pTotUncorr.Mag() );
 309
 310   Double_t correctedDca = CorrectedDCA(track);
 311
 312   Double_t cutVariable = pTotMean * correctedDca;
 313
 314   Double_t cutValue = PDCACutValue(track);
 315
 316   return ( cutVariable < cutValue );
 317 }
 318
 319 //_____________________________________________________________________________
 320 Bool_t AliAnalysisMuMuFromESD::TrackEtaCut(const AliESDMuonTrack& track) const
 321 {
 322   /// whether the track passes the eta cut
 323   return track.Eta() < -2.5 && track.Eta() > -4;
 324 }
 325
 326 //_____________________________________________________________________________
 327 UInt_t AliAnalysisMuMuFromESD::GetTrackMask(const AliESDMuonTrack& track) const
 328 {
 329   /// Get the mask of all cuts this track pass
 330
 331   UInt_t m(kAll);
 332
 333   if ( TrackRabsCut(track) ) m |= kRabs;
 334
 335   if ( TrackPtCut(track) ) m |= kPt;
 336
 337   if ( TrackMatchCut(track) ) m |= kMatched;
 338
 339   if ( TrackMatchLowCut(track) ) m |= kMatchedLow;
 340
 341   if ( TrackMatchHighCut(track) ) m |= kMatchedHigh;
 342
 343   if ( TrackEtaCut(track) ) m |= kEta;
 344
 345   if ( TrackChi2(track) ) m |= kChi2;
 346
 347   if ( TrackDCACut(track) ) m |= kDCA;
 348
 349   return m;
 350
 351 }
 352
 353 //_____________________________________________________________________________
 354 void AliAnalysisMuMuFromESD::FillHistosForTrack(const char* physics,
 355                                                 const char* triggerClassName,
 356                                                 const char* centrality,
 357                                                 const AliESDMuonTrack& track,
 358                                                 const char* /*runNumber*/)
 359 {
 360   /// Fill histograms for one track
 361
 362   if ( track.ContainTrackerData() )
 363   {
 364     TLorentzVector p;
 365
 366     track.LorentzP(p);
 367
 368     TString charge("Plus");
 369
 370     if ( track.Charge() < 0 ) charge = "Minus";
 371
 372     UInt_t mask = GetTrackMask(track);
 373
 374     Double_t theta = TMath::ATan(track.GetRAtAbsorberEnd()/AbsZEnd());
 375
 376     Double_t xdca = track.GetNonBendingCoorAtDCA();
 377     Double_t ydca = track.GetBendingCoorAtDCA();
 378     Double_t dca = TMath::Sqrt(xdca*xdca+ydca*ydca);
 379
 380     TIter next(fSingleTrackCutNames);
 381     TObjString* str;
 382
 383     while ( ( str = static_cast<TObjString*>(next()) ) )
 384     {
 385       Bool_t test = ( ( str->GetUniqueID() & mask ) == str->GetUniqueID() );
 386
 387       if ( test )
 388       {
 389
 390         TVector3 pTot(track.Px(),track.Py(),track.Pz());
 391
 392         TVector3 pTotUncorr(track.PxUncorrected(), track.PyUncorrected(), track.PzUncorrected());
 393
 394         Double_t pTotMean = 0.5 * ( pTot.Mag() + pTotUncorr.Mag() );
 395
 396         Double_t pdcacorr = pTotMean*CorrectedDCA(track);
 397         Double_t pdcacut = PDCACutValue(track);
 398
 399         Histo(physics,triggerClassName,centrality,str->String().Data(),Form("PtEtaMu%s",charge.Data()))->Fill(p.Eta(),p.Pt());
 400
 401         Histo(physics,triggerClassName,centrality,str->String().Data(),Form("PEtaMu%s",charge.Data()))->Fill(p.Eta(),p.P());
 402
 403
 404         Histo(physics,triggerClassName,centrality,str->String().Data(),Form("XYdcaMu%s",charge.Data()))->Fill(ydca,xdca);
 405
 406         Histo(physics,triggerClassName,centrality,str->String().Data(),Form("Chi2Mu%s",charge.Data()))->Fill(track.GetNormalizedChi2());
 407
 408         if ( theta >= Deg2() && theta < Deg3() )
 409         {
 410           //        Histo(physics,triggerClassName,centrality,str->String().Data(),Form("PDCAP23Mu%s",charge.Data()))->Fill(p.P(),p.P()*dca);
 411           Histo(physics,triggerClassName,centrality,str->String().Data(),Form("dcaP23Mu%s",charge.Data()))->Fill(p.P(),dca);
 412
 413           Histo(physics,triggerClassName,centrality,str->String().Data(),Form("PDCAcorrP23Mu%s",charge.Data()))->Fill(pTotMean,pdcacorr);
 414           Histo(physics,triggerClassName,centrality,str->String().Data(),Form("PDCAcutP23Mu%s",charge.Data()))->Fill(p.P(),pdcacut);
 415
 416           if ( p.Pt() > 2 )
 417           {
 418             Histo(physics,triggerClassName,centrality,str->String().Data(),Form("dcaPwPtCut23Mu%s",charge.Data()))->Fill(p.P(),dca);
 419           }
 420         }
 421         else if ( theta >= Deg3() && theta < Deg10() )
 422         {
 423           //        Histo(physics,triggerClassName,centrality,str->String().Data(),Form("PDCAP310Mu%s",charge.Data()))->Fill(p.P(),p.P()*dca);
 424           Histo(physics,triggerClassName,centrality,str->String().Data(),Form("dcaP310Mu%s",charge.Data()))->Fill(p.P(),dca);
 425           Histo(physics,triggerClassName,centrality,str->String().Data(),Form("PDCAcorrP310Mu%s",charge.Data()))->Fill(pTotMean,pdcacorr);
 426           Histo(physics,triggerClassName,centrality,str->String().Data(),Form("PDCAcutP310Mu%s",charge.Data()))->Fill(p.P(),pdcacut);
 427           if ( p.Pt() > 2 )
 428           {
 429             Histo(physics,triggerClassName,centrality,str->String().Data(),Form("dcaPwPtCut310Mu%s",charge.Data()))->Fill(p.P(),dca);
 430           }
 431         }
 432       }
 433     }
 434   }
 435 }
 436
 437 //_____________________________________________________________________________
 438 void AliAnalysisMuMuFromESD::FillHistos(const char* physics,
 439                                         const char* triggerClassName,
 440                                         const char* centrality,
 441                                         AliESDEvent& esd)
 442 {
 443   /// Fill histograms for /physics/triggerClassName/centrality
 444
 445   Histo(physics,triggerClassName,centrality,"BCX")->Fill(1.0*esd.GetBunchCrossNumber());
 446
 447   if ( fVertex )
 448   {
 449     Histo(physics,triggerClassName,centrality,"Xvertex")->Fill(fVertex->GetX());
 450     Histo(physics,triggerClassName,centrality,"Yvertex")->Fill(fVertex->GetY());
 451     Histo(physics,triggerClassName,centrality,"Zvertex")->Fill(fVertex->GetZ());
 452   }
 453
 454   AliESDVZERO* vzero = esd.GetVZEROData();
 455
 456   if ( vzero )
 457   {
 458     Histo(physics,triggerClassName,centrality,"V0Time")->Fill(vzero->GetV0ATime(),vzero->GetV0CTime());
 459
 460     Double_t v0mult = vzero->GetMTotV0A()+vzero->GetMTotV0C();
 461
 462     Histo(physics,triggerClassName,centrality,"V0Amplitude")->Fill(v0mult);
 463   }
 464
 465   // Track loop
 466
 467   TString runNumber(RunNumber(esd));
 468
 469   for (Int_t i = 0; i < esd.GetNumberOfMuonTracks(); ++i)
 470   {
 471     AliESDMuonTrack* tracki = esd.GetMuonTrack(i);
 472
 473     if ( tracki->ContainTrackerData() )
 474     {
 475       UInt_t maski = GetTrackMask(*tracki);
 476
 477       FillHistosForTrack(physics,triggerClassName,centrality,*tracki,runNumber.Data());
 478
 479       TLorentzVector pi;
 480
 481       tracki->LorentzP(pi);
 482
 483       for (Int_t j = i+1; j < esd.GetNumberOfMuonTracks(); ++j)
 484       {
 485         AliESDMuonTrack* trackj = esd.GetMuonTrack(j);
 486
 487         if ( trackj->ContainTrackerData() )
 488         {
 489           UInt_t maskj = GetTrackMask(*trackj);
 490
 491           TLorentzVector pj;
 492           trackj->LorentzP(pj);
 493
 494           pj += pi;
 495
 496           TIter next(fPairTrackCutNames);
 497           TObjString* str;
 498
 499           while ( ( str = static_cast<TObjString*>(next()) ) )
 500           {
 501             UInt_t singleTrackMask(0);
 502             UInt_t pairMask(0);
 503
 504             DecodePairCutMask(str->GetUniqueID(),singleTrackMask,pairMask);
 505
 506             Bool_t testi = ( ( maski & singleTrackMask ) == singleTrackMask ) ;
 507             Bool_t testj = ( ( maskj & singleTrackMask ) == singleTrackMask ) ;
 508             Bool_t testij(kTRUE);
 509
 510             if ( pairMask > 0 )
 511             {
 512               testij = ( (maski & pairMask) == pairMask ) &&
 513               ( (maskj & pairMask) == pairMask );
 514             }
 515
 516             if ( ( testi || testj ) && testij )
 517             {
 518               if ( tracki->Charge() != trackj->Charge() )
 519               {
 520                 Histo(physics,triggerClassName,centrality,str->String().Data(),"MinvUSPt")->Fill(pj.Pt(),pj.M());
 521               }
 522               else if ( tracki->Charge() > 0 && trackj->Charge() > 0 )
 523               {
 524                 Histo(physics,triggerClassName,centrality,str->String().Data(),"MinvPPPt")->Fill(pj.Pt(),pj.M());
 525               }
 526               else
 527               {
 528                 Histo(physics,triggerClassName,centrality,str->String().Data(),"MinvMMPt")->Fill(pj.Pt(),pj.M());
 529               }
 530             }
 531           }
 532         }
 533       }
 534     }
 535   } //track loop
 536 }
 537