MUON/DecodeRecoCocktail.C

   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 /// \ingroup macros
  19 /// \file DecodeRecoCocktail.C
  20 /// \brief add brief description
  21 ///
  22 /// \author A. De Falco, H. Woehri, INFN Cagliari, July 2006
  23 ///
  24 /// This macro reads the generation/reconstruction files in the
  25 /// input directories (recodir and simdir),
  26 /// builds a tree that contains, for each event, an array of muons and dimuons
  27 /// (TClonesArrays of AliMUONTrackLight and AliMUONPairLight objects)
  28 /// and writes the tree in an output file (outFileName)
  29 /// Note that if the path for the output file is not explicitly specified,
  30 /// it will be written in the directory containing the generation/reconstruction
  31 /// 27-Nov-2006: modified by in order to loop on files
  32 ///
  33 /// 13 Nov 2007:
  34 /// Updated this macro to work with new version of AliMUONRecoCheck. Also, we are
  35 /// now fetching reconstructed track data from ESD and not the track tree, because
  36 /// the AliMUONTrack objects for reconstructed tracks are no longer stored on disk.
  37 ///  - Artur Szostak <artursz@iafrica.com>
  38
  39 #if !defined(__CINT__) || defined(__MAKECINT__)
  40 #include <Riostream.h>
  41 #include <TClonesArray.h>
  42 #include <TFile.h>
  43 #include <TTree.h>
  44 #include <TMath.h>
  45 #include <TLorentzVector.h>
  46 #include <TParticle.h>
  47 #include <TSystem.h>
  48 #include "AliMUONRecoCheck.h"
  49 #include "AliMUONTrack.h"
  50 #include "AliMUONTrackLight.h"
  51 #include "AliMUONPairLight.h"
  52 #include "AliMUONVTrackStore.h"
  53 #include "AliESDEvent.h"
  54 #include "AliESDVertex.h"
  55 #include "AliESDMuonTrack.h"
  56 #include "AliMCEventHandler.h"
  57 #include "AliMCEvent.h"
  58 #include "AliStack.h"
  59 /*TODO: need to update this with changes made to ITS
  60 #include "AliITSVertexerPPZ.h"
  61 #include "AliITSLoader.h"
  62 */
  63 #endif
  64
  65
  66 void DecodeRecoCocktail(
  67     char* recodir=".",          // The directory containing galice.root for reconstructed data.
  68     char* simdir="generated/",  // The directory containing galice.root for simulated data.
  69     char* outFileName = "MuonLight.root" // The output filename containing AliMUONTrackLight and AliMUONPairLight objects.
  70   )
  71 {
  72 /// \param recodir      The directory containing galice.root for reconstructed data.
  73 /// \param simdir       The directory containing galice.root for simulated data.
  74 /// \param outFileName  The output filename containing AliMUONTrackLight and AliMUONPairLight objects.
  75
  76   char startingDir[200];
  77   sprintf (startingDir,"%s",gSystem->pwd());
  78   gSystem->cd(recodir);
  79
  80   TClonesArray *muonArray   = new TClonesArray("AliMUONTrackLight",10);
  81   TClonesArray *dimuonArray = new TClonesArray("AliMUONPairLight",10);
  82   TTree *treeOut = new TTree("tree","tree");
  83
  84   TFile *fout = new TFile(outFileName,"recreate");
  85   fout->cd();
  86
  87   treeOut->Branch("muons",&muonArray);
  88   treeOut->Branch("dimuons",&dimuonArray);
  89
  90   AliMUONRecoCheck *rc = new AliMUONRecoCheck("AliESDs.root", simdir);
  91
  92   /*TODO: need to update this with changes made to ITS
  93   AliITSLoader* ITSloader =  (AliITSLoader*) runLoaderSim->GetLoader("ITSLoader");
  94   AliITSVertexerPPZ *dovert = 0;
  95   if (ITSloader) {
  96     dovert = new AliITSVertexerPPZ("default",0,0);
  97     // dovert->SetDebug(0);
  98     dovert->SetDiffPhiMax(0.05);
  99     dovert->SetWindow(3.);
 100   }
 101   AliESDVertex *vert = 0;
 102   */
 103
 104   TLorentzVector v;
 105
 106   Int_t nev = rc->NumberOfEvents();
 107   for(Int_t ievent = 0; ievent < nev; ievent++){ // loop over events
 108
 109     /*TODO: need to update this with changes made to ITS
 110     runLoaderSim->GetHeader();
 111     if (ITSloader) {
 112       vert = dovert->FindVertexForCurrentEvent(ievent);
 113     }
 114     */
 115
 116     //printf ("Event %d of %d\n",ievent,nev);
 117     muonArray->Clear();     // clean muon and dimuon arrays
 118     dimuonArray->Clear();
 119
 120     AliMUONVTrackStore* recoTracks = rc->ReconstructedTracks(ievent, kFALSE);  // Use tracks from actual reconstruction.
 121     //AliMUONVTrackStore* trackRefs = rc->ReconstructibleTracks(ievent);  // Only use reconstructible reference tracks.
 122     AliMUONVTrackStore* trackRefs = rc->TrackRefs(ievent);
 123     AliStack* pstack = (const_cast<AliMCEventHandler*>(rc->GetMCEventHandler()))->MCEvent()->Stack();
 124
 125     // loop over ESD tracks
 126     Int_t nreftracks = 0;
 127     const AliESDEvent* esd = rc->GetESDEvent();
 128     Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks() ;
 129     for (Int_t iTrack = 0; iTrack <  nTracks;  iTrack++) {
 130
 131       AliESDMuonTrack* esdTrack = esd->GetMuonTrack(iTrack);
 132
 133       // skip ghosts
 134       if (!esdTrack->ContainTrackerData()) continue;
 135
 136       // find the corresponding MUON track
 137       AliMUONTrack* trackReco = (AliMUONTrack*) recoTracks->FindObject(esdTrack->GetUniqueID());
 138
 139       // try to match the reconstructed track with a simulated one
 140       Int_t nMatchClusters = 0;
 141       AliMUONTrack* matchedTrackRef = rc->FindCompatibleTrack(*trackReco, *trackRefs, nMatchClusters, kFALSE, 10.);
 142
 143       if (matchedTrackRef) {
 144
 145         //store new referenced track in the muonArray
 146         AliMUONTrackLight* muLight = new ((*muonArray)[nreftracks++]) AliMUONTrackLight();
 147
 148         // assign parameters concerning the reconstructed tracks
 149         muLight->FillFromESD(esdTrack);
 150         // muLight->FillFromAliMUONTrack(trackReco);
 151
 152         // store further information related to the simulated track
 153         muLight->SetTrackPythiaLine(matchedTrackRef->GetUniqueID());
 154         TParticle *part = pstack->Particle(matchedTrackRef->GetUniqueID());
 155         muLight->SetTrackPDGCode(part->GetPdgCode());
 156         v.SetPxPyPzE(part->Px(), part->Py(), part->Pz(), part->Energy());
 157         muLight->SetPGen(v);
 158         muLight->FillMuonHistory(pstack, part);
 159         //        muLight.PrintInfo("A");
 160
 161       }
 162
 163     } // end esd track
 164
 165     // now loop over muon pairs to build dimuons
 166     Int_t nmuons = muonArray->GetEntriesFast();
 167     Int_t ndimuons = 0;
 168     for(Int_t itrRec1 = 0; itrRec1 < nmuons-1; itrRec1++) {
 169       AliMUONTrackLight* mu1 = (AliMUONTrackLight*) muonArray->At(itrRec1);
 170       for(Int_t itrRec2 = itrRec1+1; itrRec2 < nmuons; itrRec2++){
 171         AliMUONTrackLight* mu2 = (AliMUONTrackLight*) muonArray->At(itrRec2);
 172         AliMUONPairLight dimuLight;
 173         dimuLight.SetMuons(*mu1, *mu2);
 174         TClonesArray &dimuons = *dimuonArray;
 175         new (dimuons[ndimuons++]) AliMUONPairLight(dimuLight);
 176       }
 177     }
 178
 179     treeOut->Fill();
 180   }
 181
 182   fout->cd();
 183   treeOut->Write();
 184   gSystem->cd(startingDir);
 185   fout->Close();
 186   delete muonArray;
 187   delete dimuonArray;
 188   delete treeOut;
 189   delete fout;
 190   delete rc;
 191 }
 192