[u/mrichter/AliRoot.git] / MUON / MUONdigit.C

void MUONdigit (Int_t evNumber1=0,Int_t evNumber2=0,Int_t nCathode=1) 
{
/////////////////////////////////////////////////////////////////////////
//                                                                     //
//   This macros converts the pad-hit information into digits          //
//                                                                     //
/////////////////////////////////////////////////////////////////////////

// Dynamically link some shared libs

   if (gClassTable->GetID("AliRun") < 0) {
       gSystem->Load("$ALITOP/cern.so/lib/libpdfDUMMY.so");
       gSystem->Load("$ALITOP/cern.so/lib/libPythia.so");
       gSystem->Load("$ROOTSYS/lib/libEG.so");       
       gSystem->Load("$ROOTSYS/lib/libEGPythia.so");    
       gSystem->Load("libGeant3Dummy.so");        //a dummy version of Geant3
       gSystem->Load("PHOS/libPHOSdummy.so");     //the standard Alice classes 
       gSystem->Load("libgalice.so");             //the standard Alice classes 
   }

// Connect the Root Galice file containing Geometry, Kine and Hits

   TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
   if (file) file->Close(); 
   file = new TFile("galice.root","UPDATE");
   file->ls();
//   file->Map();

   printf ("\n File loaded !");

// Get AliRun object from file or create it if not on file

   if (!gAlice) {
       gAlice = (AliRun*)file->Get("gAlice");
       if (gAlice) printf("\n AliRun object found on file\n");
       if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
   }
   printf ("\n gAlice created !\n");

   AliMUON *MUON  = gAlice->GetDetector("MUON");

   AliMUONchamber*  iChamber;
   AliMUONsegmentation*  segmentation;

   Int_t Npx[10];
   Int_t Npy[10];

   Int_t nxmax=1026;
   Int_t nymax=1026;
   AliMUONlist *elem[10][1026][1026];
   TObjArray *obj=new TObjArray;

   Int_t trk[500];
   Int_t chtrk[500];


   Int_t digits[3]; 
   TVector *trinfo;
   trinfo=new TVector(2);

//
//   Loop over events 
//

   Int_t Nh=0;
   Int_t Nh1=0;
   for (int nev=0; nev<= evNumber2; nev++) {
       Int_t nparticles = gAlice->GetEvent(nev);
       cout << "nev         " <<nev<<endl;
       cout << "nparticles  " <<nparticles<<endl;
       if (nev < evNumber1) continue;
       if (nparticles <= 0) return;

       TTree *TH = gAlice->TreeH();
       Int_t ntracks = TH->GetEntries();
       Int_t Nc=0;

       Int_t counter=0;
       //
       // loop over cathodes
       for (int icat=0; icat<nCathode; icat++) { 
	   //
	   // initialize the array of pointers 
	   printf("\n Initialize array of pointers \n");
	   for (int i=0; i<10; i++) {
	       for (int j=0; j<nymax; j++) {
		   for (int k=0; k<nxmax; k++) {	
		       elem[i][j][k]=0;
		   }
	       }
	   }
	   printf("\n Finished !\n");
     
//
//   Loop over tracks
//
	   printf("\n ntracks: %d",ntracks);
	   
	   for (Int_t track=0; track<ntracks;track++) {
	       gAlice->ResetHits();
	       Int_t nbytes += TH->GetEvent(track);
	       if (MUON)  {
		   for(AliMUONhit* mHit=(AliMUONhit*)MUON->FirstHit(-1); 
		       mHit;
		       mHit=(AliMUONhit*)MUON->NextHit()) 
		   {
		       Int_t   nch   = mHit->fChamber;  // chamber number
		       Float_t x     = mHit->fX;        // x-pos of hit
		       Float_t y     = mHit->fY;        // y-pos
		       

		       iChamber = &(MUON->Chamber(nch-1));
		       response=iChamber->GetResponseModel();
		       
		       Int_t nsec=iChamber->Nsec();
		       Int_t rmin = (Int_t)iChamber->frMin;
		       Int_t rmax = (Int_t)iChamber->frMax;
		       if (nch > 10) continue;
		       
//
//
		       for (AliMUONcluster* mPad=
				(AliMUONcluster*)MUON->FirstPad(mHit);
			    mPad;
			    mPad=(AliMUONcluster*)MUON->NextPad())
		       {
			   Int_t cathode  = mPad->fCathode;   // cathode number
			   Int_t nhit     = mPad->fHitNumber; // hit number
			   Int_t qtot     = mPad->fQ;         // charge
			   Int_t ipx      = mPad->fPadX;      // pad  X
			   Int_t ipy      = mPad->fPadY;      // pad  Y
			   Int_t iqpad    = mPad->fQpad;      // charge per pad
			   Int_t rsec     = mPad->fRSec;      // sector#
//
//
			   if (cathode != (icat+1)) continue;
			   segmentation=iChamber-> 
			       GetSegmentationModel(cathode);
			   Int_t Npx[nch-1]  = segmentation->Npx();
			   Int_t Npy[nch-1]  = segmentation->Npy();
			   
			   
			   Int_t npx=Npx[nch-1];
			   Int_t npy=Npy[nch-1];
			   Float_t thex, they;
			   segmentation->GetPadCxy(ipx,ipy,thex,they);
			   
			   Float_t rpad=TMath::Sqrt(thex*thex+they*they);
			   
			   
                	// check boundaries
			   if (rpad < rmin || iqpad ==0 || rpad > rmax) continue;

			   // fill the info array

			   trinfo(0)=(Float_t)track;
			   trinfo(1)=(Float_t)iqpad;

			   digits[0]=ipx;
			   digits[1]=ipy;
			   digits[2]=iqpad;

			   // build the list of fired pads and update the info
			   
			   AliMUONlist *pdigit=elem[nch-1][ipy+npy][ipx+npx];

			   if (pdigit==0) {
			       obj->AddAtAndExpand(
				   new AliMUONlist(rpad,digits),counter);
			       counter++;
			       Int_t last=obj->GetLast();
			       pdigit=(AliMUONlist*)obj->At(last);
			       elem[nch-1][ipy+npy][ipx+npx]=pdigit;
			       // list of tracks
			       TObjArray *trlist=
				   (TObjArray*)pdigit->TrackList();
			       trlist->Add(trinfo);
			   } else {
			       // update charge
			       (*pdigit).fSignal+=iqpad;
			       // update list of tracks
			       TObjArray* trlist=(TObjArray*)pdigit
				   ->TrackList();
			       Int_t nptracks=trlist->GetEntriesFast();
			       
			       for (Int_t tr=0;tr<nptracks;tr++) {
				   TVector *ptrk=(TVector*)trlist->At(tr);
//				   TVector &vtrk = *ptrk;
				   Int_t entries=ptrk->GetNrows();
				   trk[tr]=Int_t(ptrk(0));
				   chtrk[tr]=Int_t(ptrk(1));
				   if (trk[tr]==track) {
				       chtrk[tr]+=iqpad;
				       trlist->RemoveAt(tr);
				       trinfo(0)=trk[tr];
				       trinfo(1)=chtrk[tr];
				       trlist->AddAt(trinfo,tr);
				   } else {
				       trlist->Add(trinfo);
				   }
			       } // end loop over list of tracks for one pad
			   } // end if pdigit
		       } //end loop over clust
		   } // hit loop
	       } // if MUON
	   } // track loop

	   Int_t tracks[10];
	   Int_t charges[10];
	   cout<<"Start filling digits now !"<<endl;
	   const Int_t zero_supm = 6;
	   const Int_t adc_satm = 1024;
	   Int_t nd=0;
	   

	   // start filling the digits
	   for (Int_t id=0; id<10; id++) {
	       Int_t nd=0;
	       for (Int_t iy=0;iy<nymax;iy++) {
		   for (Int_t ix=0;ix<nxmax;ix++) {
		       AliMUONlist *address=elem[id][iy][ix];
		       if (address==0) continue; 
		       // do zero-suppression and signal truncation
		       Int_t q=address->fSignal;        
		       // 
		       if ( q <= zero_supm ) continue;
		       if ( q > adc_satm)  q=adc_satm;
		       digits[0]=address->fPadX;
		       digits[1]=address->fPadY;
		       digits[2]=address->fSignal;

		       TObjArray* trlist=(TObjArray*)address->TrackList();
		       Int_t nptracks=trlist->GetEntriesFast();
		       // this was changed to accomodate the real number of tracks
		       if (nptracks > 10) {
			   cout<<"Attention - nptracks > 10 "<<nptracks<<endl;
			   nptracks=10;
		       }
		       for (Int_t tr=0;tr<nptracks;tr++) {
			   TVector *pp=(TVector*)trlist->At(tr);
			   Int_t entries=pp->GetNrows();
			   TVector &v2 = *pp;
			   tracks[tr]=(Int_t)v2(0);
			   charges[tr]=(Int_t)v2(1);
		       } //end loop over list of tracks for one pad

		       // fill digits
		       MUON->AddDigits(id,tracks,charges,digits);
		       nd++;
		   } // end loop over ix
	       } // end loop over iy
	       cout<<"I'm out of the loops over pads"<<endl;
	       cout<<" id  nd  "<<id<<" "<<nd<<endl;
	   } //loop over chambers
	   cout<<"I'm out of the loops for digitisation"<<endl;

	   gAlice->TreeD()->Fill();
	   TTree *TD=gAlice->TreeD();
	   int ndig=TD->GetEntries();
	   cout<<"number of digits  "<<ndig<<endl;

	   TClonesArray *fDch;
	   for (int i=0;i<10;i++) {
	       fDch= MUON->DigitsAddress(i);
	       int ndig=fDch->GetEntriesFast();
	       printf (" i, ndig %d %d \n",i,ndig);
	   }


	   MUON->ResetDigits();
	   obj->Clear();
       } //end loop over cathodes

       char hname[30];
       sprintf(hname,"TreeD%d",nev);
       gAlice->TreeD()->Write(hname);
       
       file->ls();
//       file->Map();
   } // event loop 
   obj->Delete();
   delete obj;

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