X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=MUON%2FMUONTRKda.cxx;h=ff4dfd0e44d89f6bc1dac65c999858032e42286c;hb=061d7a52a64e5cfcc44a766ccb92333e4d21e843;hp=f72351004597956535e31422c8ac6a08e60f34b6;hpb=fd99693f98b46e44fea7701e280394cd469738d0;p=u%2Fmrichter%2FAliRoot.git diff --git a/MUON/MUONTRKda.cxx b/MUON/MUONTRKda.cxx index f7235100459..ff4dfd0e44d 100644 --- a/MUON/MUONTRKda.cxx +++ b/MUON/MUONTRKda.cxx @@ -1,13 +1,43 @@ /* +Contact: Jean-Luc Charvet + Link: http://aliceinfo.cern.ch/static/Offline/dimuon/muon_html/README_Mchda +Run Type: PEDESTAL, CALIBRATION + DA Type: LDC + Number of events needed: 400 events for pedestal and each calibration run + Input Files: Rawdata file (DATE format) + Output Files: local dir (not persistent) -> MUONTRKda_ped_.ped , MUONTRKda_gain_.par + FXS -> run<#>_MCH__PEDESTALS, run<#>_MCH__GAINS + Trigger types used: +*/ + +/************************************************************************** +* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +/* $Id$ */ + +/* + ------------------------------------------------------------------------- + 2008-11-14 New version: MUONTRKda.cxx,v 1.15 + ------------------------------------------------------------------------- -Version 1 for MUONTRKda MUON tracking -Working version for pedestal -Framework for gain computing -(Ch. Finck) + Version for MUONTRKda MUON tracking + (A. Baldisseri, J.-L. Charvet & Ch. Finck) -Rem: AliMUON2DMap stores all channels, even those which are not connected - if pedMean == -1, channel not connected to a pad + Rem: AliMUON2DMap stores all channels, even those which are not connected + if pedMean == -1, channel not connected to a pad */ @@ -21,26 +51,26 @@ extern "C" { #include #include #include +#include //AliRoot -#include "AliMUONRawStreamTracker.h" -#include "AliMUONDspHeader.h" -#include "AliMUONBlockHeader.h" -#include "AliMUONBusStruct.h" -#include "AliMUONDDLTracker.h" -#include "AliMUONV2DStore.h" +// #include "AliMUONLogger.h" +#include "AliMUONRawStreamTrackerHP.h" +// #include "AliMUONDspHeader.h" +// #include "AliMUONBlockHeader.h" +// #include "AliMUONBusStruct.h" +// #include "AliMUONDDLTracker.h" +#include "AliRawReader.h" +#include "AliMUONVStore.h" #include "AliMUON2DMap.h" -#include "AliMUONCalibParamNF.h" -#include "AliMUONObjectPair.h" -#include "AliMUONVDataIterator.h" -#include "AliMpDDLStore.h" +#include "AliMUONCalibParamND.h" #include "AliMpIntPair.h" - -#include "AliRawReaderDate.h" - +#include "AliMpConstants.h" +#include "AliRawDataErrorLog.h" //ROOT #include "TFile.h" +#include "TSystem.h" #include "TTree.h" #include "TH1F.h" #include "TString.h" @@ -51,373 +81,891 @@ extern "C" { #include "TF1.h" #include "TROOT.h" #include "TPluginManager.h" +#include "TFitter.h" +#include "THashTable.h" +#include +#define NFITPARAMS 4 // global variables -AliMUONV2DStore* pedestalStore = new AliMUON2DMap(kFALSE); -const Int_t kNchannels = 64; -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 +const Int_t kNChannels = AliMpConstants::ManuNofChannels(); +const Int_t kADCMax = 4095; -//________________ -Double_t fitFunc(Double_t *x, Double_t *par) -{ - //fit function for gains +AliMUONVStore* gAliPedestalStore = new AliMUON2DMap(kFALSE); - Float_t xx = x[0]; - Double_t f; +Int_t gAliNManu = 0; +Int_t gAliNChannel = 0; +UInt_t gAliRunNumber = 0; +Int_t gAliNEvents = 0; +Int_t gAliNEventsRecovered = 0; +Int_t gAliPrintLevel = 1; // global printout variable (others: 2 and 3) +Int_t gAliPlotLevel = 0; // global plot variable - if (xx < par[3]) - f = par[0] + par[1]*xx; - else - f= par[0] + par[1]*xx + par[2]*xx*xx; +Char_t gAliHistoFileName[256]; - return f; -} +// used for computing gain parameters +Int_t gAlinbpf1 = 6; // linear fit over nbf1 points +Char_t gAliHistoFileNamegain[256]="MUONTRKda_gain.data"; +Char_t gAliOutFolder[256]="."; +Char_t gAlifilename[256]; +Char_t gAlifilenam[256]="MUONTRKda_gain"; +Char_t gAliflatFile[256]=""; +ofstream gAlifilcout; -//__________ -void MakePed(Int_t busPatchId, Int_t manuId, Int_t channelId, Int_t charge) -{ +TString gAliOutputFile; +TString gAliCommand("ped"); +TTimeStamp gAlidate; - AliMUONVCalibParam* ped = - static_cast(pedestalStore->Get(busPatchId, manuId)); +class ErrorCounter : public TNamed +{ +public : + ErrorCounter(Int_t bp = 0, Int_t manu = 0, Int_t ev = 1) : busPatch(bp), manuId(manu), events(ev) {} + void Increment() {events++;} + Int_t BusPatch() const {return busPatch;} + Int_t ManuId() const {return manuId;} + Int_t Events() const {return events;} + Int_t Compare(const TObject* obj) const + { + /// Compare function + Int_t patch1, patch2, manu1, manu2; + patch1 = busPatch; + manu1 = manuId; + patch2 = ((ErrorCounter*)obj)->BusPatch(); + manu2 = ((ErrorCounter*)obj)->ManuId(); + + if (patch1 == patch2) + { + if (manu1 == manu2) + { + return 0; + } + else + return (manu1 >= manu2) ? 1 : -1; + } + else + return (patch1 >= patch2) ? 1 : -1; + } + + void Print(const Option_t* option="") const + { + TNamed::Print(option); + cout<<"bp "<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"<Iterator(); + if(x[0] <= xLim) return par[0] + par[1]*x[0]; - AliMUONObjectPair* p; - - while ( ( p = dynamic_cast(it->Next() ) ) ) - { - AliMUONVCalibParam* ped = dynamic_cast(p->Value()); - AliMpIntPair* pair = dynamic_cast(p->Key()); - busPatchId = pair->GetFirst(); - manuId = pair->GetSecond(); + Double_t yLim = par[0]+ par[1]*xLim; + return yLim + par[1]*(x[0] - xLim) + par[2]*(x[0] - xLim)*(x[0] - xLim); +} - for (channelId = 0; channelId < ped->Size() ; ++channelId) { - pedMean = ped->ValueAsFloat(channelId, 0); +//__________ +void MakePed(Int_t busPatchId, Int_t manuId, Int_t channelId, Int_t charge) +{ + /// Compute pedestals values + AliMUONVCalibParam* ped = + static_cast(gAliPedestalStore->FindObject(busPatchId, manuId)); + + if (!ped) { + gAliNManu++; + ped = new AliMUONCalibParamND(2, kNChannels,busPatchId, manuId, -1.); // put default wise -1, not connected channel + gAliPedestalStore->Add(ped); + } - if (pedMean > 0) { // connected channels + // if (gAliNEvents == 0) { + // ped->SetValueAsDouble(channelId, 0, 0.); + // ped->SetValueAsDouble(channelId, 1, 0.); + // } + + // Initialization for the first value + if (ped->ValueAsDouble(channelId, 0) == -1) ped->SetValueAsDouble(channelId, 0, 0.); + if (ped->ValueAsDouble(channelId, 1) == -1) ped->SetValueAsDouble(channelId, 1, 0.); - ped->SetValueAsFloat(channelId, 0, pedMean/(Float_t)nEvents); + Double_t pedMean = ped->ValueAsDouble(channelId, 0) + (Double_t) charge; + Double_t pedSigma = ped->ValueAsDouble(channelId, 1) + (Double_t) charge*charge; - pedMean = ped->ValueAsFloat(channelId, 0); - pedSigma = ped->ValueAsFloat(channelId, 1); + ped->SetValueAsDouble(channelId, 0, pedMean); + ped->SetValueAsDouble(channelId, 1, pedSigma); - ped->SetValueAsFloat(channelId, 1, TMath::Sqrt(TMath::Abs(pedSigma/(Float_t)nEvents - pedMean*pedMean))); +} - pedMean = ped->ValueAsFloat(channelId, 0) + 0.5 ; - pedSigma = ped->ValueAsFloat(channelId, 1); +//________________ +void MakePedStore(TString gAliOutputFile_1 = "") +{ + /// Store pedestals in ASCII files + Double_t pedMean; + Double_t pedSigma; + ofstream fileout; + Int_t busPatchId; + Int_t manuId; + Int_t channelId; + +// histo + TFile* histoFile = 0; + TTree* tree = 0; + TH1F* pedMeanHisto = 0; + TH1F* pedSigmaHisto = 0; + if (gAliCommand.CompareTo("ped") == 0) + { + sprintf(gAliHistoFileName,"%s/MUONTRKda_ped_%d.root",gAliOutFolder,gAliRunNumber); + histoFile = new TFile(gAliHistoFileName,"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 = 4096; + Int_t xmin = 0; + Int_t xmax = 4095; + pedMeanHisto = new TH1F(name,title,nx,xmin,xmax); + pedMeanHisto->SetDirectory(histoFile); + + sprintf(name,"pedsigma_allch"); + sprintf(title,"Pedestal sigma all channels"); + nx = 201; + xmin = 0; + xmax = 200; + pedSigmaHisto = new TH1F(name,title,nx,xmin,xmax); + pedSigmaHisto->SetDirectory(histoFile); + + tree = new TTree("t","Pedestal tree"); + tree->Branch("bp",&busPatchId,"bp/I"); + tree->Branch("manu",&manuId,",manu/I"); + tree->Branch("channel",&channelId,",channel/I"); + tree->Branch("pedMean",&pedMean,",pedMean/D"); + tree->Branch("pedSigma",&pedSigma,",pedSigma/D"); + } + if (!gAliOutputFile_1.IsNull()) { + fileout.open(gAliOutputFile_1.Data()); + fileout<<"//===========================================================================" << endl; + fileout<<"// Pedestal file calculated by MUONTRKda"<GetSize()) + { + fileout<<"//"<MakeIterator(); + ErrorCounter* parityerror; + while((parityerror = (ErrorCounter*) iter->Next())) + { + fileout<<"// bp "<BusPatch()<<" events used "<Events()<GetSize()) + { + cout<<"\n* Buspatches with less statistics (due to parity errors)"<MakeIterator(); + ErrorCounter* parityerror; + while((parityerror = (ErrorCounter*) iter->Next())) + { + cout<<" bp "<BusPatch()<<": events used = "<Events()<BusPatch()<<": events used = "<Events()<Fill(pedMean); - pedSigmaHisto->Fill(pedSigma); + +// iterator over pedestal + TIter next(gAliPedestalStore->CreateIterator()); + AliMUONVCalibParam* ped; - tree->Fill(); - } - } - - delete p; // cos create a new object for each iteration - } + while ( ( ped = dynamic_cast(next() ) ) ) + { + busPatchId = ped->ID0(); + manuId = ped->ID1(); + Int_t eventCounter; + + // Correct the number of events for buspatch with errors + char bpname[256]; + ErrorCounter* errorCounter; + sprintf(bpname,"bp%d",busPatchId); + if ((errorCounter = (ErrorCounter*)gAliErrorBuspatchTable->FindObject(bpname))) + { + eventCounter = gAliNEvents - errorCounter->Events(); + } + else + { + eventCounter = gAliNEvents; + } + + for (channelId = 0; channelId < ped->Size() ; ++channelId) { + pedMean = ped->ValueAsDouble(channelId, 0); + + if (pedMean > 0) { // connected channels + + ped->SetValueAsDouble(channelId, 0, pedMean/(Double_t)eventCounter); + + pedMean = ped->ValueAsDouble(channelId, 0); + pedSigma = ped->ValueAsDouble(channelId, 1); + + ped->SetValueAsDouble(channelId, 1, TMath::Sqrt(TMath::Abs(pedSigma/(Double_t)eventCounter - pedMean*pedMean))); + + pedMean = ped->ValueAsDouble(channelId, 0); + pedSigma = ped->ValueAsDouble(channelId, 1); + + + if (!gAliOutputFile_1.IsNull()) { + fileout << "\t" << busPatchId << "\t" << manuId <<"\t"<< channelId << "\t" + << pedMean <<"\t"<< pedSigma << endl; + } + + if (gAliCommand.CompareTo("ped") == 0) + { + pedMeanHisto->Fill(pedMean); + pedSigmaHisto->Fill(pedSigma); + + tree->Fill(); + } + } + } +} - // file outputs - if (!flatOutputFile.IsNull()) - fileout.close(); +// file outputs +if (!gAliOutputFile_1.IsNull()) fileout.close(); - histoFile->Write(); - histoFile->Close(); +if (gAliCommand.CompareTo("ped") == 0) +{ + 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() +void MakePedStoreForGain(Int_t injCharge) { - // store pedestal map in root file - - Int_t injCharge = 200; - - TTree* tree = 0x0; + /// Store pedestal map in root file + TTree* tree = 0x0; + + FILE *pfilew=0; + if (gAliCommand.Contains("gain") && !gAliCommand.Contains("comp")) { + if(gAliOutputFile.IsNull()) + { + sprintf(gAlifilename,"%s_%d_DAC_%d.par",gAlifilenam,gAliRunNumber,injCharge); + gAliOutputFile=gAlifilename; + } + if(!gAliOutputFile.IsNull()) + { + pfilew = fopen (gAliOutputFile.Data(),"w"); + + fprintf(pfilew,"//DUMMY FILE (to prevent Shuttle failure)\n"); + fprintf(pfilew,"//================================================\n"); + fprintf(pfilew,"// MUONTRKda: Calibration run \n"); + fprintf(pfilew,"//=================================================\n"); + fprintf(pfilew,"// * Run : %d \n",gAliRunNumber); + fprintf(pfilew,"// * Date : %s \n",gAlidate.AsString("l")); + fprintf(pfilew,"// * DAC : %d \n",injCharge); + fprintf(pfilew,"//-------------------------------------------------\n"); + fclose(pfilew); + } + } - // compute and store pedestals - MakePedStore(); + if(gAliPrintLevel==2) + { + // compute and store pedestals + sprintf(gAliflatFile,"%s/%s_%d_DAC_%d.ped",gAliOutFolder,gAlifilenam,gAliRunNumber,injCharge); + cout << "\nMUONTRKda : Flat file generated : " << gAliflatFile << "\n"; + MakePedStore(gAliflatFile); + } + else + MakePedStore(); - // store in root file - sprintf(histoFileName,"mutrkgain.root"); - - TString mode("UPDATE"); + TString mode("UPDATE"); - if (command.Contains("cre")) { - mode = "RECREATE"; - } - TFile* histoFile = new TFile(histoFileName, mode.Data(), "MUON Tracking Gains"); + if (gAliCommand.Contains("cre")) { + mode = "RECREATE"; + } + TFile* histoFile = new TFile(gAliHistoFileNamegain, 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); + // 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(gAliRunNumber, injCharge); - if (mode.CompareTo("UPDATE") == 0) { - tree = (TTree*)histoFile->Get("t"); - tree->SetBranchAddress("run",&pair); - tree->SetBranchAddress("ped",&pedestalStore); + if (mode.CompareTo("UPDATE") == 0) { + tree = (TTree*)histoFile->Get("t"); + tree->SetBranchAddress("run",&pair); + tree->SetBranchAddress("ped",&gAliPedestalStore); - } 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); + } else { + tree = new TTree("t","Pedestal tree"); + tree->Branch("run", "AliMpIntPair",&pair); + tree->Branch("ped", "AliMUON2DMap",&gAliPedestalStore); + tree->SetBranchAddress("run",&pair); + tree->SetBranchAddress("ped",&gAliPedestalStore); - } + } - tree->Fill(); - tree->Write("t", TObject::kOverwrite); // overwrite the tree - histoFile->Close(); + tree->Fill(); + tree->Write("t", TObject::kOverwrite); // overwrite the tree + histoFile->Close(); - delete pair; + delete pair; } //________________ -void MakeGainStore(TString flatOutputFile) +void MakeGainStore() { - - // 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. - - Float_t pedMean[10]; - Float_t pedSigma[10]; - Float_t injCharge[10]; - Float_t injChargeErr[10]; - - ofstream fileout; - TTimeStamp date; - - if (!flatOutputFile.IsNull()) { - fileout.open(flatOutputFile.Data()); - fileout<<"//===========================================================================" << endl; - fileout<<"// Pedestal file calculated by MUONTRKda"<Get("t"); + Int_t nEntries = tree->GetEntries(); + Int_t nbpf2 = nEntries - (nInit + gAlinbpf1) + 1; // nb pts used for 2nd order fit + + // 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); + // std::cout << map[i] << " " << run[i] << std::endl; + } + //jlc_feb_08 modif: gAliRunNumber=(UInt_t)run[0]->GetFirst(); + gAliRunNumber=(UInt_t)run[nEntries-1]->GetFirst(); + // sscanf(getenv("DATE_RUN_NUMBER"),"%d",&gAliRunNumber); + + Double_t pedMean[11]; + Double_t pedSigma[11]; + for ( Int_t i=0 ; i<11 ; i++) {pedMean[i]=0.;pedSigma[i]=1.;}; + Double_t injCharge[11]; + Double_t injChargeErr[11]; + for ( Int_t i=0 ; i<11 ; i++) {injCharge[i]=0.;injChargeErr[i]=1.;}; + + // some print + cout<<"\n ******** MUONTRKda for Gain computing (Run = " << gAliRunNumber << ")\n" << endl; + cout<<" * Date : " << gAlidate.AsString("l") << "\n" << endl; + cout << " Entries = " << nEntries << " DAC values \n" << endl; + for (Int_t i = 0; i < nEntries; ++i) { + cout<< " Run = " << run[i]->GetFirst() << " DAC = " << run[i]->GetSecond() << endl; + numrun[i] = run[i]->GetFirst(); + injCharge[i] = run[i]->GetSecond(); + injChargeErr[i] = 0.01*injCharge[i]; + if(injChargeErr[i] <= 1.) injChargeErr[i]=1.; + } + cout << "" << endl; - } + // full print out + sprintf(gAlifilename,"%s/%s_%d.log",gAliOutFolder,gAlifilenam,gAliRunNumber); + logOutputFilecomp=gAlifilename; - sprintf(histoFileName,"mutrkgain.root"); - TFile* histoFile = new TFile(histoFileName); + filcouc.open(logOutputFilecomp.Data()); + filcouc<<"//====================================================" << endl; + filcouc<<"// MUONTRKda for Gain computing (Run = " << gAliRunNumber << ")" << endl; + filcouc<<"//====================================================" << endl; + filcouc<<"// * Date : " << gAlidate.AsString("l") << "\n" << endl; - 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); + // why 2 files ? (Ch. F.) + FILE *pfilen = 0; + if(gAliPrintLevel==2) + { + sprintf(gAlifilename,"%s/%s_%d.param",gAliOutFolder,gAlifilenam,gAliRunNumber); + cout << " fit parameter file = " << gAlifilename << "\n"; + pfilen = fopen (gAlifilename,"w"); + + fprintf(pfilen,"//===================================================================\n"); + fprintf(pfilen,"// BP MANU CH. par[0] [1] [2] [3] xlim P(chi2) p1 P(chi2)2 p2\n"); + fprintf(pfilen,"//===================================================================\n"); + fprintf(pfilen,"// * Run : %d \n",gAliRunNumber); + fprintf(pfilen,"//===================================================================\n"); + } + + FILE *pfilew=0; + if(gAliOutputFile.IsNull()) + { + sprintf(gAlifilename,"%s_%d.par",gAlifilenam,gAliRunNumber); + gAliOutputFile=gAlifilename; + } + if(!gAliOutputFile.IsNull()) + { + pfilew = fopen (gAliOutputFile.Data(),"w"); + + fprintf(pfilew,"//================================================\n"); + fprintf(pfilew,"// Calibration file calculated by MUONTRKda \n"); + fprintf(pfilew,"//=================================================\n"); + fprintf(pfilew,"// * Run : %d \n",gAliRunNumber); + fprintf(pfilew,"// * Date : %s \n",gAlidate.AsString("l")); + fprintf(pfilew,"// * Statictics : %d \n",gAliNEvents); + fprintf(pfilew,"// * # of MANUS : %d \n",gAliNManu); + fprintf(pfilew,"// * # of channels : %d \n",gAliNChannel); + fprintf(pfilew,"//-------------------------------------------------\n"); + if(nInit==0) + fprintf(pfilew,"// %d DAC values fit: %d pts (1st order) %d pts (2nd order) \n",nEntries,gAlinbpf1,nbpf2); + if(nInit==1) + fprintf(pfilew,"// %d DAC values fit: %d pts (1st order) %d pts (2nd order) DAC=0 excluded\n",nEntries,gAlinbpf1,nbpf2); + fprintf(pfilew,"// RUN DAC \n"); + fprintf(pfilew,"//-----------------\n"); + for (Int_t i = 0; i < nEntries; ++i) { + tree->SetBranchAddress("run",&run[i]); + fprintf(pfilew,"// %d %5.0f \n",numrun[i],injCharge[i]); + } + fprintf(pfilew,"//=======================================\n"); + fprintf(pfilew,"// BP MANU CH. a1 a2 thres. q\n"); + fprintf(pfilew,"//=======================================\n"); + } + + FILE *pfilep = 0; + if(gAliPrintLevel==2) + { + sprintf(gAlifilename,"%s/%s_%d.peak",gAliOutFolder,gAlifilenam,gAliRunNumber); + cout << " File containing Peak mean values = " << gAlifilename << "\n"; + pfilep = fopen (gAlifilename,"w"); + + fprintf(pfilep,"//==============================================================================================================================\n"); + fprintf(pfilep,"// * Run : %d \n",gAliRunNumber); + fprintf(pfilep,"//==============================================================================================================================\n"); + fprintf(pfilep,"// BP MANU CH. Ped. <0> <1> <2> <3> <4> <5> <6> <7> <8> <9> <10> \n"); + fprintf(pfilep,"//==============================================================================================================================\n"); + fprintf(pfilep,"// DAC= %9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f fC\n",injCharge[0],injCharge[1],injCharge[2],injCharge[3],injCharge[4],injCharge[5],injCharge[6],injCharge[7],injCharge[8],injCharge[9],injCharge[10]); + fprintf(pfilep,"// %9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f\n",injChargeErr[0],injChargeErr[1],injChargeErr[2],injChargeErr[3],injChargeErr[4],injChargeErr[5],injChargeErr[6],injChargeErr[7],injChargeErr[8],injChargeErr[9],injChargeErr[10]); + fprintf(pfilep,"//==============================================================================================================================\n"); } - - // Q = f(ADC) - TF1* func = new TF1("func",fitFunc, 0., 2500., 4); - // TF1* func = new TF1("func","pol1"); - // iterates over the first pedestal run - AliMUONVDataIterator* it = map[0]->Iterator(); - AliMUONObjectPair* p; - while ( ( p = dynamic_cast(it->Next() ) ) ) + // plot out + + TFile* gainFile = 0x0; + sprintf(rootFileName,"%s/%s_%d.root",gAliOutFolder,gAlifilenam,gAliRunNumber); + gainFile = new TFile(rootFileName,"RECREATE"); + + Double_t chi2 = 0.; + Double_t chi2P2 = 0.; + Double_t prChi2 = 0; + Double_t prChi2P2 =0; + Double_t a0=0.,a1=1.,a2=0.; + Int_t busPatchId ; + Int_t manuId ; + Int_t channelId ; + Int_t threshold = 0; + Int_t q = 0; + Int_t p1 =0; + Int_t p2 =0; + Double_t gain=0; + Double_t capa=0.2; // internal capacitor (pF) + + TTree *tg = new TTree("tg","TTree avec class Manu_DiMu"); + + tg->Branch("bp",&busPatchId, "busPatchId/I"); + tg->Branch("manu",&manuId, "manuId/I"); + tg->Branch("channel",&channelId, "channelId/I"); + + tg->Branch("a0",&a0, "a0/D"); + tg->Branch("a1",&a1, "a1/D"); + tg->Branch("a2",&a2, "a2/D"); + tg->Branch("Pchi2",&prChi2, "prChi2/D"); + tg->Branch("Pchi2_2",&prChi2P2, "prChi2P2/D"); + tg->Branch("Threshold",&threshold, "threshold/I"); + tg->Branch("q",&q, "q/I"); + tg->Branch("p1",&p1, "p1/I"); + tg->Branch("p2",&p2, "p2/I"); + tg->Branch("gain",&gain, "gain/D"); + + char graphName[256]; + + // iterates over the first pedestal run + TIter next(map[0]->CreateIterator()); + AliMUONVCalibParam* p; + + Int_t nmanu = 0; + Int_t nGoodChannel = 0; + Int_t nBadChannel = 0; + Int_t noFitChannel = 0; + Int_t nplot=0; + Double_t sumProbChi2 = 0.; + Double_t sumA1 = 0.; + Double_t sumProbChi2P2 = 0.; + Double_t sumA2 = 0.; + + Double_t x[11], xErr[11], y[11], yErr[11]; + Double_t xp[11], xpErr[11], yp[11], ypErr[11]; + + Int_t uncalcountertotal=0 ; + + while ( ( p = dynamic_cast(next() ) ) ) { - ped[0] = dynamic_cast(p->Value()); - pair = dynamic_cast(p->Key()); + ped[0] = p; - Int_t busPatchId = pair->GetFirst(); - Int_t manuId = pair->GetSecond(); + busPatchId = p->ID0(); + 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(map[i]->Get(busPatchId, manuId)); + ped[i] = static_cast(map[i]->FindObject(busPatchId, manuId)); } // compute for each channel the gain parameters - for ( Int_t channelId = 0; channelId < ped[0]->Size() ; ++channelId ) { + for ( channelId = 0; channelId < ped[0]->Size() ; ++channelId ) + { - Int_t n = 0; - for (Int_t i = 0; i < nEntries; ++i) { + Int_t n = 0; + for (Int_t i = 0; i < nEntries; ++i) { - if (!ped[i]) continue; //shouldn't happen. - pedMean[i] = ped[i]->ValueAsFloat(channelId, 0); - pedSigma[i] = ped[i]->ValueAsFloat(channelId, 1); - injCharge[i] = (Float_t)run[i]->GetSecond(); - injChargeErr[i] = 1.; + if (!ped[i]) continue; //shouldn't happen. + pedMean[i] = ped[i]->ValueAsDouble(channelId, 0); + pedSigma[i] = ped[i]->ValueAsDouble(channelId, 1); - if (pedMean[i] < 0) continue; // not connected + if (pedMean[i] < 0) continue; // not connected - if (pedSigma[i] <= 0) pedSigma[i] = 1.; // should not happen. - n++; - } + if (pedSigma[i] <= 0) pedSigma[i] = 1.; // should not happen. + n++; + } + + + // print_peak_mean_values + if(gAliPrintLevel==2) + { + + fprintf(pfilep,"%4i%5i%5i%10.3f",busPatchId,manuId,channelId,0.); + fprintf(pfilep,"%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f%9.1f \n",pedMean[0],pedMean[1],pedMean[2],pedMean[3],pedMean[4],pedMean[5],pedMean[6],pedMean[7],pedMean[8],pedMean[9],pedMean[10]); + fprintf(pfilep," sig= %9.3f%9.3f%9.3f%9.3f%9.3f%9.3f%9.3f%9.3f%9.3f%9.3f%9.3f \n",pedSigma[0],pedSigma[1],pedSigma[2],pedSigma[3],pedSigma[4],pedSigma[5],pedSigma[6],pedSigma[7],pedSigma[8],pedSigma[9],pedSigma[10]); + } + + // makegain + + + // Fit Method: Linear fit over gAlinbpf1 points + parabolic fit over nbpf2 points) + // nInit=1 : 1st pt DAC=0 excluded + + // 1. - linear fit over gAlinbpf1 points + + Double_t par[4] = {0.,0.5,0.,kADCMax}; + Int_t nbs = nEntries - nInit; + if(nbs < gAlinbpf1)gAlinbpf1=nbs; - if (n > 4) { - // if (n > 1) { - //fit - TGraph *gain = new TGraphErrors(n, pedMean, pedSigma, injCharge, 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; + Int_t fitproceed=1; + for (Int_t j = 0; j < nbs; ++j) + { + Int_t k = j + nInit; + x[j] = pedMean[k]; + if(x[j]==0.)fitproceed=0; + xErr[j] = pedSigma[k]; + y[j] = injCharge[k]; + yErr[j] = injChargeErr[k]; + + } + + TGraphErrors *graphErr; + if(!fitproceed) { p1=0; p2=0; noFitChannel++;} + + if(fitproceed) + { + + TF1 *f1 = new TF1("f1",funcLin,0.,kADCMax,2); + graphErr = new TGraphErrors(gAlinbpf1, x, y, xErr, yErr); + + f1->SetParameters(0,0); + + graphErr->Fit("f1","RQ"); + + chi2 = f1->GetChisquare(); + f1->GetParameters(par); + + delete graphErr; + graphErr=0; + delete f1; + + prChi2 = TMath::Prob(chi2, gAlinbpf1 - 2); + + Double_t xLim = pedMean[nInit + gAlinbpf1 - 1]; + Double_t yLim = par[0]+par[1] * xLim; + + a0 = par[0]; + a1 = par[1]; + + // 2. - Translation : new origin (xLim, yLim) + parabolic fit over nbf2 points + + if(nbpf2 > 1) + { + for (Int_t j = 0; j < nbpf2; j++) + { + Int_t k = j + (nInit + gAlinbpf1) - 1; + xp[j] = pedMean[k] - xLim; + xpErr[j] = pedSigma[k]; + + yp[j] = injCharge[k] - yLim - par[1]*xp[j]; + ypErr[j] = injChargeErr[k]; + } + + TF1 *f2 = new TF1("f2",funcParabolic,0.,kADCMax,1); + graphErr = new TGraphErrors(nbpf2, xp, yp, xpErr, ypErr); + + graphErr->Fit(f2,"RQ"); + chi2P2 = f2->GetChisquare(); + f2->GetParameters(par); + + delete graphErr; + graphErr=0; + delete f2; + + prChi2P2 = TMath::Prob(chi2P2, nbpf2-1); + a2 = par[0]; + + // delete graphErr; + + } + + par[0] = a0; + par[1] = a1; + par[2] = a2; + par[3] = xLim; + +// p1 = TMath::Nint(ceil(prChi2*14))+1; // round up value : ceil (2.2)=3. +// p2 = TMath::Nint(ceil(prChi2P2*14))+1; + if(prChi2>0.999999)prChi2=0.999999 ; if(prChi2P2>0.999999)prChi2P2=0.9999999; // avoiding Pr(Chi2)=1 value + p1 = TMath::Nint(floor(prChi2*15))+1; // round down value : floor(2.8)=2. + p2 = TMath::Nint(floor(prChi2P2*15))+1; + q = p1*16 + p2; // fit quality + + Double_t x0 = -par[0]/par[1]; // value of x corresponding to à 0 fC + threshold = TMath::Nint(ceil(par[3]-x0)); // linear if x < threshold + + if(gAliPrintLevel==2) + { + fprintf(pfilen,"%4i %4i %2i",busPatchId,manuId,channelId); + fprintf(pfilen," %6.2f %6.4f %10.3e %4.2f %4i %8.6f %8.6f %x %8.6f %8.6f %x\n", + par[0], par[1], par[2], par[3], threshold, prChi2, floor(prChi2*15), p1, prChi2P2, floor(prChi2P2*15),p2); + } + // tests + if(par[1]< goodA1Min || par[1]> goodA1Max) p1=0; + if(par[2]< goodA2Min || par[2]> goodA2Max) p2=0; + + } // fitproceed + + if(fitproceed && p1>0 && p2>0) + { + nGoodChannel++; + sumProbChi2 += prChi2; + sumA1 += par[1]; + gain=1./(par[1]*capa); + sumProbChi2P2 += prChi2P2; + sumA2 += par[2]; + } + else // bad calibration + { + nBadChannel++; + q=0; + par[1]=0.5; a1=0.5; p1=0; + par[2]=0.; a2=0.; p2=0; + threshold=kADCMax; + + char bpmanuname[256]; + ErrorCounter* uncalcounter; + + sprintf(bpmanuname,"bp%dmanu%d",busPatchId,manuId); + if (!(uncalcounter = (ErrorCounter*)uncalBuspatchManuTable->FindObject(bpmanuname))) + { + // New buspatch_manu name + uncalcounter= new ErrorCounter (busPatchId,manuId); + uncalcounter->SetName(bpmanuname); + uncalBuspatchManuTable->Add(uncalcounter); + } + else + { + // Existing buspatch_manu name + uncalcounter->Increment(); + } + // uncalcounter->Print_uncal() + uncalcountertotal ++; + } + + if(gAliPlotLevel){ + // if(q==0 and nplot < 100) + // if(p1>1 && p2==0 and nplot < 100) + // if(p1>1 && p2>1 and nplot < 100) + // if(p1>=1 and p1<=2 and nplot < 100) + if((p1==1 || p2==1) and nplot < 100) + { + nplot++; + // cout << " nplot = " << nplot << endl; + TF1 *f2Calib = new TF1("f2Calib",funcCalib,0.,kADCMax,NFITPARAMS); + + graphErr = new TGraphErrors(nEntries,pedMean,injCharge,pedSigma,injChargeErr); + + sprintf(graphName,"BusPatch_%d_Manu_%d_Ch_%d",busPatchId, manuId,channelId); + + graphErr->SetTitle(graphName); + graphErr->SetMarkerColor(3); + graphErr->SetMarkerStyle(12); + graphErr->Write(graphName); + + sprintf(graphName,"f2_BusPatch_%d_Manu_%d_Ch_%d",busPatchId, manuId,channelId); + f2Calib->SetTitle(graphName); + f2Calib->SetLineColor(4); + f2Calib->SetParameters(par); + f2Calib->Write(graphName); + + delete graphErr; + graphErr=0; + delete f2Calib; + } } - delete gain; - } - } - delete p; + + tg->Fill(); + + if (!gAliOutputFile.IsNull()) + { + fprintf(pfilew,"%4i %5i %2i %7.4f %10.3e %4i %2x\n",busPatchId,manuId,channelId,par[1],par[2],threshold,q); + } + + } + nmanu++; + if(fmod(nmanu,500)==0)std::cout << " Nb manu = " << nmanu << std::endl; } - // file outputs for gain - if (!flatOutputFile.IsNull()) - fileout.close(); + // file outputs for gain + if (!gAliOutputFile.IsNull()) fclose(pfilew); + if(gAliPrintLevel==2){ fclose(pfilen); fclose(pfilep); } + + tg->Write(); + histoFile->Close(); + + //OutPut + if (gAliPrintLevel) + { + // print in logfile + if (uncalBuspatchManuTable->GetSize()) + { + uncalBuspatchManuTable->Sort(); // use compare + TIterator* iter = uncalBuspatchManuTable->MakeIterator(); + ErrorCounter* uncalcounter; + filcouc << "\n List of problematic BusPatch and Manu " << endl; + filcouc << " ========================================" << endl; + filcouc << " BP Manu Nb Channel " << endl ; + filcouc << " ========================================" << endl; + while((uncalcounter = (ErrorCounter*) iter->Next())) + { + filcouc << "\t" << uncalcounter->BusPatch() << "\t " << uncalcounter->ManuId() << "\t\t" << uncalcounter->Events() << endl; + } + filcouc << " ========================================" << endl; + + filcouc << " Number of bad calibrated Manu = " << uncalBuspatchManuTable->GetSize() << endl ; + filcouc << " Number of bad calibrated channel = " << uncalcountertotal << endl; + + } + + + filcouc << "\n Nb of channels in raw data = " << nmanu*64 << " (" << nmanu << " Manu)" << endl; + filcouc << " Nb of calibrated channel = " << nGoodChannel << " (" << goodA1Min << " = " << meanA1 << "\t = " << 1./(meanA1*capaManu) + << " mV/fC (capa= " << capaManu << " pF)" << endl; + filcouc << " Prob(chi2)> = " << meanProbChi2 << endl; + filcouc << "\n parabolic fit: = " << meanA2 << endl; + filcouc << " Prob(chi2)> = " << meanProbChi2P2 << "\n" << endl; + + cout << "\n = " << 1./(meanA1*capaManu) + << " mV/fC (capa= " << capaManu << " pF)" + << " Prob(chi2)> = " << meanProbChi2 << endl; + } - delete func; + filcouc.close(); + cout << "\nMUONTRKda : Output logfile : " << logOutputFilecomp << endl; + cout << "MUONTRKda : Root Histo. file : " << rootFileName << endl; + return ; } @@ -426,35 +974,42 @@ void MakeGainStore(TString flatOutputFile) // main routine int main(Int_t argc, Char_t **argv) { - - // needed for streamer application - gROOT->GetPluginManager()->AddHandler("TVirtualStreamerInfo", - "*", - "TStreamerInfo", - "RIO", - "TStreamerInfo()"); - - - Int_t printLevel = 0; // Global variable defined as extern in the others .cxx files - Int_t skipEvents = 0; - Int_t maxEvents = 1000000; - Float_t nSigma = 3; - Int_t threshold = -1; - Char_t inputFile[256]; - TString flatOutputFile; - TString crocusOutputFile; - TString crocusConfigFile; - - Int_t totalParityErr = 0; - Int_t totalGlitchErr = 0; - -// option handler - - // decode the input line + + // needed for streamer application + gROOT->GetPluginManager()->AddHandler("TVirtualStreamerInfo", + "*", + "TStreamerInfo", + "RIO", + "TStreamerInfo()"); + + TFitter *minuitFit = new TFitter(NFITPARAMS); + TVirtualFitter::SetFitter(minuitFit); + + // ofstream gAlifilcout; + + Int_t fes = 1; // by default FES is used + Int_t skipEvents = 0; + Int_t maxEvents = 1000000; + Int_t maxDateEvents = 1000000; + Int_t injCharge = 0; + Char_t inputFile[256]=""; + + Int_t nDateEvents = 0; + Int_t gGlitchErrors= 0; + Int_t gParityErrors= 0; + Int_t gPaddingErrors= 0; + Int_t recoverParityErrors = 1; + + TString fesOutputFile; + TString logOutputFile; + + // 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]) @@ -465,64 +1020,84 @@ int main(Int_t argc, Char_t **argv) break; case 'a' : i++; - flatOutputFile = argv[i]; + gAliOutputFile = argv[i]; break; - case 'o' : + case 'b' : i++; - crocusOutputFile = argv[i]; + sprintf(gAliOutFolder,argv[i]); break; case 'c' : i++; - crocusConfigFile = argv[i]; + fes=atoi(argv[i]); + break; + case 'd' : + i++; + gAliPrintLevel=atoi(argv[i]); break; case 'e' : i++; - command = argv[i]; + gAliCommand = argv[i]; break; - case 'd' : + case 'g' : i++; - printLevel=atoi(argv[i]); + gAliPlotLevel=atoi(argv[i]); + break; + case 'i' : + i++; + gAlinbpf1=atoi(argv[i]); break; case 's' : i++; skipEvents=atoi(argv[i]); break; - case 'n' : + case 'l' : i++; - sscanf(argv[i],"%d",&maxEvents); + injCharge=atoi(argv[i]); break; - case 'p' : + case 'm' : i++; - sscanf(argv[i],"%f",&nSigma); + sscanf(argv[i],"%d",&maxDateEvents); break; - case 't' : + case 'n' : i++; - sscanf(argv[i],"%d",&threshold); + sscanf(argv[i],"%d",&maxEvents); + break; + case 'r' : + i++; + sprintf(gAliHistoFileNamegain,argv[i]); + break; + case 'p' : + i++; + sscanf(argv[i],"%d",&recoverParityErrors); 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-h help (this screen)"); printf("\n"); printf("\n Input"); - printf("\n-f (default = %s)",inputFile); - printf("\n-c (default = %s)",crocusConfigFile.Data()); + printf("\n-f (default = %s)",inputFile); printf("\n"); printf("\n Output"); - printf("\n-a (default = %s)",flatOutputFile.Data()); - printf("\n-o (default = %s)",crocusOutputFile.Data()); + printf("\n-a (default = %s)",gAliOutputFile.Data()); printf("\n"); printf("\n Options"); - printf("\n-d (default = %d)",printLevel); - printf("\n-s (default = %d)",skipEvents); - printf("\n-n (default = %d)",maxEvents); - printf("\n-p (default = %f)",nSigma); - printf("\n-t (default = %d)",threshold); - printf("\n-e (default = %s)",command.Data()); + printf("\n-b (default = %s)",gAliOutFolder); + printf("\n-c (default = %d)",fes); + printf("\n-d (default = %d)",gAliPrintLevel); + printf("\n-g (default = %d)",gAliPlotLevel); + printf("\n-i (default = %d)",gAlinbpf1); + printf("\n-l (default = %d)",injCharge); + printf("\n-m (default = %d)",maxDateEvents); + printf("\n-s (default = %d)",skipEvents); + printf("\n-n (default = %d)",maxEvents); + printf("\n-r root file data for gain (default = %s)",gAliHistoFileNamegain); + printf("\n-e (default = %s)",gAliCommand.Data()); printf("\n-e make gain & create a new root file"); printf("\n-e make gain & update root file"); printf("\n-e make gain & compute gains"); + printf("\n-p (default = %d)",recoverParityErrors); printf("\n\n"); exit(-1); @@ -532,195 +1107,333 @@ int main(Int_t argc, Char_t **argv) } // end of switch } // end of for i - // set command to lower case - command.ToLower(); + // set gAliCommand to lower case + gAliCommand.ToLower(); - // decoding the events - - Int_t status; - Int_t nDateEvents = 0; - void* event; + // decoding the events - // containers - AliMUONDDLTracker* ddlTracker = 0x0; - AliMUONBlockHeader* blkHeader = 0x0; - AliMUONDspHeader* dspHeader = 0x0; - AliMUONBusStruct* busStruct = 0x0; + Int_t status=0; + // void* event; - pedMeanHisto = 0x0; - pedSigmaHisto = 0x0; + // gAliPedMeanHisto = 0x0; + // gAliPedSigmaHisto = 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 <= 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"<GetType(); - runNumber = rawReader->GetRunNumber(); - - if (eventType != PHYSICS_EVENT) - continue; // for the moment - - nEvents++; - - // decoding MUON payload - AliMUONRawStreamTracker* rawStream = new AliMUONRawStreamTracker(rawReader); + UShort_t manuId; + UChar_t channelId; + UShort_t charge; + TString key("MUONTRKda :"); - // loops over DDL - while((status = rawStream->NextDDL())) { + // AliMUONRawStreamTrackerHP* rawStream = 0; - if (printLevel) printf("\niDDL %d\n", rawStream->GetDDL()); - - ddlTracker = rawStream->GetDDLTracker(); - - // loop over block (CRT) structure - Int_t nBlock = ddlTracker->GetBlkHeaderEntries(); - for(Int_t iBlock = 0; iBlock < nBlock ;iBlock++){ - - blkHeader = ddlTracker->GetBlkHeaderEntry(iBlock); - if (printLevel) printf("Block Total length %d\n",blkHeader->GetTotalLength()); - - // loop over DSP (FRT) structure - Int_t nDsp = blkHeader->GetDspHeaderEntries(); - for(Int_t iDsp = 0; iDsp < nDsp ;iDsp++){ //DSP loop + if (gAliCommand.CompareTo("comp") != 0) + { + + // Rawdeader, RawStreamHP + AliRawReader* rawReader = AliRawReader::Create(inputFile); + AliMUONRawStreamTrackerHP* rawStream = new AliMUONRawStreamTrackerHP(rawReader); + rawStream->DisableWarnings(); + rawStream->EnabbleErrorLogger(); - dspHeader = blkHeader->GetDspHeaderEntry(iDsp); - if (printLevel) printf("Dsp length %d\n",dspHeader->GetTotalLength()); + cout << "\nMUONTRKda : Reading data from file " << inputFile << endl; - // loop over BusPatch structure - Int_t nBusPatch = dspHeader->GetBusPatchEntries(); - for(Int_t iBusPatch = 0; iBusPatch < nBusPatch; iBusPatch++) { + while (rawReader->NextEvent()) + { + if (nDateEvents >= maxDateEvents) break; + if (gAliNEvents >= maxEvents) break; + if (nDateEvents>0 && nDateEvents % 100 == 0) + cout<<"Cumulated: DATE events = " << nDateEvents << " Used events = " << gAliNEvents << endl; + + // check shutdown condition + if (daqDA_checkShutdown()) + break; + + //Skip events + while (skipEvents) + { + rawReader->NextEvent(); + skipEvents--; + } + + // starts reading + // status = monitorGetEventDynamic(&event); + // if (status < 0) { + // cout<<"EOF found"<GetType(); + gAliRunNumber = rawReader->GetRunNumber(); + + // Output log file initialisations + + if(nDateEvents==0) + { + if (gAliCommand.CompareTo("ped") == 0){ + sprintf(gAliflatFile,"%s/MUONTRKda_ped_%d.log",gAliOutFolder,gAliRunNumber); + logOutputFile=gAliflatFile; + + gAlifilcout.open(logOutputFile.Data()); + gAlifilcout<<"//=================================================" << endl; + gAlifilcout<<"// MUONTRKda for Pedestal run = " << gAliRunNumber << endl; + cout<<"\n ******** MUONTRKda for Pedestal run = " << gAliRunNumber << "\n" << endl; + } + + if (gAliCommand.Contains("gain")){ + sprintf(gAliflatFile,"%s/%s_%d_DAC_%d.log",gAliOutFolder,gAlifilenam,gAliRunNumber,injCharge); + logOutputFile=gAliflatFile; + + gAlifilcout.open(logOutputFile.Data()); + gAlifilcout<<"//=================================================" << endl; + gAlifilcout<<"// MUONTRKda for Gain run = " << gAliRunNumber << " (DAC=" << injCharge << ")" << endl; + cout<<"\n ******** MUONTRKda for Gain run = " << gAliRunNumber << " (DAC=" << injCharge << ")\n" << endl; + } + + gAlifilcout<<"//=================================================" << endl; + gAlifilcout<<"// * Date : " << gAlidate.AsString("l") << "\n" << endl; + cout<<" * Date : " << gAlidate.AsString("l") << "\n" << endl; + + } + + nDateEvents++; + + if (eventType != PHYSICS_EVENT) + continue; // for the moment + + // decoding MUON payload + // rawStream = new AliMUONRawStreamTrackerHP(rawReader); + // rawStream->DisableWarnings(); + // rawStream->EnabbleErrorLogger(); + + // First lopp over DDL's to find good events + // Error counters per event (counters in the decoding lib are for each DDL) + Bool_t eventIsErrorMessage = kFALSE; + int eventGlitchErrors = 0; + int eventParityErrors = 0; + int eventPaddingErrors = 0; + rawStream->First(); + do + { + if (rawStream->IsErrorMessage()) eventIsErrorMessage = kTRUE; + eventGlitchErrors += rawStream->GetGlitchErrors(); + eventParityErrors += rawStream->GetParityErrors(); + eventPaddingErrors += rawStream->GetPaddingErrors(); + } while(rawStream->NextDDL()); + + AliMUONRawStreamTrackerHP::AliBusPatch* busPatch; + if (!eventIsErrorMessage) + { + // Good events (no error) -> compute pedestal for all channels + rawStream->First(); + while( (busPatch = (AliMUONRawStreamTrackerHP::AliBusPatch*) rawStream->Next())) + { + for(int i = 0; i < busPatch->GetLength(); ++i) + { + if (gAliNEvents == 0) gAliNChannel++; + busPatch->GetData(i, manuId, channelId, charge); + MakePed(busPatch->GetBusPatchId(), (Int_t)manuId, (Int_t)channelId, (Int_t)charge); + } + } + gAliNEvents++; + } + else + { + // Events with errors + if (recoverParityErrors && eventParityErrors && !eventGlitchErrors&& !eventPaddingErrors) + { + // Recover parity errors -> compute pedestal for all good buspatches + if ( TEST_SYSTEM_ATTRIBUTE( rawReader->GetAttributes(), + ATTR_ORBIT_BC )) + { + gAlifilcout <<"Event recovered -> Period:"<GetEventId() ) + <<" Orbit:"<GetEventId() ) + <<" BunchCrossing:"<GetEventId() )< nbInRun:"<GetEventId() ) + <<" burstNb:"<GetEventId() ) + <<" nbInBurst:"<GetEventId() )<First(); + while( (busPatch = (AliMUONRawStreamTrackerHP::AliBusPatch*) rawStream->Next())) + { + // Check the buspatch -> if error not use it in the pedestal calculation + int errorCount = 0; + for(int i = 0; i < busPatch->GetLength(); ++i) + { + if (!busPatch->IsParityOk(i)) errorCount++; + } + if (!errorCount) + { + // Good buspatch + for(int i = 0; i < busPatch->GetLength(); ++i) + { + if (gAliNEvents == 0) gAliNChannel++; + busPatch->GetData(i, manuId, channelId, charge); + // if (busPatch->GetBusPatchId()==1719 && manuId == 1 && channelId == 0) cout <<"Recovered charge "<GetBusPatchId(), (Int_t)manuId, (Int_t)channelId, (Int_t)charge); + } + } + else + { + char bpname[256]; + ErrorCounter* errorCounter; + // Bad buspatch -> not used (just print) + gAlifilcout<<"bpId "<GetBusPatchId()<<" words "<GetLength() + <<" parity errors "<GetBusPatchId()); + if (!(errorCounter = (ErrorCounter*)gAliErrorBuspatchTable->FindObject(bpname))) + { + // New buspatch + errorCounter = new ErrorCounter(busPatch->GetBusPatchId()); + errorCounter->SetName(bpname); + gAliErrorBuspatchTable->Add(errorCounter); + } + else + { + // Existing buspatch + errorCounter->Increment(); + } + // errorCounter->Print(); + } // end of if (!errorCount) + } // end of while( (busPatch = (AliMUONRawStreamTrackerHP ... + gAliNEvents++; + gAliNEventsRecovered++; + } //end of if (recoverParityErrors && eventParityErrors && !eventGlitchErrors&& !eventPaddingErrors) + else + { + // Fatal errors reject the event + if ( TEST_SYSTEM_ATTRIBUTE( rawReader->GetAttributes(), + ATTR_ORBIT_BC )) + { + gAlifilcout <<"Event rejected -> Period:"<GetEventId() ) + <<" Orbit:"<GetEventId() ) + <<" BunchCrossing:"<GetEventId() )< nbInRun:"<GetEventId() ) + <<" burstNb:"<GetEventId() ) + <<" nbInBurst:"<GetEventId() )<GetGlitchErrors() && !rawStream->GetPaddingErrors() ... + gAlifilcout<<"Number of errors : Glitch "<IsErrorMessage()) + + if (eventGlitchErrors) gGlitchErrors++; + if (eventParityErrors) gParityErrors++; + if (eventPaddingErrors) gPaddingErrors++; + + } // while (rawReader->NextEvent()) + delete rawReader; + delete rawStream; + + + if (gAliCommand.CompareTo("ped") == 0) + { + sprintf(gAliflatFile,"MUONTRKda_ped_%d.ped",gAliRunNumber); + if(gAliOutputFile.IsNull())gAliOutputFile=gAliflatFile; + MakePedStore(gAliOutputFile); + } - busStruct = dspHeader->GetBusPatchEntry(iBusPatch); + // option gain -> update root file with pedestal results + // gain + create -> recreate root file + // gain + comp -> update root file and compute gain parameters - if (printLevel) printf("busPatchId %d", busStruct->GetBusPatchId()); - if (printLevel) printf(" BlockId %d", busStruct->GetBlockId()); - if (printLevel) printf(" DspId %d\n", busStruct->GetDspId()); + if (gAliCommand.Contains("gain")) + { + MakePedStoreForGain(injCharge); + } - Int_t busPatchId = busStruct->GetBusPatchId(); - // loop over data - Int_t dataSize = busStruct->GetLength(); - for (Int_t iData = 0; iData < dataSize; iData++) { + delete gAliPedestalStore; - if (nEvents == 1) - nChannel++; - Int_t manuId = busStruct->GetManuId(iData); - Int_t channelId = busStruct->GetChannelId(iData); - Int_t charge = busStruct->GetCharge(iData); + delete minuitFit; + TVirtualFitter::SetFitter(0); - if (printLevel) printf("manuId: %d, channelId: %d charge: %d\n", manuId, - channelId, charge); + timers.Stop(); - MakePed(busPatchId, manuId, channelId, charge); - - } // iData - } // iBusPatch - } // iDsp - } // iBlock + cout << "\nMUONTRKda : Nb of DATE events = " << nDateEvents << endl; + cout << "MUONTRKda : Nb of Glitch errors = " << gGlitchErrors << endl; + cout << "MUONTRKda : Nb of Parity errors = " << gParityErrors << endl; + cout << "MUONTRKda : Nb of Padding errors = " << gPaddingErrors << endl; + cout << "MUONTRKda : Nb of events recovered = " << gAliNEventsRecovered<< endl; + cout << "MUONTRKda : Nb of events without errors = " << gAliNEvents-gAliNEventsRecovered<< endl; + cout << "MUONTRKda : Nb of events used = " << gAliNEvents << endl; - totalParityErr += rawStream->GetPayLoad()->GetParityErrors(); - totalGlitchErr += rawStream->GetPayLoad()->GetGlitchErrors(); - } // NextDDL - - delete rawReader; - delete rawStream; + gAlifilcout << "\nMUONTRKda : Nb of DATE events = " << nDateEvents << endl; + gAlifilcout << "MUONTRKda : Nb of Glitch errors = " << gGlitchErrors << endl; + gAlifilcout << "MUONTRKda : Nb of Parity errors = " << gParityErrors << endl; + gAlifilcout << "MUONTRKda : Nb of Padding errors = " << gPaddingErrors << endl; + gAlifilcout << "MUONTRKda : Nb of events recovered = " << gAliNEventsRecovered<< endl; + gAlifilcout << "MUONTRKda : Nb of events without errors = " << gAliNEvents-gAliNEventsRecovered<< endl; + gAlifilcout << "MUONTRKda : Nb of events used = " << gAliNEvents << endl; - } // while (1) + if (gAliCommand.CompareTo("ped") == 0) + { + cout << "\nMUONTRKda : Output logfile : " << logOutputFile << endl; + cout << "MUONTRKda : Pedestal Histo file : " << gAliHistoFileName << endl; + cout << "MUONTRKda : Pedestal file (to SHUTTLE) : " << gAliOutputFile << endl; + } + else + { + cout << "\nMUONTRKda : Output logfile : " << logOutputFile << endl; + cout << "MUONTRKda : DAC data (root file) : " << gAliHistoFileNamegain << endl; + cout << "MUONTRKda : Dummy file (to SHUTTLE) : " << gAliOutputFile << endl; + } - nEvents -= totalGlitchErr; // skipped event cos of Glitch errors + } + // Compute gain parameters - cout << " Total number of parity errors in the Run: " << totalParityErr << endl; - cout << " Total number of glitch errors in the Run: " << totalGlitchErr << endl; - if (command.CompareTo("ped") == 0) - MakePedStore(flatOutputFile); + if (gAliCommand.Contains("comp")) + { + gAliOutputFile=""; - // option gain -> update root file with pedestal results - // gain + create -> recreate root file - // gain + comp -> update root file and compute gain parameters + MakeGainStore(); + cout << "MUONTRKda : Gain file (to SHUTTLE) : " << gAliOutputFile << endl; + } - if (command.Contains("gain")) - MakePedStoreForGain(); - - if (command.Contains("comp")) - MakeGainStore(flatOutputFile); - - delete pedestalStore; + if(fes) // Store IN FES + { + printf("\n ***** STORE FILE in FES ****** \n"); - timers.Stop(); + // be sure that env variable DAQDALIB_PATH is set in script file + // gSystem->Setenv("DAQDALIB_PATH", "$DATE_SITE/infoLogger"); - if (!(crocusConfigFile.IsNull())) - cout << "MUONTRKda : CROCUS command file generated : " << crocusOutputFile.Data() << endl; - else - cout << "MUONTRKda : WARNING no CROCUS command file generated" << endl; + if (!gAliOutputFile.IsNull()) + { + if (gAliCommand.CompareTo("ped") == 0) + status = daqDA_FES_storeFile(gAliOutputFile.Data(),"PEDESTALS"); + else + status = daqDA_FES_storeFile(gAliOutputFile.Data(),"GAINS"); + + if (status) + { + printf(" Failed to export file : %d\n",status); + } + else if(gAliPrintLevel) printf(" %s successfully exported to FES \n",gAliOutputFile.Data()); + } + } - 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; + gAlifilcout.close(); - printf("Execution time : R:%7.2fs C:%7.2fs\n", timers.RealTime(), timers.CpuTime()); + printf("\nExecution time : R:%7.2fs C:%7.2fs\n", timers.RealTime(), timers.CpuTime()); return status; }