[u/mrichter/AliRoot.git] / PWG / muon / AliAnalysisMuMuFromESD.cxx

/**************************************************************************
* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
*                                                                        *
* Author: The ALICE Off-line Project.                                    *
* Contributors are mentioned in the code where appropriate.              *
*                                                                        *
* Permission to use, copy, modify and distribute this software and its   *
* documentation strictly for non-commercial purposes is hereby granted   *
* without fee, provided that the above copyright notice appears in all   *
* copies and that both the copyright notice and this permission notice   *
* appear in the supporting documentation. The authors make no claims     *
* about the suitability of this software for any purpose. It is          *
* provided "as is" without express or implied warranty.                  *
**************************************************************************/

// $Id$

#include "AliAnalysisMuMuFromESD.h"

#include "AliAnalysisManager.h"
#include "AliInputEventHandler.h"
#include "AliHistogramCollection.h"
#include "AliAODEvent.h"
#include "AliESDEvent.h"
#include "AliESDMuonTrack.h"
#include "AliESDVertex.h"
#include "AliVParticle.h"
#include "Riostream.h"
#include "TCanvas.h"
#include "TDirectory.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TLorentzVector.h"
#include "TMap.h"
#include "TMath.h"
#include "TObjArray.h"
#include "TPaveText.h"
#include <cassert>
#include "AliPhysicsSelection.h"
#include "AliPhysicsSelectionTask.h"
#include "AliBackgroundSelection.h"
#include "AliCounterCollection.h"
#include "THashList.h"

#include "TGrid.h"

/// 
/// AliAnalysisMuMuFromESD
///
/// cuts must be added through AddSingleCut and AddPaircut methods
///
/// 

ClassImp(AliAnalysisMuMuFromESD)

//_____________________________________________________________________________
AliAnalysisMuMuFromESD::AliAnalysisMuMuFromESD() 
: AliAnalysisMuMu(), fVertex(0x0)
{
  // default ctor
}

//_____________________________________________________________________________
AliAnalysisMuMuFromESD::AliAnalysisMuMuFromESD(Bool_t aa)
:AliAnalysisMuMu(kTRUE,aa), fVertex(0x0)
{
  // Constructor
  Ctor();
}

//_____________________________________________________________________________
AliAnalysisMuMuFromESD::AliAnalysisMuMuFromESD(TList* triggerClasses) 
:AliAnalysisMuMu(kTRUE,triggerClasses), fVertex(0x0)
{
  // Constructor
  Ctor();
}

//_____________________________________________________________________________
void AliAnalysisMuMuFromESD::Ctor()
{
  // Common ctor

  if ( fAA ) 
  {
    AddGlobalEventSelection("SD2");
  }
  else
  {
    AddGlobalEventSelection("ALL");    
  }
  
  AddGlobalEventSelection("PS");
  
//  fBranchNames = "ESD:AliESDRun.,AliESDHeader.,MuonTracks";
}

//_____________________________________________________________________________
AliAnalysisMuMuFromESD::~AliAnalysisMuMuFromESD()
{
  /// dtor
  delete fVertex;
}

//_____________________________________________________________________________
void AliAnalysisMuMuFromESD::MuUserExec(Option_t*)
{
  // Main loop
  // Called for each event
  
  delete fVertex;
  fVertex = 0x0;
  
  AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
  
  if (esd)
  {      

    if ( fAA ) 
    {
    // consider only events with OSM2 fired
    
      UInt_t sd2 = (1<<12);
      UInt_t trigger = esd->GetHeader()->GetL0TriggerInputs();
      
      Bool_t ok = ( ( trigger & sd2 ) == sd2);


      if (!ok) return;

    }
        
    TString runNumber(RunNumber(*esd));
  
    TString triggers = esd->GetFiredTriggerClasses();
    
    if ( IsDynamicTriggerClasses() ) AddTriggerClasses(triggers.Data());

    TObjArray* a = triggers.Tokenize(" ");
    TIter next(a);
    TObjString* tname;
    
    AliInputEventHandler* inputEventHandler = static_cast<AliInputEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
    
    Bool_t isPhysicsSelected =  ( inputEventHandler->IsEventSelected() & AliVEvent::kINT7 )
    || ( inputEventHandler->IsEventSelected() & AliVEvent::kMUSH7 )
    || ( inputEventHandler->IsEventSelected() & AliVEvent::kMUL7 )
    || ( inputEventHandler->IsEventSelected() & AliVEvent::kMUS7 )
    || ( inputEventHandler->IsEventSelected() & AliVEvent::kMUU7 );

    TString centrality(DefaultCentralityName());
    
    fVertex = new AliESDVertex(*esd->GetPrimaryVertexSPD());

    if ( fAA ) 
    {
      // FIXME : should get it from centrality task...
      
      // Double_t cent(-100);
      //    Double_t v0mult = vzero->GetMTotV0A()+vzero->GetMTotV0C();
      
      //    for ( std::vector<double>::size_type i = 0 ; i < fCentralityLimits.size()-1 && cent < -10 ; ++i )
      //    {
      //      if ( v0mult > fCentralityAmplitude[i] ) 
      //      {
      //        cent = fCentralityLimits[i];
      //      }
      //    }
      //    
      // if ( cent > -1 ) 
      // {
      //   centrality = CentralityName(cent);
      // }
    }

    TString eventtype("all");
    TString et;
    
    if ( fAA ) eventtype="sd2";
    et = eventtype;
    et.ToUpper();
    
    while ( ( tname = static_cast<TObjString*>(next()) ) )
    {
      if ( !fTriggerClasses->FindObject(tname->GetName()) ) 
      {
        // skip the trigger classes we're not supposed to analyse

        continue;
      }

      fEventCounters->Count(Form("event:%s/trigger:%s/run:%d", eventtype.Data(),tname->GetName(), runNumber.Atoi()));
        
      AssertHistogramCollection(et.Data(),tname->String().Data());
        
      FillHistos(et.Data(),tname->String().Data(),centrality.Data(),*esd);    
      
      if ( isPhysicsSelected )
      {
        
        fEventCounters->Count(Form("event:ps/trigger:%s/run:%d", tname->GetName(), runNumber.Atoi()));
        
        AssertHistogramCollection("PS",tname->String().Data());
        
        FillHistos("PS",tname->String().Data(),centrality.Data(),*esd);            
      }
    }
    
    delete a;
  }
}

//_____________________________________________________________________________
const char*
AliAnalysisMuMuFromESD::RunNumber(const AliESDEvent& esd) const
{
  /// Get the current runnumber
  return Form("%09d",esd.GetRunNumber());
}

//_____________________________________________________________________________
Bool_t AliAnalysisMuMuFromESD::TrackMatchCut(const AliESDMuonTrack& track) const
{
  /// whether the track match the trigger (any pt)
  return track.GetMatchTrigger() > 0;
}

//_____________________________________________________________________________
Bool_t AliAnalysisMuMuFromESD::TrackMatchLowCut(const AliESDMuonTrack& track) const
{
  /// whether the track match the trigger (low pt)
  return track.GetMatchTrigger() > 1;
}

//_____________________________________________________________________________
Bool_t AliAnalysisMuMuFromESD::TrackMatchHighCut(const AliESDMuonTrack& track) const
{
  /// whether the track match the trigger (high pt)
  return track.GetMatchTrigger() > 2;
}

//_____________________________________________________________________________
Bool_t AliAnalysisMuMuFromESD::TrackRabsCut(const AliESDMuonTrack& track) const
{
  /// Cut between 2 and 9 degrees
  
  Double_t angle = TMath::ATan(track.GetRAtAbsorberEnd()/AbsZEnd());
  
  return ( angle > Deg2() && angle < Deg10() );
}

//_____________________________________________________________________________
Bool_t AliAnalysisMuMuFromESD::TrackPtCut(const AliESDMuonTrack& track) const
{
  /// whether the track pass the pt cut
  return track.Pt() > 1.0;
}

//_____________________________________________________________________________
Bool_t AliAnalysisMuMuFromESD::TrackChi2(const AliESDMuonTrack& track) const
{
  /// whether the track pass the chi2 cut
  return track.GetNormalizedChi2() < 4.0;
}

//_____________________________________________________________________________
Double_t AliAnalysisMuMuFromESD::CorrectedDCA(const AliESDMuonTrack& track) const
{
  /// Get corrected DCA of the track

  if (!fVertex) return -999.;
  
  TVector3 eventVertex(fVertex->GetX(),fVertex->GetY(),fVertex->GetZ());
  
  TVector3 trackDcaAtVz(track.GetNonBendingCoorAtDCA(), track.GetBendingCoorAtDCA(), fVertex->GetZ());
  
  TVector3 meanDca(-0.46, -0.92, 0.); // LHC10h1
  TVector3 dcaAtVz = trackDcaAtVz - eventVertex - meanDca;
  return dcaAtVz.Mag(); // it should also be equal to dcaAtVz.Pt().
}

//_____________________________________________________________________________
Double_t AliAnalysisMuMuFromESD::PDCACutValue(const AliESDMuonTrack& track) const
{
  /// Get P x DCA cut value of the track

  TVector3 pTot(track.Px(),track.Py(),track.Pz());

  Double_t theta = TMath::ATan(track.GetRAtAbsorberEnd()/AbsZEnd());  
  Double_t cutValue = (theta > Deg3()) ? 63. : 120.;
  
//  AliInfo(Form("theta %e cutValue %e",theta*TMath::RadToDeg(),cutValue));
  
  return 5*TMath::Sqrt(cutValue*cutValue + 0.4*0.4*pTot.Mag2());
}

//_____________________________________________________________________________
Bool_t AliAnalysisMuMuFromESD::TrackDCACut(const AliESDMuonTrack& track) const
{
  /// whether the track pass the P x DCA cut
  
  if (!fVertex) return kFALSE;
  
  TVector3 pTot(track.Px(),track.Py(),track.Pz());
  
  TVector3 pTotUncorr(track.PxUncorrected(), track.PyUncorrected(), track.PzUncorrected());

  Double_t pTotMean = 0.5 * ( pTot.Mag() + pTotUncorr.Mag() );
  
  Double_t correctedDca = CorrectedDCA(track);  
  
  Double_t cutVariable = pTotMean * correctedDca;
  
  Double_t cutValue = PDCACutValue(track);
  
  return ( cutVariable < cutValue );
}

//_____________________________________________________________________________
Bool_t AliAnalysisMuMuFromESD::TrackEtaCut(const AliESDMuonTrack& track) const
{
  /// whether the track passes the eta cut
  return track.Eta() < -2.5 && track.Eta() > -4;
}

//_____________________________________________________________________________
UInt_t AliAnalysisMuMuFromESD::GetTrackMask(const AliESDMuonTrack& track) const
{
  /// Get the mask of all cuts this track pass
  
  UInt_t m(kAll);
  
  if ( TrackRabsCut(track) ) m |= kRabs;
  
  if ( TrackPtCut(track) ) m |= kPt;  
  
  if ( TrackMatchCut(track) ) m |= kMatched;
  
  if ( TrackMatchLowCut(track) ) m |= kMatchedLow;
  
  if ( TrackMatchHighCut(track) ) m |= kMatchedHigh;
  
  if ( TrackEtaCut(track) ) m |= kEta;
  
  if ( TrackChi2(track) ) m |= kChi2;
  
  if ( TrackDCACut(track) ) m |= kDCA;
  
  return m;
  
}

//_____________________________________________________________________________
void AliAnalysisMuMuFromESD::FillHistosForTrack(const char* physics,
                                                const char* triggerClassName, 
                                                const char* centrality,
                                                const AliESDMuonTrack& track,
                                                const char* /*runNumber*/)
{
  /// Fill histograms for one track
  
  if ( track.ContainTrackerData() )
  {      
    TLorentzVector p;
    
    track.LorentzP(p);
    
    TString charge("Plus");
    
    if ( track.Charge() < 0 ) charge = "Minus";

    UInt_t mask = GetTrackMask(track);
    
    Double_t theta = TMath::ATan(track.GetRAtAbsorberEnd()/AbsZEnd());

    Double_t xdca = track.GetNonBendingCoorAtDCA();
    Double_t ydca = track.GetBendingCoorAtDCA();
    Double_t dca = TMath::Sqrt(xdca*xdca+ydca*ydca);
    
    TIter next(fSingleTrackCutNames);
    TObjString* str;
    
    while ( ( str = static_cast<TObjString*>(next()) ) )
    {
      Bool_t test = ( ( str->GetUniqueID() & mask ) == str->GetUniqueID() );

      if ( test ) 
      {
       
        TVector3 pTot(track.Px(),track.Py(),track.Pz());
        
        TVector3 pTotUncorr(track.PxUncorrected(), track.PyUncorrected(), track.PzUncorrected());
        
        Double_t pTotMean = 0.5 * ( pTot.Mag() + pTotUncorr.Mag() );
        
        Double_t pdcacorr = pTotMean*CorrectedDCA(track);
        Double_t pdcacut = PDCACutValue(track);
        
        Histo(physics,triggerClassName,centrality,str->String().Data(),Form("PtEtaMu%s",charge.Data()))->Fill(p.Eta(),p.Pt());
        
        Histo(physics,triggerClassName,centrality,str->String().Data(),Form("PEtaMu%s",charge.Data()))->Fill(p.Eta(),p.P());
        
        
        Histo(physics,triggerClassName,centrality,str->String().Data(),Form("XYdcaMu%s",charge.Data()))->Fill(ydca,xdca);
        
        Histo(physics,triggerClassName,centrality,str->String().Data(),Form("Chi2Mu%s",charge.Data()))->Fill(track.GetNormalizedChi2());
        
        if ( theta >= Deg2() && theta < Deg3() )         
        {
          //        Histo(physics,triggerClassName,centrality,str->String().Data(),Form("PDCAP23Mu%s",charge.Data()))->Fill(p.P(),p.P()*dca);
          Histo(physics,triggerClassName,centrality,str->String().Data(),Form("dcaP23Mu%s",charge.Data()))->Fill(p.P(),dca);
          
          Histo(physics,triggerClassName,centrality,str->String().Data(),Form("PDCAcorrP23Mu%s",charge.Data()))->Fill(pTotMean,pdcacorr);
          Histo(physics,triggerClassName,centrality,str->String().Data(),Form("PDCAcutP23Mu%s",charge.Data()))->Fill(p.P(),pdcacut);
          
          if ( p.Pt() > 2 )
          {
            Histo(physics,triggerClassName,centrality,str->String().Data(),Form("dcaPwPtCut23Mu%s",charge.Data()))->Fill(p.P(),dca);
          }
        }
        else if ( theta >= Deg3() && theta < Deg10() )
        {
          //        Histo(physics,triggerClassName,centrality,str->String().Data(),Form("PDCAP310Mu%s",charge.Data()))->Fill(p.P(),p.P()*dca);        
          Histo(physics,triggerClassName,centrality,str->String().Data(),Form("dcaP310Mu%s",charge.Data()))->Fill(p.P(),dca);
          Histo(physics,triggerClassName,centrality,str->String().Data(),Form("PDCAcorrP310Mu%s",charge.Data()))->Fill(pTotMean,pdcacorr);
          Histo(physics,triggerClassName,centrality,str->String().Data(),Form("PDCAcutP310Mu%s",charge.Data()))->Fill(p.P(),pdcacut);
          if ( p.Pt() > 2 )
          {
            Histo(physics,triggerClassName,centrality,str->String().Data(),Form("dcaPwPtCut310Mu%s",charge.Data()))->Fill(p.P(),dca);
          }
        }        
      }
    }
  }
}

//_____________________________________________________________________________
void AliAnalysisMuMuFromESD::FillHistos(const char* physics, 
                                        const char* triggerClassName, 
                                        const char* centrality,
                                        AliESDEvent& esd)
{
  /// Fill histograms for /physics/triggerClassName/centrality
  
  Histo(physics,triggerClassName,centrality,"BCX")->Fill(1.0*esd.GetBunchCrossNumber());

  if ( fVertex ) 
  {
    Histo(physics,triggerClassName,centrality,"Xvertex")->Fill(fVertex->GetX());
    Histo(physics,triggerClassName,centrality,"Yvertex")->Fill(fVertex->GetY());
    Histo(physics,triggerClassName,centrality,"Zvertex")->Fill(fVertex->GetZ());
  }
  
  AliESDVZERO* vzero = esd.GetVZEROData();
  
  if ( vzero ) 
  {
    Histo(physics,triggerClassName,centrality,"V0Time")->Fill(vzero->GetV0ATime(),vzero->GetV0CTime());
    
    Double_t v0mult = vzero->GetMTotV0A()+vzero->GetMTotV0C();
    
    Histo(physics,triggerClassName,centrality,"V0Amplitude")->Fill(v0mult);    
  }
  
  // Track loop
  
  TString runNumber(RunNumber(esd));
  
  for (Int_t i = 0; i < esd.GetNumberOfMuonTracks(); ++i) 
  {
    AliESDMuonTrack* tracki = esd.GetMuonTrack(i);

    if ( tracki->ContainTrackerData() ) 
    {
      UInt_t maski = GetTrackMask(*tracki);
    
      FillHistosForTrack(physics,triggerClassName,centrality,*tracki,runNumber.Data());
    
      TLorentzVector pi;

      tracki->LorentzP(pi);
    
      for (Int_t j = i+1; j < esd.GetNumberOfMuonTracks(); ++j) 
      {
        AliESDMuonTrack* trackj = esd.GetMuonTrack(j);
        
        if ( trackj->ContainTrackerData() ) 
        {
          UInt_t maskj = GetTrackMask(*trackj);

          TLorentzVector pj;
          trackj->LorentzP(pj);

          pj += pi;
        
          TIter next(fPairTrackCutNames);
          TObjString* str;
          
          while ( ( str = static_cast<TObjString*>(next()) ) )
          {
            UInt_t singleTrackMask(0);
            UInt_t pairMask(0);
            
            DecodePairCutMask(str->GetUniqueID(),singleTrackMask,pairMask);
            
            Bool_t testi = ( ( maski & singleTrackMask ) == singleTrackMask ) ;
            Bool_t testj = ( ( maskj & singleTrackMask ) == singleTrackMask ) ;
            Bool_t testij(kTRUE);
            
            if ( pairMask > 0 )
            {
              testij = ( (maski & pairMask) == pairMask ) &&
              ( (maskj & pairMask) == pairMask );
            }
            
            if ( ( testi || testj ) && testij )
            {
              if ( tracki->Charge() != trackj->Charge() )
              {
                Histo(physics,triggerClassName,centrality,str->String().Data(),"MinvUSPt")->Fill(pj.Pt(),pj.M());            
              }
              else if ( tracki->Charge() > 0 && trackj->Charge() > 0 )
              {
                Histo(physics,triggerClassName,centrality,str->String().Data(),"MinvPPPt")->Fill(pj.Pt(),pj.M());                        
              }
              else
              {
                Histo(physics,triggerClassName,centrality,str->String().Data(),"MinvMMPt")->Fill(pj.Pt(),pj.M());                        
              }
            }
          }
        }
      }
    }
  } //track loop
}
Commit	Line	Data
5eabea87	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
b0bb161e	270
b0bb161e	271	if (!fVertex) return -999.;
5eabea87	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 5TMath::Sqrt(cutValuecutValue + 0.40.4pTot.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(xdcaxdca+ydcaydca);
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,
fe0324de	441	AliESDEvent& esd)
5eabea87	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