MUON/MUONTRKda.cxx

   1 /*
   2
   3 Version 1 for MUONTRKda MUON tracking
   4 Working version for pedestal
   5 Framework for gain computing
   6 (Ch. Finck)
   7
   8
   9 Rem:  AliMUON2DMap stores all channels, even those which are not connected
  10       if pedMean == -1, channel not connected to a pad
  11
  12
  13 */
  14 extern "C" {
  15 #include <daqDA.h>
  16 }
  17
  18 #include "event.h"
  19 #include "monitor.h"
  20
  21 #include <Riostream.h>
  22 #include <stdio.h>
  23 #include <stdlib.h>
  24
  25 //AliRoot
  26 #include "AliMUONRawStreamTracker.h"
  27 #include "AliMUONDspHeader.h"
  28 #include "AliMUONBlockHeader.h"
  29 #include "AliMUONBusStruct.h"
  30 #include "AliMUONDDLTracker.h"
  31 #include "AliMUONVStore.h"
  32 #include "AliMUON2DMap.h"
  33 #include "AliMUONCalibParamND.h"
  34 #include "AliMpIntPair.h"
  35 #include "AliMpConstants.h"
  36 #include "AliRawReaderDate.h"
  37
  38 //ROOT
  39 #include "TFile.h"
  40 #include "TTree.h"
  41 #include "TH1F.h"
  42 #include "TString.h"
  43 #include "TStopwatch.h"
  44 #include "TMath.h"
  45 #include "TTimeStamp.h"
  46 #include "TGraphErrors.h"
  47 #include "TF1.h"
  48 #include "TROOT.h"
  49 #include "TPluginManager.h"
  50
  51
  52 // global variables
  53 AliMUONVStore* pedestalStore =  new AliMUON2DMap(kFALSE);
  54 const Int_t kNchannels = AliMpConstants::ManuNofChannels();
  55 Int_t nManu = 0;
  56 Int_t nChannel = 0;
  57 UInt_t runNumber = 0;
  58 Int_t nEvents = 0;
  59 TH1F* pedMeanHisto = 0x0;
  60 TH1F* pedSigmaHisto = 0x0;
  61 Char_t histoFileName[256];
  62 TString command("ped");
  63
  64 // funtions
  65
  66 //________________
  67 Double_t fitFunc(Double_t *x, Double_t *par)
  68 {
  69     //fit function for gains
  70
  71     Double_t xx = x[0];
  72     Double_t f;
  73
  74     if (xx < par[3])
  75         f = par[0] + par[1]*xx;
  76     else
  77         f= par[0] + par[1]*xx + par[2]*xx*xx;
  78
  79     return f;
  80 }
  81
  82
  83
  84 //__________
  85 void MakePed(Int_t busPatchId, Int_t manuId, Int_t channelId, Int_t charge)
  86 {
  87
  88     AliMUONVCalibParam* ped =
  89         static_cast<AliMUONVCalibParam*>(pedestalStore->FindObject(busPatchId, manuId));
  90
  91     if (!ped) {
  92       nManu++;
  93       ped = new AliMUONCalibParamND(2, kNchannels,busPatchId, manuId, -1.); // put default wise -1, not connected channel
  94       pedestalStore->Add(ped);
  95     }
  96
  97     if (nEvents == 1) {
  98       ped->SetValueAsDouble(channelId, 0, 0.);
  99       ped->SetValueAsDouble(channelId, 1, 0.);
 100     }
 101
 102     Double_t pedMean  = ped->ValueAsDouble(channelId, 0) + charge;
 103     Double_t pedSigma = ped->ValueAsDouble(channelId, 1) + charge*charge;
 104
 105     ped->SetValueAsDouble(channelId, 0, pedMean);
 106     ped->SetValueAsDouble(channelId, 1, pedSigma);
 107
 108 }
 109
 110 //________________
 111 void MakePedStore(TString flatOutputFile = "")
 112 {
 113   TTimeStamp date;
 114   Double_t pedMean;
 115   Double_t pedSigma;
 116   ofstream fileout;
 117   Int_t busPatchId;
 118   Int_t manuId;
 119   Int_t channelId;
 120
 121  // histo
 122   sprintf(histoFileName,"mutrkped-%d.root",runNumber);
 123   TFile*  histoFile = new TFile(histoFileName,"RECREATE","MUON Tracking pedestals");
 124
 125   Char_t name[255];
 126   Char_t title[255];
 127   sprintf(name,"pedmean_allch");
 128   sprintf(title,"Pedestal mean all channels");
 129   Int_t nx = 1000;
 130   Int_t xmin = 0;
 131   Int_t xmax = 1000;
 132   pedMeanHisto = new TH1F(name,title,nx,xmin,xmax);
 133   pedMeanHisto->SetDirectory(histoFile);
 134
 135   sprintf(name,"pedsigma_allch");
 136   sprintf(title,"Pedestal sigma all channels");
 137   nx = 200;
 138   xmin = 0;
 139   xmax = 50;
 140   pedSigmaHisto = new TH1F(name,title,nx,xmin,xmax);
 141   pedSigmaHisto->SetDirectory(histoFile);
 142
 143   TTree* tree = new TTree("t","Pedestal tree");
 144   tree->Branch("bp",&busPatchId,"bp/I");
 145   tree->Branch("manu",&manuId,",manu/I");
 146   tree->Branch("channel",&channelId,",channel/I");
 147
 148   if (!flatOutputFile.IsNull()) {
 149     fileout.open(flatOutputFile.Data());
 150     fileout<<"//===========================================================================" << endl;
 151     fileout<<"//                       Pedestal file calculated by MUONTRKda"<<endl;
 152     fileout<<"//===========================================================================" << endl;
 153     fileout<<"//       * Run           : " << runNumber << endl;
 154     fileout<<"//       * Date          : " << date.AsString("l") <<endl;
 155     fileout<<"//       * Statictics    : " << nEvents << endl;
 156     fileout<<"//       * # of MANUS    : " << nManu << endl;
 157     fileout<<"//       * # of channels : " << nChannel << endl;
 158     fileout<<"//"<<endl;
 159     fileout<<"//---------------------------------------------------------------------------" << endl;
 160     fileout<<"//---------------------------------------------------------------------------" << endl;
 161     fileout<<"//format : BUS_PATCH MANU_ID CHANNEL MEAN SIGMA"<<endl;
 162     fileout<<"//---------------------------------------------------------------------------" << endl;
 163
 164   }
 165
 166   // iterator over pedestal
 167   TIter next(pedestalStore->CreateIterator());
 168   AliMUONVCalibParam* ped;
 169
 170   while ( ( ped = dynamic_cast<AliMUONVCalibParam*>(next() ) ) )
 171   {
 172     busPatchId              = ped->ID0();
 173     manuId                  = ped->ID1();
 174
 175     for (channelId = 0; channelId < ped->Size() ; ++channelId) {
 176
 177       pedMean  = ped->ValueAsDouble(channelId, 0);
 178
 179       if (pedMean > 0) { // connected channels
 180
 181         ped->SetValueAsDouble(channelId, 0, pedMean/(Double_t)nEvents);
 182
 183         pedMean  = ped->ValueAsDouble(channelId, 0);
 184         pedSigma = ped->ValueAsDouble(channelId, 1);
 185
 186         ped->SetValueAsDouble(channelId, 1, TMath::Sqrt(TMath::Abs(pedSigma/(Double_t)nEvents - pedMean*pedMean)));
 187
 188         pedMean  = ped->ValueAsDouble(channelId, 0) + 0.5 ;
 189         pedSigma = ped->ValueAsDouble(channelId, 1);
 190
 191
 192         if (!flatOutputFile.IsNull()) {
 193           fileout << "\t" << busPatchId << "\t" << manuId <<"\t"<< channelId << "\t"
 194                   << pedMean <<"\t"<< pedSigma << endl;
 195         }
 196
 197         pedMeanHisto->Fill(pedMean);
 198         pedSigmaHisto->Fill(pedSigma);
 199
 200         tree->Fill();
 201       }
 202     }
 203   }
 204
 205   // file outputs
 206   if (!flatOutputFile.IsNull())
 207       fileout.close();
 208
 209   histoFile->Write();
 210   histoFile->Close();
 211
 212 //  delete tree;
 213
 214
 215 // not usable need an update of DATE ->5.28, but it compiles
 216 // uncomment when DATE version ok.
 217 // setenv DATE_FES_PATH
 218 // setenv DATE_RUN_NUMBER
 219 // setenv DATE_ROLE_NAME
 220 // setenv DATE_DETECTOR_CODE
 221
 222 //   if (!flatOutputFile.IsNull()) {
 223
 224 //     flatOutputFile.Prepend("./");
 225
 226 //     status = daqDA_FES_storeFile(flatOutputFile.Data(),"MUONTRKda_results");
 227 //     if (status) {
 228 //       printf("Failed to export file : %d\n",status);
 229 //       return -1;
 230 //     }
 231 //  }
 232
 233 }
 234
 235 //________________
 236 void MakePedStoreForGain()
 237 {
 238     // store pedestal map in root file
 239
 240     Int_t injCharge = 200;
 241
 242     TTree* tree = 0x0;
 243
 244     // compute and store pedestals
 245     MakePedStore();
 246
 247     // store in root file
 248     sprintf(histoFileName,"mutrkgain.root");
 249
 250     TString mode("UPDATE");
 251
 252     if (command.Contains("cre")) {
 253         mode = "RECREATE";
 254     }
 255     TFile* histoFile = new TFile(histoFileName, mode.Data(), "MUON Tracking Gains");
 256
 257     // second argument should be the injected charge, taken from config crocus file
 258     // put also info about run number could be usefull
 259     AliMpIntPair* pair   = new AliMpIntPair(runNumber, injCharge);
 260
 261     if (mode.CompareTo("UPDATE") == 0) {
 262       tree = (TTree*)histoFile->Get("t");
 263       tree->SetBranchAddress("run",&pair);
 264       tree->SetBranchAddress("ped",&pedestalStore);
 265
 266     } else {
 267       tree = new TTree("t","Pedestal tree");
 268       tree->Branch("run", "AliMpIntPair",&pair);
 269       tree->Branch("ped", "AliMUON2DMap",&pedestalStore);
 270       tree->SetBranchAddress("run",&pair);
 271       tree->SetBranchAddress("ped",&pedestalStore);
 272
 273     }
 274
 275     tree->Fill();
 276     tree->Write("t", TObject::kOverwrite); // overwrite the tree
 277     histoFile->Close();
 278
 279     delete pair;
 280 }
 281
 282 //________________
 283 void MakeGainStore(TString flatOutputFile)
 284 {
 285
 286     // open file mutrkgain.root
 287     // read again the pedestal for the calibration runs (9 runs ?)
 288     // need the injection charge from config file (to be done)
 289     // For each channel make a TGraphErrors (mean, sigma) vs injected charge
 290     // Fit with a polynomial fct
 291     // store the result in a flat file.
 292
 293     Double_t pedMean[10];
 294     Double_t pedSigma[10];
 295     Double_t injCharge[10];
 296     Double_t injChargeErr[10];
 297
 298     ofstream fileout;
 299     TTimeStamp date;
 300
 301     if (!flatOutputFile.IsNull()) {
 302       fileout.open(flatOutputFile.Data());
 303       fileout<<"//===========================================================================" << endl;
 304       fileout<<"//                       Pedestal file calculated by MUONTRKda"<<endl;
 305       fileout<<"//===========================================================================" << endl;
 306       fileout<<"//       * Run           : " << runNumber << endl;
 307       fileout<<"//       * Date          : " << date.AsString("l") <<endl;
 308       fileout<<"//       * Statictics    : " << nEvents << endl;
 309       fileout<<"//       * # of MANUS    : " << nManu << endl;
 310       fileout<<"//       * # of channels : " << nChannel << endl;
 311       fileout<<"//"<<endl;
 312       fileout<<"//---------------------------------------------------------------------------" << endl;
 313       fileout<<"//---------------------------------------------------------------------------" << endl;
 314       fileout<<"//format : BUS_PATCH MANU_ID CHANNEL GAIN0 GAIN1 GAIN2 XLIM CHI2"<<endl;
 315       fileout<<"//---------------------------------------------------------------------------" << endl;
 316
 317     }
 318
 319
 320     sprintf(histoFileName,"mutrkgain.root");
 321     TFile*  histoFile = new TFile(histoFileName);
 322
 323     AliMUON2DMap* map[10];
 324     AliMUONVCalibParam* ped[10];
 325     //   AliMpIntPair* pair;
 326     AliMpIntPair* run[10];
 327
 328     //read back from root file
 329     TTree* tree = (TTree*)histoFile->Get("t");
 330     Int_t nEntries = tree->GetEntries();
 331
 332     // read back info
 333     for (Int_t i = 0; i < nEntries; ++i) {
 334       map[i] = 0x0;
 335       run[i] = 0x0;
 336       tree->SetBranchAddress("ped",&map[i]);
 337       tree->SetBranchAddress("run",&run[i]);
 338       tree->GetEvent(i);
 339     }
 340
 341     // Q = f(ADC)
 342     TF1* func = new TF1("func",fitFunc, 0., 2500., 4);
 343     // TF1* func = new TF1("func","pol1");
 344
 345     // iterates over the first pedestal run
 346     TIter next(map[0]->CreateIterator());
 347     AliMUONVCalibParam* p;
 348
 349     while ( ( p = dynamic_cast<AliMUONVCalibParam*>(next() ) ) )
 350     {
 351       ped[0]  = p;
 352
 353       Int_t busPatchId = p->ID0();
 354       Int_t manuId     = p->ID1();
 355
 356       // read back pedestal from the other runs for the given (bupatch, manu)
 357       for (Int_t i = 1; i < nEntries; ++i) {
 358         ped[i] = static_cast<AliMUONVCalibParam*>(map[i]->FindObject(busPatchId, manuId));
 359       }
 360
 361       // compute for each channel the gain parameters
 362       for ( Int_t channelId = 0; channelId < ped[0]->Size() ; ++channelId ) {
 363
 364         Int_t n = 0;
 365         for (Int_t i = 0; i < nEntries; ++i) {
 366
 367           if (!ped[i]) continue; //shouldn't happen.
 368           pedMean[i]      = ped[i]->ValueAsDouble(channelId, 0);
 369           pedSigma[i]     = ped[i]->ValueAsDouble(channelId, 1);
 370           injCharge[i]    = (Double_t)run[i]->GetSecond();
 371           injChargeErr[i] = 1.;
 372
 373           if (pedMean[i] < 0) continue; // not connected
 374
 375           if (pedSigma[i] <= 0) pedSigma[i] = 1.; // should not happen.
 376           n++;
 377         }
 378
 379         if (n > 4) {
 380           // if (n > 1) {
 381           //fit
 382           TGraph *gain = new TGraphErrors(n, pedMean, injCharge, pedSigma, injChargeErr);
 383           //should set some initial parameters
 384           func->SetParameter(0,-300);  // a0
 385           func->SetParameter(1, 1.);   // a1
 386           func->SetParameter(2, 0.00001);// a2
 387           func->SetParameter(3, 1100.); // xlim in ADC
 388
 389           gain->Fit("func","q");
 390           cout  << setw(8) << func->GetParameter(0)  <<"\t"<< setw(8) << func->GetParameter(1) << endl;
 391
 392           if (!flatOutputFile.IsNull()) {
 393             fileout << "\t" << busPatchId << "\t" << manuId <<"\t"<< channelId << "\t"
 394                     << setw(8) << func->GetParameter(0)  <<"\t"<< setw(8) << func->GetParameter(1)
 395                     << setw(8) << func->GetParameter(2)  <<"\t"<< setw(5) << func->GetParameter(3)
 396                     << "\t" << func->GetChisquare() << endl;
 397           }
 398           delete gain;
 399         }
 400
 401       }
 402     }
 403
 404     // file outputs for gain
 405     if (!flatOutputFile.IsNull())
 406         fileout.close();
 407
 408     delete func;
 409
 410 }
 411
 412 //*************************************************************//
 413
 414 // main routine
 415 int main(Int_t argc, Char_t **argv)
 416 {
 417
 418     // needed for streamer application
 419     gROOT->GetPluginManager()->AddHandler("TVirtualStreamerInfo",
 420                                           "*",
 421                                           "TStreamerInfo",
 422                                           "RIO",
 423                                           "TStreamerInfo()");
 424
 425
 426     Int_t printLevel = 0;  // Global variable defined as extern in the others .cxx files
 427     Int_t skipEvents = 0;
 428     Int_t maxEvents  = 1000000;
 429     Double_t nSigma = 3;
 430     Int_t threshold = -1;
 431     Char_t inputFile[256];
 432     TString flatOutputFile;
 433     TString crocusOutputFile;
 434     TString crocusConfigFile;
 435
 436 // option handler
 437
 438    // decode the input line
 439   for (Int_t i = 1; i < argc; i++) // argument 0 is the executable name
 440   {
 441       Char_t* arg;
 442
 443       arg = argv[i];
 444       if (arg[0] != '-') continue;
 445       switch (arg[1])
 446         {
 447         case 'f' :
 448           i++;
 449           sprintf(inputFile,argv[i]);
 450           break;
 451         case 'a' :
 452           i++;
 453           flatOutputFile = argv[i];
 454           break;
 455         case 'o' :
 456           i++;
 457           crocusOutputFile = argv[i];
 458           break;
 459         case 'c' :
 460           i++;
 461           crocusConfigFile = argv[i];
 462           break;
 463         case 'e' :
 464           i++;
 465           command = argv[i];
 466           break;
 467         case 'd' :
 468           i++;
 469           printLevel=atoi(argv[i]);
 470           break;
 471         case 's' :
 472           i++;
 473           skipEvents=atoi(argv[i]);
 474           break;
 475         case 'n' :
 476           i++;
 477           sscanf(argv[i],"%d",&maxEvents);
 478           break;
 479         case 'p' :
 480           i++;
 481           sscanf(argv[i],"%lf",&nSigma);
 482           break;
 483         case 't' :
 484           i++;
 485           sscanf(argv[i],"%d",&threshold);
 486           break;
 487         case 'h' :
 488           i++;
 489           printf("\n******************* %s usage **********************",argv[0]);
 490           printf("\n%s -options, the available options are :",argv[0]);
 491           printf("\n-h help                   (this screen)");
 492           printf("\n");
 493           printf("\n Input");
 494           printf("\n-f <raw data file>        (default = %s)",inputFile);
 495           printf("\n-c <Crocus config. file>  (default = %s)",crocusConfigFile.Data());
 496           printf("\n");
 497           printf("\n Output");
 498           printf("\n-a <Flat ASCII file>      (default = %s)",flatOutputFile.Data());
 499           printf("\n-o <CROUCUS cmd file>     (default = %s)",crocusOutputFile.Data());
 500           printf("\n");
 501           printf("\n Options");
 502           printf("\n-d <print level>          (default = %d)",printLevel);
 503           printf("\n-s <skip events>          (default = %d)",skipEvents);
 504           printf("\n-n <max events>           (default = %d)",maxEvents);
 505           printf("\n-p <n sigmas>             (default = %f)",nSigma);
 506           printf("\n-t <threshold (-1 = no)>  (default = %d)",threshold);
 507           printf("\n-e <execute ped/gain>     (default = %s)",command.Data());
 508           printf("\n-e <gain create>           make gain & create a new root file");
 509           printf("\n-e <gain>                  make gain & update root file");
 510           printf("\n-e <gain compute>          make gain & compute gains");
 511
 512           printf("\n\n");
 513           exit(-1);
 514         default :
 515           printf("%s : bad argument %s (please check %s -h)\n",argv[0],argv[i],argv[0]);
 516           argc = 2; exit(-1); // exit if error
 517         } // end of switch
 518     } // end of for i
 519
 520   // set command to lower case
 521   command.ToLower();
 522
 523   // decoding the events
 524
 525   Int_t status;
 526   Int_t nDateEvents = 0;
 527
 528   void* event;
 529
 530   pedMeanHisto = 0x0;
 531   pedSigmaHisto = 0x0;
 532
 533   TStopwatch timers;
 534
 535   timers.Start(kTRUE);
 536
 537   // once we have a configuration file in db
 538   // copy locally a file from daq detector config db
 539   // The current detector is identified by detector code in variable
 540   // DATE_DETECTOR_CODE. It must be defined.
 541   // If environment variable DAQDA_TEST_DIR is defined, files are copied from DAQDA_TEST_DIR
 542   // instead of the database. The usual environment variables are not needed.
 543   if (!crocusConfigFile.IsNull()) {
 544     status = daqDA_DB_getFile("myconfig", crocusConfigFile.Data());
 545     if (status) {
 546       printf("Failed to get config file : %d\n",status);
 547       return -1;
 548     }
 549   }
 550
 551
 552   status = monitorSetDataSource(inputFile);
 553   if (status) {
 554     cerr << "ERROR : monitorSetDataSource status (hex) = " << hex << status
 555               << " " << monitorDecodeError(status) << endl;
 556     return -1;
 557   }
 558   status = monitorDeclareMp("MUON Tracking monitoring");
 559   if (status) {
 560     cerr << "ERROR : monitorDeclareMp status (hex) = " << hex << status
 561               << " " << monitorDecodeError(status) << endl;
 562     return -1;
 563   }
 564
 565   cout << "MUONTRKda : Reading data from file " << inputFile <<endl;
 566
 567   Int_t busPatchId;
 568   UShort_t manuId;
 569   UChar_t channelId;
 570   UShort_t charge;
 571
 572   while(1)
 573   {
 574     if (nEvents >= maxEvents) break;
 575     if (nEvents && nEvents % 100 == 0)
 576         cout<<"Cumulated events " << nEvents << endl;
 577
 578     // check shutdown condition
 579     if (daqDA_checkShutdown())
 580         break;
 581
 582     // Skip Events if needed
 583     while (skipEvents) {
 584       status = monitorGetEventDynamic(&event);
 585       skipEvents--;
 586     }
 587
 588     // starts reading
 589     status = monitorGetEventDynamic(&event);
 590     if (status < 0)  {
 591       cout<<"EOF found"<<endl;
 592       break;
 593     }
 594
 595     nDateEvents++;
 596
 597     // decoding rawdata headers
 598     AliRawReader *rawReader = new AliRawReaderDate(event);
 599
 600     Int_t eventType = rawReader->GetType();
 601     runNumber = rawReader->GetRunNumber();
 602
 603     if (eventType != PHYSICS_EVENT)
 604         continue; // for the moment
 605
 606     nEvents++;
 607
 608     // decoding MUON payload
 609     AliMUONRawStreamTracker* rawStream  = new AliMUONRawStreamTracker(rawReader);
 610
 611     // loops over DDL
 612     rawStream->First();
 613     while( (status = rawStream->Next(busPatchId, manuId, channelId, charge)) ) {
 614
 615       if (nEvents == 1)
 616           nChannel++;
 617
 618       if (printLevel) printf("manuId: %d, channelId: %d charge: %d\n", manuId,
 619                              channelId, charge);
 620
 621       MakePed(busPatchId, (Int_t)manuId, (Int_t)channelId, (Int_t)charge);
 622
 623     } // Next digit
 624
 625     delete rawReader;
 626     delete rawStream;
 627
 628   } // while (1)
 629
 630
 631   if (command.CompareTo("ped") == 0)
 632       MakePedStore(flatOutputFile);
 633
 634   // option gain -> update root file with pedestal results
 635   // gain + create -> recreate root file
 636   // gain + comp -> update root file and compute gain parameters
 637
 638   if (command.Contains("gain"))
 639       MakePedStoreForGain();
 640
 641   if (command.Contains("comp"))
 642       MakeGainStore(flatOutputFile);
 643
 644
 645   delete pedestalStore;
 646
 647   timers.Stop();
 648
 649   if (!(crocusConfigFile.IsNull()))
 650       cout << "MUONTRKda : CROCUS command file generated : " << crocusOutputFile.Data() << endl;
 651   else
 652       cout << "MUONTRKda : WARNING no CROCUS command file generated" << endl;
 653
 654   cout << "MUONTRKda : Flat ASCII file generated     : " << flatOutputFile << endl;
 655   cout << "MUONTRKda : Histo file generated          : " << histoFileName  << endl;
 656   cout << "MUONTRKda : Nb of DATE events     = "         << nDateEvents    << endl;
 657   cout << "MUONTRKda : Nb of events used     = "         << nEvents        << endl;
 658
 659   printf("Execution time : R:%7.2fs C:%7.2fs\n", timers.RealTime(), timers.CpuTime());
 660
 661   return status;
 662 }