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