Coding conventions
[u/mrichter/AliRoot.git] / MUON / DecodeRecoCocktail.C
CommitLineData
fa6e7866 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 directory (dirname),
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
4d309390 26// 27-Nov-2006: modified by in order to loop on files
fa6e7866 27
4d309390 28#if !defined(__CINT__) || defined(__MAKECINT__)
29#include <Riostream.h>
fa6e7866 30#include <TClonesArray.h>
31#include <TFile.h>
32#include <TTree.h>
33#include <TMath.h>
34#include <TLorentzVector.h>
35#include <TParticle.h>
36#include <TSystem.h>
59830a87 37#include <TGeoManager.h>
4d309390 38#include "AliMUONRecoCheck.h"
39#include "AliMUONTrack.h"
40#include "AliMUONTrackParam.h"
41#include "AliRunLoader.h"
fa6e7866 42#include "AliMUONTrackLight.h"
43#include "AliMUONPairLight.h"
22ccc301 44#include "AliMUONTrackExtrap.h"
4d309390 45#include "AliESD.h"
46#include "AliESDVertex.h"
47#include "AliITSVertexerPPZ.h"
48#include "AliITSLoader.h"
59830a87 49#include "AliTracker.h"
50#include "AliMagFMaps.h"
4d309390 51#endif
fa6e7866 52
4d309390 53
54void DecodeRecoCocktail(char* dirname=".", char* outFileName = "MuonLight.root", char* geoFilename = "geometry.root"){
55
56
57 char startingDir[200];
58 sprintf (startingDir,"%s",gSystem->pwd());
fa6e7866 59 gSystem->cd(dirname);
4d309390 60
61
62 TClonesArray *muonArray = new TClonesArray("AliMUONTrackLight",10);
63 TClonesArray *dimuonArray = new TClonesArray("AliMUONPairLight",10);
64 TTree *treeOut = new TTree("tree","tree");
65
fa6e7866 66 TFile *fout = new TFile(outFileName,"recreate");
4d309390 67 fout->cd();
fa6e7866 68
4d309390 69 treeOut->Branch("muons",&muonArray);
70 treeOut->Branch("dimuons",&dimuonArray);
71
72 TFile* esdFile = TFile::Open("AliESDs.root");
73 if (!esdFile || !esdFile->IsOpen()) {
21939432 74 Error("DecodeRecoCocktailNew", "opening ESD file AliESDs.root failed");
4d309390 75 return;
76 }
77
78 AliESD* esd = new AliESD();
79 TTree* treeESD = (TTree*) esdFile->Get("esdTree");
80 if (!treeESD) {
81 Error("CheckESD", "no ESD tree found");
82 return;
83 }
84 treeESD->SetBranchAddress("ESD", &esd);
85
59830a87 86 // Import TGeo geometry (needed by AliMUONTrackExtrap::ExtrapToVertex)
87 if (!gGeoManager) {
88 TGeoManager::Import(geoFilename);
89 if (!gGeoManager) {
90 Error("MUONmass_ESD", "getting geometry from file %s failed", geoFilename);
91 return;
92 }
93 }
94
95 // set mag field
96 // waiting for mag field in CDB
97 printf("Loading field map...\n");
98 AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 1, 1., 10., AliMagFMaps::k5kG);
99 AliTracker::SetFieldMap(field, kFALSE);
100 // set the magnetic field for track extrapolations
101 AliMUONTrackExtrap::SetField(AliTracker::GetFieldMap());
102
a4ee7ab9 103 AliMUONRecoCheck *rc = new AliMUONRecoCheck("galice.root", "galice_sim.root");
fa6e7866 104 AliRunLoader *runLoader = rc->GetRunLoader();
a4ee7ab9 105 AliRunLoader *runLoaderSim = rc->GetRunLoaderSim();
fa6e7866 106
a4ee7ab9 107 runLoaderSim->LoadKinematics("READ");
108 Int_t nev = runLoaderSim->GetNumberOfEvents();
4d309390 109// nevent = nevent +nev;
110// printf(" number of files and events = %d - %d \n",irun,nevent);
111
a4ee7ab9 112 AliITSLoader* ITSloader = (AliITSLoader*) runLoaderSim->GetLoader("ITSLoader");
4d309390 113 AliITSVertexerPPZ *dovert = 0;
114 if (ITSloader) {
115 dovert = new AliITSVertexerPPZ("default",0,0);
21939432 116 // dovert->SetDebug(0);
4d309390 117 dovert->SetDiffPhiMax(0.05);
118 dovert->SetWindow(3.);
119 }
120 AliESDVertex *vert = 0;
fa6e7866 121
122 TLorentzVector v;
4d309390 123
fa6e7866 124 for(Int_t ievent = 0; ievent < nev; ievent++){ // loop over events
a4ee7ab9 125 runLoaderSim->GetHeader();
4d309390 126 if (ITSloader) {
127 vert = dovert->FindVertexForCurrentEvent(ievent);
128 }
a4ee7ab9 129 //printf ("Event %d of %d\n",ievent,nev);
fa6e7866 130 muonArray->Clear(); // clean muon and dimuon arrays
131 dimuonArray->Clear();
132 runLoader->GetEvent(ievent);
a4ee7ab9 133 runLoaderSim->GetEvent(ievent);
4d309390 134 treeESD->GetEvent(ievent);
fa6e7866 135 rc->ResetTracks();
136 rc->MakeTrackRef(); // make reconstructable tracks
137
138 TClonesArray * trackRecoArray = rc->GetTrackReco();
139 TClonesArray * trackRefArray = rc->GetMuonTrackRef();
140 Int_t nTrackReco = trackRecoArray->GetEntriesFast();
4d309390 141 Int_t nTracksESD = (Int_t)esd->GetNumberOfMuonTracks() ;
142 if (nTrackReco != nTracksESD) printf ("Tracks in AliMUONTrack and in ESD do not match!\n");
fa6e7866 143 Int_t nreftracks = 0;
144 for (Int_t itrRec = 0; itrRec<nTrackReco; itrRec++) {
145 // assign parameters concerning the reconstructed tracks
4d309390 146 AliMUONTrackLight muLight;
fa6e7866 147 AliMUONTrack *trackReco = (AliMUONTrack *)trackRecoArray->At(itrRec);
a4ee7ab9 148
4d309390 149 muLight.FillFromESD(esd->GetMuonTrack(itrRec));
150 // muLight.FillFromAliMUONTrack(trackReco);
151
152 // find the reference track and store further information
a4ee7ab9 153 TParticle *part = muLight.FindRefTrack(trackReco,trackRefArray,runLoaderSim);
fa6e7866 154 if (part) {
155 v.SetPxPyPzE(part->Px(), part->Py(), part->Pz(), part->Energy());
82724a77 156 muLight.SetPGen(v);
a4ee7ab9 157 muLight.FillMuonHistory(runLoaderSim, part);
4d309390 158 // muLight.PrintInfo("A");
fa6e7866 159 //store the referenced track in the muonArray:
160 TClonesArray &muons = *muonArray;
82724a77 161 new (muons[nreftracks++]) AliMUONTrackLight(muLight);
4d309390 162 }
163 } // end reco track
fa6e7866 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);
4d309390 172 AliMUONPairLight dimuLight;
82724a77 173 dimuLight.SetMuons(*mu1, *mu2);
fa6e7866 174 TClonesArray &dimuons = *dimuonArray;
82724a77 175 new (dimuons[ndimuons++]) AliMUONPairLight(dimuLight);
fa6e7866 176 }
177 }
fa6e7866 178 treeOut->Fill();
4d309390 179 }
180
fa6e7866 181 fout->cd();
182 treeOut->Write();
183 gSystem->cd(startingDir);
82724a77 184 fout->Close();
82724a77 185 delete muonArray;
186 delete dimuonArray;
4d309390 187 delete treeOut;
188 delete fout;
189 delete rc;
190 delete esd;
fa6e7866 191}
4d309390 192