[u/mrichter/AliRoot.git] / MUON / CreateWeightedRejectList.C

/**************************************************************************
 * 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$

///
/// \ingroup macros
///
/// \file CreateWeightedRejectList.C
///
/// \brief Macro test for the creation of a RejectList object (for simulations) in the OCDB, for the MUON Tracke 
///
/// Usage:
///
/// root[0] .L CreateWeightedRejectList.C+
/// root[1] CreateWeightedRejectList("runlist.txt","local://$HOME/myLocalOCDB");
///
/// where runlist.txt has 2 integers per line = "run nevents"
/// where nevents is the number of events where there's (for instance) a CMUS1B trigger
///
///
/// Assuming the file coming from the Export feature of the logbook is runlist.logbook.txt =>
///
/// 115521;PHYSICS_1;1357;0.9
/// 115516;PHYSICS_1;944;0.9
///
/// The awk command below will output the needed format for this macro :
///
/// awk '{split ($0,a,";"); print a[1] " " a[3];}' runlist.logbook.txt =>
///
/// 115521 1357
/// 115516 944
///
/// \author Matthieu Lenhardt, Subatech
///

#include "AliCDBEntry.h"
#include "AliCDBManager.h"
#include "AliDAQ.h"
#include "AliDCSValue.h"
#include "AliMUONCDB.h"
#include "AliMUONCalibParamNI.h"
#include "AliMUONCalibrationData.h"
#include "AliMUONPadStatusMaker.h"
#include "AliMUONPadStatusMapMaker.h"
#include "AliMUONRecoParam.h"
#include "AliMUONRejectList.h"
#include "AliMUONRejectList.h"
#include "AliMUONTrackerData.h"
#include "AliMUONVCalibParam.h"
#include "AliMpBusPatch.h"
#include "AliMpCDB.h"
#include "AliMpConstants.h"
#include "AliMpDCSNamer.h"
#include "AliMpDDLStore.h"
#include "AliMpDEIterator.h"
#include "AliMpDEManager.h"
#include "AliMpDetElement.h"
#include "AliSysInfo.h"
#include "TEnv.h"
#include <Riostream.h>
#include <TFile.h>
#include <TGrid.h>
#include <TMap.h>
#include <TObjArray.h>
#include <TObject.h>
#include <TTree.h>
#include <vector>

class RunInfo
{
public:
  RunInfo(Int_t run, Long64_t nevents) : fNumber(run), fNevents(nevents) { }

  Int_t Number() const { return fNumber; }

  Long64_t Nevents() const { return fNevents; }

  void SetNevents(Long64_t n) { fNevents = n; }

private:
  Int_t fNumber;
  Long64_t fNevents;
};

std::ostream& operator<<(std::ostream& out, const RunInfo& run)
{
  out << Form("RUN %09d NEVENTS %10d",run.Number(),run.Nevents());
  return out;
}

AliMUONRejectList* CheckDE_BP_ManuPedestals(Int_t rejectMask, AliMUONPadStatusMaker& status);

Int_t AddEventsSingleRun(Int_t index, 
                         Int_t run_number, 
                         AliMUONRejectList& rejectedEvents);

std::vector<RunInfo> runs;

//______________________________________________________________________________
bool CreateWeightedRejectList(const char* runlistfile="runlist.txt", 
                              const char* rejectListPath="local://$HOME/OCDB")
{
  /// Create a weighted RejectList for the runs included in the run list
  /// The cuts are the same that the one applied in the RecoParam used to create the ESDs.
  ///
  /// \param runlistfile : an ASCII file where each line is a pair "run nevents"
  /// \param rejectListPath : set this to the folder where the RejectList will be created (must have right access to it, of course)
  
  TGrid::Connect("alien://");
  
  ifstream in(runlistfile);
  int run, nevents;
  
  while ( in >> run >> nevents ) 
  {
    runs.push_back(RunInfo(run,nevents));
  };
  
  for ( std::vector<RunInfo>::size_type i = 0; i < runs.size(); ++i )
  {
    cout << runs[i] << endl;
  }
  
  if (runs.empty()) return false;
  
  Int_t firstRun = runs[0].Number();
  Int_t lastRun = runs[runs.size()-1].Number();
  
  const char* rawOCDB = "raw://";
  
  AliCDBManager* manager = AliCDBManager::Instance();

  manager->SetCacheFlag(kFALSE);
  
  manager->SetDefaultStorage(rawOCDB);

  manager->SetSpecificStorage("MUON/Calib/RejectList", rejectListPath);
  
  AliMUONRejectList weightedRejectList;
  AliMUONRejectList rejectedEvents;
  
  Long64_t nEventsTotal = 0;
  
  for (std::vector<RunInfo>::size_type ii = 0; ii < runs.size(); ++ii)
  {
    Int_t runNumber = runs[ii].Number();
    AliSysInfo::Instance()->AddStamp(Form("RUN%d",runNumber));
    nEventsTotal += AddEventsSingleRun(ii, runNumber, rejectedEvents);
  }  
  
  AliMpDEIterator DEiter;
  
  for (DEiter.First(); !DEiter.IsDone(); DEiter.Next())
    if (DEiter.CurrentDEId() < 1100)
    {
      Int_t DEid = DEiter.CurrentDEId();
      
      Double_t nBadEventsDE = rejectedEvents.DetectionElementProbability(DEid);
      Double_t nEventsTotalDE = static_cast<Double_t>(nEventsTotal);
      Float_t probaDE = 0.0;
      if (nEventsTotalDE != 0.0)
        probaDE = nBadEventsDE / nEventsTotalDE;
      weightedRejectList.SetDetectionElementProbability(DEid, probaDE);
      
      Int_t nBusPatches = DEiter.CurrentDE()->GetNofBusPatches();
      for (Int_t ii = 0; ii < nBusPatches; ii++)
      {
        Int_t BPid = DEiter.CurrentDE()->GetBusPatchId(ii);
        
        Double_t nBadEventsBP = rejectedEvents.BusPatchProbability(BPid);
        Double_t nEventsTotalBP = nEventsTotalDE - nBadEventsDE;
        Float_t probaBP = 0.0;
        if (nEventsTotalBP != 0.0)
          probaBP = nBadEventsBP / nEventsTotalBP;
        weightedRejectList.SetBusPatchProbability(BPid, probaBP);
        
        Int_t nManus = AliMpDDLStore::Instance(kFALSE)->GetBusPatch(BPid, kFALSE)->GetNofManus();
        for (Int_t jj = 0; jj < nManus; jj++)
	      {
          Int_t manuId = AliMpDDLStore::Instance(kFALSE)->GetBusPatch(BPid, kFALSE)->GetManuId(jj);
          
          Double_t nBadEventsManu = rejectedEvents.ManuProbability(DEid, manuId);
          Double_t nEventsTotalManu = nEventsTotalBP - nBadEventsBP;
          Float_t probaManu = 0.0;
          if (nEventsTotalManu != 0.0)
            probaManu = nBadEventsManu / nEventsTotalManu;
          weightedRejectList.SetManuProbability(DEid, manuId, probaManu);
          
          for (Int_t channel = 0; channel < AliMpConstants::ManuNofChannels(); channel++)
          {
            Double_t nBadEventsChannel = rejectedEvents.ChannelProbability(DEid, manuId, channel);
            Double_t nEventsTotalChannel = nEventsTotalManu - nBadEventsManu;
            Float_t probaChannel = 0.0;
            if (nEventsTotalChannel != 0.0)
            {
              probaChannel = nBadEventsChannel / nEventsTotalChannel;
            }
            
            weightedRejectList.SetChannelProbability(DEid, manuId, channel, probaChannel);
          }
	      }
      }
    }
  
  
  AliMUONCDB::WriteToCDB(&weightedRejectList, "MUON/Calib/RejectList",
                         firstRun , lastRun, 
                         "Weighted reject List for MCH, for simulations only", 
                         "Matthieu Lenhardt");
  
  return true;
}

//____________________________________________________________________________________________
AliMUONRejectList* CheckDE_BP_ManuPedestals(Int_t rejectMask, AliMUONPadStatusMaker& status)
{
  AliMUONRejectList* nBadChannels = new AliMUONRejectList;
  
  AliMpDEIterator DEiter;
  for (DEiter.First(); !DEiter.IsDone(); DEiter.Next())
    if (DEiter.CurrentDEId() < 1100)
    {
      Int_t DEid = DEiter.CurrentDEId();
      Int_t nBusPatches = DEiter.CurrentDE()->GetNofBusPatches();
      for (Int_t ii = 0; ii < nBusPatches; ii++)
      {
        Int_t BPid = DEiter.CurrentDE()->GetBusPatchId(ii);
        
        Int_t nManus = AliMpDDLStore::Instance(kFALSE)->GetBusPatch(BPid, kFALSE)->GetNofManus();
        for (Int_t jj = 0; jj < nManus; jj++)
	      {
          Int_t manuId = AliMpDDLStore::Instance(kFALSE)->GetBusPatch(BPid, kFALSE)->GetManuId(jj);
          
          Int_t nChannels = DEiter.CurrentDE()->NofChannelsInManu(manuId);
          for (Int_t channel = 0; channel < nChannels; channel++)
          {
            Int_t padStatus = status.PadStatus(DEid, manuId, channel);
            if ((rejectMask == 0 && padStatus != 0) || (padStatus & rejectMask) != 0)
            {
              nBadChannels->SetDetectionElementProbability(DEid, nBadChannels->DetectionElementProbability(DEid) + 1.0);
              nBadChannels->SetBusPatchProbability(BPid, nBadChannels->BusPatchProbability(BPid) + 1.0);
              nBadChannels->SetManuProbability(DEid, manuId, nBadChannels->ManuProbability(DEid, manuId) + 1.0);
              nBadChannels->SetChannelProbability(DEid, manuId, channel, 1.0);
            }
          }
	      }
      }
    }
  
  return nBadChannels;
}


//____________________________________________________________________________________________
Int_t AddEventsSingleRun(Int_t index, Int_t runNumber, AliMUONRejectList& rejectedEvents)
{
  AliCDBManager::Instance()->SetRun(runNumber);
  
  AliMpCDB::LoadAll();

  Double_t maxBadChannels = 80.0;
  
  AliMUONRecoParam *recoParam = AliMUONCDB::LoadRecoParam();
  AliMUONCalibrationData calibrationData(runNumber);
  AliMUONPadStatusMaker status(calibrationData);
  
  status.SetLimits(*recoParam);
  
  Long64_t nEvents = runs[index].Nevents();
  
  Int_t rejectMask = recoParam->PadGoodnessMask();
  
  // Functions to compute the number of bad channels at the DE, BP and Manu level
  AliMUONRejectList* nBadChannels = CheckDE_BP_ManuPedestals(rejectMask, status);
  
  AliMpDEIterator DEiter;
  for (DEiter.First(); !DEiter.IsDone(); DEiter.Next())
    if (DEiter.CurrentDEId() < 1100)
    {
      Int_t DEid = DEiter.CurrentDEId();
      Int_t nBusPatches = DEiter.CurrentDE()->GetNofBusPatches();
      
      Double_t nChannelsDE = static_cast<Double_t>(AliMpDDLStore::Instance(kFALSE)->GetDetElement(DEid)->NofChannels());
      Double_t nBadChannelsDE = nBadChannels->DetectionElementProbability(DEid);
      Double_t ratioBadChannelsDE = 0.0;
      if (nChannelsDE != 0.0)
        ratioBadChannelsDE = 100.0 * nBadChannelsDE / nChannelsDE;
      
      
      Bool_t goodDE = kTRUE;
      if (ratioBadChannelsDE >= maxBadChannels)
      {
        Double_t newRejectedEventsDE = rejectedEvents.DetectionElementProbability(DEid) + static_cast<Double_t>(nEvents);
        rejectedEvents.SetDetectionElementProbability(DEid, newRejectedEventsDE);
        goodDE = kFALSE;
      }
	  	
      if (goodDE)      
        for (Int_t ii = 0; ii < nBusPatches; ii++)
        {
          Int_t BPid = DEiter.CurrentDE()->GetBusPatchId(ii);
          Int_t nManus = AliMpDDLStore::Instance(kFALSE)->GetBusPatch(BPid, kFALSE)->GetNofManus();
          
          Double_t nChannelsBP = static_cast<Double_t>(nManus * AliMpConstants::ManuNofChannels());
          Double_t nBadChannelsBP = nBadChannels->BusPatchProbability(BPid);
          Double_t ratioBadChannelsBP = 100.0 * nBadChannelsBP / nChannelsBP;
          
          Bool_t goodBP = kTRUE;
          
          if (goodBP)
            if (ratioBadChannelsBP >= maxBadChannels)
            {
              Double_t newRejectedEventsBP = rejectedEvents.BusPatchProbability(BPid) + static_cast<Double_t>(nEvents);
              rejectedEvents.SetBusPatchProbability(BPid, newRejectedEventsBP);
              goodBP = kFALSE;
            }
          
          if (goodBP)
            for (Int_t jj = 0; jj < nManus; jj++)
            {
              Int_t manuId = AliMpDDLStore::Instance(kFALSE)->GetBusPatch(BPid, kFALSE)->GetManuId(jj);
              
              Double_t nChannelsManu = DEiter.CurrentDE()->NofChannelsInManu(manuId);;
              Double_t nBadChannelsManu = nBadChannels->ManuProbability(DEid, manuId);
              Double_t ratioBadChannelsManu = 100.0 * nBadChannelsManu / nChannelsManu;
              
              Bool_t goodManu = kTRUE;
              
              if (ratioBadChannelsManu >= maxBadChannels)
              {
                Double_t newRejectedEventsManu = rejectedEvents.ManuProbability(DEid, manuId) + static_cast<Double_t>(nEvents);
                rejectedEvents.SetManuProbability(DEid, manuId, newRejectedEventsManu);
                goodManu = kFALSE;
              }
              
              
              Int_t nChannels = DEiter.CurrentDE()->NofChannelsInManu(manuId);
              
              if (goodManu)
                for (Int_t channel = 0; channel < nChannels; channel++)
                  if (nBadChannels->ChannelProbability(DEid, manuId, channel) > 0.0)
                  {
                    Double_t newRejectedEventsChannel = rejectedEvents.ChannelProbability(DEid, manuId, channel) + static_cast<Double_t>(nEvents);
                    rejectedEvents.SetChannelProbability(DEid, manuId, channel, newRejectedEventsChannel);
                  } // end of Channels loop
            } // end of MANUs loop	    
        } // end of BPs loop
    } // end of DEs loop
  
  delete nBadChannels;
  
  return nEvents;
}

//______________________________________________________________________________
//
//
// The code below might be used someday to retrieve the number of events
// with at least one muon directly from a runlist, instead of requiring
// as input both the run numbers and the number of events.
//
// This approach though is not without its problems :
//
// - even if we are using tags, as the tags are chunk-based (for the moment at
//   least), it's pretty slow to loop on all the chunks (it is becoming usual
//   to get runs with some 600 chunks or so...)
// - if part of a run is not reconstructed (for whatever reason) with the pass
//   we use to compute the number of events with at least one muon, we might
//   underestimate its weight in the final rejectlist for the next passes (
//   if for the next passes the full run is correctly reconstructed for instance)
// - a muon in the tag = either a trigger track or a tracker track, i.e. we don't
//   know for sure it's a good track...
//
// Given all those reasons, we for the moment use a dumb approach : 
//
// - go to alice-logbook.cern.ch . Select the runs you want and the trigger class
//   you want (most probably CMUS1B), and export the list as (run,nevents)
// - run this macro with this (run,nevents) list
//
//______________________________________________________________________________

#if 0

#include "TChain.h"
#include "AliDetectorTagCuts.h"
#include "AliEventTagCuts.h"
#include "AliLHCTagCuts.h"
#include "AliRunTagCuts.h"
#include "AliTagAnalysis.h"
#include "TGridResult.h"

//______________________________________________________________________________
bool CreateXMLCollectionFromRunList(const char* collectionName, 
                                    const char* runlist,
                                    const char* type,
                                    int passNumber)
{
  /// From a list of run, create a collection with only events containing at least one muon
  /// \param collectionName : output name of the collection (without extension, which will be .xml)
  /// \param runlist : text file containing one integer per line = one run number per line
  /// \param type : files to consider, either ESD or AD
  /// \param passNumber : 1 or 2 most probably (to distinguish which reco pass is used)
  
  if (!gGrid) TGrid::Connect("alien://");
  if (!gGrid) 
  {
    return 0x0;
  }
  
  TString stype(type);
  stype.ToUpper();
  
  if ( stype != "ESD" && stype != "AOD" )
  {
    cout << "Only ESD or AOD type supported" << endl;
    return false;
  }
  
  ifstream in(gSystem->ExpandPathName(runlist));
  Int_t runNumber;
  Int_t ntagfiles(0);
  
  AliTagAnalysis tagAnalysis("ESD");
  
  while ( in >> runNumber ) 
  {
    TGridResult *res = gGrid->Query("/alice/data",
                                    Form("%09d/%ss/pass%d/*%d*/Run%d.Event*.ESD.tag.root",
                                         runNumber,stype.Data(),passNumber,runNumber,runNumber));
    Int_t nFiles = res->GetEntries();
    if (!nFiles) 
    {
      continue;
    }
    
    for (Int_t i = 0; i < nFiles; ++i) 
    {
      TString filename = res->GetKey(i, "turl");
      if(filename == "") continue;
      tagAnalysis.AddTagsFile(filename.Data(),kFALSE);
      ++ntagfiles;
    }
    delete res;
  }
  
  cout << ntagfiles << " tag files added" << endl;
  
  AliRunTagCuts runCuts;
  AliEventTagCuts eventCuts;
  AliLHCTagCuts lhcCuts;
  AliDetectorTagCuts detCuts;
  
  eventCuts.SetNMuonRange(1,99999);
  
  return tagAnalysis.CreateXMLCollection(collectionName,&runCuts,&lhcCuts,&detCuts,&eventCuts);  
  
}

//______________________________________________________________________________
TChain* CreateChainFromXMLCollection(const char* collectionName, const char* type)
{
  /// Create a chain from an XML collection file.
  
  if (!gGrid) TGrid::Connect("alien://");
  if (!gGrid) 
  {
    return 0x0;
  }
  
  TString stype(type);
  stype.ToUpper();
  
  if ( stype != "ESD" && stype != "AOD" )
  {
    cout << "Only ESD or AOD type supported" << endl;
    return 0x0;
  }
  
  AliTagAnalysis tagAnalysis(stype.Data());
  
  return tagAnalysis.CreateChainFromCollection(collectionName,stype=="ESD" ? "esdTree":"aodTree");
}

#endif
Commit	Line	Data
f0938a5b	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	// $Id$
	17
	18	///
	19	/// \ingroup macros
	20	///
	21	/// \file CreateWeightedRejectList.C
	22	///
	23	/// \brief Macro test for the creation of a RejectList object (for simulations) in the OCDB, for the MUON Tracke
	24	///
	25	/// Usage:
	26	///
	27	/// root[0] .L CreateWeightedRejectList.C+
	28	/// root[1] CreateWeightedRejectList("runlist.txt","local://$HOME/myLocalOCDB");
	29	///
	30	/// where runlist.txt has 2 integers per line = "run nevents"
	31	/// where nevents is the number of events where there's (for instance) a CMUS1B trigger
	32	///
951dd287	33	///
	34	/// Assuming the file coming from the Export feature of the logbook is runlist.logbook.txt =>
	35	///
	36	/// 115521;PHYSICS_1;1357;0.9
	37	/// 115516;PHYSICS_1;944;0.9
	38	///
	39	/// The awk command below will output the needed format for this macro :
	40	///
	41	/// awk '{split ($0,a,";"); print a[1] " " a[3];}' runlist.logbook.txt =>
	42	///
	43	/// 115521 1357
	44	/// 115516 944
	45	///
f0938a5b	46	/// \author Matthieu Lenhardt, Subatech
	47	///
	48
	49	#include "AliCDBEntry.h"
	50	#include "AliCDBManager.h"
	51	#include "AliDAQ.h"
	52	#include "AliDCSValue.h"
	53	#include "AliMUONCDB.h"
	54	#include "AliMUONCalibParamNI.h"
	55	#include "AliMUONCalibrationData.h"
	56	#include "AliMUONPadStatusMaker.h"
	57	#include "AliMUONPadStatusMapMaker.h"
	58	#include "AliMUONRecoParam.h"
	59	#include "AliMUONRejectList.h"
	60	#include "AliMUONRejectList.h"
	61	#include "AliMUONTrackerData.h"
	62	#include "AliMUONVCalibParam.h"
	63	#include "AliMpBusPatch.h"
	64	#include "AliMpCDB.h"
	65	#include "AliMpConstants.h"
	66	#include "AliMpDCSNamer.h"
	67	#include "AliMpDDLStore.h"
	68	#include "AliMpDEIterator.h"
	69	#include "AliMpDEManager.h"
	70	#include "AliMpDetElement.h"
	71	#include "AliSysInfo.h"
	72	#include "TEnv.h"
	73	#include <Riostream.h>
	74	#include <TFile.h>
	75	#include <TGrid.h>
	76	#include <TMap.h>
	77	#include <TObjArray.h>
	78	#include <TObject.h>
	79	#include <TTree.h>
	80	#include <vector>
	81
	82	class RunInfo
	83	{
	84	public:
	85	RunInfo(Int_t run, Long64_t nevents) : fNumber(run), fNevents(nevents) { }
	86
	87	Int_t Number() const { return fNumber; }
	88
	89	Long64_t Nevents() const { return fNevents; }
	90
	91	void SetNevents(Long64_t n) { fNevents = n; }
	92
	93	private:
	94	Int_t fNumber;
	95	Long64_t fNevents;
	96	};
	97
	98	std::ostream& operator<<(std::ostream& out, const RunInfo& run)
	99	{
	100	out << Form("RUN %09d NEVENTS %10d",run.Number(),run.Nevents());
	101	return out;
	102	}
	103
	104	AliMUONRejectList* CheckDE_BP_ManuPedestals(Int_t rejectMask, AliMUONPadStatusMaker& status);
	105
	106	Int_t AddEventsSingleRun(Int_t index,
	107	Int_t run_number,
	108	AliMUONRejectList& rejectedEvents);
	109
110	std::vector<RunInfo> runs;
111
112	//______________________________________________________________________________
113	bool CreateWeightedRejectList(const char* runlistfile="runlist.txt",
114	const char* rejectListPath="local://$HOME/OCDB")
115	{
116	/// Create a weighted RejectList for the runs included in the run list
117	/// The cuts are the same that the one applied in the RecoParam used to create the ESDs.
118	///
119	/// \param runlistfile : an ASCII file where each line is a pair "run nevents"
120	/// \param rejectListPath : set this to the folder where the RejectList will be created (must have right access to it, of course)
121
122	TGrid::Connect("alien://");
123
124	ifstream in(runlistfile);
125	int run, nevents;
126
127	while ( in >> run >> nevents )
128	{
129	runs.push_back(RunInfo(run,nevents));
130	};
131
132	for ( std::vector<RunInfo>::size_type i = 0; i < runs.size(); ++i )
133	{
134	cout << runs[i] << endl;
135	}
136
137	if (runs.empty()) return false;
138
139	Int_t firstRun = runs[0].Number();
140	Int_t lastRun = runs[runs.size()-1].Number();
141
142	const char* rawOCDB = "raw://";
143
144	AliCDBManager* manager = AliCDBManager::Instance();
145
146	manager->SetCacheFlag(kFALSE);
147
148	manager->SetDefaultStorage(rawOCDB);
149
150	manager->SetSpecificStorage("MUON/Calib/RejectList", rejectListPath);
151
152	AliMUONRejectList weightedRejectList;
153	AliMUONRejectList rejectedEvents;
154
155	Long64_t nEventsTotal = 0;
156
157	for (std::vector<RunInfo>::size_type ii = 0; ii < runs.size(); ++ii)
158	{
159	Int_t runNumber = runs[ii].Number();
160	AliSysInfo::Instance()->AddStamp(Form("RUN%d",runNumber));
161	nEventsTotal += AddEventsSingleRun(ii, runNumber, rejectedEvents);
162	}
163
164	AliMpDEIterator DEiter;
165
166	for (DEiter.First(); !DEiter.IsDone(); DEiter.Next())
167	if (DEiter.CurrentDEId() < 1100)
168	{
169	Int_t DEid = DEiter.CurrentDEId();
170
171	Double_t nBadEventsDE = rejectedEvents.DetectionElementProbability(DEid);
172	Double_t nEventsTotalDE = static_cast<Double_t>(nEventsTotal);
173	Float_t probaDE = 0.0;
174	if (nEventsTotalDE != 0.0)
175	probaDE = nBadEventsDE / nEventsTotalDE;
176	weightedRejectList.SetDetectionElementProbability(DEid, probaDE);
177
178	Int_t nBusPatches = DEiter.CurrentDE()->GetNofBusPatches();
179	for (Int_t ii = 0; ii < nBusPatches; ii++)
180	{
181	Int_t BPid = DEiter.CurrentDE()->GetBusPatchId(ii);
182
183	Double_t nBadEventsBP = rejectedEvents.BusPatchProbability(BPid);
184	Double_t nEventsTotalBP = nEventsTotalDE - nBadEventsDE;
185	Float_t probaBP = 0.0;
186	if (nEventsTotalBP != 0.0)
187	probaBP = nBadEventsBP / nEventsTotalBP;
188	weightedRejectList.SetBusPatchProbability(BPid, probaBP);
189
190	Int_t nManus = AliMpDDLStore::Instance(kFALSE)->GetBusPatch(BPid, kFALSE)->GetNofManus();
191	for (Int_t jj = 0; jj < nManus; jj++)
192	{
193	Int_t manuId = AliMpDDLStore::Instance(kFALSE)->GetBusPatch(BPid, kFALSE)->GetManuId(jj);
194
195	Double_t nBadEventsManu = rejectedEvents.ManuProbability(DEid, manuId);
196	Double_t nEventsTotalManu = nEventsTotalBP - nBadEventsBP;
197	Float_t probaManu = 0.0;
198	if (nEventsTotalManu != 0.0)
199	probaManu = nBadEventsManu / nEventsTotalManu;
200	weightedRejectList.SetManuProbability(DEid, manuId, probaManu);
201
202	for (Int_t channel = 0; channel < AliMpConstants::ManuNofChannels(); channel++)
203	{
204	Double_t nBadEventsChannel = rejectedEvents.ChannelProbability(DEid, manuId, channel);
205	Double_t nEventsTotalChannel = nEventsTotalManu - nBadEventsManu;
206	Float_t probaChannel = 0.0;
207	if (nEventsTotalChannel != 0.0)
208	{
209	probaChannel = nBadEventsChannel / nEventsTotalChannel;
210	}
211
212	weightedRejectList.SetChannelProbability(DEid, manuId, channel, probaChannel);
213	}
214	}
215	}
216	}
217
218
219	AliMUONCDB::WriteToCDB(&weightedRejectList, "MUON/Calib/RejectList",
220	firstRun , lastRun,
221	"Weighted reject List for MCH, for simulations only",
222	"Matthieu Lenhardt");
223
224	return true;
225	}
226
227	//____________________________________________________________________________________________
228	AliMUONRejectList* CheckDE_BP_ManuPedestals(Int_t rejectMask, AliMUONPadStatusMaker& status)
229	{
230	AliMUONRejectList* nBadChannels = new AliMUONRejectList;
231
232	AliMpDEIterator DEiter;
233	for (DEiter.First(); !DEiter.IsDone(); DEiter.Next())
234	if (DEiter.CurrentDEId() < 1100)
235	{
236	Int_t DEid = DEiter.CurrentDEId();
237	Int_t nBusPatches = DEiter.CurrentDE()->GetNofBusPatches();
238	for (Int_t ii = 0; ii < nBusPatches; ii++)
239	{
240	Int_t BPid = DEiter.CurrentDE()->GetBusPatchId(ii);
241
242	Int_t nManus = AliMpDDLStore::Instance(kFALSE)->GetBusPatch(BPid, kFALSE)->GetNofManus();
243	for (Int_t jj = 0; jj < nManus; jj++)
244	{
245	Int_t manuId = AliMpDDLStore::Instance(kFALSE)->GetBusPatch(BPid, kFALSE)->GetManuId(jj);
246
247	Int_t nChannels = DEiter.CurrentDE()->NofChannelsInManu(manuId);
248	for (Int_t channel = 0; channel < nChannels; channel++)
249	{
250	Int_t padStatus = status.PadStatus(DEid, manuId, channel);
251	if ((rejectMask == 0 && padStatus != 0) \|\| (padStatus & rejectMask) != 0)
252	{
253	nBadChannels->SetDetectionElementProbability(DEid, nBadChannels->DetectionElementProbability(DEid) + 1.0);
254	nBadChannels->SetBusPatchProbability(BPid, nBadChannels->BusPatchProbability(BPid) + 1.0);
255	nBadChannels->SetManuProbability(DEid, manuId, nBadChannels->ManuProbability(DEid, manuId) + 1.0);
256	nBadChannels->SetChannelProbability(DEid, manuId, channel, 1.0);
257	}
258	}
259	}
260	}
261	}
262
263	return nBadChannels;
264	}
265
266
267	//____________________________________________________________________________________________
268	Int_t AddEventsSingleRun(Int_t index, Int_t runNumber, AliMUONRejectList& rejectedEvents)
269	{
270	AliCDBManager::Instance()->SetRun(runNumber);
271
272	AliMpCDB::LoadAll();
273
274	Double_t maxBadChannels = 80.0;
275
276	AliMUONRecoParam *recoParam = AliMUONCDB::LoadRecoParam();
277	AliMUONCalibrationData calibrationData(runNumber);
278	AliMUONPadStatusMaker status(calibrationData);
279
beffc4d5	280	status.SetLimits(*recoParam);
beffc4d5	281
f0938a5b	282	Long64_t nEvents = runs[index].Nevents();
	283
	284	Int_t rejectMask = recoParam->PadGoodnessMask();
	285
	286	// Functions to compute the number of bad channels at the DE, BP and Manu level
	287	AliMUONRejectList* nBadChannels = CheckDE_BP_ManuPedestals(rejectMask, status);
	288
	289	AliMpDEIterator DEiter;
	290	for (DEiter.First(); !DEiter.IsDone(); DEiter.Next())
	291	if (DEiter.CurrentDEId() < 1100)
	292	{
	293	Int_t DEid = DEiter.CurrentDEId();
	294	Int_t nBusPatches = DEiter.CurrentDE()->GetNofBusPatches();
	295
	296	Double_t nChannelsDE = static_cast<Double_t>(AliMpDDLStore::Instance(kFALSE)->GetDetElement(DEid)->NofChannels());
	297	Double_t nBadChannelsDE = nBadChannels->DetectionElementProbability(DEid);
	298	Double_t ratioBadChannelsDE = 0.0;
	299	if (nChannelsDE != 0.0)
	300	ratioBadChannelsDE = 100.0 * nBadChannelsDE / nChannelsDE;
	301
	302
	303	Bool_t goodDE = kTRUE;
	304	if (ratioBadChannelsDE >= maxBadChannels)
	305	{
	306	Double_t newRejectedEventsDE = rejectedEvents.DetectionElementProbability(DEid) + static_cast<Double_t>(nEvents);
	307	rejectedEvents.SetDetectionElementProbability(DEid, newRejectedEventsDE);
	308	goodDE = kFALSE;
	309	}
	310
	311	if (goodDE)
	312	for (Int_t ii = 0; ii < nBusPatches; ii++)
	313	{
	314	Int_t BPid = DEiter.CurrentDE()->GetBusPatchId(ii);
	315	Int_t nManus = AliMpDDLStore::Instance(kFALSE)->GetBusPatch(BPid, kFALSE)->GetNofManus();
	316
	317	Double_t nChannelsBP = static_cast<Double_t>(nManus * AliMpConstants::ManuNofChannels());
	318	Double_t nBadChannelsBP = nBadChannels->BusPatchProbability(BPid);
	319	Double_t ratioBadChannelsBP = 100.0 * nBadChannelsBP / nChannelsBP;
	320
	321	Bool_t goodBP = kTRUE;
	322
	323	if (goodBP)
	324	if (ratioBadChannelsBP >= maxBadChannels)
	325	{
	326	Double_t newRejectedEventsBP = rejectedEvents.BusPatchProbability(BPid) + static_cast<Double_t>(nEvents);
	327	rejectedEvents.SetBusPatchProbability(BPid, newRejectedEventsBP);
	328	goodBP = kFALSE;
	329	}
	330
	331	if (goodBP)
	332	for (Int_t jj = 0; jj < nManus; jj++)
	333	{
	334	Int_t manuId = AliMpDDLStore::Instance(kFALSE)->GetBusPatch(BPid, kFALSE)->GetManuId(jj);
	335
	336	Double_t nChannelsManu = DEiter.CurrentDE()->NofChannelsInManu(manuId);;
	337	Double_t nBadChannelsManu = nBadChannels->ManuProbability(DEid, manuId);
	338	Double_t ratioBadChannelsManu = 100.0 * nBadChannelsManu / nChannelsManu;
	339
	340	Bool_t goodManu = kTRUE;
	341
	342	if (ratioBadChannelsManu >= maxBadChannels)
	343	{
	344	Double_t newRejectedEventsManu = rejectedEvents.ManuProbability(DEid, manuId) + static_cast<Double_t>(nEvents);
	345	rejectedEvents.SetManuProbability(DEid, manuId, newRejectedEventsManu);
346	goodManu = kFALSE;
347	}
348
349
350	Int_t nChannels = DEiter.CurrentDE()->NofChannelsInManu(manuId);
351
352	if (goodManu)
353	for (Int_t channel = 0; channel < nChannels; channel++)
354	if (nBadChannels->ChannelProbability(DEid, manuId, channel) > 0.0)
355	{
356	Double_t newRejectedEventsChannel = rejectedEvents.ChannelProbability(DEid, manuId, channel) + static_cast<Double_t>(nEvents);
357	rejectedEvents.SetChannelProbability(DEid, manuId, channel, newRejectedEventsChannel);
358	} // end of Channels loop
359	} // end of MANUs loop
360	} // end of BPs loop
361	} // end of DEs loop
362
363	delete nBadChannels;
364
365	return nEvents;
366	}
367
368	//______________________________________________________________________________
369	//
370	//
371	// The code below might be used someday to retrieve the number of events
372	// with at least one muon directly from a runlist, instead of requiring
373	// as input both the run numbers and the number of events.
374	//
375	// This approach though is not without its problems :
376	//
377	// - even if we are using tags, as the tags are chunk-based (for the moment at
378	// least), it's pretty slow to loop on all the chunks (it is becoming usual
379	// to get runs with some 600 chunks or so...)
380	// - if part of a run is not reconstructed (for whatever reason) with the pass
381	// we use to compute the number of events with at least one muon, we might
382	// underestimate its weight in the final rejectlist for the next passes (
383	// if for the next passes the full run is correctly reconstructed for instance)
384	// - a muon in the tag = either a trigger track or a tracker track, i.e. we don't
385	// know for sure it's a good track...
386	//
387	// Given all those reasons, we for the moment use a dumb approach :
388	//
389	// - go to alice-logbook.cern.ch . Select the runs you want and the trigger class
390	// you want (most probably CMUS1B), and export the list as (run,nevents)
391	// - run this macro with this (run,nevents) list
392	//
393	//______________________________________________________________________________
394
395	#if 0
396
397	#include "TChain.h"
398	#include "AliDetectorTagCuts.h"
399	#include "AliEventTagCuts.h"
400	#include "AliLHCTagCuts.h"
401	#include "AliRunTagCuts.h"
402	#include "AliTagAnalysis.h"
403	#include "TGridResult.h"
404
405	//______________________________________________________________________________
406	bool CreateXMLCollectionFromRunList(const char* collectionName,
407	const char* runlist,
408	const char* type,
409	int passNumber)
410	{
411	/// From a list of run, create a collection with only events containing at least one muon
412	/// \param collectionName : output name of the collection (without extension, which will be .xml)
413	/// \param runlist : text file containing one integer per line = one run number per line
414	/// \param type : files to consider, either ESD or AD
415	/// \param passNumber : 1 or 2 most probably (to distinguish which reco pass is used)
416
417	if (!gGrid) TGrid::Connect("alien://");
418	if (!gGrid)
419	{
420	return 0x0;
421	}
422
f0938a5b	423	TString stype(type);
	424	stype.ToUpper();
	425
	426	if ( stype != "ESD" && stype != "AOD" )
	427	{
	428	cout << "Only ESD or AOD type supported" << endl;
	429	return false;
	430	}
	431
	432	ifstream in(gSystem->ExpandPathName(runlist));
	433	Int_t runNumber;
	434	Int_t ntagfiles(0);
	435
	436	AliTagAnalysis tagAnalysis("ESD");
	437
	438	while ( in >> runNumber )
	439	{
	440	TGridResult *res = gGrid->Query("/alice/data",
	441	Form("%09d/%ss/pass%d/%d/Run%d.Event*.ESD.tag.root",
	442	runNumber,stype.Data(),passNumber,runNumber,runNumber));
	443	Int_t nFiles = res->GetEntries();
	444	if (!nFiles)
	445	{
	446	continue;
	447	}
	448
	449	for (Int_t i = 0; i < nFiles; ++i)
	450	{
	451	TString filename = res->GetKey(i, "turl");
	452	if(filename == "") continue;
	453	tagAnalysis.AddTagsFile(filename.Data(),kFALSE);
	454	++ntagfiles;
	455	}
	456	delete res;
	457	}
	458
	459	cout << ntagfiles << " tag files added" << endl;
	460
	461	AliRunTagCuts runCuts;
	462	AliEventTagCuts eventCuts;
	463	AliLHCTagCuts lhcCuts;
	464	AliDetectorTagCuts detCuts;
	465
	466	eventCuts.SetNMuonRange(1,99999);
	467
	468	return tagAnalysis.CreateXMLCollection(collectionName,&runCuts,&lhcCuts,&detCuts,&eventCuts);
	469
	470	}
	471
	472	//______________________________________________________________________________
	473	TChain* CreateChainFromXMLCollection(const char* collectionName, const char* type)
	474	{
	475	/// Create a chain from an XML collection file.
	476
	477	if (!gGrid) TGrid::Connect("alien://");
	478	if (!gGrid)
	479	{
	480	return 0x0;
	481	}
	482
	483	TString stype(type);
	484	stype.ToUpper();
	485
	486	if ( stype != "ESD" && stype != "AOD" )
487	{
488	cout << "Only ESD or AOD type supported" << endl;
489	return 0x0;
490	}
491
492	AliTagAnalysis tagAnalysis(stype.Data());
493
494	return tagAnalysis.CreateChainFromCollection(collectionName,stype=="ESD" ? "esdTree":"aodTree");
495	}
496
497	#endif