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