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