4d69ae91a3b05f4700a3d774c96ac510acfb6d86
[u/mrichter/AliRoot.git] / HLT / MUON / OnlineAnalysis / AliHLTMUONTrackerCalibratorComponent.cxx
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
31 ClassImp(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>
42 using 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
77 namespace
78 {
79
80 // global variables
81 const Int_t gkNChannels = AliMpConstants::ManuNofChannels();
82 const Int_t gkADCMax    = 4095;
83
84 AliMUONVStore* gPedestalStore = NULL;
85
86 Int_t  gNManu       = 0;
87 Int_t  gNChannel    = 0;
88 UInt_t gRunNumber   = 0;
89 Int_t  gNEvents     = 0;
90 Int_t  gNDateEvents = 0;
91 Int_t  gPrintLevel  = 1;  // global printout variable
92 Int_t  gPlotLevel  = 0;  // global plot variable
93
94 TH1F*  gPedMeanHisto  = 0x0;
95 TH1F*  gPedSigmaHisto = 0x0;
96 Char_t gHistoFileName[256];
97
98 // used by makegain 
99 TString gHistoFileName_gain = "MuonTrkDA_data.root";
100 TString gRootFileName = "MuonTrkDA_gain.root";
101 Char_t filenam[256] = "MuonTrkDA_gain.param"; // if gPrintLevel  = 2
102
103 TString gCommand("ped");
104
105 }; // end of namespace
106
107
108 //________________
109 Double_t funcLin(Double_t *x, Double_t *par)
110 {
111   return par[0] + par[1]*x[0];
112 }
113
114 //________________
115 Double_t funcParabolic(Double_t *x, Double_t *par)
116 {
117   return par[0]*x[0]*x[0];
118 }
119
120 //________________
121 Double_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 //__________
132 void 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 //________________
158 void 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 //________________
288 void 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 //________________
335 void 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
718 AliHLTMUONTrackerCalibratorComponent::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
734 AliHLTMUONTrackerCalibratorComponent::~AliHLTMUONTrackerCalibratorComponent()
735 {
736         /// Default destructor.
737         
738         assert( fFitter == NULL );
739 }
740
741
742 const 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
750 void 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();
758         list.push_back( AliHLTMUONConstants::TrackingDDLRawDataType() );
759 }
760
761 AliHLTComponentDataType AliHLTMUONTrackerCalibratorComponent::GetOutputDataType()
762 {
763         /// Inherited from AliHLTComponent.
764         /// Returns the type of output block generated by this component.
765
766         //TODO: fix.
767         return AliHLTMUONConstants::TrackingDDLRawDataType();
768 }
769
770 void 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
781 AliHLTComponent* AliHLTMUONTrackerCalibratorComponent::Spawn()
782 {
783         /// Inherited from AliHLTComponent.
784         /// Creates a new instance of AliHLTMUONTrackerCalibratorComponent.
785
786         return new AliHLTMUONTrackerCalibratorComponent();
787 }
788
789
790 Int_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
906 Int_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
940 Int_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
956 Int_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.
998         iter = GetFirstInputBlock( AliHLTMUONConstants::TriggerDDLRawDataType() );
999         while (iter != NULL)
1000         {
1001                 // decoding rawdata headers
1002                 AliRawReaderMemory* rawReader = new AliRawReaderMemory(
1003                                 reinterpret_cast<UChar_t*>(iter->fPtr), iter->fSize
1004                         );
1005                 rawReader->SetEquipmentID(AliHLTMUONUtils::SpecToEquipId(iter->fSpecification));
1006                 
1007                 // decoding MUON payload
1008                 AliMUONRawStreamTracker* rawStream  = new AliMUONRawStreamTracker(rawReader);
1009                 
1010                 // loops over DDL 
1011                 rawStream->First();
1012                 while( (status = rawStream->Next(busPatchId, manuId, channelId, charge)) )
1013                 {
1014                         if (gNEvents == 1)
1015                                 gNChannel++;
1016                         
1017                         //       if (gPrintLevel) printf("manuId: %d, channelId: %d charge: %d\n", manuId, 
1018                         //                           channelId, charge);
1019                         
1020                         MakePed(busPatchId, (Int_t)manuId, (Int_t)channelId, (Int_t)charge);
1021                 } // Next digit
1022                 
1023                 delete rawReader;
1024                 delete rawStream;
1025         }
1026         
1027         if (gCommand.CompareTo("ped") == 0)
1028         {
1029                 Char_t flatFile[256];
1030                 sprintf(flatFile,"MuonTrkDA_%d_ped.ped",gRunNumber);
1031                 if(fFlatOutputFile.IsNull())fFlatOutputFile=flatFile;
1032                 MakePedStore(fFlatOutputFile);
1033         }
1034         else
1035         {
1036                 if(fFlatOutputFile.IsNull())fFlatOutputFile="MuonTrkDA_gain.par";
1037         }
1038         
1039         // option gain -> update root file with pedestal results
1040         // gain + create -> recreate root file
1041         // gain + comp -> update root file and compute gain parameters
1042         
1043         if (gCommand.Contains("gain"))
1044                 MakePedStoreForGain(fInjCharge);
1045         
1046         if (gCommand.Contains("comp"))
1047                 MakeGainStore(fFlatOutputFile);
1048         
1049         if (gCommand.CompareTo("comp") != 0)
1050         {
1051                 HLTInfo("MUONTRKda : Nb of events used     = %d", gNEvents);
1052         }
1053         if (gCommand.CompareTo("ped") == 0)
1054         {
1055                 if (!(fCrocusConfigFile.IsNull()))
1056                         HLTInfo("MUONTRKda : CROCUS command file generated : %s", fCrocusOutputFile.Data());
1057                 else
1058                         HLTInfo("MUONTRKda : WARNING no CROCUS command file generated");
1059                 HLTInfo("MUONTRKda : Histo file generated for pedestal : %s", gHistoFileName);
1060         }
1061         else
1062         {
1063                 HLTInfo("MUONTRKda : Histo file generated for gain     : %s", gHistoFileName_gain.Data());
1064                 HLTInfo("MUONTRKda : Root file generated               : %s", gRootFileName.Data());
1065         }
1066         
1067         HLTInfo("MUONTRKda : Flat ASCII file generated         : %s", fFlatOutputFile.Data());
1068         
1069         //TODO:
1070         // PushBack data to shared memory ...
1071         //PushBack(.....);
1072
1073         return 0;
1074 }
1075
1076
1077 Int_t AliHLTMUONTrackerCalibratorComponent::ShipDataToFXS(
1078                 const AliHLTComponentEventData& /*evtData*/,
1079                 AliHLTComponentTriggerData& /*trigData*/
1080         )
1081 {
1082         /// Inherited from AliHLTCalibrationProcessor.
1083         /// Push the data to the FXS to ship off to the offline CDB.
1084         
1085         //TODO:
1086         // PushBack data to FXS ...
1087         //PushToFXS( ..... ) ;
1088         
1089         return 0;
1090 }