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

void MUONdigits (Int_t evNumber1=0,Int_t evNumber2=0,Int_t nCathode=1) 
{
/////////////////////////////////////////////////////////////////////////
//   This macro is a small example of a ROOT macro
//   illustrating how to read the output of GALICE
//   and do some analysis.
//   
/////////////////////////////////////////////////////////////////////////

// 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 ("I'm after Map \n");

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

   if (!gAlice) {
       gAlice = (AliRun*)file->Get("gAlice");
       if (gAlice) printf("AliRun object found on file\n");
       if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
   }
   printf ("I'm after gAlice \n");
   
   AliMUON *MUON  = gAlice->GetDetector("MUON");
   
   AliMUONchamber*  iChamber;
   AliMUONsegmentation*  segmentation;
   
   Int_t Npx[10];
   Int_t Npy[10];
   Int_t trk[50];
   Int_t chtrk[50];


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

   Int_t digits[3]; 
   // fill the info array
   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 arrray of pointers !!!!!
	   ppe=memset(elem,0,sizeof(void*)*10*1026*1026);

	   printf("Start loop over tracks \n");     
//
//   Loop over events
//
	   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
//
//
		       if (nch >10) continue;
		   
		       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;
//
//
		       for (AliMUONcluster* mPad=(AliMUONcluster*)MUON->FirstPad(mHit);
			    mPad;
			    mPad=(AliMUONcluster*)MUON->NextPad())
		       {
			   Int_t cathode     = mPad->fCathode;   // chamber number
			   Int_t nhit     = mPad->fHitNumber; // hit number
			   Int_t qtot     = mPad->fQ;         // charge
			   Int_t ipx      = mPad->fPadX;      // pad number on X
			   Int_t ipy      = mPad->fPadY;      // pad number on Y
			   Int_t iqpad    = mPad->fQpad;      // charge per pad
			   Int_t izone    = mPad->fRSec;      // r-pos of pad
//
//
			   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);
			   
			   if (rpad < rmin || iqpad ==0 || rpad > rmax) continue;

			
			   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)*(1026*1026)+(ipy+npy)*1026+(ipx+npx)];
			   if (pdigit==0) {
			       obj->AddAtAndExpand(new AliMUONlist(nch-1,digits),counter);
			       counter++;
			       Int_t last=obj->GetLast();
			       pdigit=(AliMUONlist*)obj->At(last);
			       elem[(nch-1)*(1026*1026)+(ipy+npy)*1026+(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 last_entry=trlist->GetLast();
			       TVector *ptrk=(TVector*)trlist->At(last_entry);
			       Int_t last_track=Int_t(ptrk(0));
			       Int_t last_charge=Int_t(ptrk(1));
			       if (last_track==track) {
				   last_charge+=iqpad;
				   trlist->RemoveAt(last_entry);
				   trinfo(0)=last_track;
				   trinfo(1)=last_charge;
				   trlist->AddAt(trinfo,last_entry);
			       } else {
				   trlist->Add(trinfo);
			       }
			       // check the track list
			       Int_t nptracks=trlist->GetEntriesFast();
			       if (nptracks > 2) {
				   printf("Attention - nptracks > 2  %d \n",nptracks);
				   printf("cat,nch,ix,iy %d %d %d %d  \n",icat+1,nch,ipx,ipy);
				   for (Int_t tr=0;tr<nptracks;tr++) {
				       TVector *pptrk=(TVector*)trlist->At(tr);
				       trk[tr]=Int_t(pptrk(0));
				       chtrk[tr]=Int_t(pptrk(1));
				   }
			       } // end if nptracks
			   } //  end if pdigit
		       } //end loop over clusters
		   } // hit loop
	       } // if MUON
	   } // track loop

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


	   Int_t nentries=obj->GetEntriesFast();
	   printf(" nentries %d \n",nentries);

	   // start filling the digits
	   
	   for (Int_t nent=0;nent<nentries;nent++) {
	       AliMUONlist *address=(AliMUONlist*)obj->At(nent);
	       if (address==0) continue; 
	       // do zero-suppression and signal truncation
	       Int_t ich=address->fRpad;
	       Int_t q=address->fSignal;        
	       // if ( q <= zero_supm  || rpad > 55) continue;
	       if ( q <= zero_supm ) continue;
	       if ( q > adc_satm)  q=adc_satm;
	       digits[0]=address->fPadX;
	       digits[1]=address->fPadY;
	       digits[2]=q;

	       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 > 3 "<<nptracks<<endl;
		   nptracks=10;
	       }
               if (nptracks > 2) {
		   printf("Attention - nptracks > 2  %d \n",nptracks);
		   printf("cat,ich,ix,iy,q %d %d %d %d %d \n",icat,ich,digits[0],digits[1],q);
               }
	       for (Int_t tr=0;tr<nptracks;tr++) {
		   TVector *pp=(TVector*)trlist->At(tr);
		   Int_t entries=pp->GetNrows();
		   tracks[tr]=Int_t(pp(0));
		   charges[tr]=Int_t(pp(1));
	       }      //end loop over list of tracks for one pad
	       
	       // fill digits
	       MUON->AddDigits(ich,tracks,charges,digits);
	   }
	   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();
   } // event loop 
   file->Close();
//   cout<<"END  digitisation "<<endl;
}
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[101026*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)1010261026);
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(thexthex+theythey);
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)(10261026)+(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)(10261026)+(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	}