MUON/MUONTRKda.cxx

   1 /*
   2 /* 04/10/07
   3
   4 Version for MUONTRKda MUON tracking
   5 (A. Baldisseri, J.-L. Charvet & Ch. Finck)
   6
   7
   8 Rem:  AliMUON2DMap stores all channels, even those which are not connected
   9       if pedMean == -1, channel not connected to a pad
  10
  11
  12 */
  13 extern "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"
  30 #include "AliMUONVStore.h"
  31 #include "AliMUON2DMap.h"
  32 #include "AliMUONCalibParamND.h"
  33 #include "AliMpIntPair.h"
  34 #include "AliMpConstants.h"
  35 #include "AliRawReaderDate.h"
  36
  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"
  47 #include "TROOT.h"
  48 #include "TPluginManager.h"
  49 #include "TFitter.h"
  50
  51 #define  NFITPARAMS 4
  52
  53 // global variables
  54 const Int_t gkNChannels = AliMpConstants::ManuNofChannels();
  55 const Int_t gkADCMax    = 4095;
  56
  57 AliMUONVStore* gPedestalStore =  new AliMUON2DMap(kFALSE);
  58
  59 Int_t  gNManu       = 0;
  60 Int_t  gNChannel    = 0;
  61 UInt_t gRunNumber   = 0;
  62 Int_t  gNEvents     = 0;
  63 Int_t  gNDateEvents = 0;
  64 Int_t  gPrintLevel  = 1;  // global printout variable
  65 Int_t  gPlotLevel  = 0;  // global plot variable
  66
  67 TH1F*  gPedMeanHisto  = 0x0;
  68 TH1F*  gPedSigmaHisto = 0x0;
  69 Char_t gHistoFileName[256];
  70
  71 // used by makegain
  72 Char_t gHistoFileName_gain[256]="MuonTrkDA_data.root";
  73 Char_t gRootFileName[256]="MuonTrkDA_gain.root";
  74 Char_t filenam[256]="MuonTrkDA_gain.param"; // if gPrintLevel  = 2
  75
  76
  77 TString gCommand("ped");
  78
  79 // funtions
  80
  81
  82 //________________
  83 Double_t funcLin(Double_t *x, Double_t *par)
  84 {
  85   return par[0] + par[1]*x[0];
  86 }
  87
  88 //________________
  89 Double_t funcParabolic(Double_t *x, Double_t *par)
  90 {
  91   return par[0]*x[0]*x[0];
  92 }
  93
  94 //________________
  95 Double_t funcCalib(Double_t *x, Double_t *par)
  96 {
  97   Double_t xLim= par[3];
  98
  99   if(x[0] <= xLim) return par[0] + par[1]*x[0];
 100
 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);
 103 }
 104
 105
 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
 122
 123 //__________
 124 void MakePed(Int_t busPatchId, Int_t manuId, Int_t channelId, Int_t charge)
 125 {
 126
 127     AliMUONVCalibParam* ped =
 128         static_cast<AliMUONVCalibParam*>(gPedestalStore->FindObject(busPatchId, manuId));
 129
 130     if (!ped) {
 131       gNManu++;
 132       ped = new AliMUONCalibParamND(2, gkNChannels,busPatchId, manuId, -1.); // put default wise -1, not connected channel
 133       gPedestalStore->Add(ped);
 134     }
 135
 136     if (gNEvents == 1) {
 137       ped->SetValueAsDouble(channelId, 0, 0.);
 138       ped->SetValueAsDouble(channelId, 1, 0.);
 139     }
 140
 141     Double_t pedMean  = ped->ValueAsDouble(channelId, 0) + charge;
 142     Double_t pedSigma = ped->ValueAsDouble(channelId, 1) + charge*charge;
 143
 144     ped->SetValueAsDouble(channelId, 0, pedMean);
 145     ped->SetValueAsDouble(channelId, 1, pedSigma);
 146
 147 }
 148
 149 //________________
 150 void MakePedStore(TString flatOutputFile = "")
 151 {
 152   TTimeStamp date;
 153   Double_t pedMean;
 154   Double_t pedSigma;
 155   ofstream fileout;
 156   Int_t busPatchId;
 157   Int_t manuId;
 158   Int_t channelId;
 159
 160  // histo
 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");
 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;
 172   gPedMeanHisto = new TH1F(name,title,nx,xmin,xmax);
 173   gPedMeanHisto->SetDirectory(histoFile);
 174
 175   sprintf(name,"pedsigma_allch");
 176   sprintf(title,"Pedestal sigma all channels");
 177   nx = 200;
 178   xmin = 0;
 179   xmax = 50;
 180   gPedSigmaHisto = new TH1F(name,title,nx,xmin,xmax);
 181   gPedSigmaHisto->SetDirectory(histoFile);
 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");
 187   tree->Branch("pedMean",&pedMean,",pedMean/D");
 188   tree->Branch("pedSigma",&pedSigma,",pedSigma/D");
 189
 190   if (!flatOutputFile.IsNull()) {
 191     fileout.open(flatOutputFile.Data());
 192     fileout<<"//===========================================================================" << endl;
 193     fileout<<"//                       Pedestal file calculated by MUONTRKda"<<endl;
 194     fileout<<"//===========================================================================" << endl;
 195     fileout<<"//       * Run           : " << gRunNumber << endl;
 196     fileout<<"//       * Date          : " << date.AsString("l") <<endl;
 197     fileout<<"//       * Statictics    : " << gNEvents << endl;
 198     fileout<<"//       * # of MANUS    : " << gNManu << endl;
 199     fileout<<"//       * # of channels : " << gNChannel << endl;
 200     fileout<<"//"<<endl;
 201     fileout<<"//---------------------------------------------------------------------------" << endl;
 202     fileout<<"//---------------------------------------------------------------------------" << endl;
 203 //     fileout<<"//format : BUS_PATCH MANU_ID CHANNEL MEAN SIGMA"<<endl;
 204     fileout<<"//      BP     MANU     CH.      MEAN    SIGMA"<<endl;
 205     fileout<<"//---------------------------------------------------------------------------" << endl;
 206
 207   }
 208
 209   // iterator over pedestal
 210   TIter next(gPedestalStore->CreateIterator());
 211   AliMUONVCalibParam* ped;
 212
 213   while ( ( ped = dynamic_cast<AliMUONVCalibParam*>(next() ) ) )
 214   {
 215     busPatchId              = ped->ID0();
 216     manuId                  = ped->ID1();
 217
 218     for (channelId = 0; channelId < ped->Size() ; ++channelId) {
 219
 220       pedMean  = ped->ValueAsDouble(channelId, 0);
 221
 222       if (pedMean > 0) { // connected channels
 223
 224         ped->SetValueAsDouble(channelId, 0, pedMean/(Double_t)gNEvents);
 225
 226         pedMean  = ped->ValueAsDouble(channelId, 0);
 227         pedSigma = ped->ValueAsDouble(channelId, 1);
 228
 229         ped->SetValueAsDouble(channelId, 1, TMath::Sqrt(TMath::Abs(pedSigma/(Double_t)gNEvents - pedMean*pedMean)));
 230
 231         pedMean  = ped->ValueAsDouble(channelId, 0) + 0.5 ;
 232 //      pedMean  = ped->ValueAsDouble(channelId, 0) ;
 233         pedSigma = ped->ValueAsDouble(channelId, 1);
 234
 235
 236         if (!flatOutputFile.IsNull()) {
 237           fileout << "\t" << busPatchId << "\t" << manuId <<"\t"<< channelId << "\t"
 238                   << pedMean <<"\t"<< pedSigma << endl;
 239         }
 240
 241         gPedMeanHisto->Fill(pedMean);
 242         gPedSigmaHisto->Fill(pedSigma);
 243
 244         tree->Fill();
 245       }
 246     }
 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 //________________
 280 void MakePedStoreForGain(Int_t injCharge)
 281 {
 282     // store pedestal map in root file
 283
 284 //     Int_t injCharge = 200;
 285
 286     TTree* tree = 0x0;
 287
 288     // compute and store pedestals
 289     MakePedStore();
 290
 291     // store in root file
 292 //     sprintf(gHistoFileName,"mutrkgain.root");
 293
 294     TString mode("UPDATE");
 295
 296     if (gCommand.Contains("cre")) {
 297         mode = "RECREATE";
 298     }
 299     TFile* histoFile = new TFile(gHistoFileName_gain, mode.Data(), "MUON Tracking Gains");
 300
 301     // second argument should be the injected charge, taken from config crocus file
 302     // put also info about run number could be usefull
 303     AliMpIntPair* pair   = new AliMpIntPair(gRunNumber, injCharge);
 304
 305     if (mode.CompareTo("UPDATE") == 0) {
 306       tree = (TTree*)histoFile->Get("t");
 307       tree->SetBranchAddress("run",&pair);
 308       tree->SetBranchAddress("ped",&gPedestalStore);
 309
 310     } else {
 311       tree = new TTree("t","Pedestal tree");
 312       tree->Branch("run", "AliMpIntPair",&pair);
 313       tree->Branch("ped", "AliMUON2DMap",&gPedestalStore);
 314       tree->SetBranchAddress("run",&pair);
 315       tree->SetBranchAddress("ped",&gPedestalStore);
 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 //________________
 327 void MakeGainStore(TString flatOutputFile)
 328 {
 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
 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
 341     Double_t pedMean[11];
 342     Double_t pedSigma[11];
 343     Double_t injCharge[11];
 344     Double_t injChargeErr[11];
 345
 346     ofstream fileout;
 347     TTimeStamp date;
 348
 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");
 358
 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");
 362     }
 363
 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];
 391
 392 //  plot out
 393
 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");
 422
 423     //read back from root file
 424     TTree* tree = (TTree*)histoFile->Get("t");
 425     Int_t nEntries = tree->GetEntries();
 426 //     tree->Branch("a1",&a1, "a1/D");
 427
 428     if(gPrintLevel) cout << " nEntries = " << nEntries << " DAC values \n" << endl;
 429
 430     // read back info
 431     for (Int_t i = 0; i < nEntries; ++i) {
 432       map[i] = 0x0;
 433       run[i] = 0x0;
 434       tree->SetBranchAddress("ped",&map[i]);
 435       tree->SetBranchAddress("run",&run[i]);
 436       tree->GetEvent(i);
 437
 438 //       std::cout << map[i] << " " << run[i] << std::endl;
 439     }
 440
 441     // Q = f(ADC)
 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];
 446
 447     // iterates over the first pedestal run
 448     TIter next(map[0]->CreateIterator());
 449     AliMUONVCalibParam* p;
 450
 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
 460     while ( ( p = dynamic_cast<AliMUONVCalibParam*>(next() ) ) )
 461     {
 462       ped[0]  = p;
 463
 464 //       Int_t busPatchId = p->ID0();
 465 //       Int_t manuId     = p->ID1();
 466       busPatchId = p->ID0();
 467       manuId     = p->ID1();
 468
 469       // read back pedestal from the other runs for the given (bupatch, manu)
 470       for (Int_t i = 1; i < nEntries; ++i) {
 471         ped[i] = static_cast<AliMUONVCalibParam*>(map[i]->FindObject(busPatchId, manuId));
 472       }
 473
 474       // compute for each channel the gain parameters
 475 //       for ( Int_t channelId = 0; channelId < ped[0]->Size() ; ++channelId ) {
 476       for ( channelId = 0; channelId < ped[0]->Size() ; ++channelId ) {
 477
 478         Int_t n = 0;
 479         for (Int_t i = 0; i < nEntries; ++i) {
 480
 481           if (!ped[i]) continue; //shouldn't happen.
 482           pedMean[i]      = ped[i]->ValueAsDouble(channelId, 0);
 483           pedSigma[i]     = ped[i]->ValueAsDouble(channelId, 1);
 484           injCharge[i]    = (Double_t)run[i]->GetSecond();
 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;
 491
 492           if (pedMean[i] < 0) continue; // not connected
 493
 494           if (pedSigma[i] <= 0) pedSigma[i] = 1.; // should not happen.
 495           n++;
 496         }
 497
 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
 556           }
 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
 575         }
 576
 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         }
 673       }
 674       nmanu++;
 675     }
 676
 677     // file outputs for gain
 678     if (!flatOutputFile.IsNull())
 679       fileout.close();
 680
 681     tg->Write();
 682     histoFile->Close();
 683
 684     delete f1;
 685     delete f2;
 686
 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     }
 707 }
 708
 709 //*************************************************************//
 710
 711 // main routine
 712 int main(Int_t argc, Char_t **argv)
 713 {
 714
 715     // needed for streamer application
 716     gROOT->GetPluginManager()->AddHandler("TVirtualStreamerInfo",
 717                                           "*",
 718                                           "TStreamerInfo",
 719                                           "RIO",
 720                                           "TStreamerInfo()");
 721
 722     TFitter *minuitFit = new TFitter(NFITPARAMS);
 723     TVirtualFitter::SetFitter(minuitFit);
 724
 725     Int_t skipEvents = 0;
 726     Int_t maxEvents  = 1000000;
 727     Int_t MaxDateEvents  = 1000000;
 728     Int_t injCharge = 0;
 729     Double_t nSigma = 3;
 730     Int_t threshold = -1;
 731     Char_t inputFile[256];
 732     Char_t flatFile[256];
 733
 734     TString flatOutputFile;
 735     TString crocusOutputFile;
 736     TString crocusConfigFile;
 737
 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++;
 767           gCommand = argv[i];
 768           break;
 769         case 'd' :
 770           i++;
 771           gPrintLevel=atoi(argv[i]);
 772           break;
 773         case 'g' :
 774           i++;
 775           gPlotLevel=atoi(argv[i]);
 776           break;
 777         case 's' :
 778           i++;
 779           skipEvents=atoi(argv[i]);
 780           break;
 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;
 789         case 'n' :
 790           i++;
 791           sscanf(argv[i],"%d",&maxEvents);
 792           break;
 793         case 'p' :
 794           i++;
 795           sscanf(argv[i],"%lf",&nSigma);
 796           break;
 797         case 'r' :
 798           i++;
 799           sprintf(gHistoFileName_gain,argv[i]);
 800           break;
 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");
 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);
 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);
 827           printf("\n-r root file data for gain(default = %s)",gHistoFileName_gain);
 828           printf("\n-t <threshold (-1 = no)>  (default = %d)",threshold);
 829           printf("\n-e <execute ped/gain>     (default = %s)",gCommand.Data());
 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
 842   // set gCommand to lower case
 843   gCommand.ToLower();
 844
 845   // decoding the events
 846
 847   Int_t status;
 848   Int_t nDateEvents = 0;
 849
 850   void* event;
 851
 852   gPedMeanHisto = 0x0;
 853   gPedSigmaHisto = 0x0;
 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
 889   Int_t busPatchId;
 890   UShort_t manuId;
 891   UChar_t channelId;
 892   UShort_t charge;
 893
 894   while(1)
 895   {
 896     if (gNDateEvents >= MaxDateEvents) break;
 897     if (gNEvents >= maxEvents) break;
 898     if (gNEvents && gNEvents % 100 == 0)
 899         cout<<"Cumulated events " << gNEvents << endl;
 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
 918     gNDateEvents++;
 919
 920     // decoding rawdata headers
 921     AliRawReader *rawReader = new AliRawReaderDate(event);
 922
 923     Int_t eventType = rawReader->GetType();
 924     gRunNumber = rawReader->GetRunNumber();
 925
 926     if (eventType != PHYSICS_EVENT)
 927         continue; // for the moment
 928
 929     gNEvents++;
 930
 931     // decoding MUON payload
 932     AliMUONRawStreamTracker* rawStream  = new AliMUONRawStreamTracker(rawReader);
 933
 934     // loops over DDL
 935     rawStream->First();
 936     while( (status = rawStream->Next(busPatchId, manuId, channelId, charge)) ) {
 937
 938       if (gNEvents == 1)
 939           gNChannel++;
 940
 941 //       if (gPrintLevel) printf("manuId: %d, channelId: %d charge: %d\n", manuId,
 942 //                           channelId, charge);
 943
 944       MakePed(busPatchId, (Int_t)manuId, (Int_t)channelId, (Int_t)charge);
 945
 946     } // Next digit
 947
 948     delete rawReader;
 949     delete rawStream;
 950
 951   } // while (1)
 952
 953
 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       }
 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
 971   if (gCommand.Contains("gain"))
 972     MakePedStoreForGain(injCharge);
 973
 974   if (gCommand.Contains("comp"))
 975     MakeGainStore(flatOutputFile);
 976
 977
 978   delete gPedestalStore;
 979
 980   delete minuitFit;
 981   TVirtualFitter::SetFitter(0);
 982
 983   timers.Stop();
 984
 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     }
 998   else
 999     {
1000       cout << "\nMUONTRKda : Histo file generated for gain     : " << gHistoFileName_gain  << endl;
1001       cout << "MUONTRKda : Root file generated               : " << gRootFileName  << endl;
1002     }
1003
1004   cout << "MUONTRKda : Flat ASCII file generated         : " << flatOutputFile << endl;
1005   printf("\nExecution time : R:%7.2fs C:%7.2fs\n", timers.RealTime(), timers.CpuTime());
1006
1007   return status;
1008 }