Fixing code so that components use the new standard DDL_RAW data type.
[u/mrichter/AliRoot.git] / HLT / MUON / OnlineAnalysis / AliHLTMUONTrackerCalibratorComponent.cxx
CommitLineData
8134dd2e 1/**************************************************************************
2 * This file is property of and copyright by the ALICE HLT Project *
3 * All rights reserved. *
4 * *
5 * Primary Authors: *
6 * Artur Szostak <artursz@iafrica.com> *
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/// @file AliHLTMUONTrackerCalibratorComponent.cxx
21/// @author Artur Szostak <artursz@iafrica.com>
22/// @date
23/// @brief Implementation of the AliHLTMUONTrackerCalibratorComponent class.
24///
25
26#include "AliHLTMUONTrackerCalibratorComponent.h"
27#include "AliHLTMUONConstants.h"
28#include "AliHLTMUONUtils.h"
29#include <cassert>
30
31ClassImp(AliHLTMUONTrackerCalibratorComponent);
32
33
34///////////////////////////////////////////////////////////////////////////////
35// The code from here on was copied from MUONTRKda.cxx and addapted to the
36// HLT framework.
37//TODO: test that any of this actually works and clean up the code.
38///////////////////////////////////////////////////////////////////////////////
39
40#include <iostream>
41#include <fstream>
42using namespace std;
43
44#include <cstdio>
45#include <cstdlib>
46
47//AliRoot
48#include "AliRawReaderMemory.h"
49#include "AliMUONRawStreamTracker.h"
50#include "AliMUONDspHeader.h"
51#include "AliMUONBlockHeader.h"
52#include "AliMUONBusStruct.h"
53#include "AliMUONDDLTracker.h"
54#include "AliMUONVStore.h"
55#include "AliMUON2DMap.h"
56#include "AliMUONCalibParamND.h"
57#include "AliMpIntPair.h"
58#include "AliMpConstants.h"
59#include "AliRawReaderDate.h"
60
61//ROOT
62#include "TFile.h"
63#include "TTree.h"
64#include "TH1F.h"
65#include "TString.h"
66#include "TStopwatch.h"
67#include "TMath.h"
68#include "TTimeStamp.h"
69#include "TGraphErrors.h"
70#include "TF1.h"
71#include "TROOT.h"
72#include "TPluginManager.h"
73#include "TFitter.h"
74
75#define NFITPARAMS 4
76
77namespace
78{
79
80// global variables
81const Int_t gkNChannels = AliMpConstants::ManuNofChannels();
82const Int_t gkADCMax = 4095;
83
84AliMUONVStore* gPedestalStore = NULL;
85
86Int_t gNManu = 0;
87Int_t gNChannel = 0;
88UInt_t gRunNumber = 0;
89Int_t gNEvents = 0;
90Int_t gNDateEvents = 0;
91Int_t gPrintLevel = 1; // global printout variable
92Int_t gPlotLevel = 0; // global plot variable
93
94TH1F* gPedMeanHisto = 0x0;
95TH1F* gPedSigmaHisto = 0x0;
96Char_t gHistoFileName[256];
97
98// used by makegain
99TString gHistoFileName_gain = "MuonTrkDA_data.root";
100TString gRootFileName = "MuonTrkDA_gain.root";
101Char_t filenam[256] = "MuonTrkDA_gain.param"; // if gPrintLevel = 2
102
103TString gCommand("ped");
104
105}; // end of namespace
106
107
108//________________
109Double_t funcLin(Double_t *x, Double_t *par)
110{
111 return par[0] + par[1]*x[0];
112}
113
114//________________
115Double_t funcParabolic(Double_t *x, Double_t *par)
116{
117 return par[0]*x[0]*x[0];
118}
119
120//________________
121Double_t funcCalib(Double_t *x, Double_t *par)
122{
123 Double_t xLim= par[3];
124
125 if(x[0] <= xLim) return par[0] + par[1]*x[0];
126
127 Double_t yLim = par[0]+ par[1]*xLim;
128 return yLim + par[1]*(x[0] - xLim) + par[2]*(x[0] - xLim)*(x[0] - xLim);
129}
130
131//__________
132void MakePed(Int_t busPatchId, Int_t manuId, Int_t channelId, Int_t charge)
133{
134
135 AliMUONVCalibParam* ped =
136 static_cast<AliMUONVCalibParam*>(gPedestalStore->FindObject(busPatchId, manuId));
137
138 if (!ped) {
139 gNManu++;
140 ped = new AliMUONCalibParamND(2, gkNChannels,busPatchId, manuId, -1.); // put default wise -1, not connected channel
141 gPedestalStore->Add(ped);
142 }
143
144 if (gNEvents == 1) {
145 ped->SetValueAsDouble(channelId, 0, 0.);
146 ped->SetValueAsDouble(channelId, 1, 0.);
147 }
148
149 Double_t pedMean = ped->ValueAsDouble(channelId, 0) + charge;
150 Double_t pedSigma = ped->ValueAsDouble(channelId, 1) + charge*charge;
151
152 ped->SetValueAsDouble(channelId, 0, pedMean);
153 ped->SetValueAsDouble(channelId, 1, pedSigma);
154
155}
156
157//________________
158void MakePedStore(TString flatOutputFile = "")
159{
160 TTimeStamp date;
161 Double_t pedMean;
162 Double_t pedSigma;
163 ofstream fileout;
164 Int_t busPatchId;
165 Int_t manuId;
166 Int_t channelId;
167
168 // histo
169// sprintf(gHistoFileName,"mutrkped-%d.root",gRunNumber);
170 sprintf(gHistoFileName,"MuonTrkDA_%d_ped.root",gRunNumber);
171 TFile* histoFile = new TFile(gHistoFileName,"RECREATE","MUON Tracking pedestals");
172
173 Char_t name[255];
174 Char_t title[255];
175 sprintf(name,"pedmean_allch");
176 sprintf(title,"Pedestal mean all channels");
177 Int_t nx = 1000;
178 Int_t xmin = 0;
179 Int_t xmax = 1000;
180 gPedMeanHisto = new TH1F(name,title,nx,xmin,xmax);
181 gPedMeanHisto->SetDirectory(histoFile);
182
183 sprintf(name,"pedsigma_allch");
184 sprintf(title,"Pedestal sigma all channels");
185 nx = 200;
186 xmin = 0;
187 xmax = 50;
188 gPedSigmaHisto = new TH1F(name,title,nx,xmin,xmax);
189 gPedSigmaHisto->SetDirectory(histoFile);
190
191 TTree* tree = new TTree("t","Pedestal tree");
192 tree->Branch("bp",&busPatchId,"bp/I");
193 tree->Branch("manu",&manuId,",manu/I");
194 tree->Branch("channel",&channelId,",channel/I");
195 tree->Branch("pedMean",&pedMean,",pedMean/D");
196 tree->Branch("pedSigma",&pedSigma,",pedSigma/D");
197
198 if (!flatOutputFile.IsNull()) {
199 fileout.open(flatOutputFile.Data());
200 fileout<<"//===========================================================================" << endl;
201 fileout<<"// Pedestal file calculated by MUONTRKda"<<endl;
202 fileout<<"//===========================================================================" << endl;
203 fileout<<"// * Run : " << gRunNumber << endl;
204 fileout<<"// * Date : " << date.AsString("l") <<endl;
205 fileout<<"// * Statictics : " << gNEvents << endl;
206 fileout<<"// * # of MANUS : " << gNManu << endl;
207 fileout<<"// * # of channels : " << gNChannel << endl;
208 fileout<<"//"<<endl;
209 fileout<<"//---------------------------------------------------------------------------" << endl;
210 fileout<<"//---------------------------------------------------------------------------" << endl;
211// fileout<<"//format : BUS_PATCH MANU_ID CHANNEL MEAN SIGMA"<<endl;
212 fileout<<"// BP MANU CH. MEAN SIGMA"<<endl;
213 fileout<<"//---------------------------------------------------------------------------" << endl;
214
215 }
216
217 // iterator over pedestal
218 TIter next(gPedestalStore->CreateIterator());
219 AliMUONVCalibParam* ped;
220
221 while ( ( ped = dynamic_cast<AliMUONVCalibParam*>(next() ) ) )
222 {
223 busPatchId = ped->ID0();
224 manuId = ped->ID1();
225
226 for (channelId = 0; channelId < ped->Size() ; ++channelId) {
227
228 pedMean = ped->ValueAsDouble(channelId, 0);
229
230 if (pedMean > 0) { // connected channels
231
232 ped->SetValueAsDouble(channelId, 0, pedMean/(Double_t)gNEvents);
233
234 pedMean = ped->ValueAsDouble(channelId, 0);
235 pedSigma = ped->ValueAsDouble(channelId, 1);
236
237 ped->SetValueAsDouble(channelId, 1, TMath::Sqrt(TMath::Abs(pedSigma/(Double_t)gNEvents - pedMean*pedMean)));
238
239 pedMean = ped->ValueAsDouble(channelId, 0) + 0.5 ;
240// pedMean = ped->ValueAsDouble(channelId, 0) ;
241 pedSigma = ped->ValueAsDouble(channelId, 1);
242
243
244 if (!flatOutputFile.IsNull()) {
245 fileout << "\t" << busPatchId << "\t" << manuId <<"\t"<< channelId << "\t"
246 << pedMean <<"\t"<< pedSigma << endl;
247 }
248
249 gPedMeanHisto->Fill(pedMean);
250 gPedSigmaHisto->Fill(pedSigma);
251
252 tree->Fill();
253 }
254 }
255 }
256
257 // file outputs
258 if (!flatOutputFile.IsNull())
259 fileout.close();
260
261 histoFile->Write();
262 histoFile->Close();
263
264// delete tree;
265
266
267// not usable need an update of DATE ->5.28, but it compiles
268// uncomment when DATE version ok.
269// setenv DATE_FES_PATH
270// setenv DATE_RUN_NUMBER
271// setenv DATE_ROLE_NAME
272// setenv DATE_DETECTOR_CODE
273
274// if (!flatOutputFile.IsNull()) {
275
276// flatOutputFile.Prepend("./");
277
278// status = daqDA_FES_storeFile(flatOutputFile.Data(),"MUONTRKda_results");
279// if (status) {
280// printf("Failed to export file : %d\n",status);
281// return -1;
282// }
283// }
284
285}
286
287//________________
288void MakePedStoreForGain(Int_t injCharge)
289{
290 // store pedestal map in root file
291
292// Int_t injCharge = 200;
293
294 TTree* tree = 0x0;
295
296 // compute and store pedestals
297 MakePedStore();
298
299 // store in root file
300// sprintf(gHistoFileName,"mutrkgain.root");
301
302 TString mode("UPDATE");
303
304 if (gCommand.Contains("cre")) {
305 mode = "RECREATE";
306 }
307 TFile* histoFile = new TFile(gHistoFileName_gain, mode.Data(), "MUON Tracking Gains");
308
309 // second argument should be the injected charge, taken from config crocus file
310 // put also info about run number could be usefull
311 AliMpIntPair* pair = new AliMpIntPair(gRunNumber, injCharge);
312
313 if (mode.CompareTo("UPDATE") == 0) {
314 tree = (TTree*)histoFile->Get("t");
315 tree->SetBranchAddress("run",&pair);
316 tree->SetBranchAddress("ped",&gPedestalStore);
317
318 } else {
319 tree = new TTree("t","Pedestal tree");
320 tree->Branch("run", "AliMpIntPair",&pair);
321 tree->Branch("ped", "AliMUON2DMap",&gPedestalStore);
322 tree->SetBranchAddress("run",&pair);
323 tree->SetBranchAddress("ped",&gPedestalStore);
324
325 }
326
327 tree->Fill();
328 tree->Write("t", TObject::kOverwrite); // overwrite the tree
329 histoFile->Close();
330
331 delete pair;
332}
333
334//________________
335void MakeGainStore(TString flatOutputFile)
336{
337 Double_t goodA1Min = 0.5;
338 Double_t goodA1Max = 2.;
339 Double_t goodA2Min = -0.5E-03;
340 Double_t goodA2Max = 1.E-03;
341
342 // open file mutrkgain.root
343 // read again the pedestal for the calibration runs (9 runs ?)
344 // need the injection charge from config file (to be done)
345 // For each channel make a TGraphErrors (mean, sigma) vs injected charge
346 // Fit with a polynomial fct
347 // store the result in a flat file.
348
349 Double_t pedMean[11];
350 Double_t pedSigma[11];
351 Double_t injCharge[11];
352 Double_t injChargeErr[11];
353
354 ofstream fileout;
355 TTimeStamp date;
356
357// full print out
358
359 // why 2 files ? (Ch. F.)
360 FILE *pfilen = 0;
361 if(gPrintLevel==2)
362 {
363// Char_t filenam[256]="MuonTrkDA_gain.param";
364 cout << " fit parameter file = " << filenam << "\n";
365 pfilen = fopen (filenam,"w");
366
367 fprintf(pfilen,"//===================================================================\n");
368 fprintf(pfilen,"// BP MANU CH. a0 a1 a2 xlim P(chi2) P(chi2)2 Q\n");
369 fprintf(pfilen,"//===================================================================\n");
370 }
371
372 FILE *pfilew=0;
373 if (!flatOutputFile.IsNull())
374 {
375 pfilew = fopen (flatOutputFile.Data(),"w");
376
377 fprintf(pfilew,"//=================================================\n");
378 fprintf(pfilew,"// Calibration file calculated by MUONTRKda \n");
379 fprintf(pfilew,"//=================================================\n");
380 fprintf(pfilew,"// * Run : %d \n",gRunNumber);
381 fprintf(pfilew,"// * Date : %s \n",date.AsString("l"));
382 fprintf(pfilew,"// * Statictics : %d \n",gNEvents);
383 fprintf(pfilew,"// * # of MANUS : %d \n",gNManu);
384 fprintf(pfilew,"// * # of channels : %d \n",gNChannel);
385 fprintf(pfilew,"//-------------------------------------------------\n");
386 fprintf(pfilew,"//=======================================\n");
387 fprintf(pfilew,"// BP MANU CH. a1 a2 thres. Q\n");
388 fprintf(pfilew,"//=======================================\n");
389 }
390
391
392
393// sprintf(gHistoFileName,"mutrkgain.root");
394 TFile* histoFile = new TFile(gHistoFileName_gain);
395
396 AliMUON2DMap* map[11];
397 AliMUONVCalibParam* ped[11];
398 AliMpIntPair* run[11];
399
400// plot out
401
402 TFile* gainFile = 0x0;
403// sprintf(gRootFileName,"makegain.root");
404 gainFile = new TFile(gRootFileName,"RECREATE");
405
406 Double_t chi2 = 0.;
407 Double_t chi2P2 = 0.;
408 Double_t prChi2 = 0;
409 Double_t prChi2P2 =0;
410 Double_t a0,a1,a2;
411 Int_t busPatchId ;
412 Int_t manuId ;
413 Int_t channelId ;
414 Int_t threshold = 0;
415 Int_t Q = 0;
416
417 TTree *tg = new TTree("tg","TTree avec class Manu_DiMu");
418
419 tg->Branch("bp",&busPatchId, "busPatchId/I");
420 tg->Branch("manu",&manuId, "manuId/I");
421 tg->Branch("channel",&channelId, "channelId/I");
422
423 tg->Branch("a0",&a0, "a0/D");
424 tg->Branch("a1",&a1, "a1/D");
425 tg->Branch("a2",&a2, "a2/D");
426 tg->Branch("Pchi2",&prChi2, "prChi2/D");
427 tg->Branch("Pchi2_2",&prChi2P2, "prChi2P2/D");
428 tg->Branch("Threshold",&threshold, "threshold/I");
429 tg->Branch("Q",&Q, "Q/I");
430
431 //read back from root file
432 TTree* tree = (TTree*)histoFile->Get("t");
433 Int_t nEntries = tree->GetEntries();
434// tree->Branch("a1",&a1, "a1/D");
435
436 if(gPrintLevel) cout << " nEntries = " << nEntries << " DAC values \n" << endl;
437
438 // read back info
439 for (Int_t i = 0; i < nEntries; ++i) {
440 map[i] = 0x0;
441 run[i] = 0x0;
442 tree->SetBranchAddress("ped",&map[i]);
443 tree->SetBranchAddress("run",&run[i]);
444 tree->GetEvent(i);
445
446// std::cout << map[i] << " " << run[i] << std::endl;
447 }
448
449 // Q = f(ADC)
450 TF1 *f1 = new TF1("f1",funcLin,0.,gkADCMax,2);
451 TF1 *f2 = new TF1("f2",funcParabolic,0.,gkADCMax,1);
452
453 char graphName[256];
454
455 // iterates over the first pedestal run
456 TIter next(map[0]->CreateIterator());
457 AliMUONVCalibParam* p;
458
459 Int_t nmanu = 0;
460 Int_t nBadChannel = 0;
461 Double_t sumProbChi2 = 0.;
462 Double_t sumA1 = 0.;
463 Double_t sumProbChi2P2 = 0.;
464 Double_t sumA2 = 0.;
465
466 Double_t x[11], xErr[11], y[11], yErr[11];
467
468 while ( ( p = dynamic_cast<AliMUONVCalibParam*>(next() ) ) )
469 {
470 ped[0] = p;
471
472// Int_t busPatchId = p->ID0();
473// Int_t manuId = p->ID1();
474 busPatchId = p->ID0();
475 manuId = p->ID1();
476
477 // read back pedestal from the other runs for the given (bupatch, manu)
478 for (Int_t i = 1; i < nEntries; ++i) {
479 ped[i] = static_cast<AliMUONVCalibParam*>(map[i]->FindObject(busPatchId, manuId));
480 }
481
482 // compute for each channel the gain parameters
483// for ( Int_t channelId = 0; channelId < ped[0]->Size() ; ++channelId ) {
484 for ( channelId = 0; channelId < ped[0]->Size() ; ++channelId ) {
485
486 Int_t n = 0;
487 for (Int_t i = 0; i < nEntries; ++i) {
488
489 if (!ped[i]) continue; //shouldn't happen.
490 pedMean[i] = ped[i]->ValueAsDouble(channelId, 0);
491 pedSigma[i] = ped[i]->ValueAsDouble(channelId, 1);
492 injCharge[i] = (Double_t)run[i]->GetSecond();
493 injChargeErr[i] = 0.01*injCharge[i];
494 if(injChargeErr[i] <= 1.) injChargeErr[i]=1.;
495
496// if(n<2)cout << nEntries << " " << i << " " << injCharge[i] << endl;
497
498// cout << busPatchId << "\t" << manuId <<"\t"<< channelId << "\t" << n << " " << pedMean[i] << " " << pedSigma[i] << " " << injCharge[i] << " " << injChargeErr[i] << endl;
499
500 if (pedMean[i] < 0) continue; // not connected
501
502 if (pedSigma[i] <= 0) pedSigma[i] = 1.; // should not happen.
503 n++;
504 }
505
506 // makegain (JLC)
507
508
509 // Fit Method: Linear fit over 6 points + fit parabolic function over 3 points)
510
511 // 1. - linear fit over 6 points
512
513 Double_t par[4] = {0.,0.,0.,gkADCMax};
514
515 Int_t nInit = 1;
516 Int_t nbs = nEntries - nInit;
517 Int_t nbpf1 = 6; // linear fit over nbf1 points
518
519 for (Int_t j = 0; j < nbs; ++j)
520 {
521 Int_t k = j + nInit;
522 x[j] = pedMean[k];
523 xErr[j] = pedSigma[k];
524 y[j] = injCharge[k];
525 yErr[j] = injChargeErr[k];
526
527 }
528
529 TGraph *graphErr = new TGraphErrors(nbpf1, x, y, xErr, yErr);
530
531 f1->SetParameters(0,0);
532
533 graphErr->Fit("f1","RQ");
534
535 chi2 = f1->GetChisquare();
536 f1->GetParameters(par);
537
538 prChi2 = TMath::Prob(chi2, nbpf1 - 2);
539
540 Double_t xLim = pedMean[nInit + nbpf1 - 1];
541 Double_t yLim = par[0]+par[1] * xLim;
542
543 a0 = par[0];
544 a1 = par[1];
545
546 delete graphErr;
547
548
549 // 2. - Translation : new origin (xLim, yLim) + parabolic fit over nbf2 points
550
551 Int_t nbpf2 = nEntries - (nInit + nbpf1) + 1;
552
553 if(nbpf2 > 1)
554 {
555 for (Int_t j = 0; j < nbpf2; ++j)
556 {
557 Int_t k = j + (nInit + nbpf1) - 1;
558 x[j] = pedMean[k] - xLim;
559 xErr[j] = pedSigma[k];
560
561 y[j] = injCharge[k] - yLim - par[1]*x[j];
562 yErr[j] = injChargeErr[k];
563
564 }
565
566 TGraph *graphErr = new TGraphErrors(nbpf2, x, y, xErr, yErr);
567
568 graphErr->Fit(f2,"RQ");
569 chi2P2 = f2->GetChisquare();
570 f2->GetParameters(par);
571
572 prChi2P2 = TMath::Prob(chi2P2, nbpf2-1);
573
574 a2 = par[0];
575
576 par[0] = a0;
577 par[1] = a1;
578 par[2] = a2;
579 par[3] = xLim;
580
581 delete graphErr;
582
583 }
584
585 // Prints
586
587 Int_t p1 = TMath::Nint(ceil(prChi2*15));
588 Int_t p2 = TMath::Nint(ceil(prChi2P2*15));
589 Q = p1*16 + p2;
590
591 Double_t x0 = -par[0]/par[1]; // value of x corresponding to à 0 fC
592 threshold = TMath::Nint(ceil(par[3]-x0)); // linear if x < threshold
593
594 if(gPrintLevel==2)
595 {
596 fprintf(pfilen,"%4i %4i %2i",busPatchId,manuId,channelId);
597 fprintf(pfilen," %6.2f %6.4f %10.3e %4.2f %5.3f %5.3f %x\n",
598 par[0], par[1], par[2], par[3], prChi2, prChi2P2, Q);
599 }
600
601 // some tests
602
603 if(par[1]< goodA1Min || par[1]> goodA1Max )
604 {
605 if (gPrintLevel)
606 {
607 cout << " !!!!!!!!!!!!! Bad Calib.: BP= " << busPatchId << " Manu_Id= " << manuId <<
608 " Ch.= " << channelId << ":";
609 cout << " a1 = " << par[1] << " out of limit : [" << goodA1Min << "," << goodA1Max <<
610 "]" << endl;
611 }
612 Q=0;
613 nBadChannel++;
614 }
615 else if(par[2]< goodA2Min || par[2]> goodA2Max )
616 {
617 if (gPrintLevel)
618 {
619 cout << " !!!!!!!!!!!!! Bad Calib.: BP= " << busPatchId << " Manu_Id= " << manuId
620 << " Ch.= " << channelId << ":";
621 cout << " a2 = " << par[2] << " out of limit : [" << goodA2Min << "," << goodA2Max
622 << "]" << endl;
623 }
624 Q=0;
625 nBadChannel++;
626 }
627 else
628 {
629 sumProbChi2 += prChi2;
630 sumA1 += par[1];
631 sumProbChi2P2 += prChi2P2;
632 sumA2 += par[2];
633
634 if (gPrintLevel)
635 {
636 if(!p1)
637 {
638 cout << " ** Warning ** Bad Fit : BP= " << busPatchId << " Manu_Id= " << manuId
639 << "Ch.= " << channelId << ":";
640 cout << " a1 = " << par[1] << " P(chi2)_lin = " << prChi2 << endl ;
641 }
642 if(!p2)
643 {
644 cout << " ** Warning ** Bad Fit : BP= " << busPatchId << " Manu_Id= " << manuId
645 << " Ch.= " << channelId << ":";
646 cout << " a2 = " << par[2] << " P_(chi2)_parab = " << prChi2P2 << endl;
647 }
648 }
649 }
650
651 tg->Fill();
652
653 if (!flatOutputFile.IsNull())
654 {
655 fprintf(pfilew,"%4i %5i %2i %7.4f %10.3e %4i %2x\n",busPatchId,manuId,channelId,par[1],par[2],threshold,Q);
656 }
657
658 // Plots
659
660 if(gPlotLevel){
661 TF1 *f2Calib = new TF1("f2Calib",funcCalib,0.,gkADCMax,NFITPARAMS);
662
663 graphErr = new TGraphErrors(nEntries,pedMean,injCharge,pedSigma,injChargeErr);
664
665 sprintf(graphName,"BusPatch_%d_Manu_%d_Ch_%d",busPatchId, manuId,channelId);
666
667 graphErr->SetTitle(graphName);
668 graphErr->SetMarkerColor(3);
669 graphErr->SetMarkerStyle(12);
670 graphErr->Write(graphName);
671
672 sprintf(graphName,"f2_BusPatch_%d_Manu_%d_Ch_%d",busPatchId, manuId,channelId);
673 f2Calib->SetTitle(graphName);
674 f2Calib->SetLineColor(4);
675 f2Calib->SetParameters(par);
676 f2Calib->Write(graphName);
677
678 delete graphErr;
679 delete f2Calib;//LA
680 }
681 }
682 nmanu++;
683 }
684
685 // file outputs for gain
686 if (!flatOutputFile.IsNull())
687 fileout.close();
688
689 tg->Write();
690 histoFile->Close();
691
692 delete f1;
693 delete f2;
694
695 //OutPut
696 if (gPrintLevel)
697 {
698 cout << "\n Nb of Manu in raw data = " << nmanu << " (" << nmanu*64 << " channels)" << endl;
699 cout << "\n Nb of calibrated channels = " << nmanu*64 - nBadChannel << " (" << goodA1Min << "<a1<" << goodA1Max
700 << " and " << goodA2Min << "<a2<" << goodA2Max << ") " << endl;
701 cout << "\n Nb of UNcalibrated channels = " << nBadChannel << " (a1 or a2 out of range)\n" << endl;
702
703 Double_t meanA1 = sumA1/(nmanu*64 - nBadChannel);
704 Double_t meanProbChi2 = sumProbChi2/(nmanu*64 - nBadChannel);
705 Double_t meanA2 = sumA2/(nmanu*64 - nBadChannel);
706 Double_t meanProbChi2P2 = sumProbChi2P2/(nmanu*64 - nBadChannel);
707
708 Double_t capaManu = 0.2; // pF
709 cout << "\n linear fit : <a1> = " << meanA1 << "\t <gain> = " << 1./(meanA1*capaManu)
710 << " mV/fC (capa= " << capaManu << " pF)" << endl;
711 cout << " Prob(chi2)> = " << meanProbChi2 << endl;
712 cout << "\n parabolic fit: <a2> = " << meanA2 << endl;
713 cout << " Prob(chi2)> = " << meanProbChi2P2 << "\n" << endl;
714 }
715}
716
717
718AliHLTMUONTrackerCalibratorComponent::AliHLTMUONTrackerCalibratorComponent() :
719 AliHLTCalibrationProcessor(),
720 fFitter(NULL),
721 fSkipEvents(0),
722 fMaxEvents(1000000),
723 fFlatOutputFile(),
724 fCrocusOutputFile(),
725 fCrocusConfigFile(),
726 fInjCharge(0),
727 fNoSigma(3),
728 fThreshold(-1)
729{
730 /// Default contructor.
731}
732
733
734AliHLTMUONTrackerCalibratorComponent::~AliHLTMUONTrackerCalibratorComponent()
735{
736 /// Default destructor.
737
738 assert( fFitter == NULL );
739}
740
741
742const char* AliHLTMUONTrackerCalibratorComponent::GetComponentID()
743{
744 /// Inherited from AliHLTComponent.
745 /// Returns the component ID string for this component type.
746
747 return AliHLTMUONConstants::TrackerCalibratorId();
748}
749
750void AliHLTMUONTrackerCalibratorComponent::GetInputDataTypes(
751 vector<AliHLTComponentDataType>& list
752 )
753{
754 /// Inherited from AliHLTComponent.
755 /// Returns the list of input block types expected by this component.
756
757 list.clear();
668eee9f 758 list.push_back( AliHLTMUONConstants::DDLRawDataType() );
8134dd2e 759}
760
761AliHLTComponentDataType AliHLTMUONTrackerCalibratorComponent::GetOutputDataType()
762{
763 /// Inherited from AliHLTComponent.
764 /// Returns the type of output block generated by this component.
765
766 //TODO: fix.
668eee9f 767 return AliHLTMUONConstants::DDLRawDataType();
8134dd2e 768}
769
770void AliHLTMUONTrackerCalibratorComponent::GetOutputDataSize(
771 unsigned long& constBase, double& inputMultiplier
772 )
773{
774 /// Inherited from AliHLTComponent.
775 /// Returns an estimate of the expected output data size.
776
777 constBase = 0;
778 inputMultiplier = 2.; //TODO: is this the correct estimate.
779}
780
781AliHLTComponent* AliHLTMUONTrackerCalibratorComponent::Spawn()
782{
783 /// Inherited from AliHLTComponent.
784 /// Creates a new instance of AliHLTMUONTrackerCalibratorComponent.
785
786 return new AliHLTMUONTrackerCalibratorComponent();
787}
788
789
790Int_t AliHLTMUONTrackerCalibratorComponent::ScanArgument(int argc, const char** argv)
791{
792 /// Inherited from AliHLTCalibrationProcessor.
793 /// Parses the command line parameters.
794
795 TString arg = argv[0];
796
797 if (arg.CompareTo("-h") == 0)
798 {
799 HLTInfo("******************* usage **********************");
800 HLTInfo("The available options are :");
801 HLTInfo("-h help (this screen)");
802 HLTInfo("");
803 HLTInfo(" Input");
804 HLTInfo("-c <Crocus config. file> (default = %s)", fCrocusConfigFile.Data());
805 HLTInfo("");
806 HLTInfo(" Output");
807 HLTInfo("-a <Flat ASCII file> (default = %s)", fFlatOutputFile.Data());
808 HLTInfo("-o <CROUCUS cmd file> (default = %s)", fCrocusOutputFile.Data());
809 HLTInfo("");
810 HLTInfo(" Options");
811 HLTInfo("-d <print level> (default = %d)", gPrintLevel);
812 HLTInfo("-g <plot level> (default = %d)", gPlotLevel);
813 HLTInfo("-l <DAC level> (default = %d)", fInjCharge);
814 HLTInfo("-s <skip events> (default = %d)", fSkipEvents);
815 HLTInfo("-n <max events> (default = %d)", fMaxEvents);
816 HLTInfo("-p <n sigmas> (default = %f)", fNoSigma);
817 HLTInfo("-r root file data for gain(default = %s)", gHistoFileName_gain.Data());
818 HLTInfo("-t <threshold (-1 = no)> (default = %d)", fThreshold);
819 HLTInfo("-e <execute ped/gain> (default = %s)", gCommand.Data());
820 HLTInfo("-e <gain create> make gain & create a new root file");
821 HLTInfo("-e <gain> make gain & update root file");
822 HLTInfo("-e <gain compute> make gain & compute gains");
823 return 0; // Zero parameters parsed.
824 }
825
826 if (arg.CompareTo("-a") == 0)
827 {
828 if (argc < 2) return -EPROTO;
829 fFlatOutputFile = argv[1];
830 return 1; // 1 parameter parsed.
831 }
832 if (arg.CompareTo("-o") == 0)
833 {
834 if (argc < 2) return -EPROTO;
835 fCrocusOutputFile = argv[1];
836 return 1; // 1 parameter parsed.
837 }
838 if (arg.CompareTo("-c") == 0)
839 {
840 if (argc < 2) return -EPROTO;
841 fCrocusConfigFile = argv[1];
842 return 1; // 1 parameter parsed.
843 }
844 if (arg.CompareTo("-e") == 0)
845 {
846 if (argc < 2) return -EPROTO;
847 gCommand = argv[1];
848 gCommand.ToLower(); // set command to lower case
849 return 1; // 1 parameter parsed.
850 }
851 if (arg.CompareTo("-d") == 0)
852 {
853 if (argc < 2) return -EPROTO;
854 gPrintLevel = atoi(argv[1]);
855 return 1; // 1 parameter parsed.
856 }
857 if (arg.CompareTo("-g") == 0)
858 {
859 if (argc < 2) return -EPROTO;
860 gPlotLevel = atoi(argv[1]);
861 return 1; // 1 parameter parsed.
862 }
863 if (arg.CompareTo("-s") == 0)
864 {
865 if (argc < 2) return -EPROTO;
866 fSkipEvents = atoi(argv[1]);
867 return 1; // 1 parameter parsed.
868 }
869 if (arg.CompareTo("-l") == 0)
870 {
871 if (argc < 2) return -EPROTO;
872 fInjCharge = atoi(argv[1]);
873 return 1; // 1 parameter parsed.
874 }
875 if (arg.CompareTo("-n") == 0)
876 {
877 if (argc < 2) return -EPROTO;
878 sscanf(argv[1],"%d",&fMaxEvents);
879 return 1; // 1 parameter parsed.
880 }
881 if (arg.CompareTo("-p") == 0)
882 {
883 if (argc < 2) return -EPROTO;
884 sscanf(argv[1],"%lf",&fNoSigma);
885 return 1; // 1 parameter parsed.
886 }
887 if (arg.CompareTo("-r") == 0)
888 {
889 if (argc < 2) return -EPROTO;
890 gHistoFileName_gain = argv[1];
891 return 1; // 1 parameter parsed.
892 }
893 if (arg.CompareTo("-t") == 0)
894 {
895 if (argc < 2) return -EPROTO;
896 sscanf(argv[1],"%d",&fThreshold);
897 return 1; // 1 parameter parsed.
898 }
899
900 // Do not know what this argument is so return an error code.
901 HLTError("Bad argument %s (please check with -h)", arg.Data());
902 return -EINVAL;
903}
904
905
906Int_t AliHLTMUONTrackerCalibratorComponent::InitCalibration()
907{
908 /// Inherited from AliHLTCalibrationProcessor.
909 /// Initialise the calibration component.
910
911 TFitter* fFitter = new TFitter(NFITPARAMS);
912 TVirtualFitter::SetFitter(fFitter);
913
914 Int_t status;
915
916 gPedestalStore = new AliMUON2DMap(kFALSE);
917
918 gPedMeanHisto = 0x0;
919 gPedSigmaHisto = 0x0;
920
921 // once we have a configuration file in db
922 // copy locally a file from daq detector config db
923 // The current detector is identified by detector code in variable
924 // DATE_DETECTOR_CODE. It must be defined.
925 // If environment variable DAQDA_TEST_DIR is defined, files are copied from DAQDA_TEST_DIR
926 // instead of the database. The usual environment variables are not needed.
927 if (!fCrocusConfigFile.IsNull()) {
928 //status = daqDA_DB_getFile("myconfig", fCrocusConfigFile.Data());
929 status = 0;
930 if (status) {
931 printf("Failed to get config file : %d\n",status);
932 return -1;
933 }
934 }
935
936 return 0;
937}
938
939
940Int_t AliHLTMUONTrackerCalibratorComponent::DeinitCalibration()
941{
942 /// Inherited from AliHLTCalibrationProcessor.
943 /// Cleanup the calibration component releasing allocated memory.
944
945 delete gPedestalStore;
946 gPedestalStore = NULL;
947
948 delete fFitter;
949 fFitter = NULL;
950 TVirtualFitter::SetFitter(0);
951
952 return 0;
953}
954
955
956Int_t AliHLTMUONTrackerCalibratorComponent::ProcessCalibration(
957 const AliHLTComponentEventData& /*evtData*/,
958 AliHLTComponentTriggerData& /*trigData*/
959 )
960{
961 /// Inherited from AliHLTCalibrationProcessor.
962 /// Perform a calibration procedure on the new event.
963
964 // Skip Events if needed
965 if (fSkipEvents > 0)
966 {
967 fSkipEvents--;
968 return 0;
969 }
970
971 // Do not process more than fMaxEvents.
972 if (gNEvents >= fMaxEvents) return 0;
973 if (gNEvents && gNEvents % 100 == 0)
974 HLTInfo("Cumulated events %d", gNEvents);
975
976 gNEvents++;
977
978 gRunNumber = GetRunNo();
979
980 /*
981 Int_t eventType = GetRunType();
982 if (eventType != 7) // PHYSICS_EVENT - from event.h (DATE software)
983 {
984 HLTWarning("This event is not a physics event");
985 return 0;
986 }
987 */
988
989 Int_t status;
990 Int_t busPatchId;
991 UShort_t manuId;
992 UChar_t channelId;
993 UShort_t charge;
994
995 const AliHLTComponentBlockData* iter = NULL;
996
997 // Loop over all DDL raw data input blocks and decode the event.
668eee9f 998 iter = GetFirstInputBlock( AliHLTMUONConstants::DDLRawDataType() );
8134dd2e 999 while (iter != NULL)
1000 {
668eee9f 1001 // Make sure we have the correct muon tracker DDL type.
1002 if (not AliHLTMUONUtils::IsTrackerDDL(iter->fSpecification))
1003 {
1004 iter = GetNextInputBlock();
1005 continue;
1006 }
1007
8134dd2e 1008 // decoding rawdata headers
1009 AliRawReaderMemory* rawReader = new AliRawReaderMemory(
1010 reinterpret_cast<UChar_t*>(iter->fPtr), iter->fSize
1011 );
1012 rawReader->SetEquipmentID(AliHLTMUONUtils::SpecToEquipId(iter->fSpecification));
1013
1014 // decoding MUON payload
1015 AliMUONRawStreamTracker* rawStream = new AliMUONRawStreamTracker(rawReader);
1016
1017 // loops over DDL
1018 rawStream->First();
1019 while( (status = rawStream->Next(busPatchId, manuId, channelId, charge)) )
1020 {
1021 if (gNEvents == 1)
1022 gNChannel++;
1023
1024 // if (gPrintLevel) printf("manuId: %d, channelId: %d charge: %d\n", manuId,
1025 // channelId, charge);
1026
1027 MakePed(busPatchId, (Int_t)manuId, (Int_t)channelId, (Int_t)charge);
1028 } // Next digit
1029
1030 delete rawReader;
1031 delete rawStream;
1032 }
1033
1034 if (gCommand.CompareTo("ped") == 0)
1035 {
1036 Char_t flatFile[256];
1037 sprintf(flatFile,"MuonTrkDA_%d_ped.ped",gRunNumber);
1038 if(fFlatOutputFile.IsNull())fFlatOutputFile=flatFile;
1039 MakePedStore(fFlatOutputFile);
1040 }
1041 else
1042 {
1043 if(fFlatOutputFile.IsNull())fFlatOutputFile="MuonTrkDA_gain.par";
1044 }
1045
1046 // option gain -> update root file with pedestal results
1047 // gain + create -> recreate root file
1048 // gain + comp -> update root file and compute gain parameters
1049
1050 if (gCommand.Contains("gain"))
1051 MakePedStoreForGain(fInjCharge);
1052
1053 if (gCommand.Contains("comp"))
1054 MakeGainStore(fFlatOutputFile);
1055
1056 if (gCommand.CompareTo("comp") != 0)
1057 {
1058 HLTInfo("MUONTRKda : Nb of events used = %d", gNEvents);
1059 }
1060 if (gCommand.CompareTo("ped") == 0)
1061 {
1062 if (!(fCrocusConfigFile.IsNull()))
1063 HLTInfo("MUONTRKda : CROCUS command file generated : %s", fCrocusOutputFile.Data());
1064 else
1065 HLTInfo("MUONTRKda : WARNING no CROCUS command file generated");
1066 HLTInfo("MUONTRKda : Histo file generated for pedestal : %s", gHistoFileName);
1067 }
1068 else
1069 {
1070 HLTInfo("MUONTRKda : Histo file generated for gain : %s", gHistoFileName_gain.Data());
1071 HLTInfo("MUONTRKda : Root file generated : %s", gRootFileName.Data());
1072 }
1073
1074 HLTInfo("MUONTRKda : Flat ASCII file generated : %s", fFlatOutputFile.Data());
1075
1076 //TODO:
1077 // PushBack data to shared memory ...
1078 //PushBack(.....);
1079
1080 return 0;
1081}
1082
1083
1084Int_t AliHLTMUONTrackerCalibratorComponent::ShipDataToFXS(
1085 const AliHLTComponentEventData& /*evtData*/,
1086 AliHLTComponentTriggerData& /*trigData*/
1087 )
1088{
1089 /// Inherited from AliHLTCalibrationProcessor.
1090 /// Push the data to the FXS to ship off to the offline CDB.
1091
1092 //TODO:
1093 // PushBack data to FXS ...
1094 //PushToFXS( ..... ) ;
1095
1096 return 0;
1097}