Class AliMUONCDB:
[u/mrichter/AliRoot.git] / 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 <TGeoManager.h>
49 #include "AliMUONCDB.h"
50 #include "AliMUONRecoCheck.h"
51 #include "AliMUONTrack.h"
52 #include "AliMUONTrackLight.h"
53 #include "AliMUONPairLight.h"
54 #include "AliMUONVTrackStore.h"
55 #include "AliMUONTrackExtrap.h"
56 #include "AliESDEvent.h"
57 #include "AliESDVertex.h"
58 #include "AliMCEventHandler.h"
59 #include "AliMCEvent.h"
60 #include "AliStack.h"
61 #include "AliCDBManager.h"
62 /*TODO: need to update this with changes made to ITS
63 #include "AliITSVertexerPPZ.h"
64 #include "AliITSLoader.h"
65 */
66 #endif
67
68
69 void DecodeRecoCocktail(
70     char* recodir=".",          // The directory containing galice.root for reconstructed data.
71     char* simdir="generated/",  // The directory containing galice.root for simulated data.
72     char* outFileName = "MuonLight.root", // The output filename containing AliMUONTrackLight and AliMUONPairLight objects.
73     char* geoFilename = "geometry.root"  // The filename containing the geometry.
74   )
75 {
76 /// \param recodir      The directory containing galice.root for reconstructed data.
77 /// \param simdir       The directory containing galice.root for simulated data.
78 /// \param outFileName  The output filename containing AliMUONTrackLight and AliMUONPairLight objects.
79 /// \param geoFilename  The filename containing the geometry.
80
81   char startingDir[200]; 
82   sprintf (startingDir,"%s",gSystem->pwd()); 
83   gSystem->cd(recodir); 
84
85   TClonesArray *muonArray   = new TClonesArray("AliMUONTrackLight",10); 
86   TClonesArray *dimuonArray = new TClonesArray("AliMUONPairLight",10); 
87   TTree *treeOut = new TTree("tree","tree"); 
88
89   TFile *fout = new TFile(outFileName,"recreate"); 
90   fout->cd(); 
91
92   treeOut->Branch("muons",&muonArray); 
93   treeOut->Branch("dimuons",&dimuonArray); 
94   
95   // Import TGeo geometry (needed by AliMUONTrackExtrap::ExtrapToVertex)
96   if (!gGeoManager) {
97     TGeoManager::Import(geoFilename);
98     if (!gGeoManager) {
99       Error("MUONmass_ESD", "getting geometry from file %s failed", geoFilename);
100       return;
101     }
102   }
103   
104   // load necessary data from OCDB
105   AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
106   AliCDBManager::Instance()->SetRun(0);
107   if (!AliMUONCDB::LoadField()) return;
108   // set the magnetic field for track extrapolations
109   AliMUONTrackExtrap::SetField();
110
111   AliMUONRecoCheck *rc = new AliMUONRecoCheck("AliESDs.root", simdir);
112   Int_t nev = rc->NumberOfEvents();
113   
114   /*TODO: need to update this with changes made to ITS
115   AliITSLoader* ITSloader =  (AliITSLoader*) runLoaderSim->GetLoader("ITSLoader");
116   AliITSVertexerPPZ *dovert = 0; 
117   if (ITSloader) { 
118     dovert = new AliITSVertexerPPZ("default",0,0);
119     // dovert->SetDebug(0);
120     dovert->SetDiffPhiMax(0.05);
121     dovert->SetWindow(3.);
122   }
123   AliESDVertex *vert = 0;
124   */
125   
126   TLorentzVector v; 
127  
128   for(Int_t ievent = 0; ievent < nev; ievent++){ // loop over events 
129     
130     /*TODO: need to update this with changes made to ITS
131     runLoaderSim->GetHeader();
132     if (ITSloader) { 
133       vert = dovert->FindVertexForCurrentEvent(ievent);
134     }
135     */
136     
137     //printf ("Event %d of %d\n",ievent,nev);
138     muonArray->Clear();     // clean muon and dimuon arrays 
139     dimuonArray->Clear(); 
140     
141     AliMUONVTrackStore* recoTracks = rc->ReconstructedTracks(ievent, kFALSE);  // Use tracks from actual reconstruction.
142     //AliMUONVTrackStore* trackRefs = rc->ReconstructibleTracks(ievent);  // Only use reconstructible reference tracks.
143     AliMUONVTrackStore* trackRefs = rc->TrackRefs(ievent);
144     AliStack* pstack = (const_cast<AliMCEventHandler*>(rc->GetMCEventHandler()))->MCEvent()->Stack();
145     
146     TIter next(recoTracks->CreateIterator());
147     AliMUONTrack* trackReco = NULL;
148     
149     Int_t nTrackReco = recoTracks->GetSize();
150     Int_t nTracksESD = rc->GetESDEvent()->GetNumberOfMuonTracks();
151     if (nTrackReco != nTracksESD) printf ("Tracks in recoTracks (%d) and in ESD (%d) do not match!\n", nTrackReco, nTracksESD);
152     Int_t nreftracks = 0;
153     Int_t itrRec = 0;
154     while ( (trackReco = static_cast<AliMUONTrack*>(next())) != NULL )
155     {
156       // assign parameters concerning the reconstructed tracks
157       AliMUONTrackLight muLight;
158       
159       muLight.FillFromESD(rc->GetESDEvent()->GetMuonTrack(itrRec));
160       // muLight.FillFromAliMUONTrack(trackReco);
161       
162       // find the reference track and store further information 
163       TParticle *part = muLight.FindRefTrack(trackReco, trackRefs, pstack);
164       if (part) { 
165         v.SetPxPyPzE(part->Px(), part->Py(), part->Pz(), part->Energy());
166         muLight.SetPGen(v); 
167         muLight.FillMuonHistory(pstack, part);
168         //        muLight.PrintInfo("A");
169         //store the referenced track in the muonArray:
170         TClonesArray &muons = *muonArray;
171         new (muons[nreftracks++]) AliMUONTrackLight(muLight);
172       } 
173       
174       itrRec++;
175     }  // end reco track
176     
177     // now loop over muon pairs to build dimuons
178     Int_t nmuons = muonArray->GetEntriesFast(); 
179     Int_t ndimuons = 0; 
180     for(Int_t itrRec1 = 0; itrRec1 < nmuons-1; itrRec1++) { 
181       AliMUONTrackLight* mu1 = (AliMUONTrackLight*) muonArray->At(itrRec1); 
182       for(Int_t itrRec2 = itrRec1+1; itrRec2 < nmuons; itrRec2++){
183         AliMUONTrackLight* mu2 = (AliMUONTrackLight*) muonArray->At(itrRec2); 
184         AliMUONPairLight dimuLight;
185         dimuLight.SetMuons(*mu1, *mu2);
186         TClonesArray &dimuons = *dimuonArray;
187         new (dimuons[ndimuons++]) AliMUONPairLight(dimuLight);
188       }
189     }
190     
191     treeOut->Fill(); 
192   } 
193   
194   fout->cd(); 
195   treeOut->Write(); 
196   gSystem->cd(startingDir); 
197   fout->Close();
198   delete muonArray;
199   delete dimuonArray;
200   delete treeOut; 
201   delete fout;
202   delete rc;
203 }
204