MUON/CreateWeightedRejectList.C

   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 ///
  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 ///
  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
 280   status.SetLimits(*recoParam);
 281
 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
 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