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