[u/mrichter/AliRoot.git] / MUON / MUONTRKda.cxx

/*

Version 1 for MUONTRKda MUON tracking
Working version for pedestal
Framework for gain computing
(Ch. Finck)


Rem:  AliMUON2DMap stores all channels, even those which are not connected
      if pedMean == -1, channel not connected to a pad  


*/
extern "C" {
#include <daqDA.h>
}

#include "event.h"
#include "monitor.h"

#include <Riostream.h>
#include <stdio.h>
#include <stdlib.h>

//AliRoot
#include "AliMUONRawStreamTracker.h"
#include "AliMUONDspHeader.h"
#include "AliMUONBlockHeader.h"
#include "AliMUONBusStruct.h"
#include "AliMUONDDLTracker.h"
#include "AliMUONVStore.h"
#include "AliMUON2DMap.h"
#include "AliMUONCalibParamND.h"
#include "AliMpIntPair.h"
#include "AliMpConstants.h"
#include "AliRawReaderDate.h"

//ROOT
#include "TFile.h"
#include "TTree.h"
#include "TH1F.h"
#include "TString.h"
#include "TStopwatch.h"
#include "TMath.h"
#include "TTimeStamp.h"
#include "TGraphErrors.h"
#include "TF1.h"
#include "TROOT.h"
#include "TPluginManager.h"
#include "TFitter.h"

#define  NFITPARAMS 4

// global variables
AliMUONVStore* pedestalStore =  new AliMUON2DMap(kFALSE);
const Int_t kNchannels = AliMpConstants::ManuNofChannels();
Int_t nManu = 0;
Int_t nChannel = 0;
UInt_t runNumber = 0;
Int_t nEvents = 0;
TH1F* pedMeanHisto = 0x0;
TH1F* pedSigmaHisto = 0x0;
Char_t histoFileName[256];
TString command("ped");

// funtions

//________________
Double_t fitFunc(Double_t *x, Double_t *par)
{
    //fit function for gains

    Double_t xx = x[0];
    Double_t f;

    if (xx < par[3])
	f = par[0] + par[1]*xx;
    else
	f= par[0] + par[1]*xx + par[2]*xx*xx;

    return f;
}


//__________
void MakePed(Int_t busPatchId, Int_t manuId, Int_t channelId, Int_t charge)
{

    AliMUONVCalibParam* ped = 
	static_cast<AliMUONVCalibParam*>(pedestalStore->FindObject(busPatchId, manuId));

    if (!ped) {
      nManu++;
      ped = new AliMUONCalibParamND(2, kNchannels,busPatchId, manuId, -1.); // put default wise -1, not connected channel
      pedestalStore->Add(ped);	
    }

    if (nEvents == 1) {
      ped->SetValueAsDouble(channelId, 0, 0.);
      ped->SetValueAsDouble(channelId, 1, 0.);
    }

    Double_t pedMean  = ped->ValueAsDouble(channelId, 0) + charge;
    Double_t pedSigma = ped->ValueAsDouble(channelId, 1) + charge*charge;

    ped->SetValueAsDouble(channelId, 0, pedMean);
    ped->SetValueAsDouble(channelId, 1, pedSigma);

}

//________________
void MakePedStore(TString flatOutputFile = "")
{
  TTimeStamp date;
  Double_t pedMean;
  Double_t pedSigma;
  ofstream fileout;
  Int_t busPatchId;
  Int_t manuId;
  Int_t channelId;

 // histo
  sprintf(histoFileName,"mutrkped-%d.root",runNumber);
  TFile*  histoFile = new TFile(histoFileName,"RECREATE","MUON Tracking pedestals");

  Char_t name[255];
  Char_t title[255];
  sprintf(name,"pedmean_allch");
  sprintf(title,"Pedestal mean all channels");
  Int_t nx = 1000;
  Int_t xmin = 0;
  Int_t xmax = 1000; 
  pedMeanHisto = new TH1F(name,title,nx,xmin,xmax);
  pedMeanHisto->SetDirectory(histoFile);

  sprintf(name,"pedsigma_allch");
  sprintf(title,"Pedestal sigma all channels");
  nx = 200;
  xmin = 0;
  xmax = 50; 
  pedSigmaHisto = new TH1F(name,title,nx,xmin,xmax);
  pedSigmaHisto->SetDirectory(histoFile);
    
  TTree* tree = new TTree("t","Pedestal tree");
  tree->Branch("bp",&busPatchId,"bp/I");
  tree->Branch("manu",&manuId,",manu/I");
  tree->Branch("channel",&channelId,",channel/I");

  if (!flatOutputFile.IsNull()) {
    fileout.open(flatOutputFile.Data());
    fileout<<"//===========================================================================" << endl;
    fileout<<"//                       Pedestal file calculated by MUONTRKda"<<endl;
    fileout<<"//===========================================================================" << endl;
    fileout<<"//       * Run           : " << runNumber << endl; 
    fileout<<"//       * Date          : " << date.AsString("l") <<endl;
    fileout<<"//       * Statictics    : " << nEvents << endl;
    fileout<<"//       * # of MANUS    : " << nManu << endl;
    fileout<<"//       * # of channels : " << nChannel << endl;
    fileout<<"//"<<endl;
    fileout<<"//---------------------------------------------------------------------------" << endl;
    fileout<<"//---------------------------------------------------------------------------" << endl;
    fileout<<"//format : BUS_PATCH MANU_ID CHANNEL MEAN SIGMA"<<endl;
    fileout<<"//---------------------------------------------------------------------------" << endl;

  }

  // iterator over pedestal
  TIter next(pedestalStore->CreateIterator());
  AliMUONVCalibParam* ped;
  
  while ( ( ped = dynamic_cast<AliMUONVCalibParam*>(next() ) ) )
  {
    busPatchId              = ped->ID0();
    manuId                  = ped->ID1();

    for (channelId = 0; channelId < ped->Size() ; ++channelId) {

      pedMean  = ped->ValueAsDouble(channelId, 0);

      if (pedMean > 0) { // connected channels

	ped->SetValueAsDouble(channelId, 0, pedMean/(Double_t)nEvents);

	pedMean  = ped->ValueAsDouble(channelId, 0);
	pedSigma = ped->ValueAsDouble(channelId, 1);

	ped->SetValueAsDouble(channelId, 1, TMath::Sqrt(TMath::Abs(pedSigma/(Double_t)nEvents - pedMean*pedMean)));

	pedMean  = ped->ValueAsDouble(channelId, 0) + 0.5 ;
	pedSigma = ped->ValueAsDouble(channelId, 1);


	if (!flatOutputFile.IsNull()) {
	  fileout << "\t" << busPatchId << "\t" << manuId <<"\t"<< channelId << "\t"
		  << pedMean <<"\t"<< pedSigma << endl;
	}

	pedMeanHisto->Fill(pedMean);
	pedSigmaHisto->Fill(pedSigma);

	tree->Fill();
      }
    }
  }

  // file outputs
  if (!flatOutputFile.IsNull()) 
      fileout.close();

  histoFile->Write();
  histoFile->Close();

//  delete tree;


// not usable need an update of DATE ->5.28, but it compiles
// uncomment when DATE version ok.
// setenv DATE_FES_PATH
// setenv DATE_RUN_NUMBER
// setenv DATE_ROLE_NAME
// setenv DATE_DETECTOR_CODE

//   if (!flatOutputFile.IsNull()) {

//     flatOutputFile.Prepend("./");

//     status = daqDA_FES_storeFile(flatOutputFile.Data(),"MUONTRKda_results");
//     if (status) {
//       printf("Failed to export file : %d\n",status);
//       return -1;
//     }
//  }

}

//________________
void MakePedStoreForGain()
{
    // store pedestal map in root file

    Int_t injCharge = 200;

    TTree* tree = 0x0;

    // compute and store pedestals
    MakePedStore();

    // store in root file
    sprintf(histoFileName,"mutrkgain.root");
    
    TString mode("UPDATE");

    if (command.Contains("cre")) {
	mode = "RECREATE";
    }
    TFile* histoFile = new TFile(histoFileName, mode.Data(), "MUON Tracking Gains");

    // second argument should be the injected charge, taken from config crocus file
    // put also info about run number could be usefull
    AliMpIntPair* pair   = new AliMpIntPair(runNumber, injCharge);

    if (mode.CompareTo("UPDATE") == 0) {
      tree = (TTree*)histoFile->Get("t");
      tree->SetBranchAddress("run",&pair);
      tree->SetBranchAddress("ped",&pedestalStore);

    } else {
      tree = new TTree("t","Pedestal tree");
      tree->Branch("run", "AliMpIntPair",&pair);
      tree->Branch("ped", "AliMUON2DMap",&pedestalStore);
      tree->SetBranchAddress("run",&pair);
      tree->SetBranchAddress("ped",&pedestalStore);

    }

    tree->Fill();
    tree->Write("t", TObject::kOverwrite); // overwrite the tree
    histoFile->Close();

    delete pair;
}

//________________
void MakeGainStore(TString flatOutputFile)
{
    
    // open file mutrkgain.root
    // read again the pedestal for the calibration runs (9 runs ?)
    // need the injection charge from config file (to be done)
    // For each channel make a TGraphErrors (mean, sigma) vs injected charge
    // Fit with a polynomial fct
    // store the result in a flat file.

    Double_t pedMean[10];
    Double_t pedSigma[10];
    Double_t injCharge[10];
    Double_t injChargeErr[10];

    ofstream fileout;
    TTimeStamp date;

    if (!flatOutputFile.IsNull()) {
      fileout.open(flatOutputFile.Data());
      fileout<<"//===========================================================================" << endl;
      fileout<<"//                       Pedestal file calculated by MUONTRKda"<<endl;
      fileout<<"//===========================================================================" << endl;
      fileout<<"//       * Run           : " << runNumber << endl; 
      fileout<<"//       * Date          : " << date.AsString("l") <<endl;
      fileout<<"//       * Statictics    : " << nEvents << endl;
      fileout<<"//       * # of MANUS    : " << nManu << endl;
      fileout<<"//       * # of channels : " << nChannel << endl;
      fileout<<"//"<<endl;
      fileout<<"//---------------------------------------------------------------------------" << endl;
      fileout<<"//---------------------------------------------------------------------------" << endl;
      fileout<<"//format : BUS_PATCH MANU_ID CHANNEL GAIN0 GAIN1 GAIN2 XLIM CHI2"<<endl;
      fileout<<"//---------------------------------------------------------------------------" << endl;

    }


    sprintf(histoFileName,"mutrkgain.root");
    TFile*  histoFile = new TFile(histoFileName);

    AliMUON2DMap* map[10];
    AliMUONVCalibParam* ped[10];
    //   AliMpIntPair* pair;
    AliMpIntPair* run[10];

    //read back from root file
    TTree* tree = (TTree*)histoFile->Get("t");
    Int_t nEntries = tree->GetEntries();

    // read back info
    for (Int_t i = 0; i < nEntries; ++i) {
      map[i] = 0x0;
      run[i] = 0x0;
      tree->SetBranchAddress("ped",&map[i]);
      tree->SetBranchAddress("run",&run[i]);
      tree->GetEvent(i);
    }
      
    // Q = f(ADC)
    TF1* func = new TF1("func",fitFunc, 0., 2500., NFITPARAMS);
    // TF1* func = new TF1("func","pol1");

    // iterates over the first pedestal run
    TIter next(map[0]->CreateIterator());
    AliMUONVCalibParam* p;

    while ( ( p = dynamic_cast<AliMUONVCalibParam*>(next() ) ) )
    {
      ped[0]  = p;

      Int_t busPatchId = p->ID0();
      Int_t manuId     = p->ID1();

      // read back pedestal from the other runs for the given (bupatch, manu)
      for (Int_t i = 1; i < nEntries; ++i) {
	ped[i] = static_cast<AliMUONVCalibParam*>(map[i]->FindObject(busPatchId, manuId));
      }

      // compute for each channel the gain parameters
      for ( Int_t channelId = 0; channelId < ped[0]->Size() ; ++channelId ) {

	Int_t n = 0;
	for (Int_t i = 0; i < nEntries; ++i) {

	  if (!ped[i]) continue; //shouldn't happen.
	  pedMean[i]      = ped[i]->ValueAsDouble(channelId, 0);
	  pedSigma[i]     = ped[i]->ValueAsDouble(channelId, 1);
	  injCharge[i]    = (Double_t)run[i]->GetSecond();
	  injChargeErr[i] = 1.;

	  if (pedMean[i] < 0) continue; // not connected

	  if (pedSigma[i] <= 0) pedSigma[i] = 1.; // should not happen.
	  n++;
	}

	if (n > NFITPARAMS) {
	  // if (n > 1) {
	  //fit 
	  TGraph *gain = new TGraphErrors(n, pedMean, injCharge, pedSigma, injChargeErr);
	  //should set some initial parameters
	  func->SetParameter(0,-300);  // a0
	  func->SetParameter(1, 1.);   // a1
	  func->SetParameter(2, 0.00001);// a2
	  func->SetParameter(3, 1100.); // xlim in ADC

	  gain->Fit("func","q");
	  cout	<< setw(8) << func->GetParameter(0)  <<"\t"<< setw(8) << func->GetParameter(1) << endl;

	  if (!flatOutputFile.IsNull()) {
	    fileout << "\t" << busPatchId << "\t" << manuId <<"\t"<< channelId << "\t"
		    << setw(8) << func->GetParameter(0)  <<"\t"<< setw(8) << func->GetParameter(1) 
		    << setw(8) << func->GetParameter(2)  <<"\t"<< setw(5) << func->GetParameter(3) 
		    << "\t" << func->GetChisquare() << endl;
	  }
	  delete gain;
	}

      }
    }

    // file outputs for gain
    if (!flatOutputFile.IsNull()) 
	fileout.close();

    delete func;

}

//*************************************************************//

// main routine
int main(Int_t argc, Char_t **argv) 
{
  
    // needed for streamer application
    gROOT->GetPluginManager()->AddHandler("TVirtualStreamerInfo",
					  "*",
					  "TStreamerInfo",
					  "RIO",
					  "TStreamerInfo()"); 

    TFitter *minuitFit = new TFitter(NFITPARAMS);
    TVirtualFitter::SetFitter(minuitFit);

    Int_t printLevel = 0;  // Global variable defined as extern in the others .cxx files
    Int_t skipEvents = 0;
    Int_t maxEvents  = 1000000;
    Double_t nSigma = 3;
    Int_t threshold = -1;
    Char_t inputFile[256];
    TString flatOutputFile;
    TString crocusOutputFile;
    TString crocusConfigFile;

// option handler

   // decode the input line
  for (Int_t i = 1; i < argc; i++) // argument 0 is the executable name
  {
      Char_t* arg;
      
      arg = argv[i];
      if (arg[0] != '-') continue;
      switch (arg[1])
	{
	case 'f' : 
	  i++;
	  sprintf(inputFile,argv[i]);
	  break;
	case 'a' : 
	  i++;
	  flatOutputFile = argv[i];
	  break;
	case 'o' : 
	  i++;
	  crocusOutputFile = argv[i];
	  break;
	case 'c' : 
	  i++;
	  crocusConfigFile = argv[i];
	  break;
	case 'e' : 
	  i++;
	  command = argv[i];
	  break;
	case 'd' :
	  i++; 
	  printLevel=atoi(argv[i]);
	  break;
	case 's' :
	  i++; 
	  skipEvents=atoi(argv[i]);
	  break;
	case 'n' :
	  i++; 
	  sscanf(argv[i],"%d",&maxEvents);
	  break;
	case 'p' :
	  i++; 
	  sscanf(argv[i],"%lf",&nSigma);
	  break;
	case 't' :
	  i++; 
	  sscanf(argv[i],"%d",&threshold);
	  break;
	case 'h' :
	  i++;
	  printf("\n******************* %s usage **********************",argv[0]);
	  printf("\n%s -options, the available options are :",argv[0]);
	  printf("\n-h help                   (this screen)");
	  printf("\n");
	  printf("\n Input");
	  printf("\n-f <raw data file>        (default = %s)",inputFile); 
	  printf("\n-c <Crocus config. file>  (default = %s)",crocusConfigFile.Data()); 
	  printf("\n");
	  printf("\n Output");
	  printf("\n-a <Flat ASCII file>      (default = %s)",flatOutputFile.Data()); 
	  printf("\n-o <CROUCUS cmd file>     (default = %s)",crocusOutputFile.Data()); 
	  printf("\n");
	  printf("\n Options");
	  printf("\n-d <print level>          (default = %d)",printLevel);
	  printf("\n-s <skip events>          (default = %d)",skipEvents);
	  printf("\n-n <max events>           (default = %d)",maxEvents);
	  printf("\n-p <n sigmas>             (default = %f)",nSigma);
	  printf("\n-t <threshold (-1 = no)>  (default = %d)",threshold);
	  printf("\n-e <execute ped/gain>     (default = %s)",command.Data());
	  printf("\n-e <gain create>           make gain & create a new root file");
	  printf("\n-e <gain>                  make gain & update root file");
	  printf("\n-e <gain compute>          make gain & compute gains");

	  printf("\n\n");
	  exit(-1);
	default :
	  printf("%s : bad argument %s (please check %s -h)\n",argv[0],argv[i],argv[0]);
	  argc = 2; exit(-1); // exit if error
	} // end of switch  
    } // end of for i  

  // set command to lower case
  command.ToLower();

  // decoding the events
  
  Int_t status;
  Int_t nDateEvents = 0;

  void* event;

  pedMeanHisto = 0x0;
  pedSigmaHisto = 0x0;

  TStopwatch timers;

  timers.Start(kTRUE); 

  // once we have a configuration file in db
  // copy locally a file from daq detector config db 
  // The current detector is identified by detector code in variable
  // DATE_DETECTOR_CODE. It must be defined.
  // If environment variable DAQDA_TEST_DIR is defined, files are copied from DAQDA_TEST_DIR
  // instead of the database. The usual environment variables are not needed.
  if (!crocusConfigFile.IsNull()) {
    status = daqDA_DB_getFile("myconfig", crocusConfigFile.Data());
    if (status) {
      printf("Failed to get config file : %d\n",status);
      return -1;
    }
  }


  status = monitorSetDataSource(inputFile);
  if (status) {
    cerr << "ERROR : monitorSetDataSource status (hex) = " << hex << status
	      << " " << monitorDecodeError(status) << endl;
    return -1;
  }
  status = monitorDeclareMp("MUON Tracking monitoring");
  if (status) {
    cerr << "ERROR : monitorDeclareMp status (hex) = " << hex << status
	      << " " << monitorDecodeError(status) << endl;
    return -1;
  }

  cout << "MUONTRKda : Reading data from file " << inputFile <<endl;

  Int_t busPatchId;
  UShort_t manuId;  
  UChar_t channelId;
  UShort_t charge;

  while(1) 
  {
    if (nEvents >= maxEvents) break;
    if (nEvents && nEvents % 100 == 0) 	
	cout<<"Cumulated events " << nEvents << endl;

    // check shutdown condition 
    if (daqDA_checkShutdown()) 
	break;

    // Skip Events if needed
    while (skipEvents) {
      status = monitorGetEventDynamic(&event);
      skipEvents--;
    }

    // starts reading
    status = monitorGetEventDynamic(&event);
    if (status < 0)  {
      cout<<"EOF found"<<endl;
      break;
    }

    nDateEvents++;

    // decoding rawdata headers
    AliRawReader *rawReader = new AliRawReaderDate(event);
 
    Int_t eventType = rawReader->GetType();
    runNumber = rawReader->GetRunNumber();

    if (eventType != PHYSICS_EVENT)
	continue; // for the moment

    nEvents++;

    // decoding MUON payload
    AliMUONRawStreamTracker* rawStream  = new AliMUONRawStreamTracker(rawReader);

    // loops over DDL 
    rawStream->First();
    while( (status = rawStream->Next(busPatchId, manuId, channelId, charge)) ) {
  
      if (nEvents == 1)
	  nChannel++;
      
      if (printLevel) printf("manuId: %d, channelId: %d charge: %d\n", manuId, 
			     channelId, charge);
      
      MakePed(busPatchId, (Int_t)manuId, (Int_t)channelId, (Int_t)charge);
		  
    } // Next digit
     
    delete rawReader;
    delete rawStream;

  } // while (1)


  if (command.CompareTo("ped") == 0)
      MakePedStore(flatOutputFile);

  // option gain -> update root file with pedestal results
  // gain + create -> recreate root file
  // gain + comp -> update root file and compute gain parameters

  if (command.Contains("gain")) 
      MakePedStoreForGain();
  
  if (command.Contains("comp")) 
      MakeGainStore(flatOutputFile);
  

  delete pedestalStore;

  delete minuitFit;
  TVirtualFitter::SetFitter(0);

  timers.Stop();

  if (!(crocusConfigFile.IsNull()))
      cout << "MUONTRKda : CROCUS command file generated : " << crocusOutputFile.Data() << endl;
  else
      cout << "MUONTRKda : WARNING no CROCUS command file generated" << endl;

  cout << "MUONTRKda : Flat ASCII file generated     : " << flatOutputFile << endl;
  cout << "MUONTRKda : Histo file generated          : " << histoFileName  << endl;
  cout << "MUONTRKda : Nb of DATE events     = "         << nDateEvents    << endl;
  cout << "MUONTRKda : Nb of events used     = "         << nEvents        << endl;

  printf("Execution time : R:%7.2fs C:%7.2fs\n", timers.RealTime(), timers.CpuTime());

  return status;
}
Commit	Line	Data
3f637cfc	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
3f637cfc	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"
2bf54af2	31	#include "AliMUONVStore.h"
3f637cfc	32	#include "AliMUON2DMap.h"
0aa03544	33	#include "AliMUONCalibParamND.h"
3f637cfc	34	#include "AliMpIntPair.h"
2bf54af2	35	#include "AliMpConstants.h"
3f637cfc	36	#include "AliRawReaderDate.h"
3f637cfc	37
3f637cfc	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"
fd99693f	48	#include "TROOT.h"
fd99693f	49	#include "TPluginManager.h"
a76f513d	50	#include "TFitter.h"
3f637cfc	51
a76f513d	52	#define NFITPARAMS 4
3f637cfc	53
3f637cfc	54	// global variables
2bf54af2	55	AliMUONVStore* pedestalStore = new AliMUON2DMap(kFALSE);
2bf54af2	56	const Int_t kNchannels = AliMpConstants::ManuNofChannels();
3f637cfc	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
0aa03544	73	Double_t xx = x[0];
3f637cfc	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 =
0aa03544	91	static_cast<AliMUONVCalibParam*>(pedestalStore->FindObject(busPatchId, manuId));
3f637cfc	92
	93	if (!ped) {
	94	nManu++;
0aa03544	95	ped = new AliMUONCalibParamND(2, kNchannels,busPatchId, manuId, -1.); // put default wise -1, not connected channel
0aa03544	96	pedestalStore->Add(ped);
3f637cfc	97	}
	98
	99	if (nEvents == 1) {
0aa03544	100	ped->SetValueAsDouble(channelId, 0, 0.);
0aa03544	101	ped->SetValueAsDouble(channelId, 1, 0.);
3f637cfc	102	}
3f637cfc	103
0aa03544	104	Double_t pedMean = ped->ValueAsDouble(channelId, 0) + charge;
0aa03544	105	Double_t pedSigma = ped->ValueAsDouble(channelId, 1) + charge*charge;
3f637cfc	106
0aa03544	107	ped->SetValueAsDouble(channelId, 0, pedMean);
0aa03544	108	ped->SetValueAsDouble(channelId, 1, pedSigma);
3f637cfc	109
	110	}
	111
	112	//________________
	113	void MakePedStore(TString flatOutputFile = "")
	114	{
	115	TTimeStamp date;
0aa03544	116	Double_t pedMean;
0aa03544	117	Double_t pedSigma;
3f637cfc	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
0aa03544	169	TIter next(pedestalStore->CreateIterator());
2bf54af2	170	AliMUONVCalibParam* ped;
3f637cfc	171
2bf54af2	172	while ( ( ped = dynamic_cast<AliMUONVCalibParam*>(next() ) ) )
3f637cfc	173	{
2bf54af2	174	busPatchId = ped->ID0();
2bf54af2	175	manuId = ped->ID1();
3f637cfc	176
	177	for (channelId = 0; channelId < ped->Size() ; ++channelId) {
	178
0aa03544	179	pedMean = ped->ValueAsDouble(channelId, 0);
3f637cfc	180
	181	if (pedMean > 0) { // connected channels
	182
0aa03544	183	ped->SetValueAsDouble(channelId, 0, pedMean/(Double_t)nEvents);
3f637cfc	184
0aa03544	185	pedMean = ped->ValueAsDouble(channelId, 0);
0aa03544	186	pedSigma = ped->ValueAsDouble(channelId, 1);
3f637cfc	187
0aa03544	188	ped->SetValueAsDouble(channelId, 1, TMath::Sqrt(TMath::Abs(pedSigma/(Double_t)nEvents - pedMean*pedMean)));
3f637cfc	189
0aa03544	190	pedMean = ped->ValueAsDouble(channelId, 0) + 0.5 ;
0aa03544	191	pedSigma = ped->ValueAsDouble(channelId, 1);
3f637cfc	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	}
3f637cfc	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
0aa03544	295	Double_t pedMean[10];
	296	Double_t pedSigma[10];
	297	Double_t injCharge[10];
	298	Double_t injChargeErr[10];
3f637cfc	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];
0aa03544	327	// AliMpIntPair* pair;
3f637cfc	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
fd99693f	334	// read back info
3f637cfc	335	for (Int_t i = 0; i < nEntries; ++i) {
fd99693f	336	map[i] = 0x0;
fd99693f	337	run[i] = 0x0;
3f637cfc	338	tree->SetBranchAddress("ped",&map[i]);
	339	tree->SetBranchAddress("run",&run[i]);
	340	tree->GetEvent(i);
	341	}
	342
	343	// Q = f(ADC)
a76f513d	344	TF1* func = new TF1("func",fitFunc, 0., 2500., NFITPARAMS);
fd99693f	345	// TF1* func = new TF1("func","pol1");
3f637cfc	346
3f637cfc	347	// iterates over the first pedestal run
0aa03544	348	TIter next(map[0]->CreateIterator());
2bf54af2	349	AliMUONVCalibParam* p;
3f637cfc	350
2bf54af2	351	while ( ( p = dynamic_cast<AliMUONVCalibParam*>(next() ) ) )
3f637cfc	352	{
2bf54af2	353	ped[0] = p;
3f637cfc	354
0aa03544	355	Int_t busPatchId = p->ID0();
0aa03544	356	Int_t manuId = p->ID1();
3f637cfc	357
3f637cfc	358	// read back pedestal from the other runs for the given (bupatch, manu)
3f637cfc	359	for (Int_t i = 1; i < nEntries; ++i) {
0aa03544	360	ped[i] = static_cast<AliMUONVCalibParam*>(map[i]->FindObject(busPatchId, manuId));
3f637cfc	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.
0aa03544	370	pedMean[i] = ped[i]->ValueAsDouble(channelId, 0);
	371	pedSigma[i] = ped[i]->ValueAsDouble(channelId, 1);
	372	injCharge[i] = (Double_t)run[i]->GetSecond();
3f637cfc	373	injChargeErr[i] = 1.;
3f637cfc	374
fd99693f	375	if (pedMean[i] < 0) continue; // not connected
3f637cfc	376
	377	if (pedSigma[i] <= 0) pedSigma[i] = 1.; // should not happen.
	378	n++;
	379	}
	380
a76f513d	381	if (n > NFITPARAMS) {
fd99693f	382	// if (n > 1) {
fd99693f	383	//fit
b1611778	384	TGraph *gain = new TGraphErrors(n, pedMean, injCharge, pedSigma, injChargeErr);
3f637cfc	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");
fd99693f	392	cout << setw(8) << func->GetParameter(0) <<"\t"<< setw(8) << func->GetParameter(1) << endl;
fd99693f	393
3f637cfc	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	}
3f637cfc	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
fd99693f	420	// needed for streamer application
	421	gROOT->GetPluginManager()->AddHandler("TVirtualStreamerInfo",
	422	"*",
	423	"TStreamerInfo",
	424	"RIO",
	425	"TStreamerInfo()");
	426
a76f513d	427	TFitter *minuitFit = new TFitter(NFITPARAMS);
a76f513d	428	TVirtualFitter::SetFitter(minuitFit);
fd99693f	429
3f637cfc	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;
0aa03544	433	Double_t nSigma = 3;
3f637cfc	434	Int_t threshold = -1;
	435	Char_t inputFile[256];
	436	TString flatOutputFile;
	437	TString crocusOutputFile;
	438	TString crocusConfigFile;
	439
3f637cfc	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++;
0aa03544	485	sscanf(argv[i],"%lf",&nSigma);
3f637cfc	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
3f637cfc	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
0aa03544	571	Int_t busPatchId;
	572	UShort_t manuId;
	573	UChar_t channelId;
	574	UShort_t charge;
	575
3f637cfc	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
0aa03544	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);
3f637cfc	626
0aa03544	627	} // Next digit
3f637cfc	628
	629	delete rawReader;
	630	delete rawStream;
	631
	632	} // while (1)
	633
3f637cfc	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
a76f513d	651	delete minuitFit;
	652	TVirtualFitter::SetFitter(0);
	653
3f637cfc	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	}