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