]>
Commit | Line | Data |
---|---|---|
fe4da5cc | 1 | void MUONlist2digits (Int_t evNumber1=0,Int_t evNumber2=0,Int_t nCathode=1) |
2 | { | |
3 | ///////////////////////////////////////////////////////////////////////// | |
4 | // This macro is a small example of a ROOT macro | |
5 | // illustrating how to read the output of GALICE | |
6 | // and do some analysis. | |
7 | // | |
8 | ///////////////////////////////////////////////////////////////////////// | |
9 | ||
10 | // Dynamically link some shared libs | |
11 | ||
12 | ||
13 | if (gClassTable->GetID("AliRun") < 0) { | |
14 | gSystem->Load("$ALITOP/cern.so/lib/libpdfDUMMY.so"); | |
15 | gSystem->Load("$ALITOP/cern.so/lib/libPythia.so"); | |
16 | gSystem->Load("$ROOTSYS/lib/libEG.so"); | |
17 | gSystem->Load("$ROOTSYS/lib/libEGPythia.so"); | |
18 | gSystem->Load("libGeant3Dummy.so"); //a dummy version of Geant3 | |
19 | gSystem->Load("PHOS/libPHOSdummy.so"); //the standard Alice classes | |
20 | gSystem->Load("libgalice.so"); // the standard Alice classes | |
21 | } | |
22 | ||
23 | // Connect the Root Galice file containing Geometry, Kine and Hits | |
24 | ||
25 | TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root"); | |
26 | if (file) file->Close(); | |
27 | file = new TFile("galice.root","UPDATE"); | |
28 | file->ls(); | |
29 | file->Map(); | |
30 | ||
31 | printf ("I'm after Map \n"); | |
32 | ||
33 | // Get AliRun object from file or create it if not on file | |
34 | ||
35 | if (!gAlice) { | |
36 | gAlice = (AliRun*)file->Get("gAlice"); | |
37 | if (gAlice) printf("AliRun object found on file\n"); | |
38 | if (!gAlice) gAlice = new AliRun("gAlice","Alice test program"); | |
39 | } | |
40 | printf ("I'm after gAlice \n"); | |
41 | ||
42 | AliMUON *MUON = gAlice->GetDetector("MUON"); | |
43 | ||
44 | // AliMUONchamber& iChamber; | |
45 | AliMUONchamber* iChamber; | |
46 | AliMUONsegmentation* segmentation; | |
47 | ||
48 | const Float_t rmin[10] = {17.5,17.5,23.5,23.5,33.5,33.5,43,43,50,50}; | |
49 | const Float_t rmax[10] = {91.5,91.5,122.5,122.5,173,173,221,221,256.5,256.5}; | |
50 | ||
51 | ||
52 | Int_t Npx[10]; | |
53 | Int_t Npy[10]; | |
54 | ||
55 | ||
56 | for ( int i=0;i<10;i++) { | |
57 | Npx[i]=(Int_t)rmax[i]/0.75; | |
58 | Npy[i]=(Int_t)rmax[i]/0.5; | |
59 | } | |
60 | ||
61 | Int_t nxmax=1026; | |
62 | Int_t nymax=1026; | |
63 | AliMUONlist *elem[10][1026][1026]; | |
64 | TObjArray *obj=new TObjArray; | |
65 | ||
66 | Int_t digits[3]; | |
67 | // | |
68 | // Loop over events | |
69 | // | |
70 | ||
71 | Int_t Nh=0; | |
72 | Int_t Nh1=0; | |
73 | for (int nev=0; nev<= evNumber2; nev++) { | |
74 | Int_t nparticles = gAlice->GetEvent(nev); | |
75 | cout << "nev " << nev <<endl; | |
76 | cout << "nparticles " << nparticles <<endl; | |
77 | if (nev < evNumber1) continue; | |
78 | if (nparticles <= 0) return; | |
79 | ||
80 | TTree *TH = gAlice->TreeH(); | |
81 | Int_t ntracks = TH->GetEntries(); | |
82 | Int_t Nc=0; | |
83 | // | |
84 | Int_t counter=0; | |
85 | ||
86 | for (int icat=0;icat<nCathode;icat++) { | |
87 | ||
88 | //initialize the arrray of pointers !!!!! | |
89 | for (int i=0; i<10; i++) { | |
90 | for (int j=0; j<nymax; j++) { | |
91 | for (int k=0; k<nxmax; k++) { | |
92 | elem[i][j][k]=0; | |
93 | } | |
94 | } | |
95 | } | |
96 | ||
97 | ||
98 | // | |
99 | // Loop over events | |
100 | // | |
101 | for (Int_t track=0; track<ntracks;track++) { | |
102 | gAlice->ResetHits(); | |
103 | Int_t nbytes += TH->GetEvent(track); | |
104 | if (MUON) { | |
105 | for(AliMUONhit* mHit=(AliMUONhit*)MUON->FirstHit(-1); | |
106 | mHit; | |
107 | mHit=(AliMUONhit*)MUON->NextHit()) | |
108 | { | |
109 | Int_t nch = mHit->fChamber; // chamber number | |
110 | Float_t x = mHit->fX; // x-pos of hit | |
111 | Float_t y = mHit->fY; // y-pos | |
112 | ||
113 | iChamber = &(MUON->Chamber(nch-1)); | |
114 | response=iChamber->GetResponseModel(); | |
115 | ||
116 | Int_t nsec=iChamber->Nsec(); | |
117 | // Int_t rmin = (Int_t)iChamber->frMin; | |
118 | // Int_t rmax = (Int_t)iChamber->frMax; | |
119 | // | |
120 | // | |
121 | for (AliMUONcluster* mPad=(AliMUONcluster*)MUON->FirstPad(mHit); | |
122 | mPad; | |
123 | mPad=(AliMUONcluster*)MUON->NextPad()) | |
124 | { | |
125 | Int_t cathode = mPad->fCathode; // chamber number | |
126 | Int_t nhit = mPad->fHitNumber; // hit number | |
127 | Int_t qtot = mPad->fQ; // charge | |
128 | Int_t ipx = mPad->fPadX; // pad number on X | |
129 | Int_t ipy = mPad->fPadY; // pad number on Y | |
130 | Int_t iqpad = mPad->fQpad; // charge per pad | |
131 | Int_t rpad = mPad->fRpad; // r-pos of pad | |
132 | // | |
133 | // | |
134 | if (cathode != (icat+1)) continue; | |
135 | segmentation=iChamber->GetSegmentationModel(cathode); | |
136 | // Int_t Npx[nch-1] = segmentation->Npx(); | |
137 | // Int_t Npy[nch-1] = segmentation->Npy(); | |
138 | ||
139 | if (cathode==1) { | |
140 | Int_t npx=Npx[nch-1]; | |
141 | Int_t npy=Npy[nch-1]; | |
142 | } else { | |
143 | Int_t npx=Npy[nch-1]; | |
144 | Int_t npy=Npx[nch-1]; | |
145 | } | |
146 | ||
147 | // printf("nch, cathode, rmin, rmax, npx, npy %d %d %d %d %d %d\n", | |
148 | // nch,cathode,rmin,rmax,npx,npy); | |
149 | ||
150 | // printf("icat, iqpad, ipx, ipy %d %d %d %d \n", | |
151 | // icat, iqpad, ipx, ipy); | |
152 | ||
153 | ||
154 | // check boundaries | |
155 | if (rpad < rmin[nch-1] || iqpad ==0 || rpad > rmax[nch-1]) continue; | |
156 | // if (rpad < rmin[nch-1] || iqpad ==0 || rpad > 81.65) continue; | |
157 | // if (rpad < rmin || iqpad ==0 || rpad > rmax) continue; | |
158 | ||
159 | // fill the info array | |
160 | TVector *trinfo; | |
161 | trinfo=new TVector(2); | |
162 | ||
163 | trinfo(0)=(Float_t)track; | |
164 | trinfo(1)=(Float_t)iqpad; | |
165 | ||
166 | // Int_t digits[3]; | |
167 | digits[0]=ipx; | |
168 | digits[1]=ipy; | |
169 | digits[2]=iqpad; | |
170 | ||
171 | ||
172 | Int_t trk[50]; | |
173 | Int_t chtrk[50]; | |
174 | ||
175 | // build the list of fired pads and update the info | |
176 | AliMUONlist *pdigit=elem[nch-1][ipy+npy][ipx+npx]; | |
177 | if (pdigit==0) { | |
178 | obj->AddAtAndExpand(new AliMUONlist(rpad,digits),counter); | |
179 | counter++; | |
180 | // Int_t nentrobj=obj->GetEntriesFast(); | |
181 | Int_t last=obj->GetLast(); | |
182 | pdigit=(AliMUONlist*)obj->At(last); | |
183 | elem[nch-1][ipy+npy][ipx+npx]=pdigit; | |
184 | // Int_t q=pdigit->fSignal; | |
185 | // list of tracks | |
186 | TObjArray *trlist=(TObjArray*)pdigit->TrackList(); | |
187 | trlist->Add(trinfo); | |
188 | } else { | |
189 | // update charge | |
190 | (*pdigit).fSignal+=iqpad; | |
191 | // update list of tracks | |
192 | TObjArray* trlist=(TObjArray*)pdigit->TrackList(); | |
193 | Int_t nptracks=trlist->GetEntriesFast(); | |
194 | for (Int_t tr=0;tr<nptracks;tr++) { | |
195 | TVector *ptrk=(TVector*)trlist->At(tr); | |
196 | Int_t entries=ptrk->GetNrows(); | |
197 | TVector &vtrk = *ptrk; | |
198 | trk[tr]=(Int_t)vtrk(0); | |
199 | chtrk[tr]=(Int_t)vtrk(1); | |
200 | if (trk[tr]==track) { | |
201 | chtrk[tr]+=iqpad; | |
202 | trlist->RemoveAt(tr); | |
203 | trinfo(0)=trk[tr]; | |
204 | trinfo(1)=chtrk[tr]; | |
205 | trlist->AddAt(trinfo,tr); | |
206 | } else { | |
207 | trlist->Add(trinfo); | |
208 | } | |
209 | } //end loop over list of tracks for one pad | |
210 | } // end if pdigit | |
211 | ||
212 | ||
213 | } //end loop over clust | |
214 | ||
215 | } // hit loop | |
216 | } // if MUON | |
217 | } // track loop | |
218 | ||
219 | Int_t tracks[10]; | |
220 | Int_t charges[10]; | |
221 | cout<<"start filling digits "<<endl; | |
222 | // Int_t digits[3]; | |
223 | const Int_t zero_supm = 6; | |
224 | const Int_t adc_satm = 1024; | |
225 | Int_t nd=0; | |
226 | ||
227 | ||
228 | // start filling the digits | |
229 | for (Int_t id=0;id<10;id++) { | |
230 | Int_t nd=0; | |
231 | // Int_t npx=Npx[id]; | |
232 | // Int_t npy=Npy[id]; | |
233 | ||
234 | // Int_t nx=2*npx; | |
235 | // Int_t ny=2*npx; | |
236 | for (Int_t iy=0;iy<nymax;iy++) { | |
237 | for (Int_t ix=0;ix<nxmax;ix++) { | |
238 | AliMUONlist *address=elem[id][iy][ix]; | |
239 | if (address==0) continue; | |
240 | // do zero-suppression and signal truncation | |
241 | Int_t rpad=address->fRpad; | |
242 | Int_t q=address->fSignal; | |
243 | if ( q <= zero_supm || rpad > 55) continue; | |
244 | // if ( q <= zero_supm ) continue; | |
245 | if ( q > adc_satm) q=adc_satm; | |
246 | digits[0]=(*address).fPadX; | |
247 | digits[1]=(*address).fPadY; | |
248 | digits[2]=address->fSignal; | |
249 | ||
250 | TObjArray* trlist=(TObjArray*)address->TrackList(); | |
251 | Int_t nptracks=trlist->GetEntriesFast(); | |
252 | // this should be changed to accomodate the real number of tracks | |
253 | if (nptracks > 10) { | |
254 | cout<<"Attention - nptracks > 3 "<<nptracks<<endl; | |
255 | nptracks=10; | |
256 | } | |
257 | for (Int_t tr=0;tr<nptracks;tr++) { | |
258 | TVector *pp=(TVector*)trlist->At(tr); | |
259 | Int_t entries=pp->GetNrows(); | |
260 | TVector &v2 = *pp; | |
261 | tracks[tr]=(Int_t)v2(0); | |
262 | charges[tr]=(Int_t)v2(1); | |
263 | } //end loop over list of tracks for one pad | |
264 | ||
265 | // fill digits | |
266 | MUON->AddDigits(id+1,tracks,charges,digits); | |
267 | nd++; | |
268 | ||
269 | } // end loop over ix | |
270 | } // end loop over iy | |
271 | cout<<"I'm out of the loops over pads"<<endl; | |
272 | cout<<" id nd "<<id<<" "<<nd<<endl; | |
273 | ||
274 | // CathodeIndex[id+icat*10] = nd; // update cathodes index | |
275 | ||
276 | ||
277 | } //loop over chambers | |
278 | cout<<"I'm out of the loops for digitisation"<<endl; | |
279 | ||
280 | ||
281 | gAlice->TreeD()->Fill(); | |
282 | TTree *TD=gAlice->TreeD(); | |
283 | int ndig=TD->GetEntries(); | |
284 | cout<<"number of digits "<<ndig<<endl; | |
285 | ||
286 | int ndig1=(MUON->Dch1())->GetEntriesFast(); | |
287 | cout<<"number of digits 1 "<<ndig1<<endl; | |
288 | int ndig2=(MUON->Dch2())->GetEntriesFast(); | |
289 | cout<<"number of digits 2 "<<ndig2<<endl; | |
290 | int ndig3=(MUON->Dch3())->GetEntriesFast(); | |
291 | cout<<"number of digits 3 "<<ndig3<<endl; | |
292 | int ndig4=(MUON->Dch4())->GetEntriesFast(); | |
293 | cout<<"number of digits 4 "<<ndig4<<endl; | |
294 | int ndig5=(MUON->Dch5())->GetEntriesFast(); | |
295 | cout<<"number of digits 5 "<<ndig5<<endl; | |
296 | int ndig6=(MUON->Dch6())->GetEntriesFast(); | |
297 | cout<<"number of digits 6 "<<ndig6<<endl; | |
298 | int ndig7=(MUON->Dch7())->GetEntriesFast(); | |
299 | cout<<"number of digits 7 "<<ndig7<<endl; | |
300 | int ndig8=(MUON->Dch8())->GetEntriesFast(); | |
301 | cout<<"number of digits 8 "<<ndig8<<endl; | |
302 | int ndig9=(MUON->Dch9())->GetEntriesFast(); | |
303 | cout<<"number of digits 9 "<<ndig9<<endl; | |
304 | int ndig10=(MUON->Dch10())->GetEntriesFast(); | |
305 | cout<<"number of digits 10 "<<ndig10<<endl; | |
306 | ||
307 | MUON->ResetDigits(); | |
308 | ||
309 | ||
310 | // char hname[30]; | |
311 | // sprintf(hname,"TreeD%d",nev); | |
312 | // gAlice->TreeD()->Write(hname); | |
313 | ||
314 | } //end loop over cathodes | |
315 | ||
316 | char hname[30]; | |
317 | sprintf(hname,"TreeD%d",nev); | |
318 | gAlice->TreeD()->Write(hname); | |
319 | // MUON->ResetDigits(); | |
320 | // if (CathodeIndex) CathodeIndex->Clear(); | |
321 | ||
322 | // file->Write(); | |
323 | ||
324 | file->ls(); | |
325 | file->Map(); | |
326 | } // event loop | |
327 | ||
328 | // delete [] CathodeIndex; | |
329 | ||
330 | file->Close(); | |
331 | cout<<"END digitisation "<<endl; | |
332 | ||
333 | } |