]> git.uio.no Git - u/mrichter/AliRoot.git/blame - MUON/MUONTRKda.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
patched warnings and few bug fixes for empty events
[u/mrichter/AliRoot.git] / MUON / MUONTRKda.cxx
CommitLineData
830e3f26 1/*
2Contact: Jean-Luc Charvet <jean-luc.charvet@cea.fr>
c4ac9a95 3 Link: http://aliceinfo.cern.ch/static/Offline/dimuon/muon_html/README_Mchda
830e3f26 4Run Type: PEDESTAL, CALIBRATION
c4ac9a95 5 DA Type: LDC
6 Number of events needed: 400 events for pedestal and each calibration run
7 Input Files: Rawdata file (DATE format)
8 Output Files: local dir (not persistent) -> MUONTRKda_ped_<run#>.ped , MUONTRKda_gain_<run#>.par
9 FXS -> run<#>_MCH_<ldc>_PEDESTALS, run<#>_MCH_<ldc>_GAINS
10 Trigger types used:
830e3f26 11*/
12
c758ba48 13/**************************************************************************
c4ac9a95 14* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
15 * *
16 * Author: The ALICE Off-line Project. *
17 * Contributors are mentioned in the code where appropriate. *
18 * *
19 * Permission to use, copy, modify and distribute this software and its *
20 * documentation strictly for non-commercial purposes is hereby granted *
21 * without fee, provided that the above copyright notice appears in all *
22 * copies and that both the copyright notice and this permission notice *
23 * appear in the supporting documentation. The authors make no claims *
24 * about the suitability of this software for any purpose. It is *
25 * provided "as is" without express or implied warranty. *
26 **************************************************************************/
53f515d3 27
28/* $Id$ */
29
81028269 30/*
c4ac9a95 31 -------------------------------------------------------------------------
4a3a8cae 32 2008-10-08 New version: MUONTRKda.cxx,v 1.14
c4ac9a95 33 -------------------------------------------------------------------------
3f637cfc 34
c4ac9a95 35 Version for MUONTRKda MUON tracking
36 (A. Baldisseri, J.-L. Charvet & Ch. Finck)
3f637cfc 37
38
c4ac9a95 39 Rem: AliMUON2DMap stores all channels, even those which are not connected
40 if pedMean == -1, channel not connected to a pad
3f637cfc 41
3f637cfc 42
43*/
44extern "C" {
45#include <daqDA.h>
46}
47
48#include "event.h"
49#include "monitor.h"
50
51#include <Riostream.h>
52#include <stdio.h>
53#include <stdlib.h>
5c93a711 54#include <math.h>
3f637cfc 55
56//AliRoot
96fa832f 57#include "AliMUONLogger.h"
c4ac9a95 58#include "AliMUONRawStreamTrackerHP.h"
3f637cfc 59#include "AliMUONDspHeader.h"
60#include "AliMUONBlockHeader.h"
61#include "AliMUONBusStruct.h"
62#include "AliMUONDDLTracker.h"
2bf54af2 63#include "AliMUONVStore.h"
3f637cfc 64#include "AliMUON2DMap.h"
0aa03544 65#include "AliMUONCalibParamND.h"
3f637cfc 66#include "AliMpIntPair.h"
2bf54af2 67#include "AliMpConstants.h"
3f637cfc 68#include "AliRawReaderDate.h"
c4ac9a95 69#include "AliRawDataErrorLog.h"
3f637cfc 70
3f637cfc 71//ROOT
72#include "TFile.h"
53f515d3 73#include "TSystem.h"
3f637cfc 74#include "TTree.h"
75#include "TH1F.h"
76#include "TString.h"
77#include "TStopwatch.h"
78#include "TMath.h"
79#include "TTimeStamp.h"
80#include "TGraphErrors.h"
81#include "TF1.h"
fd99693f 82#include "TROOT.h"
83#include "TPluginManager.h"
a76f513d 84#include "TFitter.h"
c4ac9a95 85#include "THashTable.h"
4a3a8cae 86#include <THashList.h>
3f637cfc 87
a76f513d 88#define NFITPARAMS 4
3f637cfc 89
90// global variables
90293ba6 91const Int_t gkNChannels = AliMpConstants::ManuNofChannels();
92const Int_t gkADCMax = 4095;
93
94AliMUONVStore* gPedestalStore = new AliMUON2DMap(kFALSE);
95
96Int_t gNManu = 0;
97Int_t gNChannel = 0;
98UInt_t gRunNumber = 0;
99Int_t gNEvents = 0;
c4ac9a95 100Int_t gNEventsRecovered = 0;
90293ba6 101Int_t gNDateEvents = 0;
666b5407 102Int_t gPrintLevel = 1; // global printout variable (others: 2 and 3)
90293ba6 103Int_t gPlotLevel = 0; // global plot variable
830e3f26 104Int_t gFES = 1; // by default FES is used
90293ba6 105
106TH1F* gPedMeanHisto = 0x0;
107TH1F* gPedSigmaHisto = 0x0;
108Char_t gHistoFileName[256];
109
81028269 110// used for computing gain parameters
111Int_t nbpf1 = 6; // linear fit over nbf1 points
112
81028269 113Char_t gHistoFileName_gain[256]="MUONTRKda_gain.data";
53f515d3 114Char_t gRootFileName[256];
115Char_t gOutFolder[256]=".";
116Char_t filename[256];
117Char_t filenam[256]="MUONTRKda_gain";
81028269 118Char_t flatFile[256]="";
53f515d3 119
c4ac9a95 120
121ofstream filcout;
53f515d3 122
123TString flatOutputFile;
124TString logOutputFile;
81028269 125TString logOutputFile_comp;
90293ba6 126TString gCommand("ped");
53f515d3 127TTimeStamp date;
3f637cfc 128
c4ac9a95 129class ErrorCounter : public TNamed
130{
4a3a8cae 131public :
132 ErrorCounter(Int_t bp = 0, Int_t manu = 0, Int_t ev = 1) : busPatch(bp), manuId(manu), events(ev) {}
133 void Increment() {events++;}
134 Int_t BusPatch() {return busPatch;}
135 Int_t ManuId() {return manuId;}
136 Int_t Events() {return events;}
137 Int_t Compare(const TObject*) const;
138
139 void Print(Option_t* option="") const
140 {
141 TNamed::Print(option);
142 cout<<"bp "<<busPatch<<" events "<<events<<endl;
143 }
144 void Print_uncal(Option_t* option="") const
145 {
146 TNamed::Print(option);
147 cout<<"bp = "<<busPatch<< " manu = " << manuId << " uncal = "<< events <<endl;
148 }
149
150private :
151 Int_t busPatch; // Buspath ID
152 Int_t manuId; // Manu ID
153 Int_t events; // Events with error in this buspatch
154};
155
156Int_t ErrorCounter::Compare(const TObject* obj) const
157{
158 Int_t patch1, patch2, manu1, manu2;
159 patch1 = busPatch;
160 manu1 = manuId;
161 patch2 = ((ErrorCounter*)obj)->BusPatch();
162 manu2 = ((ErrorCounter*)obj)->ManuId();
163
164 if (patch1 == patch2)
165 {
166 if (manu1 == manu2)
c4ac9a95 167 {
4a3a8cae 168 return 0;
c4ac9a95 169 }
4a3a8cae 170 else
171 return (manu1 >= manu2) ? 1 : -1;
172 }
173 else
174 return (patch1 >= patch2) ? 1 : -1;
c4ac9a95 175};
4a3a8cae 176
177
178// Table for buspatches with parity errors
c4ac9a95 179THashTable* gErrorBuspatchTable = new THashTable(100,2);
180
4a3a8cae 181// Table for uncalibrated buspatches and manus
182// THashTable* gUncalBuspatchManuTable = new THashTable(1000,2);
183 THashList* gUncalBuspatchManuTable = new THashList(1000,2);
184
c4ac9a95 185
3f637cfc 186// funtions
187
90293ba6 188
189//________________
190Double_t funcLin(Double_t *x, Double_t *par)
191{
c4ac9a95 192 return par[0] + par[1]*x[0];
90293ba6 193}
194
3f637cfc 195//________________
90293ba6 196Double_t funcParabolic(Double_t *x, Double_t *par)
3f637cfc 197{
c4ac9a95 198 return par[0]*x[0]*x[0];
90293ba6 199}
3f637cfc 200
90293ba6 201//________________
202Double_t funcCalib(Double_t *x, Double_t *par)
203{
c4ac9a95 204 Double_t xLim= par[3];
3f637cfc 205
c4ac9a95 206 if(x[0] <= xLim) return par[0] + par[1]*x[0];
3f637cfc 207
c4ac9a95 208 Double_t yLim = par[0]+ par[1]*xLim;
209 return yLim + par[1]*(x[0] - xLim) + par[2]*(x[0] - xLim)*(x[0] - xLim);
3f637cfc 210}
211
212
3f637cfc 213//__________
214void MakePed(Int_t busPatchId, Int_t manuId, Int_t channelId, Int_t charge)
215{
216
c4ac9a95 217 AliMUONVCalibParam* ped =
218 static_cast<AliMUONVCalibParam*>(gPedestalStore->FindObject(busPatchId, manuId));
3f637cfc 219
c4ac9a95 220 if (!ped) {
221 gNManu++;
222 ped = new AliMUONCalibParamND(2, gkNChannels,busPatchId, manuId, -1.); // put default wise -1, not connected channel
223 gPedestalStore->Add(ped);
224 }
3f637cfc 225
c4ac9a95 226 // if (gNEvents == 0) {
227 // ped->SetValueAsDouble(channelId, 0, 0.);
228 // ped->SetValueAsDouble(channelId, 1, 0.);
229 // }
230
231 // Initialization for the first value
232 if (ped->ValueAsDouble(channelId, 0) == -1) ped->SetValueAsDouble(channelId, 0, 0.);
233 if (ped->ValueAsDouble(channelId, 1) == -1) ped->SetValueAsDouble(channelId, 1, 0.);
3f637cfc 234
c4ac9a95 235 Double_t pedMean = ped->ValueAsDouble(channelId, 0) + (Double_t) charge;
236 Double_t pedSigma = ped->ValueAsDouble(channelId, 1) + (Double_t) charge*charge;
3f637cfc 237
c4ac9a95 238 ped->SetValueAsDouble(channelId, 0, pedMean);
239 ped->SetValueAsDouble(channelId, 1, pedSigma);
3f637cfc 240
241}
242
243//________________
c4ac9a95 244void MakePedStore(TString flatOutputFile_1 = "")
245//void MakePedStore()
3f637cfc 246{
c4ac9a95 247// TTimeStamp date;
248 Double_t pedMean;
249 Double_t pedSigma;
250 ofstream fileout;
251 Int_t busPatchId;
252 Int_t manuId;
253 Int_t channelId;
254
255// histo
256 TFile* histoFile = 0;
257 TTree* tree = 0;
258 if (gCommand.CompareTo("ped") == 0)
259 {
260 sprintf(gHistoFileName,"%s/MUONTRKda_ped_%d.root",gOutFolder,gRunNumber);
261 histoFile = new TFile(gHistoFileName,"RECREATE","MUON Tracking pedestals");
262
263 Char_t name[255];
264 Char_t title[255];
265 sprintf(name,"pedmean_allch");
266 sprintf(title,"Pedestal mean all channels");
267 Int_t nx = 4096;
268 Int_t xmin = 0;
269 Int_t xmax = 4095;
270 gPedMeanHisto = new TH1F(name,title,nx,xmin,xmax);
271 gPedMeanHisto->SetDirectory(histoFile);
272
273 sprintf(name,"pedsigma_allch");
274 sprintf(title,"Pedestal sigma all channels");
275 nx = 201;
276 xmin = 0;
277 xmax = 200;
278 gPedSigmaHisto = new TH1F(name,title,nx,xmin,xmax);
279 gPedSigmaHisto->SetDirectory(histoFile);
280
281 tree = new TTree("t","Pedestal tree");
282 tree->Branch("bp",&busPatchId,"bp/I");
283 tree->Branch("manu",&manuId,",manu/I");
284 tree->Branch("channel",&channelId,",channel/I");
285 tree->Branch("pedMean",&pedMean,",pedMean/D");
286 tree->Branch("pedSigma",&pedSigma,",pedSigma/D");
3f637cfc 287 }
288
c4ac9a95 289 if (!flatOutputFile_1.IsNull()) {
290 fileout.open(flatOutputFile_1.Data());
291 fileout<<"//===========================================================================" << endl;
292 fileout<<"// Pedestal file calculated by MUONTRKda"<<endl;
293 fileout<<"//===========================================================================" << endl;
294 fileout<<"// * Run : " << gRunNumber << endl;
295 fileout<<"// * Date : " << date.AsString("l") <<endl;
296 fileout<<"// * Statictics : " << gNEvents << endl;
297 fileout<<"// * # of MANUS : " << gNManu << endl;
298 fileout<<"// * # of channels : " << gNChannel << endl;
299 if (gErrorBuspatchTable->GetSize())
300 {
301 fileout<<"//"<<endl;
302 fileout<<"// * Buspatches with less statistics (due to parity errors)"<<endl;
303 TIterator* iter = gErrorBuspatchTable->MakeIterator();
304 ErrorCounter* parityerror;
305 while((parityerror = (ErrorCounter*) iter->Next()))
306 {
307 fileout<<"// bp "<<parityerror->BusPatch()<<" events used "<<gNEvents-parityerror->Events()<<endl;
308 }
309
310 }
311 fileout<<"//"<<endl;
312 fileout<<"//---------------------------------------------------------------------------" << endl;
313 fileout<<"//---------------------------------------------------------------------------" << endl;
314 fileout<<"// BP MANU CH. MEAN SIGMA"<<endl;
315 fileout<<"//---------------------------------------------------------------------------" << endl;
3f637cfc 316
c4ac9a95 317 }
318 // print in logfile
319 if (gErrorBuspatchTable->GetSize())
320 {
321 cout<<"\n* Buspatches with less statistics (due to parity errors)"<<endl;
322 filcout<<"\n* Buspatches with less statistics (due to parity errors)"<<endl;
323 TIterator* iter = gErrorBuspatchTable->MakeIterator();
324 ErrorCounter* parityerror;
325 while((parityerror = (ErrorCounter*) iter->Next()))
326 {
327 cout<<" bp "<<parityerror->BusPatch()<<": events used = "<<gNEvents-parityerror->Events()<<endl;
328 filcout<<" bp "<<parityerror->BusPatch()<<": events used = "<<gNEvents-parityerror->Events()<<endl;
329 }
330
666b5407 331 }
3f637cfc 332
c4ac9a95 333
334// iterator over pedestal
335 TIter next(gPedestalStore->CreateIterator());
336 AliMUONVCalibParam* ped;
3f637cfc 337
c4ac9a95 338 while ( ( ped = dynamic_cast<AliMUONVCalibParam*>(next() ) ) )
339 {
340 busPatchId = ped->ID0();
341 manuId = ped->ID1();
342 Int_t eventCounter;
343
344 // Correct the number of events for buspatch with errors
345 char bpname[256];
346 ErrorCounter* errorCounter;
347 sprintf(bpname,"bp%d",busPatchId);
348 if ((errorCounter = (ErrorCounter*)gErrorBuspatchTable->FindObject(bpname)))
349 {
350 eventCounter = gNEvents - errorCounter->Events();
351 }
352 else
353 {
354 eventCounter = gNEvents;
355 }
356
357 for (channelId = 0; channelId < ped->Size() ; ++channelId) {
358 pedMean = ped->ValueAsDouble(channelId, 0);
359
360 if (pedMean > 0) { // connected channels
361
362 ped->SetValueAsDouble(channelId, 0, pedMean/(Double_t)eventCounter);
363
364 pedMean = ped->ValueAsDouble(channelId, 0);
365 pedSigma = ped->ValueAsDouble(channelId, 1);
366
367 ped->SetValueAsDouble(channelId, 1, TMath::Sqrt(TMath::Abs(pedSigma/(Double_t)eventCounter - pedMean*pedMean)));
368
4a3a8cae 369 pedMean = ped->ValueAsDouble(channelId, 0);
c4ac9a95 370 pedSigma = ped->ValueAsDouble(channelId, 1);
371
372
373 if (!flatOutputFile_1.IsNull()) {
374 fileout << "\t" << busPatchId << "\t" << manuId <<"\t"<< channelId << "\t"
375 << pedMean <<"\t"<< pedSigma << endl;
376 }
377
378 if (gCommand.CompareTo("ped") == 0)
379 {
380 gPedMeanHisto->Fill(pedMean);
381 gPedSigmaHisto->Fill(pedSigma);
382
383 tree->Fill();
384 }
385 }
386 }
387}
388
389// file outputs
390if (!flatOutputFile_1.IsNull()) fileout.close();
391
392if (gCommand.CompareTo("ped") == 0)
393{
394 histoFile->Write();
395 histoFile->Close();
396}
3f637cfc 397
398// delete tree;
399
3f637cfc 400}
401
402//________________
90293ba6 403void MakePedStoreForGain(Int_t injCharge)
3f637cfc 404{
c4ac9a95 405 // store pedestal map in root file
3f637cfc 406
90293ba6 407// Int_t injCharge = 200;
3f637cfc 408
c4ac9a95 409 TTree* tree = 0x0;
410
411 FILE *pfilew=0;
412 if (gCommand.Contains("gain") && !gCommand.Contains("comp")) {
413 if(flatOutputFile.IsNull())
414 {
415 sprintf(filename,"%s_%d_DAC_%d.par",filenam,gRunNumber,injCharge);
416 flatOutputFile=filename;
417 }
418 if(!flatOutputFile.IsNull())
419 {
420 pfilew = fopen (flatOutputFile.Data(),"w");
421
422 fprintf(pfilew,"//DUMMY FILE (to prevent Shuttle failure)\n");
423 fprintf(pfilew,"//================================================\n");
424 fprintf(pfilew,"// MUONTRKda: Calibration run \n");
425 fprintf(pfilew,"//=================================================\n");
426 fprintf(pfilew,"// * Run : %d \n",gRunNumber);
427 fprintf(pfilew,"// * Date : %s \n",date.AsString("l"));
428 fprintf(pfilew,"// * DAC : %d \n",injCharge);
429 fprintf(pfilew,"//-------------------------------------------------\n");
430 fclose(pfilew);
431 }
432 }
3f637cfc 433
c4ac9a95 434 if(gPrintLevel==2)
81028269 435 {
c4ac9a95 436 // compute and store pedestals
437 sprintf(flatFile,"%s/%s_%d_DAC_%d.ped",gOutFolder,filenam,gRunNumber,injCharge);
4a3a8cae 438 cout << "\nMUONTRKda : Flat file generated : " << flatFile << "\n";
c4ac9a95 439 MakePedStore(flatFile);
81028269 440 }
c4ac9a95 441 else
442 MakePedStore();
443
444 TString mode("UPDATE");
445
446 if (gCommand.Contains("cre")) {
447 mode = "RECREATE";
448 }
449 TFile* histoFile = new TFile(gHistoFileName_gain, mode.Data(), "MUON Tracking Gains");
450
451 // second argument should be the injected charge, taken from config crocus file
452 // put also info about run number could be usefull
453 AliMpIntPair* pair = new AliMpIntPair(gRunNumber, injCharge);
454
455 if (mode.CompareTo("UPDATE") == 0) {
456 tree = (TTree*)histoFile->Get("t");
457 tree->SetBranchAddress("run",&pair);
458 tree->SetBranchAddress("ped",&gPedestalStore);
459
460 } else {
461 tree = new TTree("t","Pedestal tree");
462 tree->Branch("run", "AliMpIntPair",&pair);
463 tree->Branch("ped", "AliMUON2DMap",&gPedestalStore);
464 tree->SetBranchAddress("run",&pair);
465 tree->SetBranchAddress("ped",&gPedestalStore);
466
81028269 467 }
c4ac9a95 468
469 tree->Fill();
470 tree->Write("t", TObject::kOverwrite); // overwrite the tree
471 histoFile->Close();
472
473 delete pair;
3f637cfc 474}
475
476//________________
53f515d3 477void MakeGainStore()
3f637cfc 478{
4a3a8cae 479 ofstream filcouc;
480
481 Int_t nInit = 1; // DAC=0 excluded from fit procedure
482 Double_t goodA1Min = 0.5;
483 Double_t goodA1Max = 2.;
484 // Double_t goodA1Min = 0.7;
485 // Double_t goodA1Max = 1.7;
486 Double_t goodA2Min = -0.5E-03;
487 Double_t goodA2Max = 1.E-03;
488
489 Int_t num_RUN[15];
490
491 // open file mutrkgain.root
492 // read again the pedestal for the calibration runs (9 runs ?)
493 // need the injection charge from config file (to be done)
494 // For each channel make a TGraphErrors (mean, sigma) vs injected charge
495 // Fit with a polynomial fct
496 // store the result in a flat file.
497
498
499 TFile* histoFile = new TFile(gHistoFileName_gain);
500
501 AliMUON2DMap* map[11];
502 AliMUONVCalibParam* ped[11];
503 AliMpIntPair* run[11];
504
505 //read back from root file
506 TTree* tree = (TTree*)histoFile->Get("t");
507 Int_t nEntries = tree->GetEntries();
508 Int_t nbpf2 = nEntries - (nInit + nbpf1) + 1; // nb pts used for 2nd order fit
509
510 // read back info
511 for (Int_t i = 0; i < nEntries; ++i) {
512 map[i] = 0x0;
513 run[i] = 0x0;
514 tree->SetBranchAddress("ped",&map[i]);
515 tree->SetBranchAddress("run",&run[i]);
516 tree->GetEvent(i);
517 // std::cout << map[i] << " " << run[i] << std::endl;
518 }
519 //jlc_feb_08 modif: gRunNumber=(UInt_t)run[0]->GetFirst();
520 gRunNumber=(UInt_t)run[nEntries-1]->GetFirst();
521 // sscanf(getenv("DATE_RUN_NUMBER"),"%d",&gRunNumber);
522
523 Double_t pedMean[11];
524 Double_t pedSigma[11];
525 for ( Int_t i=0 ; i<11 ; i++) {pedMean[i]=0.;pedSigma[i]=1.;};
526 Double_t injCharge[11];
527 Double_t injChargeErr[11];
528 for ( Int_t i=0 ; i<11 ; i++) {injCharge[i]=0.;injChargeErr[i]=1.;};
529
530 // some print
531 cout<<"\n ******** MUONTRKda for Gain computing (Run = " << gRunNumber << ")\n" << endl;
532 cout<<" * Date : " << date.AsString("l") << "\n" << endl;
533 cout << " Entries = " << nEntries << " DAC values \n" << endl;
534 for (Int_t i = 0; i < nEntries; ++i) {
535 cout<< " Run = " << run[i]->GetFirst() << " DAC = " << run[i]->GetSecond() << endl;
536 num_RUN[i] = run[i]->GetFirst();
537 injCharge[i] = run[i]->GetSecond();
538 injChargeErr[i] = 0.01*injCharge[i];
539 if(injChargeErr[i] <= 1.) injChargeErr[i]=1.;
540 }
541 cout << "" << endl;
542
543 // full print out
544
545 sprintf(filename,"%s/%s_%d.log",gOutFolder,filenam,gRunNumber);
546 logOutputFile_comp=filename;
547
548 filcouc.open(logOutputFile_comp.Data());
549 filcouc<<"//====================================================" << endl;
550 filcouc<<"// MUONTRKda for Gain computing (Run = " << gRunNumber << ")" << endl;
551 filcouc<<"//====================================================" << endl;
552 filcouc<<"// * Date : " << date.AsString("l") << "\n" << endl;
553
554
555
556 // why 2 files ? (Ch. F.)
557 FILE *pfilen = 0;
558 FILE *pfilef = 0;
559 if(gPrintLevel==2)
560 {
561 sprintf(filename,"%s/%s_%d.param",gOutFolder,filenam,gRunNumber);
562 cout << " fit parameter file = " << filename << "\n";
563 pfilen = fopen (filename,"w");
564
565 fprintf(pfilen,"//===================================================================\n");
566 fprintf(pfilen,"// BP MANU CH. a0 a1 a2 xlim P(chi2) p1 P(chi2)2 p2\n");
567 fprintf(pfilen,"//===================================================================\n");
568 fprintf(pfilen,"// * Run : %d \n",gRunNumber);
569 fprintf(pfilen,"//===================================================================\n");
570
571 sprintf(filename,"%s/%s_%d.bad",gOutFolder,filenam,gRunNumber);
572 cout << " Bad channel file = " << filename << "\n";
573 pfilef = fopen (filename,"w");
574
575 fprintf(pfilef,"//=================================================\n");
576 fprintf(pfilef,"// Bad Channel file calculated by MUONTRKda \n");
577 fprintf(pfilef,"//=================================================\n");
578 fprintf(pfilef,"// * Run : %d \n",gRunNumber);
579 fprintf(pfilef,"// * Date : %s \n",date.AsString("l"));
580 fprintf(pfilef,"//=======================================\n");
581 fprintf(pfilef,"// BP MANU CH. a1 a2 thres. Q\n");
582 fprintf(pfilef,"//=======================================\n");
583 }
584
585 FILE *pfilew=0;
586 if(flatOutputFile.IsNull())
587 {
588 sprintf(filename,"%s_%d.par",filenam,gRunNumber);
589 flatOutputFile=filename;
590 }
591 if(!flatOutputFile.IsNull())
592 {
593 pfilew = fopen (flatOutputFile.Data(),"w");
594
595 fprintf(pfilew,"//================================================\n");
596 fprintf(pfilew,"// Calibration file calculated by MUONTRKda \n");
597 fprintf(pfilew,"//=================================================\n");
598 fprintf(pfilew,"// * Run : %d \n",gRunNumber);
599 fprintf(pfilew,"// * Date : %s \n",date.AsString("l"));
600 fprintf(pfilew,"// * Statictics : %d \n",gNEvents);
601 fprintf(pfilew,"// * # of MANUS : %d \n",gNManu);
602 fprintf(pfilew,"// * # of channels : %d \n",gNChannel);
603 fprintf(pfilew,"//-------------------------------------------------\n");
604 if(nInit==0)
605 fprintf(pfilew,"// %d DAC values fit: %d pts (1st order) %d pts (2nd order) \n",nEntries,nbpf1,nbpf2);
606 if(nInit==1)
607 fprintf(pfilew,"// %d DAC values fit: %d pts (1st order) %d pts (2nd order) DAC=0 excluded\n",nEntries,nbpf1,nbpf2);
608 fprintf(pfilew,"// RUN DAC \n");
609 fprintf(pfilew,"//-----------------\n");
610 for (Int_t i = 0; i < nEntries; ++i) {
611 tree->SetBranchAddress("run",&run[i]);
612 fprintf(pfilew,"// %d %5.0f \n",num_RUN[i],injCharge[i]);
613 }
614 fprintf(pfilew,"//=======================================\n");
615 fprintf(pfilew,"// BP MANU CH. a1 a2 thres. Q\n");
616 fprintf(pfilew,"//=======================================\n");
617 }
618
619 FILE *pfilep = 0;
620 if(gPrintLevel==2)
621 {
622 sprintf(filename,"%s/%s_%d.peak",gOutFolder,filenam,gRunNumber);
623 cout << " File containing Peak mean values = " << filename << "\n";
624 pfilep = fopen (filename,"w");
625
626 fprintf(pfilep,"//==============================================================================================================================\n");
627 fprintf(pfilep,"// * Run : %d \n",gRunNumber);
628 fprintf(pfilep,"//==============================================================================================================================\n");
629 fprintf(pfilep,"// BP MANU CH. Ped. <0> <1> <2> <3> <4> <5> <6> <7> <8> <9> <10> \n");
630 fprintf(pfilep,"//==============================================================================================================================\n");
c4ac9a95 631 fprintf(pfilep,"// DAC= %9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f fC\n",injCharge[0],injCharge[1],injCharge[2],injCharge[3],injCharge[4],injCharge[5],injCharge[6],injCharge[7],injCharge[8],injCharge[9],injCharge[10]);
632 fprintf(pfilep,"// %9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f\n",injChargeErr[0],injChargeErr[1],injChargeErr[2],injChargeErr[3],injChargeErr[4],injChargeErr[5],injChargeErr[6],injChargeErr[7],injChargeErr[8],injChargeErr[9],injChargeErr[10]);
4a3a8cae 633 fprintf(pfilep,"//==============================================================================================================================\n");
634 }
635
636
637
638 // plot out
639
640 TFile* gainFile = 0x0;
641 sprintf(gRootFileName,"%s/%s_%d.root",gOutFolder,filenam,gRunNumber);
642 gainFile = new TFile(gRootFileName,"RECREATE");
643
644 Double_t chi2 = 0.;
645 Double_t chi2P2 = 0.;
646 Double_t prChi2 = 0;
647 Double_t prChi2P2 =0;
648 Double_t a0=0.,a1=1.,a2=0.;
649 Int_t busPatchId ;
650 Int_t manuId ;
651 Int_t channelId ;
652 Int_t threshold = 0;
653 Int_t Q = 0;
654 Int_t p1 =0;
655 Int_t p2 =0;
656 Double_t gain=0;
657 Double_t capa=0.2; // internal capacitor (pF)
658
659 TTree *tg = new TTree("tg","TTree avec class Manu_DiMu");
660
661 tg->Branch("bp",&busPatchId, "busPatchId/I");
662 tg->Branch("manu",&manuId, "manuId/I");
663 tg->Branch("channel",&channelId, "channelId/I");
664
665 tg->Branch("a0",&a0, "a0/D");
666 tg->Branch("a1",&a1, "a1/D");
667 tg->Branch("a2",&a2, "a2/D");
668 tg->Branch("Pchi2",&prChi2, "prChi2/D");
669 tg->Branch("Pchi2_2",&prChi2P2, "prChi2P2/D");
670 tg->Branch("Threshold",&threshold, "threshold/I");
671 tg->Branch("Q",&Q, "Q/I");
672 tg->Branch("p1",&p1, "p1/I");
673 tg->Branch("p2",&p2, "p2/I");
674 tg->Branch("gain",&gain, "gain/D");
675
676 char graphName[256];
677
678 // iterates over the first pedestal run
679 TIter next(map[0]->CreateIterator());
680 AliMUONVCalibParam* p;
681
682 Int_t nmanu = 0;
683 Int_t nGoodChannel = 0;
684 Int_t nBadChannel = 0;
685 Int_t nBadChannel_a1 = 0;
686 Int_t nBadChannel_a2 = 0;
687 Int_t noFitChannel = 0;
688 Int_t nplot=0;
689 Double_t sumProbChi2 = 0.;
690 Double_t sumA1 = 0.;
691 Double_t sumProbChi2P2 = 0.;
692 Double_t sumA2 = 0.;
693
694 Double_t x[11], xErr[11], y[11], yErr[11];
695 Double_t xp[11], xpErr[11], yp[11], ypErr[11];
696
697 Int_t uncalcountertotal=0 ;
698
699 while ( ( p = dynamic_cast<AliMUONVCalibParam*>(next() ) ) )
700 {
701 ped[0] = p;
702
703 busPatchId = p->ID0();
704 manuId = p->ID1();
705
706 // read back pedestal from the other runs for the given (bupatch, manu)
707 for (Int_t i = 1; i < nEntries; ++i) {
708 ped[i] = static_cast<AliMUONVCalibParam*>(map[i]->FindObject(busPatchId, manuId));
709 }
710
711 // compute for each channel the gain parameters
712 for ( channelId = 0; channelId < ped[0]->Size() ; ++channelId )
90293ba6 713 {
90293ba6 714
4a3a8cae 715 Int_t n = 0;
716 for (Int_t i = 0; i < nEntries; ++i) {
90293ba6 717
4a3a8cae 718 if (!ped[i]) continue; //shouldn't happen.
719 pedMean[i] = ped[i]->ValueAsDouble(channelId, 0);
720 pedSigma[i] = ped[i]->ValueAsDouble(channelId, 1);
90293ba6 721
4a3a8cae 722 if (pedMean[i] < 0) continue; // not connected
90293ba6 723
4a3a8cae 724 if (pedSigma[i] <= 0) pedSigma[i] = 1.; // should not happen.
725 n++;
726 }
90293ba6 727
90293ba6 728
4a3a8cae 729 // print_peak_mean_values
730 if(gPrintLevel==2)
731 {
90293ba6 732
4a3a8cae 733 fprintf(pfilep,"%4i%5i%5i%10.3f",busPatchId,manuId,channelId,0.);
734 fprintf(pfilep,"%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f \n",pedMean[0],pedMean[1],pedMean[2],pedMean[3],pedMean[4],pedMean[5],pedMean[6],pedMean[7],pedMean[8],pedMean[9],pedMean[10]);
735 fprintf(pfilep," sig= %9.3f%9.3f%9.3f%9.3f%9.3f%9.3f%9.3f%9.3f%9.3f%9.3f%9.3f \n",pedSigma[0],pedSigma[1],pedSigma[2],pedSigma[3],pedSigma[4],pedSigma[5],pedSigma[6],pedSigma[7],pedSigma[8],pedSigma[9],pedSigma[10]);
736 }
90293ba6 737
4a3a8cae 738 // makegain
90293ba6 739
90293ba6 740
4a3a8cae 741 // Fit Method: Linear fit over nbpf1 points + parabolic fit over nbpf2 points)
742 // nInit=1 : 1st pt DAC=0 excluded
90293ba6 743
4a3a8cae 744 // 1. - linear fit over nbpf1 points
90293ba6 745
4a3a8cae 746 Double_t par[4] = {0.,0.5,0.,gkADCMax};
747 Int_t nbs = nEntries - nInit;
748 if(nbs < nbpf1)nbpf1=nbs;
90293ba6 749
4a3a8cae 750 Int_t FitProceed=1;
751 for (Int_t j = 0; j < nbs; ++j)
752 {
753 Int_t k = j + nInit;
754 x[j] = pedMean[k];
755 if(x[j]==0.)FitProceed=0;
756 xErr[j] = pedSigma[k];
757 y[j] = injCharge[k];
758 yErr[j] = injChargeErr[k];
90293ba6 759
4a3a8cae 760 }
53f515d3 761
4a3a8cae 762 TGraphErrors *graphErr;
763 if(!FitProceed) { p1=0; p2=0; noFitChannel++;}
90293ba6 764
4a3a8cae 765 if(FitProceed)
766 {
c4ac9a95 767
4a3a8cae 768 TF1 *f1 = new TF1("f1",funcLin,0.,gkADCMax,2);
769 graphErr = new TGraphErrors(nbpf1, x, y, xErr, yErr);
53f515d3 770
4a3a8cae 771 f1->SetParameters(0,0);
53f515d3 772
4a3a8cae 773 graphErr->Fit("f1","RQ");
53f515d3 774
4a3a8cae 775 chi2 = f1->GetChisquare();
776 f1->GetParameters(par);
53f515d3 777
4a3a8cae 778 delete graphErr;
779 graphErr=0;
780 delete f1;
53f515d3 781
4a3a8cae 782 prChi2 = TMath::Prob(chi2, nbpf1 - 2);
53f515d3 783
4a3a8cae 784 Double_t xLim = pedMean[nInit + nbpf1 - 1];
785 Double_t yLim = par[0]+par[1] * xLim;
90293ba6 786
4a3a8cae 787 a0 = par[0];
788 a1 = par[1];
81028269 789
4a3a8cae 790 // 2. - Translation : new origin (xLim, yLim) + parabolic fit over nbf2 points
81028269 791
4a3a8cae 792 if(nbpf2 > 1)
793 {
794 for (Int_t j = 0; j < nbpf2; j++)
795 {
796 Int_t k = j + (nInit + nbpf1) - 1;
797 xp[j] = pedMean[k] - xLim;
798 xpErr[j] = pedSigma[k];
90293ba6 799
4a3a8cae 800 yp[j] = injCharge[k] - yLim - par[1]*xp[j];
801 ypErr[j] = injChargeErr[k];
3f637cfc 802
4a3a8cae 803 }
90293ba6 804
4a3a8cae 805 TF1 *f2 = new TF1("f2",funcParabolic,0.,gkADCMax,1);
806 graphErr = new TGraphErrors(nbpf2, xp, yp, xpErr, ypErr);
90293ba6 807
4a3a8cae 808 graphErr->Fit(f2,"RQ");
809 chi2P2 = f2->GetChisquare();
810 f2->GetParameters(par);
90293ba6 811
4a3a8cae 812 delete graphErr;
813 graphErr=0;
814 delete f2;
90293ba6 815
4a3a8cae 816 prChi2P2 = TMath::Prob(chi2P2, nbpf2-1);
817 a2 = par[0];
53f515d3 818
4a3a8cae 819 // delete graphErr;
53f515d3 820
4a3a8cae 821 }
90293ba6 822
4a3a8cae 823 par[0] = a0;
824 par[1] = a1;
825 par[2] = a2;
826 par[3] = xLim;
53f515d3 827
828
4a3a8cae 829 // Prints
53f515d3 830
4a3a8cae 831 p1 = TMath::Nint(ceil(prChi2*14))+1;
832 p2 = TMath::Nint(ceil(prChi2P2*14))+1;
833 Q = p1*16 + p2; // fit quality
90293ba6 834
4a3a8cae 835 Double_t x0 = -par[0]/par[1]; // value of x corresponding to à 0 fC
836 threshold = TMath::Nint(ceil(par[3]-x0)); // linear if x < threshold
90293ba6 837
4a3a8cae 838 if(gPrintLevel==2)
839 {
840 fprintf(pfilen,"%4i %4i %2i",busPatchId,manuId,channelId);
841 fprintf(pfilen," %6.2f %6.4f %10.3e %4.2f %5.3f %x %5.3f %x\n",
842 par[0], par[1], par[2], par[3], prChi2, p1, prChi2P2, p2);
843 }
90293ba6 844
4a3a8cae 845 // some tests
846
847 if(par[1]< goodA1Min || par[1]> goodA1Max)
848 {
849 p1=0;
850 nBadChannel_a1++;
851 if (gPrintLevel && nBadChannel_a1 < 1)
852 {
853 cout << " !!!!! " << nBadChannel_a1 << " !!!!!!!! Bad Calib.: BP= " << busPatchId << " Manu_Id= " << manuId <<
854 " Ch.= " << channelId << ":";
855 cout << " a1 = " << par[1] << " out of limit : [" << goodA1Min << "," << goodA1Max <<
856 "]" << endl;
857 }
858 }
90293ba6 859
53f515d3 860
4a3a8cae 861 if(par[2]< goodA2Min || par[2]> goodA2Max)
862 {
863 p2=0;
864 nBadChannel_a2++;
865 if (gPrintLevel && nBadChannel_a2 < 1)
866 {
867 cout << " !!!!! " << nBadChannel_a2 << " !!!!!!!! Bad Calib.: BP= " << busPatchId << " Manu_Id= " << manuId
868 << " Ch.= " << channelId << ":";
869 cout << " a2 = " << par[2] << " out of limit : [" << goodA2Min << "," << goodA2Max
870 << "]" << endl;
c4ac9a95 871
4a3a8cae 872 for (Int_t j = 0; j < nbpf2; ++j)
873 {cout << j << " " << x[j] << " " << xErr[j] << " " << y[j] << " " << yErr[j] << endl;}
874 }
875 }
876 } // FitProceed
877
878 if(FitProceed && p1>0 && p2>0)
879 {
880 nGoodChannel++;
881 sumProbChi2 += prChi2;
882 sumA1 += par[1];
883 gain=1./(par[1]*capa);
884 sumProbChi2P2 += prChi2P2;
885 sumA2 += par[2];
886 }
887 else // bad calibration
888 {
889 nBadChannel++;
890 Q=0;
891 if(gPrintLevel==2)fprintf(pfilef,"%4i %5i %2i %7.4f %10.3e %4i %2x\n",busPatchId,manuId,channelId,par[1],par[2],threshold,Q);
892 par[1]=0.5; a1=0.5; p1=0;
893 par[2]=0.; a2=0.; p2=0;
894 threshold=gkADCMax;
895
896 char bpmanuname[256];
897 ErrorCounter* uncalcounter;
898
899 sprintf(bpmanuname,"bp%dmanu%d",busPatchId,manuId);
900 if (!(uncalcounter = (ErrorCounter*)gUncalBuspatchManuTable->FindObject(bpmanuname)))
901 {
902 // New buspatch_manu name
903 uncalcounter= new ErrorCounter (busPatchId,manuId);
904 uncalcounter->SetName(bpmanuname);
905 gUncalBuspatchManuTable->Add(uncalcounter);
906 }
907 else
908 {
909 // Existing buspatch_manu name
910 uncalcounter->Increment();
911 }
912 // uncalcounter->Print_uncal()
913 uncalcountertotal ++;
914 }
915
916 if(gPlotLevel){
917 // if(Q==0 and nplot < 100)
918 // if(p1>1 && p2==0 and nplot < 100)
919 if(p1>1 && p2>1 and nplot < 100)
920 // if(p1>=1 and p1<=2 and nplot < 100)
921 {
922 nplot++;
923 // cout << " nplot = " << nplot << endl;
924 TF1 *f2Calib = new TF1("f2Calib",funcCalib,0.,gkADCMax,NFITPARAMS);
c4ac9a95 925
4a3a8cae 926 graphErr = new TGraphErrors(nEntries,pedMean,injCharge,pedSigma,injChargeErr);
c4ac9a95 927
4a3a8cae 928 sprintf(graphName,"BusPatch_%d_Manu_%d_Ch_%d",busPatchId, manuId,channelId);
c4ac9a95 929
4a3a8cae 930 graphErr->SetTitle(graphName);
931 graphErr->SetMarkerColor(3);
932 graphErr->SetMarkerStyle(12);
933 graphErr->Write(graphName);
c4ac9a95 934
4a3a8cae 935 sprintf(graphName,"f2_BusPatch_%d_Manu_%d_Ch_%d",busPatchId, manuId,channelId);
936 f2Calib->SetTitle(graphName);
937 f2Calib->SetLineColor(4);
938 f2Calib->SetParameters(par);
939 f2Calib->Write(graphName);
53f515d3 940
4a3a8cae 941 delete graphErr;
942 graphErr=0;
943 delete f2Calib;
944 }
945 }
53f515d3 946
c4ac9a95 947
4a3a8cae 948 tg->Fill();
c4ac9a95 949
4a3a8cae 950 if (!flatOutputFile.IsNull())
951 {
952 fprintf(pfilew,"%4i %5i %2i %7.4f %10.3e %4i %2x\n",busPatchId,manuId,channelId,par[1],par[2],threshold,Q);
953 }
c4ac9a95 954
4a3a8cae 955 }
956 nmanu++;
957 if(fmod(nmanu,100)==0)std::cout << " Nb manu = " << nmanu << std::endl;
958 }
959
960 // file outputs for gain
961 if (!flatOutputFile.IsNull()) fclose(pfilew);
962 if(gPrintLevel==2){ fclose(pfilen); fclose(pfilef); fclose(pfilep); }
963
964 tg->Write();
965 histoFile->Close();
966
967 //OutPut
968 if (gPrintLevel)
969 {
970 // print in logfile
971 if (gUncalBuspatchManuTable->GetSize())
972 {
973 gUncalBuspatchManuTable->Sort(); // use compare
974 TIterator* iter = gUncalBuspatchManuTable->MakeIterator();
975 ErrorCounter* uncalcounter;
976 filcouc << "\n List of problematic BusPatch and Manu " << endl;
977 filcouc << " ========================================" << endl;
978 filcouc << " BP Manu Nb Channel " << endl ;
979 filcouc << " ========================================" << endl;
980 while((uncalcounter = (ErrorCounter*) iter->Next()))
981 {
982 filcouc << "\t" << uncalcounter->BusPatch() << "\t " << uncalcounter->ManuId() << "\t\t" << uncalcounter->Events() << endl;
983 }
984 filcouc << " ========================================" << endl;
985
986 filcouc << " Number of bad calibrated Manu = " << gUncalBuspatchManuTable->GetSize() << endl ;
987 filcouc << " Number of bad calibrated channel = " << uncalcountertotal << endl;
988
989 }
c4ac9a95 990
c4ac9a95 991
4a3a8cae 992 filcouc << "\n Nb of channels in raw data = " << nmanu*64 << " (" << nmanu << " Manu)" << endl;
993 filcouc << " Nb of calibrated channel = " << nGoodChannel << " (" << goodA1Min << "<a1<" << goodA1Max
994 << " and " << goodA2Min << "<a2<" << goodA2Max << ") " << endl;
995 filcouc << " Nb of uncalibrated channel = " << nBadChannel << " (" << noFitChannel << " unfitted)" << endl;
c4ac9a95 996
4a3a8cae 997 cout << "\n Nb of channels in raw data = " << nmanu*64 << " (" << nmanu << " Manu)" << endl;
998 cout << " Nb of calibrated channel = " << nGoodChannel << " (" << goodA1Min << "<a1<" << goodA1Max
999 << " and " << goodA2Min << "<a2<" << goodA2Max << ") " << endl;
1000 cout << " Nb of uncalibrated channel = " << nBadChannel << " (" << noFitChannel << " unfitted)" << endl;
c4ac9a95 1001
4a3a8cae 1002 Double_t meanA1 = sumA1/(nGoodChannel);
1003 Double_t meanProbChi2 = sumProbChi2/(nGoodChannel);
1004 Double_t meanA2 = sumA2/(nGoodChannel);
1005 Double_t meanProbChi2P2 = sumProbChi2P2/(nGoodChannel);
c4ac9a95 1006
4a3a8cae 1007 Double_t capaManu = 0.2; // pF
1008 filcouc << "\n linear fit : <a1> = " << meanA1 << "\t <gain> = " << 1./(meanA1*capaManu)
1009 << " mV/fC (capa= " << capaManu << " pF)" << endl;
1010 filcouc << " Prob(chi2)> = " << meanProbChi2 << endl;
1011 filcouc << "\n parabolic fit: <a2> = " << meanA2 << endl;
1012 filcouc << " Prob(chi2)> = " << meanProbChi2P2 << "\n" << endl;
c4ac9a95 1013
4a3a8cae 1014 cout << "\n <gain> = " << 1./(meanA1*capaManu)
1015 << " mV/fC (capa= " << capaManu << " pF)"
1016 << " Prob(chi2)> = " << meanProbChi2 << endl;
1017 }
53f515d3 1018
4a3a8cae 1019 filcouc.close();
53f515d3 1020
4a3a8cae 1021 return ;
53f515d3 1022
3f637cfc 1023}
1024
1025//*************************************************************//
1026
1027// main routine
1028int main(Int_t argc, Char_t **argv)
1029{
fd99693f 1030
4a3a8cae 1031 // needed for streamer application
1032 gROOT->GetPluginManager()->AddHandler("TVirtualStreamerInfo",
1033 "*",
1034 "TStreamerInfo",
1035 "RIO",
1036 "TStreamerInfo()");
fd99693f 1037
4a3a8cae 1038 TFitter *minuitFit = new TFitter(NFITPARAMS);
1039 TVirtualFitter::SetFitter(minuitFit);
81028269 1040
4a3a8cae 1041 // ofstream filcout;
90293ba6 1042
4a3a8cae 1043 Int_t skipEvents = 0;
1044 Int_t maxEvents = 1000000;
1045 Int_t MaxDateEvents = 1000000;
1046 Int_t injCharge = 0;
1047 Char_t inputFile[256]="";
1048 Char_t inputFile1[256]="";
1049 Char_t inputFile2[256]="";
666b5407 1050
4a3a8cae 1051 Int_t gGlitchErrors= 0;
1052 Int_t gParityErrors= 0;
1053 Int_t gPaddingErrors= 0;
1054 Int_t recoverParityErrors = 1;
c4ac9a95 1055
4a3a8cae 1056 TString fesOutputFile;
3f637cfc 1057
4a3a8cae 1058 // option handler
3f637cfc 1059
4a3a8cae 1060 // decode the input line
1061 for (Int_t i = 1; i < argc; i++) // argument 0 is the executable name
1062 {
1063 Char_t* arg;
c4ac9a95 1064
4a3a8cae 1065 arg = argv[i];
1066 if (arg[0] != '-') continue;
1067 switch (arg[1])
53f515d3 1068 {
4a3a8cae 1069 case 'f' :
1070 i++;
1071 sprintf(inputFile1,argv[i]);
1072 break;
1073 case 'z' :
1074 i++;
1075 sprintf(inputFile2,argv[i]);
1076 break;
1077 case 'a' :
1078 i++;
1079 flatOutputFile = argv[i];
1080 break;
1081 case 'b' :
1082 i++;
1083 sprintf(gOutFolder,argv[i]);
1084 break;
1085 case 'c' :
1086 i++;
1087 gFES=atoi(argv[i]);
1088 break;
1089 case 'd' :
1090 i++;
1091 gPrintLevel=atoi(argv[i]);
1092 break;
1093 case 'e' :
1094 i++;
1095 gCommand = argv[i];
1096 break;
1097 case 'g' :
1098 i++;
1099 gPlotLevel=atoi(argv[i]);
1100 break;
1101 case 'i' :
1102 i++;
1103 nbpf1=atoi(argv[i]);
1104 break;
1105 case 's' :
1106 i++;
1107 skipEvents=atoi(argv[i]);
1108 break;
1109 case 'l' :
1110 i++;
1111 injCharge=atoi(argv[i]);
1112 break;
1113 case 'm' :
1114 i++;
1115 sscanf(argv[i],"%d",&MaxDateEvents);
1116 break;
1117 case 'n' :
1118 i++;
1119 sscanf(argv[i],"%d",&maxEvents);
1120 break;
1121 case 'r' :
1122 i++;
1123 sprintf(gHistoFileName_gain,argv[i]);
1124 break;
1125 case 'p' :
1126 i++;
1127 sscanf(argv[i],"%d",&recoverParityErrors);
1128 break;
1129 case 'h' :
1130 i++;
1131 printf("\n******************* %s usage **********************",argv[0]);
1132 printf("\n%s -options, the available options are :",argv[0]);
1133 printf("\n-h help (this screen)");
1134 printf("\n");
1135 printf("\n Input");
1136 printf("\n-f <raw data file> (default = %s)",inputFile1);
1137 printf("\n-z <raw data file2> (default = %s)",inputFile2);
1138 printf("\n");
1139 printf("\n Output");
1140 printf("\n-a <Flat ASCII file> (default = %s)",flatOutputFile.Data());
1141 printf("\n");
1142 printf("\n Options");
1143 printf("\n-b <output directory> (default = %s)",gOutFolder);
1144 printf("\n-c <FES switch> (default = %d)",gFES);
1145 printf("\n-d <print level> (default = %d)",gPrintLevel);
1146 printf("\n-g <plot level> (default = %d)",gPlotLevel);
1147 printf("\n-i <nb linear points> (default = %d)",nbpf1);
1148 printf("\n-l <DAC level> (default = %d)",injCharge);
1149 printf("\n-m <max date events> (default = %d)",MaxDateEvents);
1150 printf("\n-s <skip events> (default = %d)",skipEvents);
1151 printf("\n-n <max events> (default = %d)",maxEvents);
1152 printf("\n-r root file data for gain (default = %s)",gHistoFileName_gain);
1153 printf("\n-e <execute ped/gain> (default = %s)",gCommand.Data());
1154 printf("\n-e <gain create> make gain & create a new root file");
1155 printf("\n-e <gain> make gain & update root file");
1156 printf("\n-e <gain compute> make gain & compute gains");
1157 printf("\n-p <Recover parity errors> (default = %d)",recoverParityErrors);
1158
1159 printf("\n\n");
1160 exit(-1);
1161 default :
1162 printf("%s : bad argument %s (please check %s -h)\n",argv[0],argv[i],argv[0]);
1163 argc = 2; exit(-1); // exit if error
1164 } // end of switch
1165 } // end of for i
1166
1167 // set gCommand to lower case
1168 gCommand.ToLower();
1169
1170
1171 // decoding the events
1172
1173 Int_t status=0;
1174 void* event;
1175
1176 gPedMeanHisto = 0x0;
1177 gPedSigmaHisto = 0x0;
1178
1179 TStopwatch timers;
1180
1181 timers.Start(kTRUE);
1182
1183 UShort_t manuId;
1184 UChar_t channelId;
1185 UShort_t charge;
1186 TString key("MUONTRKda :");
1187
1188 AliMUONRawStreamTrackerHP* rawStream = 0;
1189
1190 int nrawdata=1;
1191 if (strlen(inputFile2) > 0) nrawdata=2; // 2 raw data files to read (offline case)
1192
1193 if (gCommand.CompareTo("comp") != 0)
1194 {
1195
1196 for(int iraw = 0; iraw < nrawdata; ++iraw)
1197 {
1198 if(iraw==0) strcpy(inputFile,inputFile1);
1199 if(iraw==1) strcpy(inputFile,inputFile2);
1200
1201 status = monitorSetDataSource(inputFile);
1202
1203 if (status) {
1204 cerr << "ERROR : monitorSetDataSource status (hex) = " << hex << status
1205 << " " << monitorDecodeError(status) << endl;
1206 return -1;
1207 }
1208 status = monitorDeclareMp("MUON Tracking monitoring");
1209 if (status) {
1210 cerr << "ERROR : monitorDeclareMp status (hex) = " << hex << status
1211 << " " << monitorDecodeError(status) << endl;
1212 return -1;
1213 }
53f515d3 1214
c4ac9a95 1215
4a3a8cae 1216 cout << "\nMUONTRKda : Reading data from file " << inputFile << " gNEvents = " << gNEvents << endl;
c4ac9a95 1217
4a3a8cae 1218 while(1)
1219 {
1220 if (gNDateEvents >= MaxDateEvents) break;
1221 if (gNEvents >= maxEvents) break;
1222 if (gNDateEvents>0 && gNDateEvents % 100 == 0)
1223 cout<<"Cumulated: DATE events = " << gNDateEvents << " Used events = " << gNEvents << endl;
c4ac9a95 1224
4a3a8cae 1225 // check shutdown condition
1226 if (daqDA_checkShutdown())
1227 break;
c4ac9a95 1228
4a3a8cae 1229 // Skip Events if needed
1230 while (skipEvents) {
1231 status = monitorGetEventDynamic(&event);
1232 skipEvents--;
1233 }
c4ac9a95 1234
4a3a8cae 1235 // starts reading
1236 status = monitorGetEventDynamic(&event);
1237 if (status < 0) {
1238 cout<<"EOF found"<<endl;
1239 break;
1240 }
c4ac9a95 1241
4a3a8cae 1242 // decoding rawdata headers
1243 AliRawReader *rawReader = new AliRawReaderDate(event);
c4ac9a95 1244
4a3a8cae 1245 Int_t eventType = rawReader->GetType();
1246 gRunNumber = rawReader->GetRunNumber();
c4ac9a95 1247
4a3a8cae 1248 // Output log file initialisations
c4ac9a95 1249
4a3a8cae 1250 if(gNDateEvents==0)
1251 {
1252 if (gCommand.CompareTo("ped") == 0){
1253 sprintf(flatFile,"%s/MUONTRKda_ped_%d.log",gOutFolder,gRunNumber);
1254 logOutputFile=flatFile;
1255
1256 filcout.open(logOutputFile.Data());
1257 filcout<<"//=================================================" << endl;
1258 filcout<<"// MUONTRKda for Pedestal run = " << gRunNumber << endl;
1259 cout<<"\n ******** MUONTRKda for Pedestal run = " << gRunNumber << "\n" << endl;
1260 }
c4ac9a95 1261
4a3a8cae 1262 if (gCommand.Contains("gain")){
1263 sprintf(flatFile,"%s/%s_%d_DAC_%d.log",gOutFolder,filenam,gRunNumber,injCharge);
1264 logOutputFile=flatFile;
1265
1266 filcout.open(logOutputFile.Data());
1267 filcout<<"//=================================================" << endl;
1268 filcout<<"// MUONTRKda for Gain run = " << gRunNumber << " (DAC=" << injCharge << ")" << endl;
1269 cout<<"\n ******** MUONTRKda for Gain run = " << gRunNumber << " (DAC=" << injCharge << ")\n" << endl;
1270 }
c4ac9a95 1271
4a3a8cae 1272 filcout<<"//=================================================" << endl;
1273 filcout<<"// * Date : " << date.AsString("l") << "\n" << endl;
1274 cout<<" * Date : " << date.AsString("l") << "\n" << endl;
c4ac9a95 1275
4a3a8cae 1276 }
c4ac9a95 1277
4a3a8cae 1278 gNDateEvents++;
c4ac9a95 1279
4a3a8cae 1280 if (eventType != PHYSICS_EVENT)
1281 continue; // for the moment
c4ac9a95 1282
4a3a8cae 1283 // decoding MUON payload
1284 rawStream = new AliMUONRawStreamTrackerHP(rawReader);
1285 rawStream->DisableWarnings();
1286 rawStream->EnabbleErrorLogger();
c4ac9a95 1287
4a3a8cae 1288 // First lopp over DDL's to find good events
1289 // Error counters per event (counters in the decoding lib are for each DDL)
1290 Bool_t eventIsErrorMessage = kFALSE;
1291 int eventGlitchErrors = 0;
1292 int eventParityErrors = 0;
1293 int eventPaddingErrors = 0;
1294 rawStream->First();
1295 do
1296 {
1297 if (rawStream->IsErrorMessage()) eventIsErrorMessage = kTRUE;
1298 eventGlitchErrors += rawStream->GetGlitchErrors();
1299 eventParityErrors += rawStream->GetParityErrors();
1300 eventPaddingErrors += rawStream->GetPaddingErrors();
1301 } while(rawStream->NextDDL());
1302
1303 AliMUONRawStreamTrackerHP::AliBusPatch* busPatch;
1304 if (!eventIsErrorMessage)
1305 {
1306 // Good events (no error) -> compute pedestal for all channels
1307 rawStream->First();
1308 while( (busPatch = (AliMUONRawStreamTrackerHP::AliBusPatch*) rawStream->Next()))
1309 {
1310 for(int i = 0; i < busPatch->GetLength(); ++i)
1311 {
1312 if (gNEvents == 0) gNChannel++;
1313 busPatch->GetData(i, manuId, channelId, charge);
1314 MakePed(busPatch->GetBusPatchId(), (Int_t)manuId, (Int_t)channelId, (Int_t)charge);
1315 }
1316 }
1317 gNEvents++;
1318 }
1319 else
1320 {
1321 // Events with errors
1322 if (recoverParityErrors && eventParityErrors && !eventGlitchErrors&& !eventPaddingErrors)
1323 {
1324 // Recover parity errors -> compute pedestal for all good buspatches
1325 if ( TEST_SYSTEM_ATTRIBUTE( rawReader->GetAttributes(),
1326 ATTR_ORBIT_BC ))
c4ac9a95 1327 {
4a3a8cae 1328 filcout <<"Event recovered -> Period:"<<EVENT_ID_GET_PERIOD( rawReader->GetEventId() )
1329 <<" Orbit:"<<EVENT_ID_GET_ORBIT( rawReader->GetEventId() )
1330 <<" BunchCrossing:"<<EVENT_ID_GET_BUNCH_CROSSING( rawReader->GetEventId() )<<endl;
1331 }
1332 else
c4ac9a95 1333 {
4a3a8cae 1334 filcout <<"Event recovered -> nbInRun:"<<EVENT_ID_GET_NB_IN_RUN( rawReader->GetEventId() )
1335 <<" burstNb:"<<EVENT_ID_GET_BURST_NB( rawReader->GetEventId() )
1336 <<" nbInBurst:"<<EVENT_ID_GET_NB_IN_BURST( rawReader->GetEventId() )<<endl;
c4ac9a95 1337 }
4a3a8cae 1338 rawStream->First();
1339 while( (busPatch = (AliMUONRawStreamTrackerHP::AliBusPatch*) rawStream->Next()))
c4ac9a95 1340 {
4a3a8cae 1341 // Check the buspatch -> if error not use it in the pedestal calculation
1342 int errorCount = 0;
1343 for(int i = 0; i < busPatch->GetLength(); ++i)
1344 {
1345 if (!busPatch->IsParityOk(i)) errorCount++;
1346 }
1347 if (!errorCount)
1348 {
1349 // Good buspatch
1350 for(int i = 0; i < busPatch->GetLength(); ++i)
c4ac9a95 1351 {
4a3a8cae 1352 if (gNEvents == 0) gNChannel++;
1353 busPatch->GetData(i, manuId, channelId, charge);
1354 // if (busPatch->GetBusPatchId()==1719 && manuId == 1 && channelId == 0) cout <<"Recovered charge "<<charge<<endl;
1355 MakePed(busPatch->GetBusPatchId(), (Int_t)manuId, (Int_t)channelId, (Int_t)charge);
1356 }
1357 }
1358 else
1359 {
1360 char bpname[256];
1361 ErrorCounter* errorCounter;
1362 // Bad buspatch -> not used (just print)
1363 filcout<<"bpId "<<busPatch->GetBusPatchId()<<" words "<<busPatch->GetLength()
1364 <<" parity errors "<<errorCount<<endl;
1365 // Number of events where this buspatch is missing
1366 sprintf(bpname,"bp%d",busPatch->GetBusPatchId());
1367 if (!(errorCounter = (ErrorCounter*)gErrorBuspatchTable->FindObject(bpname)))
c4ac9a95 1368 {
4a3a8cae 1369 // New buspatch
1370 errorCounter = new ErrorCounter(busPatch->GetBusPatchId());
1371 errorCounter->SetName(bpname);
1372 gErrorBuspatchTable->Add(errorCounter);
1373 }
1374 else
1375 {
1376 // Existing buspatch
1377 errorCounter->Increment();
1378 }
1379 // errorCounter->Print();
1380 } // end of if (!errorCount)
1381 } // end of while( (busPatch = (AliMUONRawStreamTrackerHP ...
1382 gNEvents++;
1383 gNEventsRecovered++;
1384 } //end of if (recoverParityErrors && eventParityErrors && !eventGlitchErrors&& !eventPaddingErrors)
1385 else
1386 {
1387 // Fatal errors reject the event
1388 if ( TEST_SYSTEM_ATTRIBUTE( rawReader->GetAttributes(),
1389 ATTR_ORBIT_BC ))
1390 {
1391 filcout <<"Event rejected -> Period:"<<EVENT_ID_GET_PERIOD( rawReader->GetEventId() )
1392 <<" Orbit:"<<EVENT_ID_GET_ORBIT( rawReader->GetEventId() )
1393 <<" BunchCrossing:"<<EVENT_ID_GET_BUNCH_CROSSING( rawReader->GetEventId() )<<endl;
1394 }
1395 else
1396 {
1397 filcout <<"Event rejected -> nbInRun:"<<EVENT_ID_GET_NB_IN_RUN( rawReader->GetEventId() )
1398 <<" burstNb:"<<EVENT_ID_GET_BURST_NB( rawReader->GetEventId() )
1399 <<" nbInBurst:"<<EVENT_ID_GET_NB_IN_BURST( rawReader->GetEventId() )<<endl;
c4ac9a95 1400
4a3a8cae 1401 }
1402 } // end of if (!rawStream->GetGlitchErrors() && !rawStream->GetPaddingErrors() ...
1403 filcout<<"Number of errors : Glitch "<<eventGlitchErrors
1404 <<" Parity "<<eventParityErrors
1405 <<" Padding "<<eventPaddingErrors<<endl;
1406 filcout<<endl;
1407 } // end of if (!rawStream->IsErrorMessage())
c4ac9a95 1408
4a3a8cae 1409 if (eventGlitchErrors) gGlitchErrors++;
1410 if (eventParityErrors) gParityErrors++;
1411 if (eventPaddingErrors) gPaddingErrors++;
c4ac9a95 1412
4a3a8cae 1413 delete rawReader;
1414 delete rawStream;
c4ac9a95 1415
4a3a8cae 1416 } // while (1)
1417 }
c4ac9a95 1418
c4ac9a95 1419
c4ac9a95 1420
4a3a8cae 1421 if (gCommand.CompareTo("ped") == 0)
1422 {
1423 sprintf(flatFile,"MUONTRKda_ped_%d.ped",gRunNumber);
1424 if(flatOutputFile.IsNull())flatOutputFile=flatFile;
1425 MakePedStore(flatOutputFile);
1426 }
c4ac9a95 1427
4a3a8cae 1428 // option gain -> update root file with pedestal results
1429 // gain + create -> recreate root file
1430 // gain + comp -> update root file and compute gain parameters
c4ac9a95 1431
4a3a8cae 1432 if (gCommand.Contains("gain"))
1433 {
1434 MakePedStoreForGain(injCharge);
1435 }
c4ac9a95 1436
c4ac9a95 1437
4a3a8cae 1438 delete gPedestalStore;
53f515d3 1439
4a3a8cae 1440 delete minuitFit;
1441 TVirtualFitter::SetFitter(0);
1442
1443 timers.Stop();
3f637cfc 1444
4a3a8cae 1445 cout << "\nMUONTRKda : Nb of DATE events = " << gNDateEvents << endl;
1446 cout << "MUONTRKda : Nb of Glitch errors = " << gGlitchErrors << endl;
1447 cout << "MUONTRKda : Nb of Parity errors = " << gParityErrors << endl;
1448 cout << "MUONTRKda : Nb of Padding errors = " << gPaddingErrors << endl;
1449 cout << "MUONTRKda : Nb of events recovered = " << gNEventsRecovered<< endl;
1450 cout << "MUONTRKda : Nb of events without errors = " << gNEvents-gNEventsRecovered<< endl;
1451 cout << "MUONTRKda : Nb of events used = " << gNEvents << endl;
53f515d3 1452
4a3a8cae 1453 filcout << "\nMUONTRKda : Nb of DATE events = " << gNDateEvents << endl;
1454 filcout << "MUONTRKda : Nb of Glitch errors = " << gGlitchErrors << endl;
1455 filcout << "MUONTRKda : Nb of Parity errors = " << gParityErrors << endl;
1456 filcout << "MUONTRKda : Nb of Padding errors = " << gPaddingErrors << endl;
1457 filcout << "MUONTRKda : Nb of events recovered = " << gNEventsRecovered<< endl;
1458 filcout << "MUONTRKda : Nb of events without errors = " << gNEvents-gNEventsRecovered<< endl;
1459 filcout << "MUONTRKda : Nb of events used = " << gNEvents << endl;
53f515d3 1460
4a3a8cae 1461 if (gCommand.CompareTo("ped") == 0)
1462 {
1463 cout << "\nMUONTRKda : Output logfile : " << logOutputFile << endl;
1464 cout << "MUONTRKda : Pedestal Histo file : " << gHistoFileName << endl;
1465 cout << "MUONTRKda : Pedestal file (to SHUTTLE) : " << flatOutputFile << endl;
1466 }
1467 else
81028269 1468 {
4a3a8cae 1469 cout << "\nMUONTRKda : Output logfile : " << logOutputFile << endl;
1470 cout << "MUONTRKda : DAC data (root file) : " << gHistoFileName_gain << endl;
1471 cout << "MUONTRKda : Dummy file (to SHUTTLE) : " << flatOutputFile << endl;
1472 }
53f515d3 1473
4a3a8cae 1474 }
81028269 1475
4a3a8cae 1476 // Compute gain parameters
53f515d3 1477
53f515d3 1478
4a3a8cae 1479 if (gCommand.Contains("comp"))
1480 {
1481 flatOutputFile="";
c4ac9a95 1482
4a3a8cae 1483 MakeGainStore();
c4ac9a95 1484
4a3a8cae 1485 cout << "\nMUONTRKda : Output logfile : " << logOutputFile_comp << endl;
1486 cout << "MUONTRKda : Root Histo. file : " << gRootFileName << endl;
1487 cout << "MUONTRKda : Gain file (to SHUTTLE) : " << flatOutputFile << endl;
1488 }
c4ac9a95 1489
4a3a8cae 1490
1491 if(gFES) // Store IN FES
1492 {
1493 printf("\n ***** STORE FILE in FES ****** \n");
1494
1495 // be sure that env variable DAQDALIB_PATH is set in script file
1496 // gSystem->Setenv("DAQDALIB_PATH", "$DATE_SITE/infoLogger");
1497
1498 if (!flatOutputFile.IsNull())
1499 {
1500 if (gCommand.CompareTo("ped") == 0)
1501 status = daqDA_FES_storeFile(flatOutputFile.Data(),"PEDESTALS");
1502 else
1503 status = daqDA_FES_storeFile(flatOutputFile.Data(),"GAINS");
1504
1505 if (status)
1506 {
1507 printf(" Failed to export file : %d\n",status);
1508 }
1509 else if(gPrintLevel) printf(" %s successfully exported to FES \n",flatOutputFile.Data());
53f515d3 1510 }
4a3a8cae 1511 }
3f637cfc 1512
4a3a8cae 1513 filcout.close();
81028269 1514
4a3a8cae 1515 printf("\nExecution time : R:%7.2fs C:%7.2fs\n", timers.RealTime(), timers.CpuTime());
3f637cfc 1516
4a3a8cae 1517 return status;
c758ba48 1518}