Possibility to propagate tracks to the DCA to the primary vertex at the AOD level...
[u/mrichter/AliRoot.git] / TOF / TOFphysicsda.cxx
1 /*
2
3 TOF DA for online calibration
4
5 Version: 1.0
6 Contact: Roberto.Preghenella@bo.infn.it
7 Reference Run: 115401
8 Run Type: PHYSICS
9 DA Type: MON
10 Number of events needed:
11 Input Files: no input
12 Output Files: TOFdaHits.root TOFdaReadout.root
13 Event types used: PHYSICS_EVENT
14
15 */
16
17 #define FILE_HITS "TOFdaHits.root"
18 #define FILE_READOUT "TOFdaReadout.root"
19
20 #define EVENT_DETECTOR_TOF 5
21 #define READOUT_INFO_HISTO 1
22
23 // DATE
24 #include "event.h"
25 #include "monitor.h"
26 #include "daqDA.h"
27
28 #include <stdio.h>
29 #include <stdlib.h>
30 #include <errno.h>
31
32 //AliRoot
33 #include "TROOT.h"
34 #include "AliTOFRawStream.h"
35 #include "AliRawReaderDate.h"
36 #include "AliRawReader.h"
37 #include "AliDAQ.h"
38 #include "AliTOFHitData.h"
39 #include "AliTOFHitDataBuffer.h"
40 #include "AliTOFDaConfigHandler.h"
41 #include "AliTOFHitField.h"
42 #include "AliLog.h"
43 #include "AliTOFGeometry.h"
44 #include "AliTOFDecoderV2.h"
45 #include "AliTOFDecoderSummaryData.h"
46 #include "AliTOFDRMSummaryData.h"
47 #include "AliTOFTRMSummaryData.h"
48 #include "AliTOFChainSummaryData.h"
49 #include "AliTOFTDCHitBuffer.h"
50 #include "AliTOFTDCHit.h"
51 #include "AliTOFTDCErrorBuffer.h"
52 #include "AliTOFTDCError.h"
53
54 //ROOT
55 #include "TFile.h"
56 #include "TKey.h"
57 #include "TH2S.h"
58 #include "TObject.h"
59 #include "TMath.h"
60 #include "TSystem.h"
61 #include "TROOT.h"
62 #include "TPluginManager.h"
63 #include "TSAXParser.h"
64 #include "TTree.h"
65
66 /* Main routine
67       Arguments: 
68       1- monitoring data source
69 */
70 int 
71 main(int argc, char **argv) 
72 {
73   
74   /* magic line from Rene */
75   gROOT->GetPluginManager()->AddHandler("TVirtualStreamerInfo",
76                                         "*",
77                                         "TStreamerInfo",
78                                         "RIO",
79                                         "TStreamerInfo()");
80   
81
82   /* log start of process */
83   printf("TOF DA started\n");  
84   
85   /* check that we got some arguments = list of files */
86   if (argc!=2) {
87     printf("Wrong number of arguments\n");
88     return -1;
89   }
90
91   /*
92    * CONFIG
93    */
94   
95   /* retrieve config file */
96   int getConfigFile = daqDA_DB_getFile("TOFPhysicsConfig.xml","TOFPhysicsConfig.xml");
97   if (getConfigFile != 0){
98     printf("Failed to retrieve config file from DB! returning...\n");
99     return -1;
100   }
101   /* parse config file */
102   AliTOFDaConfigHandler* tofHandler = new AliTOFDaConfigHandler();
103   TSAXParser *parser = new TSAXParser();
104   parser->ConnectToHandler("AliTOFDaConfigHandler", tofHandler);  
105   if (parser->ParseFile("./TOFPhysicsConfig.xml") != 0) {
106     printf("Failed parsing config file! retunring... \n");
107     return -1;
108   }
109   /* setup config params */
110   Int_t meanMultiplicity = tofHandler->GetMeanMultiplicity(); /* average expected TOF multiplicity */
111   Int_t maxHits = tofHandler->GetMaxHits(); /* max number of hits to be collected */
112   printf("current settings:\n");
113   printf(" - meanMultiplicity = %d\n", meanMultiplicity);
114   printf(" - maxHits          = %d\n", maxHits);
115   /* constants */
116   const Int_t nChannels = 157248;
117   Int_t noiseCheckTrigger = 10; /* first noise check after 10 events */
118   Float_t meanChannelRate = (Float_t)meanMultiplicity / (Float_t)nChannels; /* average expected channel rate (hits/event) */
119   Float_t noiseThreshold = 10. * meanChannelRate; /* noise threshold (hits/event) */
120   Int_t minNoiseHits = 10; /* min number of channel hits to check noise */
121   /* counters and flags */
122   Int_t nPhysicsEvents, nCollectedPhysicsEvents, totHits;
123   Int_t nChHits[nChannels];
124   Bool_t inhibitCollection;
125   Bool_t noiseFlag[nChannels];
126   /* variables */
127   Int_t ddl, slot, trm, chain, tdc, channel, index, timebin, totbin, deltaBC, l0l1latency, det[5], dummy;
128   Float_t noiseHitThreshold;
129
130   /*
131    * INIT 
132    */
133
134   /* init counters and flags */
135   nPhysicsEvents = 0;
136   nCollectedPhysicsEvents = 0;
137   totHits = 0;
138   inhibitCollection = kFALSE;
139   for (Int_t ich = 0; ich < nChannels; ich++) {
140     nChHits[ich] = 0;
141     noiseFlag[ich] = kFALSE;
142   }
143
144   /* TOF raw data handling */
145   AliTOFRawStream *rawStream = new AliTOFRawStream();
146   AliTOFDecoderV2 *decoder = rawStream->GetDecoderV2();
147   AliTOFDecoderSummaryData *decodersd;
148   AliTOFDRMSummaryData *drmsd;
149   //  AliTOFLTMSummaryData *ltmsd;
150   AliTOFTRMSummaryData *trmsd;  
151   AliTOFChainSummaryData *chainsd;
152   AliTOFTDCHitBuffer *hitBuffer;
153   AliTOFTDCHit *hit;
154   AliTOFTDCErrorBuffer *errorBuffer;
155   AliTOFTDCError *error;
156   UShort_t errorFlags;
157   UChar_t *data = 0x0;
158   Int_t dataSize;
159   Int_t dataWords;
160   Int_t currentDDL;
161   const AliRawDataHeader *currentCDH;
162   Int_t currentMiniEventID;
163   Int_t currentEventID1;
164   Int_t currentL0BCID ;
165   Int_t currentBunchID;
166   Bool_t skipTRM, skipChain;
167   Int_t trmIndex, chainIndex, tdcIndex;
168   Double_t chainEff;
169   
170   /* open HITS output file */
171   TFile *fileOutHits = new TFile(FILE_HITS, "RECREATE"); 
172   /* create hit field data structure */
173   AliTOFHitField *hitField = new AliTOFHitField();
174   /* create temporary tree */
175   TTree *tempTree = new TTree("tempTree", "temporary tree");
176   tempTree->Branch("hit", "AliTOFHitField", &hitField);
177   /* create output tree */
178   TTree *outTree = new TTree("hitTree", "hit tree");
179   outTree->Branch("hit", "AliTOFHitField", &hitField);
180
181   /* open READOUT output file */
182   TFile *fileOutReadout = new TFile(FILE_READOUT, "RECREATE"); 
183   /* create chain readout efficiency histo */
184   TH1F *hChainEfficiency = new TH1F("hChainEfficiency", "Chain efficiency;chain;efficiency", 1440, 0., 1440.);
185
186 #if READOUT_INFO_HISTO
187   /* create TRM data histo */
188   TH1F *hTRMData = new TH1F("hTRMData", "TRM data;TRM;frequency", 720, 0., 720.);
189   /* create TRM empty event frequency histo */
190   TH1F *hTRMEmptyEvent = new TH1F("hTRMEmptyEvent", "TRM empty event error;TRM;frequency", 720, 0., 720.);
191   /* create TRM bad event counter frequency histo */
192   TH1F *hTRMBadEventCounter = new TH1F("hTRMBadEventCounter", "TRM bad event counter;TRM;frequency", 720, 0., 720.);
193   /* create TRM bad CRC frequency histo */
194   TH1F *hTRMBadCRC = new TH1F("hTRMBadCRC", "TRM bad CRC;TRM;frequency", 720, 0., 720.);
195
196   /* create chain data histo */
197   TH1F *hChainData = new TH1F("hChainData", "Chain data;chain;frequency", 1440, 0., 1440.);
198   /* create chain bad status frequency histo */
199   TH1F *hChainBadStatus = new TH1F("hChainBadStatus", "Chain bad status;chain;frequency", 1440, 0., 1440.);
200   /* create chain bad event counter frequency histo */
201   TH1F *hChainBadEventCounter = new TH1F("hChainBadEventCounter", "Chain bad event counter;chain;status;frequency", 1440, 0., 1440.);
202
203   /* create TDC error frequency histo */
204   TH1F *hTDCError = new TH1F("hTDCError", "TDC error;TDC;frequency", 21600, 0., 21600.);
205   /* create TDC error flags frequency histo */
206   TH2F *hTDCErrorFlags = new TH2F("hTDCErrorFlags", "TDC error flags;TDC;error flag;frequency", 21600, 0., 21600., 15, 0., 15);
207 #endif
208   
209
210   /*
211    * ONLINE MONITOR
212    */
213
214   AliLog::SetGlobalLogLevel(AliLog::kFatal);
215   struct eventHeaderStruct *event;
216   int ret;
217   /* define monitoring table */
218   char *monTable[5] = {
219     "ALL", "no",
220     "PHY", "yes",
221     NULL
222   };
223   ret = monitorDeclareTable(monTable);
224   if (ret != 0) {
225     printf("monitorDeclareTable() failed: %s\n", monitorDecodeError(ret));
226     return -1;
227   }
228   /* define data source : this is argument 1 */  
229   ret = monitorSetDataSource(argv[1]);
230   if (ret != 0) {
231     printf("monitorSetDataSource() failed : %s\n",monitorDecodeError(ret));
232     return -1;
233   }
234   /* declare monitoring program */
235   ret = monitorDeclareMp("TOFdaPhysics");
236   if (ret != 0) {
237     printf("monitorDeclareMp() failed : %s\n",monitorDecodeError(ret));
238     return -1;
239   }
240   /* define wait event timeout - 1s max */
241   monitorSetNowait();
242   monitorSetNoWaitNetworkTimeout(1000);
243
244   /* variables */
245   
246   /* loop over events */
247   while (1) {
248     
249     /* check shutdown condition */
250     if (daqDA_checkShutdown()) break;
251     
252     /*
253      * NOISE CHECK
254      */
255
256     /* check inhibit collection */
257     if (!inhibitCollection) {
258       /* check number of events and check noise */
259       if (nCollectedPhysicsEvents >= noiseCheckTrigger || totHits >= maxHits) {
260         noiseHitThreshold = noiseThreshold * nCollectedPhysicsEvents;
261         printf("noise check triggered after %d events: threshold is %f hits\n", nCollectedPhysicsEvents, noiseHitThreshold);
262         /* loop over all channels */
263         for (Int_t ich = 0; ich < nChannels; ich++) {
264           /* check */
265           if (nChHits[ich] < minNoiseHits || noiseFlag[ich] || nChHits[ich] < noiseHitThreshold) continue;
266           printf("channel %d tagged as noisy (%d hits): disabled\n", ich, nChHits[ich]);
267           noiseFlag[ich] = kTRUE;
268           totHits -= nChHits[ich];
269         } /* end of loop over all channels */
270         /* set new noise check trigger value */
271         noiseCheckTrigger *= 10;
272       } /* end of noise check */    
273       
274       /* inhibit hit collection when maximum number of hits exceeded */
275       if (totHits >= maxHits) {
276         printf("maximum number of hits exceeded (%d): inhibit hit collection\n", maxHits);
277         inhibitCollection = kTRUE;
278       }
279     }
280     
281     /*
282      * GET EVENT
283      */
284
285     /* get next event (blocking call until timeout) */
286     ret = monitorGetEventDynamic((void **)&event);
287     if (ret == MON_ERR_EOF) {
288       printf ("End of File detected\n");
289       break; /* end of monitoring file has been reached */
290     }
291     if (ret != 0) {
292       printf("monitorGetEventDynamic() failed (ret=%d errno=%d): %s\n", ret, errno, monitorDecodeError(ret));
293       break;
294     }
295     /* retry if got no event */
296     if (event==NULL) continue;
297     /* check TOF in partecipating detectors */
298     if (!TEST_DETECTOR_IN_PATTERN(event->eventDetectorPattern, EVENT_DETECTOR_TOF)) {
299       free(event);
300       continue;
301     }
302     /* check event type */
303     if (event->eventType != PHYSICS_EVENT) {
304       printf("not a physics event: %d\n", event->eventType);
305       free(event);
306       continue;
307     }
308     /* increment number of physics events */
309     nPhysicsEvents++;
310     if (!inhibitCollection) nCollectedPhysicsEvents++;
311
312     /*
313      * DECODE EVENT
314      */
315
316     /* create and setup raw reader */
317     AliRawReader *rawReader = new AliRawReaderDate((void *)event);
318     rawReader->Reset();
319     rawReader->Select("TOF", 0, AliDAQ::NumberOfDdls("TOF") - 1);
320     /* setup raw stream */
321     rawStream->SetRawReader(rawReader);
322
323     /* loop over DDLs - rawReader->ReadHeader() */
324     while (rawReader->ReadHeader()) {
325
326       /* read equipment data */
327       dataSize = rawReader->GetDataSize();
328       data = new UChar_t[dataSize];
329       if (!rawReader->ReadNext(data, dataSize)){
330         delete [] data;
331         continue;
332       }
333
334       /* decode data */
335       dataWords = dataSize / 4;
336       decoder->Decode((UInt_t *)data, dataWords);
337       delete [] data;
338
339       /* read equipment info */
340       currentDDL = rawReader->GetDDLID();
341       currentCDH = rawReader->GetDataHeader();
342       currentMiniEventID = currentCDH->GetMiniEventID();
343       currentEventID1 = currentCDH->GetEventID1();
344
345       /* read decoder summary data */
346       decodersd = decoder->GetDecoderSummaryData();
347       
348       /* check DRM header/trailer */
349       drmsd = decodersd->GetDRMSummaryData();
350       if (!drmsd->GetHeader() || !drmsd->GetTrailer()) continue;
351
352       /* get DRM data */
353       currentL0BCID = drmsd->GetL0BCID();
354       
355       /* loop over TRM to get hits */
356       for (Int_t itrm = 0; itrm < 10; itrm++) {
357         trmsd = drmsd->GetTRMSummaryData(itrm);
358         trmIndex = itrm + 10 * currentDDL;
359         skipTRM = kFALSE;
360         
361         /* check header/trailer */
362         if (!trmsd->GetHeader() || !trmsd->GetTrailer()) continue;
363         
364 #if READOUT_INFO_HISTO
365         /* fill TRM data */
366         hTRMData->Fill(trmIndex);
367         /* fill TRM empty event */
368         if (trmsd->GetEBit() != 0) {
369           hTRMEmptyEvent->Fill(trmIndex);
370           skipTRM = kTRUE;
371         }
372         /* fill TRM bad event counter */
373         if (trmsd->GetEventCounter() != drmsd->GetLocalEventCounter()) {
374           hTRMBadEventCounter->Fill(trmIndex);
375           skipTRM = kTRUE;
376         }
377         /* fill TRM bad CRC */
378         if (trmsd->GetEventCRC() != trmsd->GetDecoderCRC()) {
379           hTRMBadCRC->Fill(trmIndex);
380           skipTRM = kTRUE;
381         }
382 #else
383         /* check bad condition and skip TRM */
384         if ( trmsd->GetEBit() != 0 ||
385              trmsd->GetEventCounter() != drmsd->GetLocalEventCounter() ||
386              trmsd->GetEventCRC() != trmsd->GetDecoderCRC() ) continue;
387 #endif
388         
389         /* loop over chains */
390         for (Int_t ichain = 0; ichain < 2; ichain++) {
391           chainsd = trmsd->GetChainSummaryData(ichain);
392           chainIndex = ichain + 2 * itrm + 20 * currentDDL;
393           skipChain = kFALSE;
394           
395           /* check header/trailer */
396           if (!chainsd->GetHeader() || !chainsd->GetTrailer()) continue;
397           currentBunchID = chainsd->GetBunchID();
398           hitBuffer = chainsd->GetTDCPackedHitBuffer();
399           errorBuffer = chainsd->GetTDCErrorBuffer();
400
401 #if READOUT_INFO_HISTO
402           /* fill chain data */
403           hChainData->Fill(chainIndex);
404           /* check chain bad status */
405           if (chainsd->GetStatus() != 0) {
406             hChainBadStatus->Fill(chainIndex);
407             skipChain = kTRUE;
408           }
409           /* check chain bad event counter */
410           if (chainsd->GetEventCounter() != drmsd->GetLocalEventCounter()) {    
411             hChainBadEventCounter->Fill(chainIndex);
412             skipChain = kTRUE;
413           }
414           /* fill TDC error frequency histo */
415           for (Int_t ierr = 0; ierr < errorBuffer->GetEntries(); ierr++) {
416             error = errorBuffer->GetError(ierr);
417             tdc = error->GetTDCID();
418             tdcIndex = tdc + 15 * ichain + 30 * itrm + 300 * currentDDL;
419             hTDCError->Fill(tdcIndex); 
420             errorFlags = error->GetErrorFlags();
421             for (Int_t ierflg = 0; ierflg < 15; ierflg++)
422               if (errorFlags & (1 << ierflg))
423                 hTDCErrorFlags->Fill(tdcIndex, ierflg); 
424           }
425 #else
426           /* check bad condition and skip chain */
427           if ( chainsd->GetStatus() != 0 ||
428                chainsd->GetEventCounter() != drmsd->GetLocalEventCounter() ) continue;
429 #endif
430           
431           /*
432            * CHAIN READOUT EFFICIENCY
433            */
434           
435           /* compute number of available channels removing TDCs in error */
436           chainEff = (120. - 8. * errorBuffer->GetEntries()) / 120.;
437           /* fill chain readout efficiency histo */
438           if (!skipTRM && !skipChain)
439             hChainEfficiency->Fill(chainIndex, chainEff);
440           
441           /*
442            * HIT MANIPULATION
443            */
444           
445           /* check inhibit collection */
446           if (inhibitCollection) continue;
447
448           /* loop over hits in buffer */
449           for (Int_t ihit = 0; ihit < hitBuffer->GetEntries(); ihit++) {
450             
451             /* get hit */
452             hit = hitBuffer->GetHit(ihit);
453             
454             /* get channel info */
455             ddl = currentDDL;
456             slot = trmsd->GetSlotID();
457             trm = slot - 3;
458             chain = chainsd->GetChain();
459             tdc = hit->GetTDCID();
460             channel = hit->GetChan();
461             /* get index */
462             rawStream->EquipmentId2VolumeId(ddl, slot, chain, tdc, channel, det);
463             dummy = det[4];
464             det[4] = det[3];
465             det[3] = dummy;
466             /* check valid index */
467             if (det[0] < 0 || det[0] > 17 ||
468                 det[1] < 0 || det[1] > 5 ||
469                 det[2] < 0 || det[2] > 18 ||
470                 det[3] < 0 || det[3] > 1 ||
471                 det[4] < 0 || det[4] > 47) continue;
472             index = AliTOFGeometry::GetIndex(det);
473             
474             /* check noise flag */
475             if (noiseFlag[index]) continue;
476             /* increment number of channel hits and total hits */
477             nChHits[index]++;
478             totHits++;
479             /* get signal info */
480             timebin = hit->GetHitTime();
481             totbin = hit->GetTOTWidth();
482             deltaBC = currentBunchID - currentEventID1;
483             l0l1latency = currentMiniEventID - currentL0BCID;
484             /* set hit field data */
485             hitField->SetIndex(index);
486             hitField->SetTimeBin(timebin);
487             hitField->SetTOTBin(totbin);
488             hitField->SetDeltaBC(deltaBC);
489             hitField->SetL0L1Latency(l0l1latency);
490             /* fill temp tree */
491             tempTree->Fill();
492
493           } /* end of loop over hits in buffer */
494         } /* end of loop over chains */
495       } /* end of loop over TRMs */
496     } /* end of loop over DDLs - rawReader->ReadHeader() */
497     
498     /* delete raw reader */
499     delete rawReader;
500     /* free event */
501     free(event);
502     
503   } /* end of loop over events */
504   
505   /* final noise check */
506   noiseHitThreshold = noiseThreshold * nCollectedPhysicsEvents;
507   printf("final noise check after collectiong %d events: threshold is %f hits\n", nCollectedPhysicsEvents, noiseHitThreshold);
508   /* loop over all channels */
509   for (Int_t ich = 0; ich < nChannels; ich++) {
510     /* check */
511     if (nChHits[ich] < minNoiseHits || noiseFlag[ich] || nChHits[ich] < noiseHitThreshold) continue;
512     printf("channel %d tagged as noisy (%d hits): disabled\n", ich, nChHits[ich]);
513     noiseFlag[ich] = kTRUE;
514     totHits -= nChHits[ich];
515   } /* end of loop over all channels */
516   
517   /* copy hits into output tree from temp tree */
518   printf("copy hits from temporary tree into output tree\n");
519   printf("temporary tree contains %d hits\n", (Int_t)tempTree->GetEntries());
520   for (Int_t ihit = 0; ihit < tempTree->GetEntries(); ihit++) {
521     /* get entry */
522     tempTree->GetEntry(ihit);
523     /* check noise flag */
524     if (noiseFlag[hitField->GetIndex()]) continue;
525     /* fill output tree */
526     outTree->Fill();
527   } /* end of copy hits into output tree from temp tree */
528   printf("output tree contains %d hits\n", (Int_t)outTree->GetEntries());
529
530   /* write output tree on HITS file */
531   fileOutHits->cd();
532   outTree->Write();
533   fileOutHits->Close();
534   /* export file to FXS */
535   if (daqDA_FES_storeFile(FILE_HITS, "HITS"))
536     return -2;
537
538 #if READOUT_INFO_HISTO
539   hTRMData->Sumw2();
540   hTRMEmptyEvent->Sumw2();
541   hTRMBadEventCounter->Sumw2();
542   hTRMBadCRC->Sumw2();
543
544   hChainData->Sumw2();
545   hChainBadStatus->Sumw2();
546   hChainBadEventCounter->Sumw2();
547
548   hTDCError->Sumw2();
549   hTDCErrorFlags->Sumw2();
550
551   /* divide histos */
552   if (nPhysicsEvents > 0) {
553     hTRMEmptyEvent->Divide(hTRMData);
554     hTRMBadEventCounter->Divide(hTRMData);
555     hTRMBadCRC->Divide(hTRMData);
556     hTRMData->Scale(1. / nPhysicsEvents);
557
558     hChainBadStatus->Divide(hChainData);
559     hChainBadEventCounter->Divide(hChainData);
560     hChainData->Scale(1. / nPhysicsEvents);
561
562     hTDCError->Scale(1. / nPhysicsEvents);
563     hTDCErrorFlags->Scale(1. / nPhysicsEvents);
564   }
565 #endif
566   
567   /* scale chain efficiency by number of physics events */
568   printf("found %d physics events\n", nPhysicsEvents);
569   hChainEfficiency->Sumw2();
570   if (nPhysicsEvents > 0) hChainEfficiency->Scale(1. / (nPhysicsEvents));
571   
572   /* write efficiency histo on READOUT file */
573   fileOutReadout->cd();
574 #if READOUT_INFO_HISTO
575   hTRMData->Write();
576   hTRMEmptyEvent->Write();
577   hTRMBadEventCounter->Write();
578   hTRMBadCRC->Write();
579   
580   hChainData->Write();
581   hChainBadStatus->Write();
582   hChainBadEventCounter->Write();
583   
584   hTDCError->Write();
585   hTDCErrorFlags->Write();
586 #endif
587   hChainEfficiency->Write();
588   fileOutReadout->Close();
589   /* export file to FXS */
590   if (daqDA_FES_storeFile(FILE_READOUT, "READOUT"))
591     return -2;
592
593   return 0;
594 }