PWG/muon/AliAnalysisTaskESDMuonFilter.cxx

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