Applying fixes and updates (Indra)
[u/mrichter/AliRoot.git] / HLT / MUON / OnlineAnalysis / AliHLTMUONTriggerReconstructor.cxx
1 /**************************************************************************
2  * This file is property of and copyright by the ALICE HLT Project        * 
3  * All rights reserved.                                                   *
4  *                                                                        *
5  * Primary Authors:                                                       *
6  *   Indranil Das <indra.das@saha.ac.in>                                  *
7  *                                                                        *
8  * Permission to use, copy, modify and distribute this software and its   *
9  * documentation strictly for non-commercial purposes is hereby granted   *
10  * without fee, provided that the above copyright notice appears in all   *
11  * copies and that both the copyright notice and this permission notice   *
12  * appear in the supporting documentation. The authors make no claims     *
13  * about the suitability of this software for any purpose. It is          * 
14  * provided "as is" without express or implied warranty.                  *
15  **************************************************************************/
16
17 /* $Id$ */
18
19 /**********************************************************************
20  Created on : 16/05/2007
21  Purpose    : This class is supposed to read the tracker DDL files and 
22               give the output AliMUONCoreTriggerRecord
23  Author     : Indranil Das, HEP Division, SINP
24  Email      : indra.das@saha.ac.in | indra.ehep@gmail.com
25 **********************************************************************/
26
27 ///*
28 //
29 //  The trigger reconstructor class is designed to deal the rawdata inputfiles
30 //  to findout the the reconstructed hits at the trigger DDL. The output is send
31 //  to the output block for further processing.
32 //
33 //  Author : Indranil Das ( indra.das@saha.ac.in || indra.ehep@gmail.com )
34 // 
35 //*/
36
37 #include "AliHLTMUONTriggerReconstructor.h"
38
39 const int AliHLTMUONTriggerReconstructor::fgkDetectorId = 0xB00;
40 const int AliHLTMUONTriggerReconstructor::fgkDDLOffSet = 20 ;
41 const int AliHLTMUONTriggerReconstructor::fgkNofDDL = 2 ;
42
43 const int AliHLTMUONTriggerReconstructor::fgkDDLHeaderSize = 8;
44 const int AliHLTMUONTriggerReconstructor::fgkEvenLutSize = 2602351 + 1; 
45 const int AliHLTMUONTriggerReconstructor::fgkOddLutSize = 2528735 + 1;
46
47 const int AliHLTMUONTriggerReconstructor::fgkLutLine = 10496;
48
49 const int AliHLTMUONTriggerReconstructor::fgkMinIdManuChannel[2] = {819616, 862288};
50 const int AliHLTMUONTriggerReconstructor::fgkMaxIdManuChannel[2] = {3421966, 3391022};
51
52 const float AliHLTMUONTriggerReconstructor::fgkHalfPadSizeXB[3] = {8.5, 17.0, 25.5};
53 const float AliHLTMUONTriggerReconstructor::fgkHalfPadSizeYNB[2] = {25.5, 34.0};
54
55 const int AliHLTMUONTriggerReconstructor::fgkDetElem = 9*4 ; // 9 detele per half chamber
56
57
58 AliHLTMUONTriggerReconstructor::AliHLTMUONTriggerReconstructor()
59   :
60   fPadData(NULL),
61   fLookUpTableData(NULL),
62   fRecPoints(NULL),
63   fRecPointsCount(NULL),
64   fMaxRecPointsCount(0),
65   fMaxFiredPerDetElem(),
66   fDetElemToDataId(),
67   fDDLId(0),
68   fIdOffSet(0)
69 {
70   // ctor 
71   
72   if(AliHLTMUONTriggerReconstructor::fgkEvenLutSize > AliHLTMUONTriggerReconstructor::fgkOddLutSize){
73     fPadData = new AliHLTMUONHitReconstructor::DHLTPad[AliHLTMUONTriggerReconstructor::fgkEvenLutSize];
74   }
75   else{
76     fPadData = new AliHLTMUONHitReconstructor::DHLTPad[AliHLTMUONTriggerReconstructor::fgkOddLutSize];
77   }
78
79   bzero(fGetIdTotalData,104*64*2*sizeof(int));
80 }
81
82
83 AliHLTMUONTriggerReconstructor::~AliHLTMUONTriggerReconstructor()
84 {
85   // dtor
86   delete []fPadData;
87   delete []fLookUpTableData;
88 }
89
90 bool AliHLTMUONTriggerReconstructor::SetRegToLocCardMap(RegToLoc* regToLoc)
91 {
92   if(!memcpy(fRegToLocCard,regToLoc,128*sizeof(RegToLoc)))
93     return false;
94
95   for(int i=0;i<128;i++){
96     HLTDebug("DDL : %d, reg : %d, loc : %d",fRegToLocCard[i].fTrigDDL,
97             fRegToLocCard[i].fRegId,fRegToLocCard[i].fLocId);
98   }
99
100   return true;
101 }
102
103 bool AliHLTMUONTriggerReconstructor::LoadLookUpTable(AliHLTMUONHitReconstructor::DHLTLut* lookUpTableData, int lookUpTableId)
104 {
105   if(lookUpTableId<fgkDDLOffSet || lookUpTableId>= fgkDDLOffSet + fgkNofDDL){
106     HLTError("DDL number is out of range (must be %d<=iDDL<%d)",fgkDDLOffSet,fgkDDLOffSet+fgkNofDDL);
107     return false;
108   }
109   
110   fDDLId = lookUpTableId;
111
112   int lutSize = ((lookUpTableId%2)==0) ? fgkEvenLutSize : fgkOddLutSize ;
113   int nofLutLine = fgkLutLine ;
114   fIdOffSet = fgkMinIdManuChannel[lookUpTableId%2];
115
116   int detManuChannelId;
117
118   fLookUpTableData = new AliHLTMUONHitReconstructor::DHLTLut[lutSize];
119
120   memset(fLookUpTableData,-1,lutSize*sizeof(AliHLTMUONHitReconstructor::DHLTLut));
121
122   for(int i=0; i<nofLutLine; i++){
123
124     detManuChannelId = lookUpTableData[i].fIdManuChannel - fIdOffSet + 1;
125     fLookUpTableData[detManuChannelId].fIdManuChannel = lookUpTableData[i].fIdManuChannel - fIdOffSet;
126     fLookUpTableData[detManuChannelId].fIX = lookUpTableData[i].fIX ;
127     fLookUpTableData[detManuChannelId].fIY = lookUpTableData[i].fIY ;
128     fLookUpTableData[detManuChannelId].fRealX = lookUpTableData[i].fRealX ;
129     fLookUpTableData[detManuChannelId].fRealY = lookUpTableData[i].fRealY ;
130     fLookUpTableData[detManuChannelId].fRealZ = lookUpTableData[i].fRealZ ;
131     fLookUpTableData[detManuChannelId].fPcbZone = lookUpTableData[i].fPcbZone ;
132     fLookUpTableData[detManuChannelId].fPlane = lookUpTableData[i].fPlane ;
133   }
134   return true;
135
136   return true;
137 }
138
139 bool AliHLTMUONTriggerReconstructor::Run(int *rawData, int *rawDataSize, AliHLTMUONTriggerRecordStruct trigRecord[], int *nofTrigRec)
140 {
141
142   fRecPoints = &trigRecord[0];
143   fMaxRecPointsCount = *nofTrigRec;
144   fRecPointsCount = nofTrigRec;
145   *fRecPointsCount = 0;
146   fMaxFiredPerDetElem.clear();
147   fDetElemToDataId.clear();
148
149   fPadData[0].fDetElemId = 0;
150   fPadData[0].fBuspatchId = 0;
151   fPadData[0].fIdManuChannel = 0;
152   fPadData[0].fIX = 0 ;
153   fPadData[0].fIY = 0 ;
154   fPadData[0].fRealX = 0.0 ;
155   fPadData[0].fRealY = 0.0 ;
156   fPadData[0].fRealZ = 0.0 ;
157   fPadData[0].fPlane = -1 ;
158   fPadData[0].fPcbZone = -1 ;
159   fPadData[0].fCharge = 0 ;
160   
161   if(!ReadDDL(rawData,rawDataSize)){
162     HLTError("Failed to read the complete DDL file\n");
163     return false;
164   }
165
166   if(!FindTrigHits()){
167     HLTError("Failed to generate RecHits\n");
168     return false;
169   }
170     
171   return true;
172 }
173
174
175 bool AliHLTMUONTriggerReconstructor::ReadDDL(int *rawData, int *rawDataSize)
176 {
177
178   int idManuChannel ;
179   
180   int index = 0;
181   int dataCount = 0;
182   int detElemId = 0 ;
183   int reg_output,reg_phys_trig_occur ;
184   int iLocIndex,loc,locDec,triggY,sign,loDev,triggX;
185   int iRegLoc, locId ;
186   short pattern[2][4]; // 2 stands for two cathode planes and 4 stands for 4 chambers
187
188   Int_t offset,ithSwitch,secondLocation,idetElemId;
189
190   int shiftIndex = 10 - 6 - 1; // the one comes due to indexing from zero
191
192   DataIdIndex dataIndex;
193 #ifdef DEBUG
194   int globalcard_data_occurance = (rawData[index]>>10)&0x1; //Set to 1 if global info present in DDL else set to 0 
195   int version = (rawData[index]>>12)&0xFF; // software version
196   int serial_number =  (rawData[index]>>20)&0xF; // serial number set to 0xF 
197 #endif
198   int phys_trig_occur = (rawData[index]>>30)&0x1; // 1 for physics trigger, 0 for software trigger
199   
200   HLTDebug("globalcard_data_occurance  %d, version  %d, serial_number  %d, phys_trig_occur  %d",
201          globalcard_data_occurance,version,serial_number,phys_trig_occur);
202
203   if(!phys_trig_occur) // for software trigger
204     index += 8 ;// corresponding to scalar words
205   
206   index += 1 ; // To skip the separator 0xDEADFACE
207   
208   index += 4 ; // corresponding to global input
209   
210   index += 1 ; // reaches to global output
211   
212   if((fDDLId - AliHLTMUONTriggerReconstructor::fgkDDLOffSet) == 0){ //if globalData is present in DDL 0 (presummed may be changed)
213 #ifdef DEBUG
214     int singleLpt = rawData[index] & 0x1;
215     int singleHpt = (rawData[index] >> 1) & 0x1;
216     
217     int pairUnlikeLpt = (rawData[index] >> 4)  & 0x1;
218     int pairUnlikeHpt = (rawData[index] >> 5)  & 0x1;
219     
220     int pairLikeLpt = (rawData[index] >> 2)  & 0x1;
221     int pairLikeHpt = (rawData[index] >> 3)  & 0x1;
222 #endif
223     HLTDebug("singleLpt : %x, singleHpt : %x, pairUnlikeLpt : %x, pairUnlikeHpt : %x, pairLikeLpt : %x, pairLikeHpt : %x",
224              singleLpt,singleHpt,pairUnlikeLpt,pairUnlikeHpt,pairLikeLpt,pairLikeHpt);
225   }
226
227   if(!phys_trig_occur)
228     index += 10 ;// corresponds to scalar words
229
230   index += 1; // separator 0xDEADBEEF 
231
232   for (int iReg = 0; iReg < 8; iReg++) {
233     index += 1; // DARC Status Word
234     index += 1; // Regeional Word
235     reg_output = rawData[index] & 0xFF;
236     reg_phys_trig_occur = ( rawData[index] >> 31) & 0x1;
237     
238     index += 2; // 2 words for regional input
239     
240     index += 1; // L0 counter
241     
242     if(!reg_phys_trig_occur)
243       index += 10;
244     
245     index += 1 ; // end of Regeonal header 0xBEEFFACE
246
247     for(int iLoc = 0; iLoc < 16 ; iLoc++){
248
249       iLocIndex = index ;      
250
251       loc = (rawData[index+5] >> 19) &  0xF ;
252       
253       locDec = (rawData[index+5] >> 15) & 0xF;
254       triggY = (rawData[index+5] >> 14) & 0x1 ;
255       sign = (rawData[index+5] >> 9) & 0x1;
256       loDev = (rawData[index+5] >> 5) & 0xF ;
257       triggX = (loDev >> 4 & 0x1 ) && !(loDev & 0xF) ;
258
259       if( locDec != 0x9 ){ // check for Dec
260         
261         iRegLoc = iReg*16 + iLoc;
262         locId = fRegToLocCard[iRegLoc].fLocId ; 
263         
264         if(locId<=234){ // to avoid the copy locCards
265           
266           index += 1;
267           pattern[0][0] = rawData[index] & 0xFFFF; // x-strip pattern for chaber 0 
268           pattern[0][1] = (rawData[index] >> 16) & 0xFFFF; // x-strip pattern for chaber 1
269           index += 1; 
270           pattern[0][2] = rawData[index] & 0xFFFF; 
271           pattern[0][3] = (rawData[index] >> 16) & 0xFFFF; 
272           
273           index += 1;
274           pattern[1][0] = rawData[index] & 0xFFFF; // y-strip pattern for chaber 0
275           pattern[1][1] = (rawData[index] >> 16) & 0xFFFF; // y-strip pattern for chaber 0 
276           index += 1; 
277           pattern[1][2] = rawData[index] & 0xFFFF; 
278           pattern[1][3] = (rawData[index] >> 16) & 0xFFFF; 
279           
280           if(pattern[0][0] || pattern[0][1] || pattern[0][2] || pattern[0][3]
281              || pattern[1][0] || pattern[1][1] || pattern[1][2] || pattern[1][3]
282              ){
283
284             HLTDebug("iReg: %d, iLoc :%d, locId : %d,X : %x, %x, %x, %x ...Y : %x, %x, %x, %x",
285                     iReg,iLoc,locId,pattern[0][0],pattern[0][1],pattern[0][2],pattern[0][3],
286                     pattern[1][0],pattern[1][1],pattern[1][2],pattern[1][3]);
287
288             for(int iChamber = 0; iChamber < 4 ; iChamber++){ //4 chambers per DDL 
289               for(int iPlane = 0; iPlane < 2 ; iPlane++){// 2 cathode plane
290                 if(pattern[iPlane][iChamber]){
291                   detElemId = fRegToLocCard[iRegLoc].fDetElemId[iChamber];
292                   HLTDebug("\tdetElemId : %d\n",detElemId);
293                   for (Int_t ibitxy = 0; ibitxy < 16; ++ibitxy) {
294                     if ((pattern[iPlane][iChamber] >> ibitxy) & 0x1) {
295           
296                       // not quite sure about this
297                       offset = 0;
298                       ithSwitch = (fRegToLocCard[iRegLoc].fSwitch >> shiftIndex) & 0x1;
299                       if (iPlane && ithSwitch) offset = -8;
300                       
301                       secondLocation = ibitxy + offset;
302                       
303                       idetElemId = detElemId%1000;
304
305                       idetElemId &= 0x1FF ;
306                       iPlane &= 0x1 ;
307                       locId &= 0xFF ;
308                       secondLocation &= 0xF ;
309                       
310                       idManuChannel &= 0x0;
311                       idManuChannel = (idManuChannel|idetElemId)<<1;  
312                       idManuChannel = (idManuChannel|iPlane)<<8;  
313                       idManuChannel = (idManuChannel|locId)<<4 ;
314                       idManuChannel |= secondLocation  ;
315
316                       idManuChannel -= fIdOffSet ;
317                       
318                       if(fLookUpTableData[idManuChannel+1].fIdManuChannel == -1) //skip uninitialized values
319                         continue;
320
321                       fPadData[idManuChannel].fDetElemId = detElemId;
322                       fPadData[idManuChannel].fIdManuChannel = idManuChannel;
323                       fPadData[idManuChannel].fIX = fLookUpTableData[idManuChannel+1].fIX;
324                       fPadData[idManuChannel].fIY = fLookUpTableData[idManuChannel+1].fIY;
325                       fPadData[idManuChannel].fRealX = fLookUpTableData[idManuChannel+1].fRealX;
326                       fPadData[idManuChannel].fRealY = fLookUpTableData[idManuChannel+1].fRealY;
327                       fPadData[idManuChannel].fRealZ = fLookUpTableData[idManuChannel+1].fRealZ;
328                       fPadData[idManuChannel].fPcbZone = fLookUpTableData[idManuChannel+1].fPcbZone;
329                       fPadData[idManuChannel].fPlane = fLookUpTableData[idManuChannel+1].fPlane;
330                       HLTDebug("\t Hit Found fo ich : %d, iPlane : %d, detelem %d, id : %d, at (%lf, %lf, %lf) cm"
331                               ,iChamber,fLookUpTableData[idManuChannel+1].fPlane,detElemId,fLookUpTableData[idManuChannel+1].fIdManuChannel,
332                               fPadData[idManuChannel].fRealX,
333                               fPadData[idManuChannel].fRealY,fPadData[idManuChannel].fRealZ);
334                               
335                       if(fMaxFiredPerDetElem[detElemId] == 0){
336                         DataIdIndex first;
337                         first.push_back(idManuChannel);
338                         fDetElemToDataId[detElemId] = first;
339                       }else{
340                         dataIndex =  fDetElemToDataId[detElemId];
341                         dataIndex.push_back(idManuChannel);
342                         fDetElemToDataId[detElemId] = dataIndex;
343                       }
344
345                       fMaxFiredPerDetElem[detElemId] = fMaxFiredPerDetElem[detElemId] + 1;
346
347                       dataCount ++;
348                       
349                     }//pattern maching is found 
350                   }// loop of ibitxy
351                 }// if pattern
352               }// iplane
353             }// ichamber
354             
355           }// if any non zero pattern found
356
357
358           index += 1 ; // skipping the last word though it is important
359           
360         }// if locId <=234
361       }// Dec Condn
362         
363       
364       if(!reg_phys_trig_occur)
365          index += 45;
366         
367       index += 1; // end of local Data 0xCAFEFADE
368
369       HLTDebug("iReg %d, iLoc %d, locId : %d, trigY %x, triggX %x, loDev %x, dec %x, sign %x,rawData : %x",
370                iReg,iLoc,locId,triggY,triggX,loDev,dec,sign, rawData[index]);
371
372       index = iLocIndex + 6 ; //important to reset the index counter for fake locids like 235 
373       
374      }// iLoc loop
375      
376   }// iReg Loop
377
378   return true;
379 }
380
381 bool AliHLTMUONTriggerReconstructor::FindTrigHits() 
382 {
383
384   map<int,DataIdIndex>::iterator it;
385
386   for(it = fDetElemToDataId.begin(); it != fDetElemToDataId.end(); it++){
387     HLTDebug("Nof data found in Detelem : %d = %d",it->first,(it->second).size());
388     if(!MergeTrigHits(it->second))
389       return false;
390   }// loop over detection element
391
392   DataIdIndex dataIndex;
393   for(it = fDetElemToDataId.begin(); it != fDetElemToDataId.end(); it++){
394     dataIndex = it->second;
395     for(size_t i=0;i<dataIndex.size();i++){
396       fPadData[dataIndex.at(i)].fDetElemId = 0;
397       fPadData[dataIndex.at(i)].fBuspatchId = 0;
398       fPadData[dataIndex.at(i)].fIdManuChannel = 0;
399       fPadData[dataIndex.at(i)].fIX = 0 ;
400       fPadData[dataIndex.at(i)].fIY = 0 ;
401       fPadData[dataIndex.at(i)].fRealX = 0.0 ;
402       fPadData[dataIndex.at(i)].fRealY = 0.0 ;
403       fPadData[dataIndex.at(i)].fRealZ = 0.0 ;
404       fPadData[dataIndex.at(i)].fPlane = -1 ;
405       fPadData[dataIndex.at(i)].fPcbZone = -1 ;
406       fPadData[dataIndex.at(i)].fCharge = 0 ;
407     }// data per detelem loop  
408   }//detelem loop
409   
410   return true;
411 }
412
413 bool AliHLTMUONTriggerReconstructor::MergeTrigHits(DataIdIndex& dataIndex)
414 {
415   int idManuChannelB, idManuChannelNB;
416   float halfPadLengthX,halfPadLengthY;
417   float diffX,diffY;
418
419   HLTDebug("\tThe bending plane hits are :");
420   for(size_t iPad=0;iPad<dataIndex.size();iPad++){
421     idManuChannelB   = dataIndex.at(iPad);
422     if(fPadData[idManuChannelB].fPlane == 0){
423       HLTDebug("\t detelem :%d, pcbzone : %d, (%f, %f, %f) cm",fPadData[idManuChannelB].fDetElemId,fPadData[idManuChannelB].fPcbZone,fPadData[idManuChannelB].fRealX,
424               fPadData[idManuChannelB].fRealY,fPadData[idManuChannelB].fRealZ);
425     }
426   }
427
428   HLTDebug("\tThe non-bending plane hits are :");
429   for(size_t jPad=0;jPad<dataIndex.size();jPad++){
430     idManuChannelNB   = dataIndex.at(jPad);
431     if(fPadData[idManuChannelNB].fPlane == 1){
432       HLTDebug("\t detelem :%d, pcbzone : %d,(%f, %f, %f) cm",fPadData[idManuChannelNB].fDetElemId,fPadData[idManuChannelNB].fPcbZone,fPadData[idManuChannelNB].fRealX,
433               fPadData[idManuChannelNB].fRealY,fPadData[idManuChannelNB].fRealZ);
434     }
435   }
436
437
438   for(size_t iPad=0;iPad<dataIndex.size();iPad++){
439     idManuChannelB   = dataIndex.at(iPad);
440     if(fPadData[idManuChannelB].fPlane == 0){
441       
442       halfPadLengthX = AliHLTMUONTriggerReconstructor::fgkHalfPadSizeXB[fPadData[idManuChannelB].fPcbZone] ;
443
444       for(size_t jPad=0;jPad<dataIndex.size();jPad++){
445         idManuChannelNB   = dataIndex.at(jPad);;
446         if(fPadData[idManuChannelNB].fPlane == 1){
447           
448           halfPadLengthY = AliHLTMUONTriggerReconstructor::fgkHalfPadSizeYNB[fPadData[idManuChannelNB].fPcbZone] ;
449           
450           if(fabsf(fPadData[idManuChannelNB].fRealX) > fabsf(fPadData[idManuChannelB].fRealX))
451             diffX = fabsf(fPadData[idManuChannelNB].fRealX) - fabsf(fPadData[idManuChannelB].fRealX);
452           else
453             diffX = fabsf(fPadData[idManuChannelB].fRealX) - fabsf(fPadData[idManuChannelNB].fRealX) ;
454
455           
456           if(fabsf(fPadData[idManuChannelNB].fRealY) > fabsf(fPadData[idManuChannelB].fRealY))
457             diffY = fabsf(fPadData[idManuChannelNB].fRealY) - fabsf(fPadData[idManuChannelB].fRealY);
458           else
459             diffY = fabsf(fPadData[idManuChannelB].fRealY) - fabsf(fPadData[idManuChannelNB].fRealY) ;
460           HLTDebug("\tdiffX %f,  halfPadLengthX %f,  diffY  %f, halfPadLengthY  %f\n",diffX,halfPadLengthX,diffY,halfPadLengthY);
461
462           if(diffX < halfPadLengthX + 1.0 && diffY < halfPadLengthY + 1.0 ){// added redundancy of 1.0 cm due to the pb of geometrical segmentation 
463
464             AliHLTMUONRecHitStruct hit;
465             hit.fX = fPadData[idManuChannelNB].fRealX;
466             hit.fY = fPadData[idManuChannelB].fRealY;
467             hit.fZ = fPadData[idManuChannelNB].fRealZ;
468
469             fRecPoints[(*fRecPointsCount)].fHit[0] = hit;
470             fRecPoints[(*fRecPointsCount)].fId = fPadData[idManuChannelB].fDetElemId ;
471             
472             (*fRecPointsCount)++;
473             if((*fRecPointsCount) == fMaxRecPointsCount){
474               HLTFatal("Nof RecHit (i.e. %d) exceeds the max nof RecHit limit %d\n",(*fRecPointsCount),fMaxRecPointsCount);
475               return false;
476             }
477
478             HLTDebug("\t\t\tdetelem : %d, x %f, y %f, z %f\n",fPadData[idManuChannelB].fDetElemId,fPadData[idManuChannelNB].fRealX,
479                    fPadData[idManuChannelB].fRealY,fPadData[idManuChannelB].fRealZ);
480           }
481           
482         }//condn for non-bending plane
483       }//for loop for non-bending plane
484
485     }// condn for bending plane
486   }// for loop for bending plane
487
488   return true;
489 }
490