af9df62e |
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 | |
e54bf126 |
18 | /// \ingroup macros |
19 | /// \file MUONAlignment.C |
20 | /// \brief Macro for MUON alignment using physics tracks. |
21 | /// |
22 | /// The macro uses the AliMUONAlignment class to calculate the alignment parameters. |
23 | /// An array for the alignment parameters is created and can be filled with |
24 | /// initial values that will be used as starting values by the alignment |
25 | /// algorithm. |
26 | /// |
27 | /// By default the macro run over galice.root in the working directory. If a file list |
28 | /// of galice.root is provided as third argument the macro will run over all files. |
29 | /// The macro loop over the files, events and tracks. For each track |
30 | /// AliMUONAlignment::ProcessTrack(AliMUONTrack * track) and then |
31 | /// AliMUONAlignment::LocalFit(Int_t iTrack, Double_t *lTrackParam, Int_t |
32 | /// lSingleFit) are called. After all tracks have been procesed and fitted |
33 | /// AliMUONAlignment::GlobalFit(Double_t *parameters,Double_t *errors,Double_t *pulls) |
34 | /// is called. The array parameters contains the obatained misalignement parameters. |
35 | /// A realigned geometry is generated in a local CDB. |
36 | /// |
37 | /// \author: B. Becker and J. Castillo |
af9df62e |
38 | |
2688ea9b |
39 | #if !defined(__CINT__) || defined(__MAKECINT__) |
40 | |
41 | #include "AliMUONAlignment.h" |
42 | #include "AliMUONTrack.h" |
43 | #include "AliMUONTrackExtrap.h" |
44 | #include "AliMUONTrackParam.h" |
45 | #include "AliMUONGeometryTransformer.h" |
103e6575 |
46 | #include "AliMUONESDInterface.h" |
47 | |
61fed964 |
48 | #include "AliESDEvent.h" |
49 | #include "AliESDMuonTrack.h" |
2688ea9b |
50 | #include "AliMagFMaps.h" |
51 | #include "AliTracker.h" |
52 | #include "AliCDBManager.h" |
53 | #include "AliCDBMetaData.h" |
54 | #include "AliCDBId.h" |
a02fe82d |
55 | #include "AliGeomManager.h" |
2688ea9b |
56 | |
57 | #include <TString.h> |
2688ea9b |
58 | #include <TError.h> |
59 | #include <TH1.h> |
60 | #include <TGraphErrors.h> |
61 | #include <TFile.h> |
61fed964 |
62 | #include <TTree.h> |
2688ea9b |
63 | #include <TClonesArray.h> |
64 | #include <Riostream.h> |
65 | |
66 | #include <fstream> |
67 | |
68 | #endif |
69 | |
61fed964 |
70 | void MUONAlignment(Int_t nEvents = 100000, char* geoFilename = "geometry.root", TString esdFileName = "AliESDs.root", TString fileList = "") |
af9df62e |
71 | { |
72 | |
59830a87 |
73 | // Import TGeo geometry (needed by AliMUONTrackExtrap::ExtrapToVertex) |
a02fe82d |
74 | if ( ! AliGeomManager::GetGeometry() ) { |
75 | AliGeomManager::LoadGeometry(geoFilename); |
76 | if (! AliGeomManager::GetGeometry() ) { |
77 | Error("MUONAlignment", "getting geometry from file %s failed", geoFilename); |
59830a87 |
78 | return; |
79 | } |
80 | } |
81 | |
82 | // set mag field |
83 | // waiting for mag field in CDB |
84 | printf("Loading field map...\n"); |
85 | AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 1, 1., 10., AliMagFMaps::k5kG); |
86 | AliTracker::SetFieldMap(field, kFALSE); |
87 | // set the magnetic field for track extrapolations |
88 | AliMUONTrackExtrap::SetField(AliTracker::GetFieldMap()); |
89 | |
af9df62e |
90 | Double_t parameters[3*156]; |
91 | Double_t errors[3*156]; |
92 | Double_t pulls[3*156]; |
93 | for(Int_t k=0;k<3*156;k++) { |
94 | parameters[k]=0.; |
95 | errors[k]=0.; |
96 | pulls[k]=0.; |
97 | } |
98 | |
99 | Double_t trackParams[8] = {0.,0.,0.,0.,0.,0.,0.,0.}; |
100 | |
af9df62e |
101 | // Set initial values here, good guess may help convergence |
102 | // St 1 |
2688ea9b |
103 | // Int_t iPar = 0; |
104 | // parameters[iPar++] = 0.010300 ; parameters[iPar++] = 0.010600 ; parameters[iPar++] = 0.000396 ; |
af9df62e |
105 | |
106 | bool bLoop = kFALSE; |
107 | ifstream sFileList; |
108 | if (fileList.Contains(".list")) { |
109 | cout << "Reading file list: " << fileList.Data() << endl; |
110 | bLoop = kTRUE; |
111 | |
112 | TString fullListName("."); |
113 | fullListName +="/"; |
114 | fullListName +=fileList; |
115 | sFileList.open(fileList.Data()); |
116 | } |
117 | |
118 | TH1F *fInvBenMom = new TH1F("fInvBenMom","fInvBenMom",200,-0.1,0.1); |
119 | TH1F *fBenMom = new TH1F("fBenMom","fBenMom",200,-40,40); |
120 | |
121 | AliMUONAlignment* alig = new AliMUONAlignment(); |
122 | alig->InitGlobalParameters(parameters); |
123 | |
a02fe82d |
124 | AliMUONGeometryTransformer *transform = new AliMUONGeometryTransformer(); |
125 | transform->LoadGeometryData(); |
af9df62e |
126 | alig->SetGeometryTransformer(transform); |
7f8c1308 |
127 | |
128 | // Do alignment with magnetic field off |
129 | alig->SetBFieldOn(kFALSE); |
587965bd |
130 | |
131 | // Set tracking station to use |
132 | Bool_t bStOnOff[5] = {kTRUE,kTRUE,kTRUE,kTRUE,kTRUE}; |
133 | |
134 | // Fix parameters or add constraints here |
135 | for (Int_t iSt=0; iSt<5; iSt++) |
136 | if (!bStOnOff[iSt]) alig->FixStation(iSt+1); |
137 | |
138 | // Left and right sides of the detector are independent, one can choose to align |
139 | // only one side |
140 | Bool_t bSpecLROnOff[2] = {kTRUE,kTRUE}; |
141 | alig->FixHalfSpectrometer(bStOnOff,bSpecLROnOff); |
142 | |
143 | // Set predifined global constrains: X, Y, P, XvsZ, YvsZ, PvsZ, XvsY, YvsY, PvsY |
144 | Bool_t bVarXYT[9] = {kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE}; |
145 | Bool_t bDetTLBR[4] = {kFALSE,kTRUE,kFALSE,kTRUE}; |
146 | alig->AddConstraints(bStOnOff,bVarXYT,bDetTLBR,bSpecLROnOff); |
147 | |
af9df62e |
148 | |
149 | char cFileName[100]; |
7f8c1308 |
150 | |
af9df62e |
151 | Int_t lMaxFile = 1000; |
152 | Int_t iFile = 0; |
587965bd |
153 | Int_t iEvent = 0; |
af9df62e |
154 | bool bKeepLoop = kTRUE; |
587965bd |
155 | Int_t iTrackTot=0; |
156 | Int_t iTrackOk=0; |
157 | |
af9df62e |
158 | while(bKeepLoop && iFile<lMaxFile){ |
159 | iFile++; |
160 | if (bLoop) { |
161 | sFileList.getline(cFileName,100); |
162 | if (sFileList.eof()) bKeepLoop = kFALSE; |
163 | } |
164 | else { |
61fed964 |
165 | sprintf(cFileName,esdFileName.Data()); |
af9df62e |
166 | bKeepLoop = kFALSE; |
167 | } |
61fed964 |
168 | if (!strstr(cFileName,"AliESDs.root")) continue; |
af9df62e |
169 | cout << "Using file: " << cFileName << endl; |
61fed964 |
170 | |
171 | // load ESD event |
172 | TFile* esdFile = TFile::Open(cFileName); // open the file |
173 | if (!esdFile || !esdFile->IsOpen()) { |
174 | cout << "opening ESD file " << cFileName << "failed" << endl; |
175 | continue; |
176 | } |
177 | TTree* esdTree = (TTree*) esdFile->Get("esdTree"); // get the tree |
178 | if (!esdTree) { |
179 | cout << "no ESD tree found" << endl; |
180 | esdFile->Close(); |
181 | continue; |
182 | } |
183 | AliESDEvent* esdEvent = new AliESDEvent(); // link ESD event to the tree |
184 | esdEvent->ReadFromTree(esdTree); |
7f8c1308 |
185 | |
61fed964 |
186 | Int_t nevents = esdTree->GetEntries(); |
587965bd |
187 | cout << "... with " << nevents << endl; |
af9df62e |
188 | for(Int_t event = 0; event < nevents; event++) { |
587965bd |
189 | if (iEvent >= nEvents){ |
190 | bKeepLoop = kFALSE; |
191 | break; |
192 | } |
193 | iEvent++; |
194 | |
61fed964 |
195 | if (esdTree->GetEvent(event) <= 0) { |
196 | cout << "fails to read ESD object for event " << event << endl; |
197 | continue; |
198 | } |
7f8c1308 |
199 | |
61fed964 |
200 | Int_t nTracks = Int_t(esdEvent->GetNumberOfMuonTracks()); |
201 | if (!event%100) cout << " there are " << nTracks << " tracks in event " << event << endl; |
202 | for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) { |
203 | AliESDMuonTrack* esdTrack = esdEvent->GetMuonTrack(iTrack); |
204 | if (!esdTrack->ClustersStored()) continue; |
205 | Double_t invBenMom = esdTrack->GetInverseBendingMomentum(); |
af9df62e |
206 | fInvBenMom->Fill(invBenMom); |
207 | fBenMom->Fill(1./invBenMom); |
208 | if (TMath::Abs(invBenMom)<=1.04) { |
103e6575 |
209 | AliMUONTrack track; |
210 | AliMUONESDInterface::ESDToMUON(*esdTrack, track); |
61fed964 |
211 | alig->ProcessTrack(&track); |
af9df62e |
212 | alig->LocalFit(iTrackOk++,trackParams,0); |
213 | } |
587965bd |
214 | iTrackTot++; |
af9df62e |
215 | } |
216 | } |
61fed964 |
217 | delete esdEvent; |
218 | esdFile->Close(); |
587965bd |
219 | cout << "Processed " << iTrackTot << " Tracks so far." << endl; |
af9df62e |
220 | } |
221 | alig->GlobalFit(parameters,errors,pulls); |
222 | |
223 | cout << "Done with GlobalFit " << endl; |
224 | |
225 | // Store results |
226 | Double_t DEid[156] = {0}; |
227 | Double_t MSDEx[156] = {0}; |
228 | Double_t MSDEy[156] = {0}; |
229 | Double_t MSDExt[156] = {0}; |
230 | Double_t MSDEyt[156] = {0}; |
231 | Double_t DEidErr[156] = {0}; |
232 | Double_t MSDExErr[156] = {0}; |
233 | Double_t MSDEyErr[156] = {0}; |
234 | Double_t MSDExtErr[156] = {0}; |
235 | Double_t MSDEytErr[156] = {0}; |
236 | Int_t lNDetElem = 4*2+4*2+18*2+26*2+26*2; |
237 | Int_t lNDetElemCh[10] = {4,4,4,4,18,18,26,26,26,26}; |
2688ea9b |
238 | // Int_t lSNDetElemCh[10] = {4,8,12,16,34,52,78,104,130,156}; |
af9df62e |
239 | Int_t idOffset = 0; // 400 |
240 | Int_t lSDetElemCh = 0; |
241 | for(Int_t iDE=0; iDE<lNDetElem; iDE++){ |
242 | DEidErr[iDE] = 0.; |
243 | DEid[iDE] = idOffset+100; |
244 | DEid[iDE] += iDE; |
245 | lSDetElemCh = 0; |
246 | for(Int_t iCh=0; iCh<9; iCh++){ |
247 | lSDetElemCh += lNDetElemCh[iCh]; |
248 | if (iDE>=lSDetElemCh) { |
249 | DEid[iDE] += 100; |
250 | DEid[iDE] -= lNDetElemCh[iCh]; |
251 | } |
252 | } |
253 | MSDEx[iDE]=parameters[3*iDE+0]; |
254 | MSDEy[iDE]=parameters[3*iDE+1]; |
255 | MSDExt[iDE]=parameters[3*iDE+2]; |
256 | MSDEyt[iDE]=parameters[3*iDE+2]; |
257 | MSDExErr[iDE]=(Double_t)alig->GetParError(3*iDE+0); |
258 | MSDEyErr[iDE]=(Double_t)alig->GetParError(3*iDE+1); |
259 | MSDExtErr[iDE]=(Double_t)alig->GetParError(3*iDE+2); |
260 | MSDEytErr[iDE]=(Double_t)alig->GetParError(3*iDE+2); |
261 | } |
262 | |
263 | cout << "Let's create graphs ... " << endl; |
264 | |
265 | TGraphErrors *gMSDEx = new TGraphErrors(lNDetElem,DEid,MSDEx,DEidErr,MSDExErr); |
266 | TGraphErrors *gMSDEy = new TGraphErrors(lNDetElem,DEid,MSDEy,DEidErr,MSDEyErr); |
267 | TGraphErrors *gMSDExt = new TGraphErrors(lNDetElem,DEid,MSDExt,DEidErr,MSDExtErr); |
268 | TGraphErrors *gMSDEyt = new TGraphErrors(lNDetElem,DEid,MSDEyt,DEidErr,MSDEytErr); |
269 | |
270 | cout << "... graphs created, open file ... " << endl; |
271 | |
272 | TFile *hFile = new TFile("measShifts.root","RECREATE"); |
273 | |
274 | cout << "... file opened ... " << endl; |
275 | |
276 | gMSDEx->Write("gMSDEx"); |
277 | gMSDEy->Write("gMSDEy"); |
278 | gMSDExt->Write("gMSDExt"); |
279 | gMSDEyt->Write("gMSDEyt"); |
280 | fInvBenMom->Write(); |
281 | fBenMom->Write(); |
282 | hFile->Close(); |
283 | |
284 | cout << "... and closed!" << endl; |
285 | // Re Align |
286 | AliMUONGeometryTransformer *newTransform = alig->ReAlign(transform,parameters,true); |
287 | newTransform->WriteTransformations("transform2ReAlign.dat"); |
288 | |
289 | // Generate realigned data in local cdb |
2688ea9b |
290 | const TClonesArray* array = newTransform->GetMisAlignmentData(); |
af9df62e |
291 | |
292 | // CDB manager |
293 | AliCDBManager* cdbManager = AliCDBManager::Instance(); |
294 | cdbManager->SetDefaultStorage("local://ReAlignCDB"); |
295 | |
296 | AliCDBMetaData* cdbData = new AliCDBMetaData(); |
297 | cdbData->SetResponsible("Dimuon Offline project"); |
298 | cdbData->SetComment("MUON alignment objects with residual misalignment"); |
299 | AliCDBId id("MUON/Align/Data", 0, 0); |
2688ea9b |
300 | cdbManager->Put(const_cast<TClonesArray*>(array), id, cdbData); |
af9df62e |
301 | |
61fed964 |
302 | } |
303 | |