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