In alignment code:
[u/mrichter/AliRoot.git] / MUON / MUONAlignment.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 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
38
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"
46 #include "AliMUONESDInterface.h"
47
48 #include "AliESDEvent.h"
49 #include "AliESDMuonTrack.h"
50 #include "AliMagF.h"
51 #include "AliTracker.h"
52 #include "AliCDBManager.h"
53 #include "AliCDBMetaData.h"
54 #include "AliCDBId.h"
55 #include "AliGeomManager.h"
56
57 #include <TString.h>
58 #include <TError.h>
59 #include <TH1.h>
60 #include <TGraphErrors.h>
61 #include <TFile.h>
62 #include <TTree.h>
63 #include <TClonesArray.h>
64 #include <Riostream.h>
65
66 #include <fstream>
67
68 #endif
69
70 void MUONAlignment(Int_t nEvents = 100000, char* geoFilename = "geometry.root", TString esdFileName = "AliESDs.root", TString fileList = "")
71 {
72  
73   // Import TGeo geometry (needed by AliMUONTrackExtrap::ExtrapToVertex)
74   if ( ! AliGeomManager::GetGeometry() ) {
75     AliGeomManager::LoadGeometry(geoFilename);
76     if (! AliGeomManager::GetGeometry() ) {
77       Error("MUONAlignment", "getting geometry from file %s failed", geoFilename);
78       return;
79     }
80   }
81   
82   // set  mag field 
83   // waiting for mag field in CDB 
84   if (!TGeoGlobalMagField::Instance()->GetField()) {
85     printf("Loading field map...\n");
86     //    AliMagF* field = new AliMagF("Maps","Maps",2,1.,1., 10.,AliMagF::k5kG);
87     AliMagF* field = new AliMagF("Maps","Maps",2,0.,0., 10.,AliMagF::k5kG);
88     TGeoGlobalMagField::Instance()->SetField(field);
89   }
90   // set the magnetic field for track extrapolations
91   AliMUONTrackExtrap::SetField();
92
93   Double_t parameters[4*156];
94   Double_t errors[4*156];
95   Double_t pulls[4*156];
96   for(Int_t k=0;k<4*156;k++) {
97     parameters[k]=0.;
98     errors[k]=0.;
99     pulls[k]=0.;
100   }
101
102   Double_t trackParams[8] = {0.,0.,0.,0.,0.,0.,0.,0.};
103
104   // Set initial values here, good guess may help convergence
105   // St 1 
106   //  Int_t iPar = 0;
107   //  parameters[iPar++] =  0.010300 ;  parameters[iPar++] =  0.010600 ;  parameters[iPar++] =  0.000396 ;  
108
109   bool bLoop = kFALSE;
110   ifstream sFileList;
111   if (fileList.Contains(".list")) {
112     cout << "Reading file list: " << fileList.Data() << endl;
113     bLoop = kTRUE;
114
115     TString fullListName(".");
116     fullListName +="/";
117     fullListName +=fileList;
118     sFileList.open(fileList.Data());
119   }
120   
121   TH1F *fInvBenMom = new TH1F("fInvBenMom","fInvBenMom",200,-0.1,0.1); 
122   TH1F *fBenMom = new TH1F("fBenMom","fBenMom",200,-40,40); 
123
124   AliMUONAlignment* alig = new AliMUONAlignment();
125   alig->InitGlobalParameters(parameters);
126
127   AliMUONGeometryTransformer *transform = new AliMUONGeometryTransformer();
128   transform->LoadGeometryData();
129   alig->SetGeometryTransformer(transform);
130
131   // Do alignment with magnetic field off
132   alig->SetBFieldOn(kFALSE);
133   
134   // Set tracking station to use
135   Bool_t bStOnOff[5] = {kTRUE,kTRUE,kTRUE,kTRUE,kTRUE};
136   Bool_t bChOnOff[10] = {kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE};
137
138   // Set degrees of freedom to align (see AliMUONAlignment)
139   alig->AllowVariations(bChOnOff);
140
141   // Fix parameters or add constraints here
142 //   for (Int_t iSt=0; iSt<5; iSt++)
143 //     if (!bStOnOff[iSt]) alig->FixStation(iSt+1);
144   for (Int_t iCh=0; iCh<10; iCh++)
145     if (!bChOnOff[iCh]) alig->FixChamber(iCh+1);
146
147   // Left and right sides of the detector are independent, one can choose to align 
148   // only one side
149   Bool_t bSpecLROnOff[2] = {kTRUE,kTRUE};
150   alig->FixHalfSpectrometer(bChOnOff,bSpecLROnOff);
151
152   alig->SetChOnOff(bChOnOff);
153   alig->SetSpecLROnOff(bChOnOff);
154
155   // Here we can fix some detection elements
156   alig->FixDetElem(908);
157   alig->FixDetElem(1020);
158
159   // Set predifined global constrains: X, Y, P, XvsZ, YvsZ, PvsZ, XvsY, YvsY, PvsY
160   Bool_t bVarXYT[9] = {kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE};
161   Bool_t bDetTLBR[4] = {kFALSE,kTRUE,kFALSE,kTRUE};
162   //  alig->AddConstraints(bChOnOff,bVarXYT,bDetTLBR,bSpecLROnOff);
163
164
165   char cFileName[100];  
166     
167   Int_t lMaxFile = 1000;
168   Int_t iFile = 0;
169   Int_t iEvent = 0;
170   bool bKeepLoop = kTRUE;
171   Int_t iTrackTot=0;
172   Int_t iTrackOk=0;
173
174   while(bKeepLoop && iFile<lMaxFile){
175     iFile++;
176     if (bLoop) {
177       sFileList.getline(cFileName,100);
178       if (sFileList.eof()) bKeepLoop = kFALSE;
179     }
180     else {
181       sprintf(cFileName,esdFileName.Data());
182       bKeepLoop = kFALSE;
183     }
184     if (!strstr(cFileName,"AliESDs")) continue;      
185     cout << "Using file: " << cFileName << endl;
186     
187     // load ESD event
188     TFile* esdFile = TFile::Open(cFileName); // open the file
189     if (!esdFile || !esdFile->IsOpen()) {
190       cout << "opening ESD file " << cFileName << "failed" << endl;
191       continue;
192     }
193     TTree* esdTree = (TTree*) esdFile->Get("esdTree"); // get the tree
194     if (!esdTree) {
195       cout << "no ESD tree found" << endl;
196       esdFile->Close();
197       continue;
198     }
199     AliESDEvent* esdEvent = new AliESDEvent(); // link ESD event to the tree
200     esdEvent->ReadFromTree(esdTree);
201
202     Int_t nevents = esdTree->GetEntries();
203     cout << "... with " << nevents << endl;
204     for(Int_t event = 0; event < nevents; event++) {
205       if (iEvent >= nEvents){
206         bKeepLoop = kFALSE;
207         break;
208       }
209       iEvent++;
210
211       if (esdTree->GetEvent(event) <= 0) {
212         cout << "fails to read ESD object for event " << event << endl;
213         continue;
214       }
215
216       Int_t nTracks = Int_t(esdEvent->GetNumberOfMuonTracks());
217       if (!event%100) cout << " there are " << nTracks << " tracks in event " << event << endl;
218       for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
219         AliESDMuonTrack* esdTrack = esdEvent->GetMuonTrack(iTrack);
220         if (!esdTrack->ClustersStored()) continue;
221         Double_t invBenMom = esdTrack->GetInverseBendingMomentum();
222         fInvBenMom->Fill(invBenMom);
223         fBenMom->Fill(1./invBenMom);
224         if (TMath::Abs(invBenMom)<=1.04) {
225           AliMUONTrack track;
226           AliMUONESDInterface::ESDToMUON(*esdTrack, track);
227           alig->ProcessTrack(&track);
228           alig->LocalFit(iTrackOk++,trackParams,0);
229         }
230         iTrackTot++;
231       }
232     }
233     delete esdEvent;
234     esdFile->Close();
235     cout << "Processed " << iTrackTot << " Tracks so far." << endl;
236   }
237   alig->GlobalFit(parameters,errors,pulls);
238
239   cout << "Done with GlobalFit " << endl;
240
241   // Store results
242   Double_t DEid[156] = {0};
243   Double_t MSDEx[156] = {0};
244   Double_t MSDEy[156] = {0};
245   Double_t MSDEz[156] = {0};
246   Double_t MSDEp[156] = {0};
247   Double_t DEidErr[156] = {0};
248   Double_t MSDExErr[156] = {0};
249   Double_t MSDEyErr[156] = {0};
250   Double_t MSDEzErr[156] = {0};
251   Double_t MSDEpErr[156] = {0};
252   Int_t lNDetElem = 4*2+4*2+18*2+26*2+26*2;
253   Int_t lNDetElemCh[10] = {4,4,4,4,18,18,26,26,26,26};
254   // Int_t lSNDetElemCh[10] = {4,8,12,16,34,52,78,104,130,156};
255   Int_t idOffset = 0; // 400
256   Int_t lSDetElemCh = 0;
257   for(Int_t iDE=0; iDE<lNDetElem; iDE++){
258     DEidErr[iDE] = 0.;
259     DEid[iDE] = idOffset+100;
260     DEid[iDE] += iDE; 
261     lSDetElemCh = 0;
262     for(Int_t iCh=0; iCh<9; iCh++){
263       lSDetElemCh += lNDetElemCh[iCh];
264       if (iDE>=lSDetElemCh) {
265         DEid[iDE] += 100;
266         DEid[iDE] -= lNDetElemCh[iCh];
267       }
268     }
269     MSDEx[iDE]=parameters[4*iDE+0];
270     MSDEy[iDE]=parameters[4*iDE+1];
271     MSDEz[iDE]=parameters[4*iDE+3];
272     MSDEp[iDE]=parameters[4*iDE+2];
273     MSDExErr[iDE]=(Double_t)alig->GetParError(4*iDE+0);
274     MSDEyErr[iDE]=(Double_t)alig->GetParError(4*iDE+1);
275     MSDEzErr[iDE]=(Double_t)alig->GetParError(4*iDE+3);
276     MSDEpErr[iDE]=(Double_t)alig->GetParError(4*iDE+2);
277   }
278
279   cout << "Let's create graphs ...  " << endl;
280
281   TGraphErrors *gMSDEx = new TGraphErrors(lNDetElem,DEid,MSDEx,DEidErr,MSDExErr); 
282   TGraphErrors *gMSDEy = new TGraphErrors(lNDetElem,DEid,MSDEy,DEidErr,MSDEyErr); 
283   TGraphErrors *gMSDEz = new TGraphErrors(lNDetElem,DEid,MSDEz,DEidErr,MSDEzErr); 
284   TGraphErrors *gMSDEp = new TGraphErrors(lNDetElem,DEid,MSDEp,DEidErr,MSDEpErr); 
285
286   cout << "... graphs created, open file ...  " << endl;
287
288   TFile *hFile = new TFile("measShifts.root","RECREATE");
289
290   cout << "... file opened ...  " << endl;
291
292   gMSDEx->Write("gMSDEx");
293   gMSDEy->Write("gMSDEy");
294   gMSDEz->Write("gMSDEz");
295   gMSDEp->Write("gMSDEp");
296   fInvBenMom->Write();
297   fBenMom->Write();
298   hFile->Close();
299   
300   cout << "... and closed!" << endl;
301   // Re Align
302   AliMUONGeometryTransformer *newTransform = alig->ReAlign(transform,parameters,true); 
303   newTransform->WriteTransformations("transform2ReAlign.dat");
304   
305   // Generate realigned data in local cdb
306   const TClonesArray* array = newTransform->GetMisAlignmentData();
307
308   // 100 mum residual resolution for chamber misalignments?
309   alig->SetAlignmentResolution(array,-1,0.01,0.01,0.004,0.003);
310    
311   // CDB manager
312   AliCDBManager* cdbManager = AliCDBManager::Instance();
313   cdbManager->SetDefaultStorage("local://ReAlignCDB");
314   
315   AliCDBMetaData* cdbData = new AliCDBMetaData();
316   cdbData->SetResponsible("Dimuon Offline project");
317   cdbData->SetComment("MUON alignment objects with residual misalignment");
318   AliCDBId id("MUON/Align/Data", 0, AliCDBRunRange::Infinity()); 
319   cdbManager->Put(const_cast<TClonesArray*>(array), id, cdbData);
320
321 }
322