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 // A. De Falco, H. Woehri, INFN Cagliari, July 2006
  19 // This macro reads the generation/reconstruction files in the
  20 // input directories (recodir and simdir),
  21 // builds a tree that contains, for each event, an array of muons and dimuons
  22 // (TClonesArrays of AliMUONTrackLight and AliMUONPairLight objects)
  23 // and writes the tree in an output file (outFileName)
  24 // Note that if the path for the output file is not explicitly specified,
  25 // it will be written in the directory containing the generation/reconstruction
  26 // 27-Nov-2006: modified by in order to loop on files
  27
  28 // 13 Nov 2007:
  29 // Updated this macro to work with new version of AliMUONRecoCheck. Also, we are
  30 // now fetching reconstructed track data from ESD and not the track tree, because
  31 // the AliMUONTrack objects for reconstructed tracks are no longer stored on disk.
  32 //  - Artur Szostak <artursz@iafrica.com>
  33
  34 #if !defined(__CINT__) || defined(__MAKECINT__)
  35 #include <Riostream.h>
  36 #include <TClonesArray.h>
  37 #include <TFile.h>
  38 #include <TTree.h>
  39 #include <TMath.h>
  40 #include <TLorentzVector.h>
  41 #include <TParticle.h>
  42 #include <TSystem.h>
  43 #include <TGeoManager.h>
  44 #include "AliMUONRecoCheck.h"
  45 #include "AliMUONTrack.h"
  46 #include "AliMUONTrackParam.h"
  47 #include "AliMUONTrackLight.h"
  48 #include "AliMUONPairLight.h"
  49 #include "AliMUONVTrackStore.h"
  50 #include "AliMUONTrackExtrap.h"
  51 #include "AliESDEvent.h"
  52 #include "AliESDVertex.h"
  53 #include "AliMCEventHandler.h"
  54 #include "AliMCEvent.h"
  55 /*TODO: need to update this with changes made to ITS
  56 #include "AliITSVertexerPPZ.h"
  57 #include "AliITSLoader.h"
  58 */
  59 #include "AliTracker.h"
  60 #include "AliMagFMaps.h"
  61 #endif
  62
  63
  64 void DecodeRecoCocktail(
  65     char* recodir=".",          // The directory containing galice.root for reconstructed data.
  66     char* simdir="generated/",  // The directory containing galice.root for simulated data.
  67     char* outFileName = "MuonLight.root", // The output filename containing AliMUONTrackLight and AliMUONPairLight objects.
  68     char* geoFilename = "generated/geometry.root"  // The filename containing the geometry.
  69   )
  70 {
  71   char startingDir[200];
  72   sprintf (startingDir,"%s",gSystem->pwd());
  73   gSystem->cd(recodir);
  74
  75   TClonesArray *muonArray   = new TClonesArray("AliMUONTrackLight",10);
  76   TClonesArray *dimuonArray = new TClonesArray("AliMUONPairLight",10);
  77   TTree *treeOut = new TTree("tree","tree");
  78
  79   TFile *fout = new TFile(outFileName,"recreate");
  80   fout->cd();
  81
  82   treeOut->Branch("muons",&muonArray);
  83   treeOut->Branch("dimuons",&dimuonArray);
  84
  85   TFile* esdFile = TFile::Open("AliESDs.root");
  86   if (!esdFile || !esdFile->IsOpen()) {
  87     Error("DecodeRecoCocktailNew", "opening ESD file AliESDs.root failed");
  88     return;
  89   }
  90
  91   AliESDEvent* esd = new AliESDEvent();
  92   TTree* treeESD = (TTree*) esdFile->Get("esdTree");
  93   if (!treeESD) {
  94     Error("CheckESD", "no ESD tree found");
  95     return;
  96   }
  97   esd->ReadFromTree(treeESD);
  98
  99   // Import TGeo geometry (needed by AliMUONTrackExtrap::ExtrapToVertex)
 100   if (!gGeoManager) {
 101     TGeoManager::Import(geoFilename);
 102     if (!gGeoManager) {
 103       Error("MUONmass_ESD", "getting geometry from file %s failed", geoFilename);
 104       return;
 105     }
 106   }
 107
 108   // set  mag field
 109   // waiting for mag field in CDB
 110   printf("Loading field map...\n");
 111   AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 1, 1., 10., AliMagFMaps::k5kG);
 112   AliTracker::SetFieldMap(field, kFALSE);
 113   // set the magnetic field for track extrapolations
 114   AliMUONTrackExtrap::SetField(AliTracker::GetFieldMap());
 115
 116   AliMUONRecoCheck *rc = new AliMUONRecoCheck("galice.root", simdir);
 117   Int_t nev = rc->NumberOfEvents();
 118
 119   /*TODO: need to update this with changes made to ITS
 120   AliITSLoader* ITSloader =  (AliITSLoader*) runLoaderSim->GetLoader("ITSLoader");
 121   AliITSVertexerPPZ *dovert = 0;
 122   if (ITSloader) {
 123     dovert = new AliITSVertexerPPZ("default",0,0);
 124     // dovert->SetDebug(0);
 125     dovert->SetDiffPhiMax(0.05);
 126     dovert->SetWindow(3.);
 127   }
 128   AliESDVertex *vert = 0;
 129   */
 130
 131   TLorentzVector v;
 132
 133   for(Int_t ievent = 0; ievent < nev; ievent++){ // loop over events
 134     treeESD->GetEvent(ievent);
 135     rc->GetMCEventHandler()->GetEvent(ievent);
 136     AliStack* pstack = rc->GetMCEventHandler()->MCEvent()->Stack();
 137
 138     /*TODO: need to update this with changes made to ITS
 139     runLoaderSim->GetHeader();
 140     if (ITSloader) {
 141       vert = dovert->FindVertexForCurrentEvent(ievent);
 142     }
 143     */
 144
 145     //printf ("Event %d of %d\n",ievent,nev);
 146     muonArray->Clear();     // clean muon and dimuon arrays
 147     dimuonArray->Clear();
 148
 149     //AliMUONVTrackStore* recoTracks = rc->ReconstructedTracks(ievent);  // Use tracks from actual reconstruction, but this does not work anymore.
 150     AliMUONVTrackStore* recoTracks = rc->ReconstructedTracks(esd);
 151     //AliMUONVTrackStore* trackRefs = rc->ReconstructibleTracks(ievent);  // Only use reconstructible reference tracks.
 152     AliMUONVTrackStore* trackRefs = rc->TrackRefs(ievent);
 153
 154     TIter next(recoTracks->CreateIterator());
 155     AliMUONTrack* trackReco = NULL;
 156
 157     Int_t nTrackReco = recoTracks->GetSize();
 158     Int_t nTracksESD = (Int_t)esd->GetNumberOfMuonTracks();
 159     if (nTrackReco != nTracksESD) printf ("Tracks in recoTracks (%d) and in ESD (%d) do not match!\n", nTrackReco, nTracksESD);
 160     Int_t nreftracks = 0;
 161     Int_t itrRec = 0;
 162     while ( (trackReco = static_cast<AliMUONTrack*>(next())) != NULL )
 163     {
 164       // assign parameters concerning the reconstructed tracks
 165       AliMUONTrackLight muLight;
 166
 167       muLight.FillFromESD(esd->GetMuonTrack(itrRec));
 168       // muLight.FillFromAliMUONTrack(trackReco);
 169
 170       // find the reference track and store further information
 171       TParticle *part = muLight.FindRefTrack(trackReco, trackRefs, pstack, kTRUE);
 172       if (part) {
 173         v.SetPxPyPzE(part->Px(), part->Py(), part->Pz(), part->Energy());
 174         muLight.SetPGen(v);
 175         muLight.FillMuonHistory(pstack, part);
 176         //        muLight.PrintInfo("A");
 177         //store the referenced track in the muonArray:
 178         TClonesArray &muons = *muonArray;
 179         new (muons[nreftracks++]) AliMUONTrackLight(muLight);
 180       }
 181
 182       itrRec++;
 183     }  // end reco track
 184
 185     // now loop over muon pairs to build dimuons
 186     Int_t nmuons = muonArray->GetEntriesFast();
 187     Int_t ndimuons = 0;
 188     for(Int_t itrRec1 = 0; itrRec1 < nmuons-1; itrRec1++) {
 189       AliMUONTrackLight* mu1 = (AliMUONTrackLight*) muonArray->At(itrRec1);
 190       for(Int_t itrRec2 = itrRec1+1; itrRec2 < nmuons; itrRec2++){
 191         AliMUONTrackLight* mu2 = (AliMUONTrackLight*) muonArray->At(itrRec2);
 192         AliMUONPairLight dimuLight;
 193         dimuLight.SetMuons(*mu1, *mu2);
 194         TClonesArray &dimuons = *dimuonArray;
 195         new (dimuons[ndimuons++]) AliMUONPairLight(dimuLight);
 196       }
 197     }
 198
 199     treeOut->Fill();
 200
 201     delete recoTracks;
 202   }
 203
 204   fout->cd();
 205   treeOut->Write();
 206   gSystem->cd(startingDir);
 207   fout->Close();
 208   delete muonArray;
 209   delete dimuonArray;
 210   delete treeOut;
 211   delete fout;
 212   delete rc;
 213   delete esd;
 214 }
 215