]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/AlignITSmacro3Rphi.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Improved cluster finder algorithm for SDD
[u/mrichter/AliRoot.git] / ITS / AlignITSmacro3Rphi.C
CommitLineData
e8189707 1//////////////////////////////////////////////////////////////////////////
2// Alice ITS first detector alignment program. //
3// //
4// version: 0.0.0 Draft. //
5// Date: April 18 1999 //
6// By: Bjorn S. Nilsen //
7// //
8//////////////////////////////////////////////////////////////////////////
9
10#include <fstream.h>
11
12// Data structure to hold averaged clusts.
13struct ClustAl_sl{
14 Int_t lay,lad,det;
15 Float_t xg,yg,zg,xl,yl,zl;
16};
17struct ClustAl_tl{
18 Int_t track,nclust; // track number and number of data points.
19 ClustAl_sl *clust; // data points to fit.
20 Float_t a,b,c,d,a0,b0,c0,d0,qual; // fit parameters and fit quality.
21 Float_t px,py,pz,p,pt;
22 // x=a+bz and y=c+dz;
23 // x=a0+b0*y and z=c0+d0*y in coordinate system of clust[0].lay,lad,det
24};
25
26//
27void AlignITSmacro3Rphi(const char *Rfilename="galice_ITS_B0.root",
28 const char *sfile="Align_ITS_B0",
29 Float_t sigma1=0.0,Float_t sigma2=0.0,Float_t sigma3=0.0,
30 Int_t evNumber=0) {
31/////////////////////////////////////////////////////////////////////////
32// This macro is a small example of a ROOT macro illustrating how to
33// read the output of GALICE and fill some histograms.
34//
35// Root > .L AlignITSmacro2.C // this loads the macro in memory
36// Root > AlignITSmacro2(); // by default process first event
37// Root > AlignITSmacro2("galice.root"); // process file galice.root
38// Root > AlignITSmacro2("galice.root",3); // process file galice.root
39// the first 4 events
40// or
41// Root > .x AlignITSmacro2.C; //by default process first event
42// Root > .x AlignITSmacro2.C("galice.root"); // process file galice.root
43// Root > .x AlignITSmacro2.C("galice.root",3); // process file galice.root
44// the first 4 events
45/////////////////////////////////////////////////////////////////////////
46//
47 gROOT->Reset(); // Reset root to it's default state
48// Dynamically link some shared libs
49 if (gClassTable->GetID("AliRun") < 0) {
50 gROOT->LoadMacro("loadlibs.C");
51 loadlibs();
52 } // end if gClassTable...
53
54 // Connect the Root Galice file containing Geometry, Kine and Clusts
55 TFile *Rfile = (TFile*)gROOT->GetListOfFiles()->FindObject(Rfilename);
56 if(!Rfile) Rfile = new TFile(Rfilename);
57 printf("reading from file %s\n",Rfilename);
58
59 // Get AliRun object from file or create it if not on file
60 if(!gAlice) {
61 gAlice = (AliRun*)Rfile->Get("gAlice");
62 if( gAlice) printf("AliRun object found on file\n");
63 if(!gAlice) gAlice = new AliRun("gAlice","Alice test program");
64 } /* end if gAlice */
65
66 Float_t v0[3] = {0.0,0.0,0.0};
67 Float_t tran[3] = {0.0,0.0,0.0},rot[3] = {0.0,0.0,0.0};
68 Float_t trans[15] ={0.0E-0,1.0E-4,4.0E-4,7.0E-4,1.0E-3,
69 2.0E-3,4.0E-3,6.0E-3,8.0E-3,1.0E-2,
70 2.0E-2,3.0E-2,5.0E-2,7.5E-2,1.0E-1}; // cm
71 Float_t rots[15] ={0.0E-0,1.0E-4,4.0E-4,7.0E-4,1.0E-3,
72 2.0E-3,4.0E-3,6.0E-3,8.0E-3,1.0E-2,
73 2.0E-2,3.0E-2,5.0E-2,7.5E-2,1.0E-1}; // rad
74 Int_t nparticles,Ntrkp,ntrk,Itimes,Isigmas;
75 AliITS *ITS = 0;
76 TTree *TH = 0;
77 AliITSgeom *gm,gm2;
78 char Hfilename[80];
79// char Gfilename[80];
80 TFile *Hfile = 0;
81 ClustAl_tl *trk = 0;
82 Float_t Rdta[6];
83 Float_t Fdta[24],Fdr[12],Fdta0[24],Fdr0[12];
84 Int_t Ndta[12];
85 FILE *fp1,*fp2;
86
87 for(Int_t evnt=0;evnt<=evNumber;evnt++){
88 // define some variables for later use.
89 nparticles = gAlice->GetEvent(evnt);
90 printf("nparticles %d\n",nparticles);
91 if (nparticles <= 0) continue; /* get next event */
92
93 // Get pointers to Alice detectors and Clusts containers
94 ITS = (AliITS*)gAlice->GetDetector("ITS");
95 if(!ITS) return; /* error no ITS data exit */
96 TH = gAlice->TreeH();
97 Ntrkp = TH->GetEntries();
98 gm = ITS->GetITSgeom();
99
100 // Array (stucture) of clusts for the first and second layer
101 // this should be replaced with either clusters or digits
102 // when they are proporly defined.
103 trk = new ClustAl_tl[Ntrkp];
104
105 printf("Ntrkp=%d\n",Ntrkp);
106
107 HitsToClustAl(trk,ntrk,Ntrkp,TH,ITS,1.0);
108 printf("Filled data structures ntrk=%d fitting lines next\n",ntrk);
109//
110 for(Itimes=1;Itimes<2;Itimes++){
111 if(Itimes==0){
112 fp1 = fopen("RvariationsR.csv","w");
113 fp2 = fopen("RvariationsN.csv","w");
114 fprintf(fp1,"Rr,Rs,Rrx1,Rerx1,Rsrx1,Rserx1,Rrz1,Rerz1,Rsrz1,Rserz1,"
115 "Rrd1,Rerd1,Rsrd1,Resrd1,"
116 "Rrx2,Rerx2,Rsrx2,Rserx2,Rrz2,Rerz2,Rsrz2,Rserz2,"
117 "Rrd2,Rerd2,Rsrd2,Resrd2,"
118 "Rrx3,Rerx3,Rsrx3,Rserx3,Rrz3,Rerz3,Rsrz3,Rserz3,"
119 "Rrd3,Rerd3,Rsrd3,Resrd3,"
120 "Rrx4,Rerx4,Rsrx4,Rserx4,Rrz4,Rerz4,Rsrz4,Rserz4,"
121 "Rrd4,Rerd4,Rsrd4,Resrd4,"
122 "Rrx5,Rerx5,Rsrx5,Rserx5,Rrz5,Rerz5,Rsrz5,Rserz5,"
123 "Rrd5,Rerd5,Rsrd5,Resrd5,"
124 "Rrx6,Rerx6,Rsrx6,Rserx6,Rrz6,Rerz6,Rsrz6,Rserz6,"
125 "Rrd6,Rerd6,Rsrd6,Resrd6");
126 fprintf(fp2,"Rr,Rs,RN1,RN2,RN3,RN4,RN5,RN6,RN7,RN8,RN9,RN10,"
127 "RN50,Rcut");
128 }else if(Itimes==1){
129 fp1 = fopen("RPhivariationsR.csv","w");
130 fp2 = fopen("RPhivariationsN.csv","w");
131 fprintf(fp1,"Fr,Fs,Frx1,Ferx1,Fsrx1,Fserx1,Frz1,Ferz1,Fsrz1,Fserz1,"
132 "Frd1,Ferd1,Fsrd1,Fesrd1,"
133 "Frx2,Ferx2,Fsrx2,Fserx2,Frz2,Ferz2,Fsrz2,Fserz2,"
134 "Frd2,Ferd2,Fsrd2,Fesrd2,"
135 "Frx3,Ferx3,Fsrx3,Fserx3,Frz3,Ferz3,Fsrz3,Fserz3,"
136 "Frd3,Ferd3,Fsrd3,Fesrd3,"
137 "Frx4,Ferx4,Fsrx4,Fserx4,Frz4,Ferz4,Fsrz4,Fserz4,"
138 "Frd4,Ferd4,Fsrd4,Fesrd4,"
139 "Frx5,Ferx5,Fsrx5,Fserx5,Frz5,Ferz5,Fsrz5,Fserz5,"
140 "Frd5,Ferd5,Fsrd5,Fesrd5,"
141 "Frx6,Ferx6,Fsrx6,Fserx6,Frz6,Ferz6,Fsrz6,Fserz6,"
142 "Frd6,Ferd6,Fsrd6,Fesrd6");
143 fprintf(fp2,"Fr,Fs,FN1,FN2,FN3,FN4,FN5,FN6,FN7,FN8,FN9,FN10,"
144 "FN50,Fcut");
145 }else if(Itimes==2){
146 fp1 = fopen("ZvariationsR.csv","w");
147 fp2 = fopen("ZvariationsN.csv","w");
148 fprintf(fp1,"Zr,Zs,Zrx1,Zerx1,Zsrx1,Zserx1,Zrz1,Zerz1,Zsrz1,Zserz1,"
149 "Zrd1,Zerd1,Zsrd1,Zesrd1,"
150 "Zrx2,Zerx2,Zsrx2,Zserx2,Zrz2,Zerz2,Zsrz2,Zserz2,"
151 "Zrd2,Zerd2,Zsrd2,Zesrd2,"
152 "Zrx3,Zerx3,Zsrx3,Zserx3,Zrz3,Zerz3,Zsrz3,Zserz3,"
153 "Zrd3,Zerd3,Zsrd3,Zesrd3,"
154 "Zrx4,Zerx4,Zsrx4,Zserx4,Zrz4,Zerz4,Zsrz4,Zserz4,"
155 "Zrd4,Zerd4,Zsrd4,Zesrd4,"
156 "Zrx5,Zerx5,Zsrx5,Zserx5,Zrz5,Zerz5,Zsrz5,Zserz5,"
157 "Zrd5,Zerd5,Zsrd5,Zesrd5,"
158 "Zrx6,Zerx6,Zsrx6,Zserx6,Zrz6,Zerz6,Zsrz6,Zserz6,"
159 "Zrd6,Zerd6,Zsrd6,Zesrd6");
160 fprintf(fp2,"Zr,Zs,ZN1,ZN2,ZN3,ZN4,ZN5,ZN6,ZN7,ZN8,ZN9,ZN10,"
161 "ZN50,Zcut");
162 }else if(Itimes==3){
163 fp1 = fopen("AvariationsR.csv","w");
164 fp2 = fopen("AvariationsN.csv","w");
165 fprintf(fp1,"Ar,As,Arx1,Aerx1,Asrx1,Aserx1,Arz1,Aerz1,Asrz1,Aserz1,"
166 "Ard1,Aerd1,Asrd1,Aesrd1,"
167 "Arx2,Aerx2,Asrx2,Aserx2,Arz2,Aerz2,Asrz2,Aserz2,"
168 "Ard2,Aerd2,Asrd2,Aesrd2,"
169 "Arx3,Aerx3,Asrx3,Aserx3,Arz3,Aerz3,Asrz3,Aserz3,"
170 "Ard3,Aerd3,Asrd3,Aesrd3,"
171 "Arx4,Aerx4,Asrx4,Aserx4,Arz4,Aerz4,Asrz4,Aserz4,"
172 "Ard4,Aerd4,Asrd4,Aesrd4,"
173 "Arx5,Aerx5,Asrx5,Aserx5,Arz5,Aerz5,Asrz5,Aserz5,"
174 "Ard5,Aerd5,Asrd5,Aesrd5,"
175 "Arx6,Aerx6,Asrx6,Aserx6,Arz6,Aerz6,Asrz6,Aserz6,"
176 "Ard6,Aerd6,Asrd6,Aesrd6");
177 fprintf(fp2,"Ar,As,AN1,AN2,AN3,AN4,AN5,AN6,AN7,AN8,AN9,AN10,"
178 "AN50,Acut");
179 }else if(Itimes==4){
180 fp1 = fopen("BvariationsR.csv","w");
181 fp2 = fopen("BvariationsN.csv","w");
182 fprintf(fp1,"Br,Bs,Brx1,Berx1,Bsrx1,Bserx1,Brz1,Berz1,Bsrz1,Bserz1,"
183 "Brd1,Berd1,Bsrd1,Besrd1,"
184 "Brx2,Berx2,Bsrx2,Bserx2,Brz2,Berz2,Bsrz2,Bserz2,"
185 "Brd2,Berd2,Bsrd2,Besrd2,"
186 "Brx3,Berx3,Bsrx3,Bserx3,Brz3,Berz3,Bsrz3,Bserz3,"
187 "Brd3,Berd3,Bsrd3,Besrd3,"
188 "Brx4,Berx4,Bsrx4,Bserx4,Brz4,Berz4,Bsrz4,Bserz4,"
189 "Brd4,Berd4,Bsrd4,Besrd4,"
190 "Brx5,Berx5,Bsrx5,Bserx5,Brz5,Berz5,Bsrz5,Bserz5,"
191 "Brd5,Berd5,Bsrd5,Besrd5,"
192 "Brx6,Berx6,Bsrx6,Bserx6,Brz6,Berz6,Bsrz6,Bserz6,"
193 "Brd6,Berd6,Bsrd6,Besrd6");
194 fprintf(fp2,"Br,Bs,BN1,BN2,BN3,BN4,BN5,BN6,BN7,BN8,BN9,BN10,"
195 "BN50,Bcut");
196 }else if(Itimes==5){
197 fp1 = fopen("CvariationsR.csv","w");
198 fp2 = fopen("CvariationsN.csv","w");
199 fprintf(fp1,"Cr,Cs,Crx1,Cerx1,Csrx1,Cserx1,Crz1,Cerz1,Csrz1,Cserz1,"
200 "Crd1,Cerd1,Csrd1,Cesrd1,"
201 "Crx2,Cerx2,Csrx2,Cserx2,Crz2,Cerz2,Csrz2,Cserz2,"
202 "Crd2,Cerd2,Csrd2,Cesrd2,"
203 "Crx3,Cerx3,Csrx3,Cserx3,Crz3,Cerz3,Csrz3,Cserz3,"
204 "Crd3,Cerd3,Csrd3,Cesrd3,"
205 "Crx4,Cerx4,Csrx4,Cserx4,Crz4,Cerz4,Csrz4,Cserz4,"
206 "Crd4,Cerd4,Csrd4,Cesrd4,"
207 "Crx5,Cerx5,Csrx5,Cserx5,Crz5,Cerz5,Csrz5,Cserz5,"
208 "Crd5,Cerd5,Csrd5,Cesrd5,"
209 "Crx6,Cerx6,Csrx6,Cserx6,Crz6,Cerz6,Csrz6,Cserz6,"
210 "Crd6,Cerd6,Csrd6,Cesrd6");
211 fprintf(fp2,"Cr,Cs,CN1,CN2,CN3,CN4,CN5,CN6,CN7,CN8,CN9,CN10,"
212 "CN50,Ccut");
213 } // end if
214 fprintf(fp1,"\n");
215 fprintf(fp2,"\n");
216 for(Isigmas=0;Isigmas<15;Isigmas++){
217//
218// tran[0] = sigma1;
219// tran[1] = sigma2;
220// tran[2] = sigma3;
221 if(Itimes==0){ tran[0] = trans[Isigmas];
222 }else tran[0] = 0.0;
223 if(Itimes==1){ tran[1] = trans[Isigmas];
224 }else tran[1] = 0.0;
225 if(Itimes==2){ tran[2] = trans[Isigmas];
226 }else tran[2] = 0.0;
227 if(Itimes==3){ rot[0] = rots[Isigmas];
228 }else rot[0] = 0.0;
229 if(Itimes==4){ rot[1] = rots[Isigmas];
230 }else rot[1] = 0.0;
231 if(Itimes==5){ rot[2] = rots[Isigmas];
232 }else rot[2] = 0.0;
233 printf("tran= %e %e %e (cm), rot=%e %e %e (rad)\n",
234 tran[0],tran[1],tran[2],rot[0],rot[1],rot[2]);
235//
236 gm2 = *gm;
237 gm2.RandomCylindericalChange(tran,rot);
238
239 FillGlobalPositions(trk,ntrk,&gm2);
240 // setup to save all created histograms.
241 sprintf(Hfilename,"%s_%04.0fr%04.0fp%04.0fz%04.0fx%04.0fy%04.0fz.root",
242 sfile,
243 10000.*tran[0],10000.*tran[1],10000.*tran[2],
244 10000.* rot[0],10000.* rot[1],10000.* rot[2]);
245 Hfile = (TFile*)gROOT->GetListOfFiles()->FindObject(Hfilename);
246 if(Hfile) Hfile->Close();
247 Hfile = new TFile(Hfilename,"RECREATE","Histograms from AlignITS.C");
248 printf("Histograms saved to file %s\n",Hfilename);
249 //
250 PlotGeomChanges(gm,&gm2,Hfile,Rdta);
251 //
252 // fit all tracks and do a track quality hist.
253 FitAllTracks(trk,ntrk,v0,&gm2,sfile,Hfile,Fdta,Ndta);
254//
255 for(Int_t l=0;l<6;l++){
256 Fdr[2*l+0] = sqrt(Fdta[4*l+0]*Fdta[4*l+0] +
257 Fdta[4*l+2]*Fdta[4*l+2]);
258 Fdr[2*l+1] = sqrt(Fdta[4*l+0]*Fdta[4*l+0]*
259 Fdta[4*l+1]*Fdta[4*l+1] +
260 Fdta[4*l+2]*Fdta[4*l+2]*
261 Fdta[4*l+3]*Fdta[4*l+3])/Fdr[2*l+0];
262 } // end for l
263 if(Isigmas==0){
264 for(Int_t fp=0;fp<24;fp++){
265 Fdta0[fp] = Fdta[fp];
266 if(fp<12) Fdr0[fp] = Fdr[fp];
267 } // end for fp
268 } // end if Itimes==0&&Isigmas==0
269 if(Itimes==0) fprintf(fp1,"%e,%e,",tran[0],Rdta[0]);
270 else if(Itimes==1) fprintf(fp1,"%e,%e,",tran[1],Rdta[1]);
271 else if(Itimes==2) fprintf(fp1,"%e,%e,",tran[2],Rdta[2]);
272 else if(Itimes==3) fprintf(fp1,"%e,%e,",rot[0],Rdta[3]);
273 else if(Itimes==4) fprintf(fp1,"%e,%e,",rot[1],Rdta[4]);
274 else if(Itimes==5) fprintf(fp1,"%e,%e,",rot[2],Rdta[5]);
275 if(Itimes==0) fprintf(fp2,"%e,%e,",tran[0],Rdta[0]);
276 else if(Itimes==1) fprintf(fp2,"%e,%e,",tran[1],Rdta[1]);
277 else if(Itimes==2) fprintf(fp2,"%e,%e,",tran[2],Rdta[2]);
278 else if(Itimes==3) fprintf(fp2,"%e,%e,",rot[0],Rdta[3]);
279 else if(Itimes==4) fprintf(fp2,"%e,%e,",rot[1],Rdta[4]);
280 else if(Itimes==5) fprintf(fp2,"%e,%e,",rot[2],Rdta[5]);
281 fprintf(fp1,"%e,%e,%e,%e,",Fdta[0],Fdta[1],Fdta[0]/Fdta0[0],Fdta[1]/Fdta0[0]);
282 fprintf(fp1,"%e,%e,%e,%e,",Fdta[2],Fdta[3],Fdta[2]/Fdta0[2],Fdta[3]/Fdta0[2]);
283 fprintf(fp1,"%e,%e,%e,%e,",Fdr[0],Fdr[1],Fdr[0]/Fdr0[0],Fdr[1]/Fdr0[0]);
284 fprintf(fp1,"%e,%e,%e,%e,",Fdta[4],Fdta[5],Fdta[4]/Fdta0[4],Fdta[5]/Fdta0[4]);
285 fprintf(fp1,"%e,%e,%e,%e,",Fdta[6],Fdta[7],Fdta[6]/Fdta0[6],Fdta[7]/Fdta0[6]);
286 fprintf(fp1,"%e,%e,%e,%e,",Fdr[2],Fdr[3],Fdr[2]/Fdr0[2],Fdr[3]/Fdr0[2]);
287 fprintf(fp1,"%e,%e,%e,%e,",Fdta[8],Fdta[9],Fdta[8]/Fdta0[8],Fdta[9]/Fdta0[8]);
288 fprintf(fp1,"%e,%e,%e,%e,",Fdta[10],Fdta[11],Fdta[10]/Fdta0[10],Fdta[11]/Fdta0[10]);
289 fprintf(fp1,"%e,%e,%e,%e,",Fdr[4],Fdr[5],Fdr[4]/Fdr0[4],Fdr[5]/Fdr0[4]);
290 fprintf(fp1,"%e,%e,%e,%e,",Fdta[12],Fdta[12],Fdta[12]/Fdta0[12],Fdta[13]/Fdta0[12]);
291 fprintf(fp1,"%e,%e,%e,%e,",Fdta[14],Fdta[13],Fdta[14]/Fdta0[14],Fdta[15]/Fdta0[14]);
292 fprintf(fp1,"%e,%e,%e,%e,",Fdr[6],Fdr[7],Fdr[6]/Fdr0[6],Fdr[7]/Fdr0[6]);
293 fprintf(fp1,"%e,%e,%e,%e,",Fdta[16],Fdta[15],Fdta[16]/Fdta0[16],Fdta[17]/Fdta0[16]);
294 fprintf(fp1,"%e,%e,%e,%e,",Fdta[18],Fdta[17],Fdta[18]/Fdta0[18],Fdta[19]/Fdta0[18]);
295 fprintf(fp1,"%e,%e,%e,%e,",Fdr[8],Fdr[9],Fdr[8]/Fdr0[8],Fdr[9]/Fdr0[8]);
296 fprintf(fp1,"%e,%e,%e,%e,",Fdta[20],Fdta[19],Fdta[20]/Fdta0[20],Fdta[21]/Fdta0[20]);
297 fprintf(fp1,"%e,%e,%e,%e,",Fdta[22],Fdta[21],Fdta[22]/Fdta0[22],Fdta[23]/Fdta0[22]);
298 fprintf(fp1,"%e,%e,%e,%e\n",Fdr[10],Fdr[11],Fdr[10]/Fdr0[10],Fdr[11]/Fdr0[10]);
299 fprintf(fp2,"%d,%d,%d,%d,",Ndta[0],Ndta[1],Ndta[2],Ndta[3]);
300 fprintf(fp2,"%d,%d,%d,%d,",Ndta[4],Ndta[5],Ndta[6],Ndta[7]);
301 fprintf(fp2,"%d,%d,%d,%d\n",Ndta[8],Ndta[9],Ndta[10],Ndta[11]);
302// FitVertexAll(trk,ntrk,sfile,Hfile);
303 // Find all 2 track vertecies and hist. them.
304 Hfile->Close();
305 //
306 } // end for Isigmas
307 fclose(fp1);
308 fclose(fp2);
309 } // end for Itimes
310//
311 printf("Event %d done\n",evnt);
312//
313 deleteClustAl(trk,ntrk); // subrotine to delet memory allocated
314 // inside HitsToclustAl.
315 delete[] trk; // now delet memory allocated above.
316 } // end for evnt
317 Rfile->Close();
318 return;
319}