store also difference in local Y
[u/mrichter/AliRoot.git] / ITS / ITSAlignMille.C
CommitLineData
d3603b4d 1#if !defined(__CINT__) || defined(__MAKECINT__)
2#include <Riostream.h>
3#include <TChain.h>
4#include <TClonesArray.h>
5#include <TFile.h>
6#include <TMath.h>
7#include <TStopwatch.h>
8#include "AliAlignObjParams.h"
9#include "AliTrackPointArray.h"
10#include "AliITSAlignMille.h"
11#endif
12
13//********************************************************************
14// Example to run ITS alignment using Millepede
15//
16// Read tracks from AliTrackPoints.N.root and fill Millepede.
17// Millepede configuration is taken from AliITSAlignMille.conf
18// Results are written as AliAlignObjParams in outfile.
19//
20// Origin: M. Lunardon
21//********************************************************************
22
23
24int ITSAlignMille(int fromev=1, int toev=1, int fromentry=0, int nentries=-1, char *outfile="ITSAlignMille.root", char *confile="AliITSAlignMille.conf",int nreqpts=3) {
25
26 //AliLog::SetGlobalLogLevel(6);
27
28 TFile *fout=new TFile(outfile,"recreate");
29 if (!fout->IsOpen()) {
30 printf("\nCannot open output file!\n");
31 return -1;
32 }
33
34 TChain *chain=new TChain("spTree");
35 char dir[100]="AliTrackPoints";
36 char st[200];
37
38 for (int iev=fromev; iev<=toev; iev++) {
39 sprintf(st,"%s/AliTrackPoints.%d.root",dir,iev);
40 chain->Add(st);
41 }
42
43 int toentry=chain->GetEntries();
44 printf("There are %d entries in chain\n",toentry);
45
46 if (nentries>0 && nentries<(toentry-fromentry)) toentry=fromentry+nentries;
47
48 AliTrackPointArray *tpa = 0;
49 chain->SetBranchAddress("SP", &tpa);
50
51 ////////////////////////////////////////////
52
53 AliITSAlignMille *alig = new AliITSAlignMille(confile);
54 if (!alig->IsConfigured()) return -3;
55
56 Int_t nmod=alig->GetNModules();
57 alig->SetMinNPtsPerTrack(nreqpts);
58
59 // require 4 points in SPD (one per layer, up and down)
60 if (nreqpts>3) {
61 alig->SetRequiredPoint("LAYER",1,1,1);
62 alig->SetRequiredPoint("LAYER",1,-1,1);
63 alig->SetRequiredPoint("LAYER",2,1,1);
64 alig->SetRequiredPoint("LAYER",2,-1,1);
65 }
66
67 // correction for SSD bug (1) : inverisione sensor 18 e 19
68 // alig->SetBug(1);
69
70 Double_t *parameters = new Double_t[nmod*6];
71 Double_t *errors = new Double_t[nmod*6];
72 Double_t *pulls = new Double_t[nmod*6];
73
74 for(Int_t k=0;k<nmod*6;k++) {
75 parameters[k]=0.;
76 errors[k]=0.;
77 pulls[k]=0.;
78 }
79
80 // need array with size fNLocal*2
81 Double_t trackParams[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
82 alig->InitGlobalParameters(parameters);
83 alig->Print();
84
85 Double_t sigmatra=alig->GetParSigTranslations();
86 Double_t sigmarot=alig->GetParSigRotations();
87 ////////////////////////////////////////////////////////////////////
88
89 TStopwatch stimer;
90 stimer.Start();
91
92 /////////////////////
93
94 int iTrackOk=0; // number of good passed track
95 // loop on spTree entries
96 // one entry = one track
97 for (int i=fromentry; i<toentry; i++) {
98
99 chain->GetEvent(i);
100 if (!alig->ProcessTrack(tpa)) {
101 if (!(iTrackOk%500))
102 printf("Calling LocalFit n. %d\n",iTrackOk);
103 alig->LocalFit(iTrackOk++,trackParams,0);
104 }
105 }
106 printf("\nMillepede fed with %d tracks\n",iTrackOk);
107 printf("Total number of rejected points because of bad EqLoc : %d\n\n",alig->GetTotBadLocEqPoints());
108
109 stimer.Print();
110 stimer.Continue();
111
112 // make global fit
113 alig->GlobalFit(parameters,errors,pulls);
114
115 cout << "Done with GlobalFit " << endl;
116
117 ////////////////////////////////////////////////////////////
118 // output
119
120 printf("\nProcessed points statistics:\n");
121 Int_t maxstat=0;
122
123 printf("\nOutput values and point statistics:\n\n");
124 for (int i=0; i<nmod; i++) {
125 Int_t statm=alig->GetProcessedPoints()[i];
126 if (statm>maxstat) maxstat=statm;
127 printf("index: %-4d stat: %-7d pars: %f %f %f %f %f %f\n",alig->GetModuleIndexArray()[i], statm,
128 parameters[i*6+0],
129 parameters[i*6+1],
130 parameters[i*6+2],
131 parameters[i*6+3],
132 parameters[i*6+4],
133 parameters[i*6+5]);
134
135 }
136 FILE *pf=fopen("ITSAlignMille.out","w");
137 if (pf) {
138 fprintf(pf,"# param dx dy dz dpsi dtheta dphi\n");
139 for (int i=0; i<nmod; i++) {
140 fprintf(pf," %-5d %-10f %-10f %-10f %-10f %-10f %-10f\n",i,
141 parameters[i*6+0],
142 parameters[i*6+1],
143 parameters[i*6+2],
144 parameters[i*6+3],
145 parameters[i*6+4],
146 parameters[i*6+5]);
147
148 }
149 fclose(pf);
150 }
151
152 printf("Max statistics = %d\n",maxstat);
153 if (maxstat<1) {
154 printf("No points for alignment! quitting now...\n");
155 return -1;
156 }
157
158 TClonesArray *array = new TClonesArray("AliAlignObjParams",4000);
159 TClonesArray &alobj = *array;
160
161 // first object: ITS
162 new(alobj[0]) AliAlignObjParams("ITS", 0, 0., 0., 0., 0., 0., 0., kTRUE);
163
164 UShort_t volid;
165 Char_t *symname;
166 Double_t dx,dy,dz,dpsi,dtheta,dphi;
167 Double_t corr[21];
168 Double_t deltalocal[6];
169 Double_t t[3],r[3];
170
171 AliAlignObjParams *tmpa=NULL;
172
173 // quality factor: 0 = NOT ALIGNED or OFF
174 // N = ALIGNED with statistics = N
175 UInt_t QF;
176
177 // all ITS sensitive modules
178 for (int idx=0; idx<2198; idx++) {
179 volid=AliITSAlignMilleModule::GetVolumeIDFromIndex(idx);
180 symname = AliGeomManager::SymName(volid);
181 // default null misalignment
182 for (int jj=0;jj<21;jj++) corr[jj]=0.0;
183 for (int jj=0;jj<3;jj++) {t[jj]=0.0;r[jj]=0.0;}
184 QF=0;
185
186 if (alig->IsContained(volid)>=0) { // defined modules or inside a supermodule
187 alig->SetCurrentModule(idx); // set the supermodule that contains idx
188 int iidx=alig->GetCurrentModuleInternalIndex();
189
190 // check if all 6 parameters have been evaluated by millepede
191 Bool_t isoff=0;
192 Double_t parsig=0;
193 for (int kk=0; kk<6; kk++) {
194 parsig=sigmatra;
195 if (kk>2) parsig=sigmarot;
196 if (pulls[iidx*6+kk]==0.0 && errors[iidx*6+kk]==parsig) isoff=1;
197 }
198
199 // check if module was fixed
200 Bool_t isfixed=1;
201 for (int kk=0; kk<6; kk++) {
202 if (parameters[iidx*6+kk]!=0.0 || pulls[iidx*6+kk]!=0.0 || errors[iidx*6+kk]!=0.0) isfixed=0;
203 }
204
205 if (!isoff && !isfixed) { // OK, has been evaluated
206 deltalocal[0] = parameters[iidx*6+0];
207 deltalocal[1] = parameters[iidx*6+1];
208 deltalocal[2] = parameters[iidx*6+2];
209 deltalocal[3] = parameters[iidx*6+3];
210 deltalocal[4] = parameters[iidx*6+4];
211 deltalocal[5] = parameters[iidx*6+5];
212 tmpa = alig->GetCurrentModule()->GetSensitiveVolumeMisalignment(volid,deltalocal);
213 if (!tmpa) {
214 printf("error transforming millepede parameters for module %d\n",idx);
215 continue;
216 }
217 tmpa->GetPars(t,r);
218 // at the moment sigma of supermodule is given to sensitive modules. to be fixed...
219 corr[0] = errors[iidx*6+0]*errors[iidx*6+0];
220 corr[2] = errors[iidx*6+1]*errors[iidx*6+1];
221 corr[5] = errors[iidx*6+2]*errors[iidx*6+2];
222 corr[9] = errors[iidx*6+3]*errors[iidx*6+3];
223 corr[14]= errors[iidx*6+4]*errors[iidx*6+4];
224 corr[20]= errors[iidx*6+5]*errors[iidx*6+5];
225 Int_t statm=alig->GetProcessedPoints()[iidx];
226 QF=statm;
227 }
228 else {
229 if (isoff) {
230 printf("Module %d is OFF!\n",idx);
231 QF=0;
232 }
233 if (isfixed) {
234 printf("Module %d was FIXED!\n",idx);
235 if (alig->GetPreAlignmentQualityFactor(idx)>0)
236 QF=alig->GetPreAlignmentQualityFactor(idx);
237 else
238 QF=0;
239 }
240 }
241
242 }
243
244 new(alobj[idx+1]) AliAlignObjParams(symname,volid,t[0],t[1],t[2],r[0],r[1],r[2],kFALSE);
245 AliAlignObjParams* alo = (AliAlignObjParams*) array->UncheckedAt(idx+1);
246 alo->SetCorrMatrix(corr);
247 alo->SetUniqueID(QF);
248 }
249
250 fout->WriteObject(array,"ITSAlignObjs","kSingleKey");
251 fout->Close();
252
253 stimer.Stop();
254 stimer.Print();
255
256 return 0;
257}