1 void MUONdigits (Int_t evNumber1=0,Int_t evNumber2=0,Int_t nCathode=1)
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.
8 /////////////////////////////////////////////////////////////////////////
10 // Dynamically link some shared libs
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
23 // Connect the Root Galice file containing Geometry, Kine and Hits
25 TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
26 if (file) file->Close();
27 file = new TFile("galice.root","UPDATE");
31 printf ("I'm after Map \n");
33 // Get AliRun object from file or create it if not on file
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");
40 printf ("I'm after gAlice \n");
42 AliMUON *MUON = gAlice->GetDetector("MUON");
44 AliMUONchamber* iChamber;
45 AliMUONsegmentation* segmentation;
55 AliMUONlist *elem[10*1026*1026];
56 AliMUONlist **ppe=elem;
57 TObjArray *obj=new TObjArray;
60 // fill the info array
62 trinfo=new TVector(2);
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;
77 TTree *TH = gAlice->TreeH();
78 Int_t ntracks = TH->GetEntries();
84 for (int icat=0;icat<nCathode;icat++) {
86 //initialize the arrray of pointers !!!!!
87 ppe=memset(elem,0,sizeof(void*)*10*1026*1026);
89 printf("Start loop over tracks \n");
93 for (Int_t track=0; track<ntracks;track++) {
95 Int_t nbytes += TH->GetEvent(track);
97 for(AliMUONhit* mHit=(AliMUONhit*)MUON->FirstHit(-1);
99 mHit=(AliMUONhit*)MUON->NextHit())
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
106 if (nch >10) continue;
108 iChamber = &(MUON->Chamber(nch-1));
109 response=iChamber->GetResponseModel();
111 Int_t nsec=iChamber->Nsec();
112 Int_t rmin = (Int_t)iChamber->frMin;
113 Int_t rmax = (Int_t)iChamber->frMax;
116 for (AliMUONcluster* mPad=(AliMUONcluster*)MUON->FirstPad(mHit);
118 mPad=(AliMUONcluster*)MUON->NextPad())
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
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();
134 Int_t npx=Npx[nch-1];
135 Int_t npy=Npy[nch-1];
138 segmentation->GetPadCxy(ipx,ipy,thex,they);
139 Float_t rpad=TMath::Sqrt(thex*thex+they*they);
141 if (rpad < rmin || iqpad ==0 || rpad > rmax) continue;
144 trinfo(0)=(Float_t)track;
145 trinfo(1)=(Float_t)iqpad;
151 // build the list of fired pads and update the info
153 elem[(nch-1)*(1026*1026)+(ipy+npy)*1026+(ipx+npx)];
155 obj->AddAtAndExpand(new AliMUONlist(nch-1,digits),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;
161 TObjArray *trlist=(TObjArray*)pdigit->TrackList();
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) {
174 trlist->RemoveAt(last_entry);
175 trinfo(0)=last_track;
176 trinfo(1)=last_charge;
177 trlist->AddAt(trinfo,last_entry);
181 // check the track list
182 Int_t nptracks=trlist->GetEntriesFast();
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));
193 } //end loop over clusters
200 cout<<"start filling digits "<<endl;
201 const Int_t zero_supm = 6;
202 const Int_t adc_satm = 1024;
206 Int_t nentries=obj->GetEntriesFast();
207 printf(" nentries %d \n",nentries);
209 // start filling the digits
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;
224 TObjArray* trlist=(TObjArray*)address->TrackList();
225 Int_t nptracks=trlist->GetEntriesFast();
226 // this was changed to accomodate the real number of tracks
228 cout<<"Attention - nptracks > 3 "<<nptracks<<endl;
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);
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
243 MUON->AddDigits(ich,tracks,charges,digits);
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;
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);
258 } //end loop over cathodes
260 sprintf(hname,"TreeD%d",nev);
261 gAlice->TreeD()->Write(hname);
265 // cout<<"END digitisation "<<endl;