[u/mrichter/AliRoot.git] / PWG / muon / AliAnalysisTaskESDMuonFilter.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$ */

//
// Add the muon tracks to the generic AOD track branch during the 
// filtering of the ESD. 
//
// Authors: R. Arnaldi 5/5/08 and L. Aphecetche January 2011
//
// Note that we :
//   - completely disable all the branches that are not required by (most) the muon analyses,
//     e.g. cascades, v0s, kinks, jets, etc...
//   - filter severely the tracks (keep only muon tracks) and vertices (keep only primary -including
//     pile-up - vertices) branches 
// 
// (see AddFilteredAOD method)
//

#include "AliAnalysisTaskESDMuonFilter.h"

#include "AliAODDimuon.h"
#include "AliAODEvent.h"
#include "AliAODHandler.h"
#include "AliAODExtension.h"
#include "AliAODMCParticle.h"
#include "AliAODMuonReplicator.h"
#include "AliAODVertex.h"
#include "AliAnalysisFilter.h"
#include "AliAnalysisManager.h"
#include "AliCodeTimer.h"
#include "AliESDEvent.h"
#include "AliESDInputHandler.h"
#include "AliESDMuonTrack.h"
#include "AliESDVertex.h"
#include "AliESDtrack.h"
#include "AliLog.h"
#include "AliMCEvent.h"
#include "AliMCEventHandler.h"
#include "AliMultiplicity.h"
#include <TChain.h>
#include <TFile.h>
#include <TParticle.h>

using std::cout;
using std::endl;
ClassImp(AliAnalysisTaskESDMuonFilter)
ClassImp(AliAnalysisNonMuonTrackCuts)

////////////////////////////////////////////////////////////////////////

AliAnalysisNonMuonTrackCuts::AliAnalysisNonMuonTrackCuts()
{
  // default ctor 
}

Bool_t AliAnalysisNonMuonTrackCuts::IsSelected(TObject* obj)
{
  // Returns true if the object is a muon track
  AliAODTrack* track = dynamic_cast<AliAODTrack*>(obj);
  if (track && track->IsMuonTrack()) return kTRUE;
  return kFALSE;
}

AliAnalysisNonPrimaryVertices::AliAnalysisNonPrimaryVertices()
{
  // default ctor   
}

Bool_t AliAnalysisNonPrimaryVertices::IsSelected(TObject* obj)
{
  // Returns true if the object is a primary vertex

  AliAODVertex* vertex = dynamic_cast<AliAODVertex*>(obj);
  if (vertex)
  {
    if ( vertex->GetType() == AliAODVertex::kPrimary ||
        vertex->GetType() == AliAODVertex::kMainSPD ||
        vertex->GetType() == AliAODVertex::kPileupSPD ||
        vertex->GetType() == AliAODVertex::kPileupTracks ||
        vertex->GetType() == AliAODVertex::kMainTPC )
    {
      return kTRUE;
    }
  }
  
//  enum AODVtx_t {kUndef=-1, kPrimary, kKink, kV0, kCascade, kMulti, kMainSPD, kPileupSPD, kPileupTracks,kMainTPC};

  return kFALSE;
  
}

AliAnalysisTaskESDMuonFilter::AliAnalysisTaskESDMuonFilter(Bool_t onlyMuon, Bool_t keepAllEvents, Int_t mcMode):
  AliAnalysisTaskSE(),
  fTrackFilter(0x0),
  fEnableMuonAOD(kFALSE),
  fEnableDimuonAOD(kFALSE),
  fOnlyMuon(onlyMuon),
  fKeepAllEvents(keepAllEvents),
  fMCMode(mcMode)
{
  // Default constructor
}

AliAnalysisTaskESDMuonFilter::AliAnalysisTaskESDMuonFilter(const char* name, Bool_t onlyMuon, Bool_t keepAllEvents, Int_t mcMode):
  AliAnalysisTaskSE(name),
  fTrackFilter(0x0),
  fEnableMuonAOD(kFALSE),
  fEnableDimuonAOD(kFALSE),
  fOnlyMuon(onlyMuon),
  fKeepAllEvents(keepAllEvents),
  fMCMode(mcMode)
{
  // Constructor
}

//______________________________________________________________________________
void AliAnalysisTaskESDMuonFilter::UserCreateOutputObjects()
{
  // Create the output container
  if (fTrackFilter) OutputTree()->GetUserInfo()->Add(fTrackFilter);
}

//______________________________________________________________________________
void AliAnalysisTaskESDMuonFilter::PrintTask(Option_t *option, Int_t indent) const
{
  // Specify how we are configured
  
  AliAnalysisTaskSE::PrintTask(option,indent);
  
  TString spaces(' ',indent+3);
  
  if ( fOnlyMuon ) 
  {
    cout << spaces.Data() << "Keep only muon information " << endl;        
  }
  else 
  {
    cout << spaces.Data() << "Keep all information from standard AOD" << endl;
  }

  if ( fKeepAllEvents ) 
  {
    cout << spaces.Data() << "Keep all events, regardless of number of muons" << endl;    
  }
  else 
  {
    cout << spaces.Data() << "Keep only events with at least one muon" << endl;
  }
  
  if ( fMCMode > 0 ) 
  {
    cout << spaces.Data() << "Assuming work on MC data (i.e. will transmit MC branches)" << endl;
  }
}

//______________________________________________________________________________
void AliAnalysisTaskESDMuonFilter::AddFilteredAOD(const char* aodfilename, const char* title)
{
  // Add an output filtered and replicated aod
  
  AliAODHandler *aodH = (AliAODHandler*)((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
  if (!aodH) Fatal("UserCreateOutputObjects", "No AOD handler");

  AliAODExtension* ext = aodH->AddFilteredAOD(aodfilename,title);

  if (!ext) return;
  
  if ( fOnlyMuon ) 
  {    
    
    AliAODMuonReplicator* murep = new AliAODMuonReplicator("MuonReplicator",
                                                           "remove non muon tracks and non primary or pileup vertices",
                                                           new AliAnalysisNonMuonTrackCuts,
                                                           new AliAnalysisNonPrimaryVertices,
                                                           fMCMode);
    
    ext->DropUnspecifiedBranches(); // all branches not part of a FilterBranch call (below) will be dropped
    
    ext->FilterBranch("tracks",murep);    
    ext->FilterBranch("vertices",murep);  
    ext->FilterBranch("dimuons",murep);
    ext->FilterBranch("AliAODVZERO",murep);
    ext->FilterBranch("AliAODTZERO",murep);
    
    if ( fMCMode > 0 ) 
    {
      // MC branches will be copied (if present), as they are, but only
      // for events with at least one muon. 
      // For events w/o muon, mcparticles array will be empty and mcheader will be dummy
      // (e.g. strlen(GetGeneratorName())==0)
      
      ext->FilterBranch("mcparticles",murep);
      ext->FilterBranch("mcHeader",murep);
    }
  }  
}

//______________________________________________________________________________
void AliAnalysisTaskESDMuonFilter::Init()
{
  // Initialization
  if(fEnableMuonAOD) AddFilteredAOD("AliAOD.Muons.root", "MuonEvents");
  if(fEnableDimuonAOD) AddFilteredAOD("AliAOD.Dimuons.root", "DimuonEvents");    
}


//______________________________________________________________________________
void AliAnalysisTaskESDMuonFilter::UserExec(Option_t */*option*/)
{
  // Execute analysis for current event					    
  
  Long64_t ientry = Entry();
  if(fDebug)printf("Muon Filter: Analysing event # %5d\n", (Int_t) ientry);
  
  ConvertESDtoAOD();
}

//______________________________________________________________________________
void AliAnalysisTaskESDMuonFilter::ConvertESDtoAOD() 
{
  // ESD Muon Filter analysis task executed for each event
  
  AliCodeTimerAuto("",0);
  
  AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
  if (!esd) 
  {
    AliError("Could not get input ESD event");
    return;    
  }
  
  AliMCEventHandler *mcH = static_cast<AliMCEventHandler*>((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
    
  // Define arrays for muons
  Double_t pos[3];
  Double_t p[3];
  Double_t pid[10];
  
  // has to be changed once the muon pid is provided by the ESD
  for (Int_t i = 0; i < 10; pid[i++] = 0.) {}
  pid[AliAODTrack::kMuon]=1.;
  
  AliAODHeader* header = AODEvent()->GetHeader();
  AliAODTrack *aodTrack = 0x0;
  AliESDMuonTrack *esdMuTrack = 0x0;
  
  // Access to the AOD container of tracks
  TClonesArray &tracks = *(AODEvent()->GetTracks());
  Int_t jTracks = header->GetRefMultiplicity();
  
  // Read primary vertex from AOD event 
  AliAODVertex *primary = AODEvent()->GetPrimaryVertex();
  if (fDebug && primary) primary->Print();
  
  // Loop on muon tracks to fill the AOD track branch
  Int_t nMuTracks = esd->GetNumberOfMuonTracks();

  for (Int_t iTrack=0; iTrack<nMuTracks; ++iTrack) esd->GetMuonTrack(iTrack)->SetESDEvent(esd);
  
  // Update number of positive and negative tracks from AOD event (M.G.)
  Int_t nPosTracks = header->GetRefMultiplicityPos();
  Int_t nNegTracks = header->GetRefMultiplicityNeg();
  
  // Access to the AOD container of dimuons
  TClonesArray &dimuons = *(AODEvent()->GetDimuons());
  
  Int_t nMuons=0;
  Int_t nDimuons=0;
  Int_t jDimuons=0;
  Int_t nMuonTrack[100];
  UChar_t itsClusMap(0);
  
  for(int imuon=0;imuon<100;imuon++) nMuonTrack[imuon]=0;
  
  for (Int_t nMuTrack = 0; nMuTrack < nMuTracks; ++nMuTrack)
  {
    esdMuTrack = esd->GetMuonTrack(nMuTrack);
    
    if (!esdMuTrack->ContainTrackerData()) continue;
    
    UInt_t selectInfo(0);
    
    // Track selection
    if (fTrackFilter) {
     	selectInfo = fTrackFilter->IsSelected(esdMuTrack);
     	if (!selectInfo) {
     	  continue;
     	}  
    }
    
    p[0] = esdMuTrack->Px(); 
    p[1] = esdMuTrack->Py(); 
    p[2] = esdMuTrack->Pz();
    
    pos[0] = esdMuTrack->GetNonBendingCoor(); 
    pos[1] = esdMuTrack->GetBendingCoor(); 
    pos[2] = esdMuTrack->GetZ();
    
    if (mcH) mcH->SelectParticle(esdMuTrack->GetLabel()); // to insure that particle's ancestors will be in output MC branches
    
    aodTrack = new(tracks[jTracks++]) AliAODTrack(esdMuTrack->GetUniqueID(), // ID
                                                  esdMuTrack->GetLabel(), // label
                                                  p, // momentum
                                                  kTRUE, // cartesian coordinate system
                                                  pos, // position
                                                  kFALSE, // isDCA
                                                  0x0, // covariance matrix
                                                  esdMuTrack->Charge(), // charge
                                                  itsClusMap, // ITSClusterMap
                                                  pid, // pid
                                                  primary, // primary vertex
                                                  kFALSE, // used for vertex fit?
                                                  kFALSE, // used for primary vertex fit?
                                                  AliAODTrack::kPrimary,// track type
                                                  selectInfo); 
    
    aodTrack->SetXYAtDCA(esdMuTrack->GetNonBendingCoorAtDCA(), esdMuTrack->GetBendingCoorAtDCA());
    aodTrack->SetPxPyPzAtDCA(esdMuTrack->PxAtDCA(), esdMuTrack->PyAtDCA(), esdMuTrack->PzAtDCA());
    aodTrack->SetRAtAbsorberEnd(esdMuTrack->GetRAtAbsorberEnd());
    aodTrack->ConvertAliPIDtoAODPID();
    aodTrack->SetChi2perNDF(esdMuTrack->GetChi2() / (2.*esdMuTrack->GetNHit() - 5.));
    aodTrack->SetChi2MatchTrigger(esdMuTrack->GetChi2MatchTrigger());
    aodTrack->SetHitsPatternInTrigCh(esdMuTrack->GetHitsPatternInTrigCh());
    UInt_t pattern = esdMuTrack->GetHitsPatternInTrigCh();
    AliESDMuonTrack::AddEffInfo(pattern, 0, esdMuTrack->LoCircuit(), (AliESDMuonTrack::EAliTriggerChPatternFlag)0);
    aodTrack->SetMUONtrigHitsMapTrg(pattern);
    aodTrack->SetMUONtrigHitsMapTrk(esdMuTrack->GetHitsPatternInTrigChTrk());
    aodTrack->SetMuonClusterMap(esdMuTrack->GetMuonClusterMap());
    aodTrack->SetMatchTrigger(esdMuTrack->GetMatchTrigger());
    aodTrack->Connected(esdMuTrack->IsConnected());
    primary->AddDaughter(aodTrack);
    
    if (esdMuTrack->Charge() > 0) nPosTracks++;
    else nNegTracks++;
    
    nMuonTrack[nMuons]= jTracks-1;
    ++nMuons;
  }
  
  if(nMuons>=2) 
  { 
    for(int i=0;i<nMuons;i++){
      Int_t index0 = nMuonTrack[i];
      for(int j=i+1;j<nMuons;j++){
        Int_t index1 = nMuonTrack[j];
        new(dimuons[jDimuons++]) AliAODDimuon(tracks.At(index0),tracks.At(index1));
        ++nDimuons;
        if (fDebug > 1){
          AliAODDimuon *dimuon0 = (AliAODDimuon*)dimuons.At(jDimuons-1);
          printf("Dimuon: mass = %f, px=%f, py=%f, pz=%f\n",dimuon0->M(),dimuon0->Px(),dimuon0->Py(),dimuon0->Pz());  
          AliAODTrack  *mu0 = (AliAODTrack*) dimuon0->GetMu(0);
          AliAODTrack  *mu1 = (AliAODTrack*) dimuon0->GetMu(1);
          printf("Muon0 px=%f py=%f pz=%f\n",mu0->Px(),mu0->Py(),mu0->Pz());
          printf("Muon1 px=%f py=%f pz=%f\n",mu1->Px(),mu1->Py(),mu1->Pz());
        }  
      }
    }
  }
  
  
  header->SetRefMultiplicity(jTracks); 
  header->SetRefMultiplicityPos(nPosTracks);
  header->SetRefMultiplicityNeg(nNegTracks);
  header->SetNumberOfMuons(nMuons);
  header->SetNumberOfDimuons(nDimuons);
  
  AliAODHandler* handler = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
  
  if ( handler && fEnableMuonAOD && ( (nMuons>0) || fKeepAllEvents ) )
  {
    AliAODExtension *extMuons = handler->GetFilteredAOD("AliAOD.Muons.root");
    if ( extMuons ) extMuons->SelectEvent();
  }
  
  if ( handler && fEnableDimuonAOD && ( (nMuons>1) || fKeepAllEvents )  )
  {
    AliAODExtension *extDimuons = handler->GetFilteredAOD("AliAOD.Dimuons.root");
    if ( extDimuons ) extDimuons->SelectEvent();
  }
  
}

void AliAnalysisTaskESDMuonFilter::Terminate(Option_t */*option*/)
{
  // Terminate analysis
  //
  if (fDebug > 1) printf("AnalysisESDfilter: Terminate() \n");
}
Commit	Line	Data
1d1ea3ce	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
27de2dfb	16	/* $Id$ */
27de2dfb	17
26ba01d4	18	//
1d1ea3ce	19	// Add the muon tracks to the generic AOD track branch during the
26ba01d4	20	// filtering of the ESD.
	21	//
	22	// Authors: R. Arnaldi 5/5/08 and L. Aphecetche January 2011
	23	//
	24	// Note that we :
	25	// - completely disable all the branches that are not required by (most) the muon analyses,
	26	// e.g. cascades, v0s, kinks, jets, etc...
	27	// - filter severely the tracks (keep only muon tracks) and vertices (keep only primary -including
	28	// pile-up - vertices) branches
	29	//
	30	// (see AddFilteredAOD method)
	31	//
aba11173	32
1d1ea3ce	33	#include "AliAnalysisTaskESDMuonFilter.h"
dcc1f876	34
dcc1f876	35	#include "AliAODDimuon.h"
1d1ea3ce	36	#include "AliAODEvent.h"
1d1ea3ce	37	#include "AliAODHandler.h"
14d6fad5	38	#include "AliAODExtension.h"
dcc1f876	39	#include "AliAODMCParticle.h"
	40	#include "AliAODMuonReplicator.h"
	41	#include "AliAODVertex.h"
1d1ea3ce	42	#include "AliAnalysisFilter.h"
dcc1f876	43	#include "AliAnalysisManager.h"
	44	#include "AliCodeTimer.h"
	45	#include "AliESDEvent.h"
	46	#include "AliESDInputHandler.h"
1d1ea3ce	47	#include "AliESDMuonTrack.h"
1d1ea3ce	48	#include "AliESDVertex.h"
dcc1f876	49	#include "AliESDtrack.h"
1d1ea3ce	50	#include "AliLog.h"
fc3a4c45	51	#include "AliMCEvent.h"
fc3a4c45	52	#include "AliMCEventHandler.h"
dcc1f876	53	#include "AliMultiplicity.h"
dcc1f876	54	#include <TChain.h>
	55	#include <TFile.h>
	56	#include <TParticle.h>
1d1ea3ce	57
3a7af7bd	58	using std::cout;
3a7af7bd	59	using std::endl;
1d1ea3ce	60	ClassImp(AliAnalysisTaskESDMuonFilter)
26ba01d4	61	ClassImp(AliAnalysisNonMuonTrackCuts)
1d1ea3ce	62
	63	////////////////////////////////////////////////////////////////////////
	64
26ba01d4	65	AliAnalysisNonMuonTrackCuts::AliAnalysisNonMuonTrackCuts()
	66	{
	67	// default ctor
	68	}
	69
	70	Bool_t AliAnalysisNonMuonTrackCuts::IsSelected(TObject* obj)
	71	{
	72	// Returns true if the object is a muon track
	73	AliAODTrack* track = dynamic_cast<AliAODTrack*>(obj);
	74	if (track && track->IsMuonTrack()) return kTRUE;
	75	return kFALSE;
	76	}
	77
	78	AliAnalysisNonPrimaryVertices::AliAnalysisNonPrimaryVertices()
	79	{
	80	// default ctor
	81	}
	82
	83	Bool_t AliAnalysisNonPrimaryVertices::IsSelected(TObject* obj)
	84	{
	85	// Returns true if the object is a primary vertex
	86
	87	AliAODVertex* vertex = dynamic_cast<AliAODVertex*>(obj);
	88	if (vertex)
	89	{
	90	if ( vertex->GetType() == AliAODVertex::kPrimary \|\|
	91	vertex->GetType() == AliAODVertex::kMainSPD \|\|
	92	vertex->GetType() == AliAODVertex::kPileupSPD \|\|
	93	vertex->GetType() == AliAODVertex::kPileupTracks \|\|
	94	vertex->GetType() == AliAODVertex::kMainTPC )
	95	{
	96	return kTRUE;
	97	}
	98	}
	99
	100	// enum AODVtx_t {kUndef=-1, kPrimary, kKink, kV0, kCascade, kMulti, kMainSPD, kPileupSPD, kPileupTracks,kMainTPC};
	101
	102	return kFALSE;
	103
	104	}
	105
14d6fad5	106	AliAnalysisTaskESDMuonFilter::AliAnalysisTaskESDMuonFilter(Bool_t onlyMuon, Bool_t keepAllEvents, Int_t mcMode):
aba11173	107	AliAnalysisTaskSE(),
f785fc59	108	fTrackFilter(0x0),
4a497116	109	fEnableMuonAOD(kFALSE),
26ba01d4	110	fEnableDimuonAOD(kFALSE),
26ba01d4	111	fOnlyMuon(onlyMuon),
14d6fad5	112	fKeepAllEvents(keepAllEvents),
14d6fad5	113	fMCMode(mcMode)
1d1ea3ce	114	{
	115	// Default constructor
	116	}
	117
14d6fad5	118	AliAnalysisTaskESDMuonFilter::AliAnalysisTaskESDMuonFilter(const char* name, Bool_t onlyMuon, Bool_t keepAllEvents, Int_t mcMode):
aba11173	119	AliAnalysisTaskSE(name),
f785fc59	120	fTrackFilter(0x0),
4a497116	121	fEnableMuonAOD(kFALSE),
26ba01d4	122	fEnableDimuonAOD(kFALSE),
26ba01d4	123	fOnlyMuon(onlyMuon),
14d6fad5	124	fKeepAllEvents(keepAllEvents),
14d6fad5	125	fMCMode(mcMode)
1d1ea3ce	126	{
	127	// Constructor
	128	}
	129
3d270be0	130	//______________________________________________________________________________
1d1ea3ce	131	void AliAnalysisTaskESDMuonFilter::UserCreateOutputObjects()
	132	{
	133	// Create the output container
aba11173	134	if (fTrackFilter) OutputTree()->GetUserInfo()->Add(fTrackFilter);
1d1ea3ce	135	}
1d1ea3ce	136
3d270be0	137	//______________________________________________________________________________
	138	void AliAnalysisTaskESDMuonFilter::PrintTask(Option_t *option, Int_t indent) const
	139	{
	140	// Specify how we are configured
	141
	142	AliAnalysisTaskSE::PrintTask(option,indent);
	143
	144	TString spaces(' ',indent+3);
	145
	146	if ( fOnlyMuon )
	147	{
	148	cout << spaces.Data() << "Keep only muon information " << endl;
	149	}
	150	else
	151	{
	152	cout << spaces.Data() << "Keep all information from standard AOD" << endl;
	153	}
	154
	155	if ( fKeepAllEvents )
	156	{
	157	cout << spaces.Data() << "Keep all events, regardless of number of muons" << endl;
	158	}
	159	else
	160	{
	161	cout << spaces.Data() << "Keep only events with at least one muon" << endl;
	162	}
14d6fad5	163
	164	if ( fMCMode > 0 )
	165	{
	166	cout << spaces.Data() << "Assuming work on MC data (i.e. will transmit MC branches)" << endl;
	167	}
3d270be0	168	}
	169
	170	//______________________________________________________________________________
26ba01d4	171	void AliAnalysisTaskESDMuonFilter::AddFilteredAOD(const char* aodfilename, const char* title)
	172	{
	173	// Add an output filtered and replicated aod
	174
	175	AliAODHandler aodH = (AliAODHandler)((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
	176	if (!aodH) Fatal("UserCreateOutputObjects", "No AOD handler");
	177
	178	AliAODExtension* ext = aodH->AddFilteredAOD(aodfilename,title);
	179
3d270be0	180	if (!ext) return;
	181
	182	if ( fOnlyMuon )
	183	{
14d6fad5	184
	185	AliAODMuonReplicator* murep = new AliAODMuonReplicator("MuonReplicator",
	186	"remove non muon tracks and non primary or pileup vertices",
	187	new AliAnalysisNonMuonTrackCuts,
	188	new AliAnalysisNonPrimaryVertices,
	189	fMCMode);
dcc1f876	190
	191	ext->DropUnspecifiedBranches(); // all branches not part of a FilterBranch call (below) will be dropped
	192
26ba01d4	193	ext->FilterBranch("tracks",murep);
3d270be0	194	ext->FilterBranch("vertices",murep);
3d270be0	195	ext->FilterBranch("dimuons",murep);
5393f2c0	196	ext->FilterBranch("AliAODVZERO",murep);
3493cd3f	197	ext->FilterBranch("AliAODTZERO",murep);
5393f2c0	198
14d6fad5	199	if ( fMCMode > 0 )
	200	{
	201	// MC branches will be copied (if present), as they are, but only
	202	// for events with at least one muon.
	203	// For events w/o muon, mcparticles array will be empty and mcheader will be dummy
	204	// (e.g. strlen(GetGeneratorName())==0)
	205
	206	ext->FilterBranch("mcparticles",murep);
	207	ext->FilterBranch("mcHeader",murep);
	208	}
26ba01d4	209	}
	210	}
	211
3d270be0	212	//______________________________________________________________________________
1d1ea3ce	213	void AliAnalysisTaskESDMuonFilter::Init()
	214	{
	215	// Initialization
26ba01d4	216	if(fEnableMuonAOD) AddFilteredAOD("AliAOD.Muons.root", "MuonEvents");
26ba01d4	217	if(fEnableDimuonAOD) AddFilteredAOD("AliAOD.Dimuons.root", "DimuonEvents");
1d1ea3ce	218	}
	219
	220
3d270be0	221	//______________________________________________________________________________
1d1ea3ce	222	void AliAnalysisTaskESDMuonFilter::UserExec(Option_t /option*/)
	223	{
	224	// Execute analysis for current event
26ba01d4	225
1d1ea3ce	226	Long64_t ientry = Entry();
edee4062	227	if(fDebug)printf("Muon Filter: Analysing event # %5d\n", (Int_t) ientry);
aba11173	228
1d1ea3ce	229	ConvertESDtoAOD();
	230	}
	231
3d270be0	232	//______________________________________________________________________________
1d1ea3ce	233	void AliAnalysisTaskESDMuonFilter::ConvertESDtoAOD()
	234	{
	235	// ESD Muon Filter analysis task executed for each event
dcc1f876	236
	237	AliCodeTimerAuto("",0);
	238
1d1ea3ce	239	AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
6fa3e5c7	240	if (!esd)
	241	{
	242	AliError("Could not get input ESD event");
	243	return;
fc3a4c45	244	}
6fa3e5c7	245
6fa3e5c7	246	AliMCEventHandler mcH = static_cast<AliMCEventHandler>((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
fc3a4c45	247
1d1ea3ce	248	// Define arrays for muons
	249	Double_t pos[3];
	250	Double_t p[3];
	251	Double_t pid[10];
aba11173	252
aba11173	253	// has to be changed once the muon pid is provided by the ESD
a6e0ebfe	254	for (Int_t i = 0; i < 10; pid[i++] = 0.) {}
aba11173	255	pid[AliAODTrack::kMuon]=1.;
aba11173	256
1d1ea3ce	257	AliAODHeader* header = AODEvent()->GetHeader();
1d1ea3ce	258	AliAODTrack *aodTrack = 0x0;
aba11173	259	AliESDMuonTrack *esdMuTrack = 0x0;
aba11173	260
1d1ea3ce	261	// Access to the AOD container of tracks
	262	TClonesArray &tracks = *(AODEvent()->GetTracks());
	263	Int_t jTracks = header->GetRefMultiplicity();
	264
	265	// Read primary vertex from AOD event
aba11173	266	AliAODVertex *primary = AODEvent()->GetPrimaryVertex();
6fa3e5c7	267	if (fDebug && primary) primary->Print();
aba11173	268
1d1ea3ce	269	// Loop on muon tracks to fill the AOD track branch
1d1ea3ce	270	Int_t nMuTracks = esd->GetNumberOfMuonTracks();
866d8d78	271
daf100e1	272	for (Int_t iTrack=0; iTrack<nMuTracks; ++iTrack) esd->GetMuonTrack(iTrack)->SetESDEvent(esd);
aba11173	273
aba11173	274	// Update number of positive and negative tracks from AOD event (M.G.)
1d1ea3ce	275	Int_t nPosTracks = header->GetRefMultiplicityPos();
1d1ea3ce	276	Int_t nNegTracks = header->GetRefMultiplicityNeg();
aba11173	277
866d8d78	278	// Access to the AOD container of dimuons
866d8d78	279	TClonesArray &dimuons = *(AODEvent()->GetDimuons());
866d8d78	280
866d8d78	281	Int_t nMuons=0;
29e7e51b	282	Int_t nDimuons=0;
866d8d78	283	Int_t jDimuons=0;
397081d5	284	Int_t nMuonTrack[100];
6fa3e5c7	285	UChar_t itsClusMap(0);
87ac315c	286
397081d5	287	for(int imuon=0;imuon<100;imuon++) nMuonTrack[imuon]=0;
26ba01d4	288
	289	for (Int_t nMuTrack = 0; nMuTrack < nMuTracks; ++nMuTrack)
	290	{
aba11173	291	esdMuTrack = esd->GetMuonTrack(nMuTrack);
	292
	293	if (!esdMuTrack->ContainTrackerData()) continue;
daf100e1	294
6fa3e5c7	295	UInt_t selectInfo(0);
6fa3e5c7	296
aba11173	297	// Track selection
	298	if (fTrackFilter) {
	299	selectInfo = fTrackFilter->IsSelected(esdMuTrack);
	300	if (!selectInfo) {
	301	continue;
	302	}
87ac315c	303	}
26ba01d4	304
aba11173	305	p[0] = esdMuTrack->Px();
	306	p[1] = esdMuTrack->Py();
	307	p[2] = esdMuTrack->Pz();
	308
	309	pos[0] = esdMuTrack->GetNonBendingCoor();
	310	pos[1] = esdMuTrack->GetBendingCoor();
	311	pos[2] = esdMuTrack->GetZ();
	312
6fa3e5c7	313	if (mcH) mcH->SelectParticle(esdMuTrack->GetLabel()); // to insure that particle's ancestors will be in output MC branches
fc3a4c45	314
aba11173	315	aodTrack = new(tracks[jTracks++]) AliAODTrack(esdMuTrack->GetUniqueID(), // ID
26ba01d4	316	esdMuTrack->GetLabel(), // label
	317	p, // momentum
	318	kTRUE, // cartesian coordinate system
	319	pos, // position
	320	kFALSE, // isDCA
	321	0x0, // covariance matrix
	322	esdMuTrack->Charge(), // charge
6fa3e5c7	323	itsClusMap, // ITSClusterMap
26ba01d4	324	pid, // pid
	325	primary, // primary vertex
	326	kFALSE, // used for vertex fit?
	327	kFALSE, // used for primary vertex fit?
	328	AliAODTrack::kPrimary,// track type
	329	selectInfo);
aba11173	330
	331	aodTrack->SetXYAtDCA(esdMuTrack->GetNonBendingCoorAtDCA(), esdMuTrack->GetBendingCoorAtDCA());
	332	aodTrack->SetPxPyPzAtDCA(esdMuTrack->PxAtDCA(), esdMuTrack->PyAtDCA(), esdMuTrack->PzAtDCA());
f43586f0	333	aodTrack->SetRAtAbsorberEnd(esdMuTrack->GetRAtAbsorberEnd());
aba11173	334	aodTrack->ConvertAliPIDtoAODPID();
	335	aodTrack->SetChi2perNDF(esdMuTrack->GetChi2() / (2.*esdMuTrack->GetNHit() - 5.));
	336	aodTrack->SetChi2MatchTrigger(esdMuTrack->GetChi2MatchTrigger());
	337	aodTrack->SetHitsPatternInTrigCh(esdMuTrack->GetHitsPatternInTrigCh());
0a2dcc83	338	UInt_t pattern = esdMuTrack->GetHitsPatternInTrigCh();
	339	AliESDMuonTrack::AddEffInfo(pattern, 0, esdMuTrack->LoCircuit(), (AliESDMuonTrack::EAliTriggerChPatternFlag)0);
	340	aodTrack->SetMUONtrigHitsMapTrg(pattern);
	341	aodTrack->SetMUONtrigHitsMapTrk(esdMuTrack->GetHitsPatternInTrigChTrk());
aba11173	342	aodTrack->SetMuonClusterMap(esdMuTrack->GetMuonClusterMap());
aba11173	343	aodTrack->SetMatchTrigger(esdMuTrack->GetMatchTrigger());
5c15a68b	344	aodTrack->Connected(esdMuTrack->IsConnected());
aba11173	345	primary->AddDaughter(aodTrack);
	346
	347	if (esdMuTrack->Charge() > 0) nPosTracks++;
	348	else nNegTracks++;
866d8d78	349
ca3ff602	350	nMuonTrack[nMuons]= jTracks-1;
26ba01d4	351	++nMuons;
87ac315c	352	}
26ba01d4	353
	354	if(nMuons>=2)
	355	{
	356	for(int i=0;i<nMuons;i++){
	357	Int_t index0 = nMuonTrack[i];
	358	for(int j=i+1;j<nMuons;j++){
	359	Int_t index1 = nMuonTrack[j];
76b26647	360	new(dimuons[jDimuons++]) AliAODDimuon(tracks.At(index0),tracks.At(index1));
26ba01d4	361	++nDimuons;
	362	if (fDebug > 1){
	363	AliAODDimuon dimuon0 = (AliAODDimuon)dimuons.At(jDimuons-1);
	364	printf("Dimuon: mass = %f, px=%f, py=%f, pz=%f\n",dimuon0->M(),dimuon0->Px(),dimuon0->Py(),dimuon0->Pz());
	365	AliAODTrack mu0 = (AliAODTrack) dimuon0->GetMu(0);
	366	AliAODTrack mu1 = (AliAODTrack) dimuon0->GetMu(1);
	367	printf("Muon0 px=%f py=%f pz=%f\n",mu0->Px(),mu0->Py(),mu0->Pz());
	368	printf("Muon1 px=%f py=%f pz=%f\n",mu1->Px(),mu1->Py(),mu1->Pz());
	369	}
87ac315c	370	}
866d8d78	371	}
26ba01d4	372	}
26ba01d4	373
aba11173	374
1d1ea3ce	375	header->SetRefMultiplicity(jTracks);
	376	header->SetRefMultiplicityPos(nPosTracks);
	377	header->SetRefMultiplicityNeg(nNegTracks);
29e7e51b	378	header->SetNumberOfMuons(nMuons);
29e7e51b	379	header->SetNumberOfDimuons(nDimuons);
866d8d78	380
6fa3e5c7	381	AliAODHandler* handler = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
	382
	383	if ( handler && fEnableMuonAOD && ( (nMuons>0) \|\| fKeepAllEvents ) )
26ba01d4	384	{
6fa3e5c7	385	AliAODExtension *extMuons = handler->GetFilteredAOD("AliAOD.Muons.root");
6fa3e5c7	386	if ( extMuons ) extMuons->SelectEvent();
866d8d78	387	}
26ba01d4	388
6fa3e5c7	389	if ( handler && fEnableDimuonAOD && ( (nMuons>1) \|\| fKeepAllEvents ) )
26ba01d4	390	{
6fa3e5c7	391	AliAODExtension *extDimuons = handler->GetFilteredAOD("AliAOD.Dimuons.root");
6fa3e5c7	392	if ( extDimuons ) extDimuons->SelectEvent();
f785fc59	393	}
26ba01d4	394
1d1ea3ce	395	}
	396
	397	void AliAnalysisTaskESDMuonFilter::Terminate(Option_t /option*/)
	398	{
aba11173	399	// Terminate analysis
	400	//
	401	if (fDebug > 1) printf("AnalysisESDfilter: Terminate() \n");
1d1ea3ce	402	}