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 using std::cout;
  59 using std::endl;
  60 ClassImp(AliAnalysisTaskESDMuonFilter)
  61 ClassImp(AliAnalysisNonMuonTrackCuts)
  62
  63 ////////////////////////////////////////////////////////////////////////
  64
  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
 106 AliAnalysisTaskESDMuonFilter::AliAnalysisTaskESDMuonFilter(Bool_t onlyMuon, Bool_t keepAllEvents, Int_t mcMode):
 107   AliAnalysisTaskSE(),
 108   fTrackFilter(0x0),
 109   fEnableMuonAOD(kFALSE),
 110   fEnableDimuonAOD(kFALSE),
 111   fOnlyMuon(onlyMuon),
 112   fKeepAllEvents(keepAllEvents),
 113   fMCMode(mcMode)
 114 {
 115   // Default constructor
 116 }
 117
 118 AliAnalysisTaskESDMuonFilter::AliAnalysisTaskESDMuonFilter(const char* name, Bool_t onlyMuon, Bool_t keepAllEvents, Int_t mcMode):
 119   AliAnalysisTaskSE(name),
 120   fTrackFilter(0x0),
 121   fEnableMuonAOD(kFALSE),
 122   fEnableDimuonAOD(kFALSE),
 123   fOnlyMuon(onlyMuon),
 124   fKeepAllEvents(keepAllEvents),
 125   fMCMode(mcMode)
 126 {
 127   // Constructor
 128 }
 129
 130 //______________________________________________________________________________
 131 void AliAnalysisTaskESDMuonFilter::UserCreateOutputObjects()
 132 {
 133   // Create the output container
 134   if (fTrackFilter) OutputTree()->GetUserInfo()->Add(fTrackFilter);
 135 }
 136
 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   }
 163
 164   if ( fMCMode > 0 )
 165   {
 166     cout << spaces.Data() << "Assuming work on MC data (i.e. will transmit MC branches)" << endl;
 167   }
 168 }
 169
 170 //______________________________________________________________________________
 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
 180   if (!ext) return;
 181
 182   if ( fOnlyMuon )
 183   {
 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);
 190
 191     ext->DropUnspecifiedBranches(); // all branches not part of a FilterBranch call (below) will be dropped
 192
 193     ext->FilterBranch("tracks",murep);
 194     ext->FilterBranch("vertices",murep);
 195     ext->FilterBranch("dimuons",murep);
 196     ext->FilterBranch("AliAODVZERO",murep);
 197     ext->FilterBranch("AliAODTZERO",murep);
 198
 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     }
 209   }
 210 }
 211
 212 //______________________________________________________________________________
 213 void AliAnalysisTaskESDMuonFilter::Init()
 214 {
 215   // Initialization
 216   if(fEnableMuonAOD) AddFilteredAOD("AliAOD.Muons.root", "MuonEvents");
 217   if(fEnableDimuonAOD) AddFilteredAOD("AliAOD.Dimuons.root", "DimuonEvents");
 218 }
 219
 220
 221 //______________________________________________________________________________
 222 void AliAnalysisTaskESDMuonFilter::UserExec(Option_t */*option*/)
 223 {
 224   // Execute analysis for current event
 225
 226   Long64_t ientry = Entry();
 227   if(fDebug)printf("Muon Filter: Analysing event # %5d\n", (Int_t) ientry);
 228
 229   ConvertESDtoAOD();
 230 }
 231
 232 //______________________________________________________________________________
 233 void AliAnalysisTaskESDMuonFilter::ConvertESDtoAOD()
 234 {
 235   // ESD Muon Filter analysis task executed for each event
 236
 237   AliCodeTimerAuto("",0);
 238
 239   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
 240   if (!esd)
 241   {
 242     AliError("Could not get input ESD event");
 243     return;
 244   }
 245
 246   AliMCEventHandler *mcH = static_cast<AliMCEventHandler*>((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
 247
 248   // Define arrays for muons
 249   Double_t pos[3];
 250   Double_t p[3];
 251   Double_t pid[10];
 252
 253   // has to be changed once the muon pid is provided by the ESD
 254   for (Int_t i = 0; i < 10; pid[i++] = 0.) {}
 255   pid[AliAODTrack::kMuon]=1.;
 256
 257   AliAODHeader* header = AODEvent()->GetHeader();
 258   AliAODTrack *aodTrack = 0x0;
 259   AliESDMuonTrack *esdMuTrack = 0x0;
 260
 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
 266   AliAODVertex *primary = AODEvent()->GetPrimaryVertex();
 267   if (fDebug && primary) primary->Print();
 268
 269   // Loop on muon tracks to fill the AOD track branch
 270   Int_t nMuTracks = esd->GetNumberOfMuonTracks();
 271
 272   for (Int_t iTrack=0; iTrack<nMuTracks; ++iTrack) esd->GetMuonTrack(iTrack)->SetESDEvent(esd);
 273
 274   // Update number of positive and negative tracks from AOD event (M.G.)
 275   Int_t nPosTracks = header->GetRefMultiplicityPos();
 276   Int_t nNegTracks = header->GetRefMultiplicityNeg();
 277
 278   // Access to the AOD container of dimuons
 279   TClonesArray &dimuons = *(AODEvent()->GetDimuons());
 280
 281   Int_t nMuons=0;
 282   Int_t nDimuons=0;
 283   Int_t jDimuons=0;
 284   Int_t nMuonTrack[100];
 285   UChar_t itsClusMap(0);
 286
 287   for(int imuon=0;imuon<100;imuon++) nMuonTrack[imuon]=0;
 288
 289   for (Int_t nMuTrack = 0; nMuTrack < nMuTracks; ++nMuTrack)
 290   {
 291     esdMuTrack = esd->GetMuonTrack(nMuTrack);
 292
 293     if (!esdMuTrack->ContainTrackerData()) continue;
 294
 295     UInt_t selectInfo(0);
 296
 297     // Track selection
 298     if (fTrackFilter) {
 299         selectInfo = fTrackFilter->IsSelected(esdMuTrack);
 300         if (!selectInfo) {
 301           continue;
 302         }
 303     }
 304
 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
 313     if (mcH) mcH->SelectParticle(esdMuTrack->GetLabel()); // to insure that particle's ancestors will be in output MC branches
 314
 315     aodTrack = new(tracks[jTracks++]) AliAODTrack(esdMuTrack->GetUniqueID(), // ID
 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
 323                                                   itsClusMap, // ITSClusterMap
 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);
 330
 331     aodTrack->SetXYAtDCA(esdMuTrack->GetNonBendingCoorAtDCA(), esdMuTrack->GetBendingCoorAtDCA());
 332     aodTrack->SetPxPyPzAtDCA(esdMuTrack->PxAtDCA(), esdMuTrack->PyAtDCA(), esdMuTrack->PzAtDCA());
 333     aodTrack->SetRAtAbsorberEnd(esdMuTrack->GetRAtAbsorberEnd());
 334     aodTrack->ConvertAliPIDtoAODPID();
 335     aodTrack->SetChi2perNDF(esdMuTrack->GetChi2() / (2.*esdMuTrack->GetNHit() - 5.));
 336     aodTrack->SetChi2MatchTrigger(esdMuTrack->GetChi2MatchTrigger());
 337     aodTrack->SetHitsPatternInTrigCh(esdMuTrack->GetHitsPatternInTrigCh());
 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());
 342     aodTrack->SetMuonClusterMap(esdMuTrack->GetMuonClusterMap());
 343     aodTrack->SetMatchTrigger(esdMuTrack->GetMatchTrigger());
 344     aodTrack->Connected(esdMuTrack->IsConnected());
 345     primary->AddDaughter(aodTrack);
 346
 347     if (esdMuTrack->Charge() > 0) nPosTracks++;
 348     else nNegTracks++;
 349
 350     nMuonTrack[nMuons]= jTracks-1;
 351     ++nMuons;
 352   }
 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];
 360         new(dimuons[jDimuons++]) AliAODDimuon(tracks.At(index0),tracks.At(index1));
 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         }
 370       }
 371     }
 372   }
 373
 374
 375   header->SetRefMultiplicity(jTracks);
 376   header->SetRefMultiplicityPos(nPosTracks);
 377   header->SetRefMultiplicityNeg(nNegTracks);
 378   header->SetNumberOfMuons(nMuons);
 379   header->SetNumberOfDimuons(nDimuons);
 380
 381   AliAODHandler* handler = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 382
 383   if ( handler && fEnableMuonAOD && ( (nMuons>0) || fKeepAllEvents ) )
 384   {
 385     AliAODExtension *extMuons = handler->GetFilteredAOD("AliAOD.Muons.root");
 386     if ( extMuons ) extMuons->SelectEvent();
 387   }
 388
 389   if ( handler && fEnableDimuonAOD && ( (nMuons>1) || fKeepAllEvents )  )
 390   {
 391     AliAODExtension *extDimuons = handler->GetFilteredAOD("AliAOD.Dimuons.root");
 392     if ( extDimuons ) extDimuons->SelectEvent();
 393   }
 394
 395 }
 396
 397 void AliAnalysisTaskESDMuonFilter::Terminate(Option_t */*option*/)
 398 {
 399   // Terminate analysis
 400   //
 401   if (fDebug > 1) printf("AnalysisESDfilter: Terminate() \n");
 402 }