]>
Commit | Line | Data |
---|---|---|
fe4da5cc | 1 | void MUONdigits (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 | if (gClassTable->GetID("AliRun") < 0) { | |
13 | gSystem->Load("$ALITOP/cern.so/lib/libpdfDUMMY.so"); | |
14 | gSystem->Load("$ALITOP/cern.so/lib/libPythia.so"); | |
15 | gSystem->Load("$ROOTSYS/lib/libEG.so"); | |
16 | gSystem->Load("$ROOTSYS/lib/libEGPythia.so"); | |
17 | gSystem->Load("libGeant3Dummy.so"); //a dummy version of Geant3 | |
18 | gSystem->Load("PHOS/libPHOSdummy.so"); //the standard Alice classes | |
19 | gSystem->Load("libgalice.so"); //the standard Alice classes | |
20 | } | |
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 | AliMUONsegmentation* segmentation; | |
46 | ||
47 | Int_t Npx[10]; | |
48 | Int_t Npy[10]; | |
49 | Int_t trk[50]; | |
50 | Int_t chtrk[50]; | |
51 | ||
52 | ||
53 | Int_t nxmax=1026; | |
54 | Int_t nymax=1026; | |
55 | AliMUONlist *elem[10*1026*1026]; | |
56 | AliMUONlist **ppe=elem; | |
57 | TObjArray *obj=new TObjArray; | |
58 | ||
59 | Int_t digits[3]; | |
60 | // fill the info array | |
61 | TVector *trinfo; | |
62 | trinfo=new TVector(2); | |
63 | ||
64 | // | |
65 | // Loop over events | |
66 | // | |
67 | ||
68 | Int_t Nh=0; | |
69 | Int_t Nh1=0; | |
70 | for (int nev=0; nev<= evNumber2; nev++) { | |
71 | Int_t nparticles = gAlice->GetEvent(nev); | |
72 | cout << "nev " <<nev<<endl; | |
73 | cout << "nparticles " <<nparticles<<endl; | |
74 | if (nev < evNumber1) continue; | |
75 | if (nparticles <= 0) return; | |
76 | ||
77 | TTree *TH = gAlice->TreeH(); | |
78 | Int_t ntracks = TH->GetEntries(); | |
79 | Int_t Nc=0; | |
80 | // | |
81 | Int_t counter=0; | |
82 | ||
83 | // loop over cathodes | |
84 | for (int icat=0;icat<nCathode;icat++) { | |
85 | ||
86 | //initialize the arrray of pointers !!!!! | |
87 | ppe=memset(elem,0,sizeof(void*)*10*1026*1026); | |
88 | ||
89 | printf("Start loop over tracks \n"); | |
90 | // | |
91 | // Loop over events | |
92 | // | |
93 | for (Int_t track=0; track<ntracks;track++) { | |
94 | gAlice->ResetHits(); | |
95 | Int_t nbytes += TH->GetEvent(track); | |
96 | if (MUON) { | |
97 | for(AliMUONhit* mHit=(AliMUONhit*)MUON->FirstHit(-1); | |
98 | mHit; | |
99 | mHit=(AliMUONhit*)MUON->NextHit()) | |
100 | { | |
101 | Int_t nch = mHit->fChamber; // chamber number | |
102 | Float_t x = mHit->fX; // x-pos of hit | |
103 | Float_t y = mHit->fY; // y-pos | |
104 | // | |
105 | // | |
106 | if (nch >10) continue; | |
107 | ||
108 | iChamber = &(MUON->Chamber(nch-1)); | |
109 | response=iChamber->GetResponseModel(); | |
110 | ||
111 | Int_t nsec=iChamber->Nsec(); | |
112 | Int_t rmin = (Int_t)iChamber->frMin; | |
113 | Int_t rmax = (Int_t)iChamber->frMax; | |
114 | // | |
115 | // | |
116 | for (AliMUONcluster* mPad=(AliMUONcluster*)MUON->FirstPad(mHit); | |
117 | mPad; | |
118 | mPad=(AliMUONcluster*)MUON->NextPad()) | |
119 | { | |
120 | Int_t cathode = mPad->fCathode; // chamber number | |
121 | Int_t nhit = mPad->fHitNumber; // hit number | |
122 | Int_t qtot = mPad->fQ; // charge | |
123 | Int_t ipx = mPad->fPadX; // pad number on X | |
124 | Int_t ipy = mPad->fPadY; // pad number on Y | |
125 | Int_t iqpad = mPad->fQpad; // charge per pad | |
126 | Int_t izone = mPad->fRSec; // r-pos of pad | |
127 | // | |
128 | // | |
129 | if (cathode != (icat+1)) continue; | |
130 | segmentation=iChamber->GetSegmentationModel(cathode); | |
131 | Int_t Npx[nch-1] = segmentation->Npx(); | |
132 | Int_t Npy[nch-1] = segmentation->Npy(); | |
133 | ||
134 | Int_t npx=Npx[nch-1]; | |
135 | Int_t npy=Npy[nch-1]; | |
136 | Float_t thex, they; | |
137 | ||
138 | segmentation->GetPadCxy(ipx,ipy,thex,they); | |
139 | Float_t rpad=TMath::Sqrt(thex*thex+they*they); | |
140 | ||
141 | if (rpad < rmin || iqpad ==0 || rpad > rmax) continue; | |
142 | ||
143 | ||
144 | trinfo(0)=(Float_t)track; | |
145 | trinfo(1)=(Float_t)iqpad; | |
146 | ||
147 | digits[0]=ipx; | |
148 | digits[1]=ipy; | |
149 | digits[2]=iqpad; | |
150 | ||
151 | // build the list of fired pads and update the info | |
152 | AliMUONlist *pdigit= | |
153 | elem[(nch-1)*(1026*1026)+(ipy+npy)*1026+(ipx+npx)]; | |
154 | if (pdigit==0) { | |
155 | obj->AddAtAndExpand(new AliMUONlist(nch-1,digits),counter); | |
156 | counter++; | |
157 | Int_t last=obj->GetLast(); | |
158 | pdigit=(AliMUONlist*)obj->At(last); | |
159 | elem[(nch-1)*(1026*1026)+(ipy+npy)*1026+(ipx+npx)]=pdigit; | |
160 | // list of tracks | |
161 | TObjArray *trlist=(TObjArray*)pdigit->TrackList(); | |
162 | trlist->Add(trinfo); | |
163 | } else { | |
164 | // update charge | |
165 | (*pdigit).fSignal+=iqpad; | |
166 | // update list of tracks | |
167 | TObjArray* trlist=(TObjArray*)pdigit->TrackList(); | |
168 | Int_t last_entry=trlist->GetLast(); | |
169 | TVector *ptrk=(TVector*)trlist->At(last_entry); | |
170 | Int_t last_track=Int_t(ptrk(0)); | |
171 | Int_t last_charge=Int_t(ptrk(1)); | |
172 | if (last_track==track) { | |
173 | last_charge+=iqpad; | |
174 | trlist->RemoveAt(last_entry); | |
175 | trinfo(0)=last_track; | |
176 | trinfo(1)=last_charge; | |
177 | trlist->AddAt(trinfo,last_entry); | |
178 | } else { | |
179 | trlist->Add(trinfo); | |
180 | } | |
181 | // check the track list | |
182 | Int_t nptracks=trlist->GetEntriesFast(); | |
183 | if (nptracks > 2) { | |
184 | printf("Attention - nptracks > 2 %d \n",nptracks); | |
185 | printf("cat,nch,ix,iy %d %d %d %d \n",icat+1,nch,ipx,ipy); | |
186 | for (Int_t tr=0;tr<nptracks;tr++) { | |
187 | TVector *pptrk=(TVector*)trlist->At(tr); | |
188 | trk[tr]=Int_t(pptrk(0)); | |
189 | chtrk[tr]=Int_t(pptrk(1)); | |
190 | } | |
191 | } // end if nptracks | |
192 | } // end if pdigit | |
193 | } //end loop over clusters | |
194 | } // hit loop | |
195 | } // if MUON | |
196 | } // track loop | |
197 | ||
198 | Int_t tracks[10]; | |
199 | Int_t charges[10]; | |
200 | cout<<"start filling digits "<<endl; | |
201 | const Int_t zero_supm = 6; | |
202 | const Int_t adc_satm = 1024; | |
203 | Int_t nd=0; | |
204 | ||
205 | ||
206 | Int_t nentries=obj->GetEntriesFast(); | |
207 | printf(" nentries %d \n",nentries); | |
208 | ||
209 | // start filling the digits | |
210 | ||
211 | for (Int_t nent=0;nent<nentries;nent++) { | |
212 | AliMUONlist *address=(AliMUONlist*)obj->At(nent); | |
213 | if (address==0) continue; | |
214 | // do zero-suppression and signal truncation | |
215 | Int_t ich=address->fRpad; | |
216 | Int_t q=address->fSignal; | |
217 | // if ( q <= zero_supm || rpad > 55) continue; | |
218 | if ( q <= zero_supm ) continue; | |
219 | if ( q > adc_satm) q=adc_satm; | |
220 | digits[0]=address->fPadX; | |
221 | digits[1]=address->fPadY; | |
222 | digits[2]=q; | |
223 | ||
224 | TObjArray* trlist=(TObjArray*)address->TrackList(); | |
225 | Int_t nptracks=trlist->GetEntriesFast(); | |
226 | // this was changed to accomodate the real number of tracks | |
227 | if (nptracks > 10) { | |
228 | cout<<"Attention - nptracks > 3 "<<nptracks<<endl; | |
229 | nptracks=10; | |
230 | } | |
231 | if (nptracks > 2) { | |
232 | printf("Attention - nptracks > 2 %d \n",nptracks); | |
233 | printf("cat,ich,ix,iy,q %d %d %d %d %d \n",icat,ich,digits[0],digits[1],q); | |
234 | } | |
235 | for (Int_t tr=0;tr<nptracks;tr++) { | |
236 | TVector *pp=(TVector*)trlist->At(tr); | |
237 | Int_t entries=pp->GetNrows(); | |
238 | tracks[tr]=Int_t(pp(0)); | |
239 | charges[tr]=Int_t(pp(1)); | |
240 | } //end loop over list of tracks for one pad | |
241 | ||
242 | // fill digits | |
243 | MUON->AddDigits(ich,tracks,charges,digits); | |
244 | } | |
245 | cout<<"I'm out of the loops for digitisation"<<endl; | |
246 | gAlice->TreeD()->Fill(); | |
247 | TTree *TD=gAlice->TreeD(); | |
248 | int ndig=TD->GetEntries(); | |
249 | cout<<"number of digits "<<ndig<<endl; | |
250 | TClonesArray *fDch; | |
251 | for (int i=0;i<10;i++) { | |
252 | fDch= MUON->DigitsAddress(i); | |
253 | int ndig=fDch->GetEntriesFast(); | |
254 | printf (" i, ndig %d %d \n",i,ndig); | |
255 | } | |
256 | MUON->ResetDigits(); | |
257 | obj->Clear(); | |
258 | } //end loop over cathodes | |
259 | char hname[30]; | |
260 | sprintf(hname,"TreeD%d",nev); | |
261 | gAlice->TreeD()->Write(hname); | |
262 | file->ls(); | |
263 | } // event loop | |
264 | file->Close(); | |
265 | // cout<<"END digitisation "<<endl; | |
266 | } |