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