PWG3/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 // Add the muon tracks to the generic AOD track branch during the
  17 // filtering of the ESD - R. Arnaldi 5/5/08
  18
  19 #include <TChain.h>
  20 #include <TFile.h>
  21 #include <TParticle.h>
  22
  23 #include "AliAnalysisTaskESDMuonFilter.h"
  24 #include "AliAnalysisManager.h"
  25 #include "AliESDEvent.h"
  26 #include "AliAODEvent.h"
  27 #include "AliESDInputHandler.h"
  28 #include "AliAODHandler.h"
  29 #include "AliAnalysisFilter.h"
  30 #include "AliESDtrack.h"
  31 #include "AliESDMuonTrack.h"
  32 #include "AliESDVertex.h"
  33 #include "AliMultiplicity.h"
  34 #include "AliLog.h"
  35 #include "AliStack.h"
  36 #include "AliMCEvent.h"
  37 #include "AliMCEventHandler.h"
  38 #include "AliAODMCParticle.h"
  39 #include "AliAODDimuon.h"
  40
  41 ClassImp(AliAnalysisTaskESDMuonFilter)
  42
  43 ////////////////////////////////////////////////////////////////////////
  44
  45 AliAnalysisTaskESDMuonFilter::AliAnalysisTaskESDMuonFilter():
  46   AliAnalysisTaskSE(),
  47   fTrackFilter(0x0),
  48   fEnableMuonAOD(kFALSE),
  49   fEnableDimuonAOD(kFALSE)
  50 {
  51   // Default constructor
  52 }
  53
  54 AliAnalysisTaskESDMuonFilter::AliAnalysisTaskESDMuonFilter(const char* name):
  55   AliAnalysisTaskSE(name),
  56   fTrackFilter(0x0),
  57   fEnableMuonAOD(kFALSE),
  58   fEnableDimuonAOD(kFALSE)
  59 {
  60   // Constructor
  61 }
  62
  63 void AliAnalysisTaskESDMuonFilter::UserCreateOutputObjects()
  64 {
  65   // Create the output container
  66   if (fTrackFilter) OutputTree()->GetUserInfo()->Add(fTrackFilter);
  67 }
  68
  69 void AliAnalysisTaskESDMuonFilter::Init()
  70 {
  71   // Initialization
  72   if (fDebug > 1) AliInfo("Init() \n");
  73   // From Andrei
  74   AliAODHandler *aodH = (AliAODHandler*)((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
  75   if (!aodH) Fatal("UserCreateOutputObjects", "No AOD handler. Aborting.");
  76   if(fEnableMuonAOD)aodH->AddFilteredAOD("AliAOD.Muons.root", "MuonEvents");
  77   if(fEnableDimuonAOD)aodH->AddFilteredAOD("AliAOD.Dimuons.root", "DimuonEvents");
  78 }
  79
  80
  81 void AliAnalysisTaskESDMuonFilter::UserExec(Option_t */*option*/)
  82 {
  83   // Execute analysis for current event
  84   Long64_t ientry = Entry();
  85   if(fDebug)printf("Muon Filter: Analysing event # %5d\n", (Int_t) ientry);
  86
  87   ConvertESDtoAOD();
  88 }
  89
  90 void AliAnalysisTaskESDMuonFilter::ConvertESDtoAOD()
  91 {
  92   // ESD Muon Filter analysis task executed for each event
  93   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
  94
  95   // Fetch Stack for debuggging if available
  96   AliStack *pStack = 0;
  97   AliMCEventHandler *mcH = 0;
  98   if(MCEvent()){
  99     pStack = MCEvent()->Stack();
 100     mcH = (AliMCEventHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
 101   }
 102
 103   // Define arrays for muons
 104   Double_t pos[3];
 105   Double_t p[3];
 106   Double_t pid[10];
 107
 108   // has to be changed once the muon pid is provided by the ESD
 109   for (Int_t i = 0; i < 10; pid[i++] = 0.) {}
 110   pid[AliAODTrack::kMuon]=1.;
 111
 112   AliAODHeader* header = AODEvent()->GetHeader();
 113   AliAODTrack *aodTrack = 0x0;
 114   AliESDMuonTrack *esdMuTrack = 0x0;
 115
 116   // Access to the AOD container of tracks
 117   TClonesArray &tracks = *(AODEvent()->GetTracks());
 118   Int_t jTracks = header->GetRefMultiplicity();
 119
 120   // Read primary vertex from AOD event
 121   AliAODVertex *primary = AODEvent()->GetPrimaryVertex();
 122   if(fDebug)primary->Print();
 123
 124   // Loop on muon tracks to fill the AOD track branch
 125   Int_t nMuTracks = esd->GetNumberOfMuonTracks();
 126
 127   for (Int_t iTrack=0; iTrack<nMuTracks; ++iTrack) esd->GetMuonTrack(iTrack)->SetESDEvent(esd);
 128
 129   // Update number of positive and negative tracks from AOD event (M.G.)
 130   Int_t nPosTracks = header->GetRefMultiplicityPos();
 131   Int_t nNegTracks = header->GetRefMultiplicityNeg();
 132
 133   // Access to the AOD container of dimuons
 134   TClonesArray &dimuons = *(AODEvent()->GetDimuons());
 135   AliAODDimuon *aodDimuon = 0x0;
 136
 137   Bool_t MuonsExist = kFALSE;
 138   Bool_t DimuonsExist = kFALSE;
 139   Int_t nMuons=0;
 140   Int_t nDimuons=0;
 141   Int_t jDimuons=0;
 142   Int_t nMuonTrack[100];
 143
 144   for(int imuon=0;imuon<100;imuon++) nMuonTrack[imuon]=0;
 145
 146   for (Int_t nMuTrack = 0; nMuTrack < nMuTracks; ++nMuTrack) {
 147     esdMuTrack = esd->GetMuonTrack(nMuTrack);
 148
 149     if (!esdMuTrack->ContainTrackerData()) continue;
 150
 151     MuonsExist = kTRUE;
 152
 153     UInt_t selectInfo = 0;
 154     // Track selection
 155     if (fTrackFilter) {
 156         selectInfo = fTrackFilter->IsSelected(esdMuTrack);
 157         if (!selectInfo) {
 158           continue;
 159         }
 160     }
 161
 162     p[0] = esdMuTrack->Px();
 163     p[1] = esdMuTrack->Py();
 164     p[2] = esdMuTrack->Pz();
 165
 166     pos[0] = esdMuTrack->GetNonBendingCoor();
 167     pos[1] = esdMuTrack->GetBendingCoor();
 168     pos[2] = esdMuTrack->GetZ();
 169
 170     if(mcH)mcH->SelectParticle(esdMuTrack->GetLabel());
 171
 172     aodTrack = new(tracks[jTracks++]) AliAODTrack(esdMuTrack->GetUniqueID(), // ID
 173                                                   esdMuTrack->GetLabel(), // label
 174                                                   p, // momentum
 175                                                   kTRUE, // cartesian coordinate system
 176                                                   pos, // position
 177                                                   kFALSE, // isDCA
 178                                                   0x0, // covariance matrix
 179                                                   esdMuTrack->Charge(), // charge
 180                                                   0, // ITSClusterMap
 181                                                   pid, // pid
 182                                                   primary, // primary vertex
 183                                                   kFALSE, // used for vertex fit?
 184                                                   kFALSE, // used for primary vertex fit?
 185                                                   AliAODTrack::kPrimary,// track type
 186                                                   selectInfo);
 187
 188     aodTrack->SetXYAtDCA(esdMuTrack->GetNonBendingCoorAtDCA(), esdMuTrack->GetBendingCoorAtDCA());
 189     aodTrack->SetPxPyPzAtDCA(esdMuTrack->PxAtDCA(), esdMuTrack->PyAtDCA(), esdMuTrack->PzAtDCA());
 190     aodTrack->SetRAtAbsorberEnd(esdMuTrack->GetRAtAbsorberEnd());
 191     aodTrack->ConvertAliPIDtoAODPID();
 192     aodTrack->SetChi2perNDF(esdMuTrack->GetChi2() / (2.*esdMuTrack->GetNHit() - 5.));
 193     aodTrack->SetChi2MatchTrigger(esdMuTrack->GetChi2MatchTrigger());
 194     aodTrack->SetHitsPatternInTrigCh(esdMuTrack->GetHitsPatternInTrigCh());
 195     aodTrack->SetMuonClusterMap(esdMuTrack->GetMuonClusterMap());
 196     aodTrack->SetMatchTrigger(esdMuTrack->GetMatchTrigger());
 197     aodTrack->Connected(esdMuTrack->IsConnected());
 198     primary->AddDaughter(aodTrack);
 199
 200     if (esdMuTrack->Charge() > 0) nPosTracks++;
 201     else nNegTracks++;
 202
 203     nMuonTrack[nMuons]= jTracks-1.;
 204     nMuons++;
 205   }
 206
 207     if(nMuons>=2) DimuonsExist = kTRUE;
 208     if(DimuonsExist) {
 209       for(int i=0;i<nMuons;i++){
 210         Int_t index0 = nMuonTrack[i];
 211         for(int j=i+1;j<nMuons;j++){
 212           Int_t index1 = nMuonTrack[j];
 213           aodDimuon = new(dimuons[jDimuons++]) AliAODDimuon(tracks.At(index0),tracks.At(index1));
 214           nDimuons++;
 215           if (fDebug > 1){
 216             AliAODDimuon *dimuon0 = (AliAODDimuon*)dimuons.At(jDimuons-1);
 217             printf("Dimuon: mass = %f, px=%f, py=%f, pz=%f\n",dimuon0->M(),dimuon0->Px(),dimuon0->Py(),dimuon0->Pz());
 218             AliAODTrack  *mu0 = (AliAODTrack*) dimuon0->GetMu(0);
 219             AliAODTrack  *mu1 = (AliAODTrack*) dimuon0->GetMu(1);
 220             printf("Muon0 px=%f py=%f pz=%f\n",mu0->Px(),mu0->Py(),mu0->Pz());
 221             printf("Muon1 px=%f py=%f pz=%f\n",mu1->Px(),mu1->Py(),mu1->Pz());
 222           }
 223         }
 224       }
 225     }
 226
 227
 228   header->SetRefMultiplicity(jTracks);
 229   header->SetRefMultiplicityPos(nPosTracks);
 230   header->SetRefMultiplicityNeg(nNegTracks);
 231   header->SetNumberOfMuons(nMuons);
 232   header->SetNumberOfDimuons(nDimuons);
 233
 234   // From Andrei
 235   if(fEnableMuonAOD && MuonsExist){
 236     AliAODExtension *extMuons = dynamic_cast<AliAODHandler*>
 237     ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler())->GetFilteredAOD("AliAOD.Muons.root");
 238     extMuons->SelectEvent();
 239   }
 240
 241   if(fEnableDimuonAOD && DimuonsExist){
 242     AliAODExtension *extDimuons = dynamic_cast<AliAODHandler*>
 243    ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler())->GetFilteredAOD("AliAOD.Dimuons.root");
 244     extDimuons->SelectEvent();
 245   }
 246
 247 }
 248
 249 void AliAnalysisTaskESDMuonFilter::Terminate(Option_t */*option*/)
 250 {
 251   // Terminate analysis
 252   //
 253   if (fDebug > 1) printf("AnalysisESDfilter: Terminate() \n");
 254 }
 255
 256 void  AliAnalysisTaskESDMuonFilter::PrintMCInfo(AliStack *pStack,Int_t label){
 257   if(!pStack)return;
 258   label = TMath::Abs(label);
 259   TParticle *part = pStack->Particle(label);
 260   Printf("########################");
 261   Printf("%s:%d %d UniqueID %d PDG %d P %3.3f",(char*)__FILE__,__LINE__,label,part->GetUniqueID(),part->GetPdgCode(),part->P());
 262   part->Print();
 263   TParticle* mother = part;
 264   Int_t imo = part->GetFirstMother();
 265   Int_t nprim = pStack->GetNprimary();
 266   //  while((imo >= nprim) && (mother->GetUniqueID() == 4)) {
 267   while((imo >= nprim)) {
 268     mother =  pStack->Particle(imo);
 269     Printf("Mother %s:%d Label %d UniqueID %d PDG %d P %3.3f",(char*)__FILE__,__LINE__,imo,mother->GetUniqueID(),mother->GetPdgCode(),mother->P());
 270     mother->Print();
 271     imo =  mother->GetFirstMother();
 272   }
 273   Printf("########################");
 274 }