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