Macros to produce bad modules map from ESD

diff --git a/PHOS/macros/BadModules/Readme b/PHOS/macros/BadModules/Readme
new file mode 100644 (file)
index 0000000..116367c
--- /dev/null
+++ b/PHOS/macros/BadModules/Readme
@@ -0,0 +1,9 @@
+This dir contains compilable macros for selecting bad modules
+using ESD. This is done in 2 steps:
+1. Make scan of ESD data using macro scanESD.C. It produces 
+   set of spectra of (all) clusters with centers in different channels.
+2. Analyze produced spectra with macro makeBadModESD.C. It produces
+   bad map in the form of TH2S and CDB entry (locally).
+
+
+ 

diff --git a/PHOS/macros/BadModules/makeBadModESD.C b/PHOS/macros/BadModules/makeBadModESD.C
new file mode 100644 (file)
index 0000000..601f85b
--- /dev/null
+++ b/PHOS/macros/BadModules/makeBadModESD.C
@@ -0,0 +1,169 @@
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include "TH3.h"
+#include "TF1.h"
+#include "TFile.h"
+#include "TLine.h"
+#include "TCanvas.h"
+#include "AliPHOSEmcBadChannelsMap.h"
+#include "AliCDBId.h"
+#include "AliCDBMetaData.h"
+#include "AliCDBEntry.h"
+#include "AliCDBManager.h"
+#endif
+void makeBadModESD(Int_t mod=3){
+  //Here we use 3D distributions of pedestals and spectra 
+  //produced by the macro scanESD.C 
+  //Having these distributions we calculate distributions
+  //over number of clusters in "Soft" and "Hard" regions and 
+  //exlude cells, havig considerable deviations from avrage
+  //Result is stored in TH2S in root file and as CDB root file
+  //
+  //Currect version is valid for one PHOS module (e.g. cosmic run)
+  //Author D.Peressounko 
+  
+  TFile * f = new TFile("scan.root") ;
+  TH3D * hSp = (TH3D*)f->Get("hEsd") ;
+
+  //Boundary of "Soft" and "Hard" regions
+  //Soft region used to exclude dead or almost dead cells
+  //Hard region used to exclude "singers"
+  //All energies in GeV
+  Double_t aSoftMin=0.22 ;  
+  Double_t aSoftMax=0.5 ;
+  Double_t aHardMin=1. ;
+  Double_t aHardMax=2. ;
+
+  //Final bad map
+  TH2D * hBadMap = new TH2D("hBadMap","Bad Modules map",64,0.,64.,56,0.,56.) ;
+
+  TH2D * hSoft = new TH2D("hSoft","Number of soft clusters",64,0.,64.,56,0.,56.) ;
+  TH2D * hHard = new TH2D("hHard","Number of hard clusters",64,0.,64.,56,0.,56.) ;
+
+  //The range of these histos depends on statistics, one probaly will have to increase it  
+  TH1D * hAllSoft = new TH1D("hallSoft","Summary of Soft",200,0.,200.) ;
+  TH1D * hAllHard = new TH1D("hallHard","Summary of Hard",100,0.,100.) ;
+
+  TAxis * ax = hSp->GetZaxis() ;
+  Int_t iSoftMin=ax->FindBin(aSoftMin) ;
+  Int_t iSoftMax=ax->FindBin(aSoftMax) ;
+  Int_t iHardMin=ax->FindBin(aHardMin) ;
+  Int_t iHardMax=ax->FindBin(aHardMax) ;
+
+  //Fill histograms with number of counts per cell
+  for(Int_t i=1; i<=hSp->GetXaxis()->GetNbins(); i++){
+    for(Int_t j=1; j<=hSp->GetYaxis()->GetNbins(); j++){
+      Double_t soft=0 ;
+      for(Int_t k=iSoftMin; k<=iSoftMax;k++){
+         soft+=hSp->GetBinContent(i,j,k) ;
+      }
+      Double_t hard=0. ;
+      for(Int_t k=iHardMin; k<=iHardMax;k++){
+         hard+=hSp->GetBinContent(i,j,k) ;
+      }
+      hSoft->SetBinContent(i,j,soft) ;
+      hHard->SetBinContent(i,j,hard) ;
+      hAllSoft->Fill(soft) ;
+      hAllHard->Fill(hard) ;
+    }
+  }
+
+
+  //=====================================================================================
+  //
+  //Usually the part below needs some manual intervention to properly select good modules region
+  //
+  //Fit overall distribution with Gauss and calculate limits of "goodness" 
+  //usually +-3-4 sigma for soft and ~4-5 sigma for hard (less statistics)
+  TF1 * gaus = new TF1("gaus","[0]*exp(-(x-[1])*(x-[1])/2./[2]/[2])",0.,1000.) ;
+  gaus->SetParameters(100.,10.,20.) ;
+  gaus->FixParameter(1,0.) ;
+  hAllSoft->Fit(gaus,"","",10,80) ;
+  Double_t m = gaus->GetParameter(1) ;
+  Double_t s = TMath::Abs(gaus->GetParameter(2)) ;
+  //one can use either sigma or hardwired parameterization
+  //of limits of goodness
+  Double_t cutSoftMin= 8. ; //TMath::Max(0.,m-4.*s) ;
+  Double_t cutSoftMax= 80. ; //m+4.*s ;
+  Double_t softH= gaus->GetParameter(0) ;
+
+  //repeat same procedure for hard
+  gaus->SetParameters(100.,2.,10.) ;
+  gaus->FixParameter(1,0.) ;
+  hAllHard->Fit(gaus,"","",1.,20.) ;
+  m = gaus->GetParameter(1) ;
+  s = TMath::Abs(gaus->GetParameter(2)) ;
+  Double_t cutHardMin= TMath::Max(0.,m-5.*s) ;
+  Double_t cutHardMax= m+5.*s ;
+  Double_t hardH= gaus->GetParameter(0) ;
+
+  //=====================================================================================
+
+  //Draw  results
+  TCanvas * cm = new TCanvas("Soft","Soft",15,22,500,300) ;
+  hSoft->Draw("colz") ;
+
+  TCanvas * cr = new TCanvas("Hard","Hard",520,22,500,300) ;
+  hHard->Draw("colz") ;
+
+  TCanvas * cma = new TCanvas("AllSoft","Soft Multiplicity",15,330,500,300) ;
+  cma->SetLogy() ;
+  hAllSoft->Draw() ;
+  TLine * l1 = new TLine(cutSoftMin,0.1,cutSoftMin,2.*softH) ;
+  l1->SetLineColor(2) ;
+  l1->Draw() ;
+  TLine * l2 = new TLine(cutSoftMax,0.1,cutSoftMax,2.*softH) ;
+  l2->SetLineColor(2) ;
+  l2->Draw() ;
+
+  TCanvas * cra = new TCanvas("AllHard","Hard Multiplicity",520,330,500,300) ;
+  cra->SetLogy() ;
+  hAllHard->Draw() ;
+  printf("Min=%f, max=%f \n",cutHardMin,cutHardMax) ;
+  TLine * l3 = new TLine(cutHardMin,0.1,cutHardMin,2.*hardH) ;
+  l3->SetLineColor(2) ;
+  l3->Draw() ;
+  TLine * l4 = new TLine(cutHardMax,0.1,cutHardMax,2.*hardH) ;
+  l4->SetLineColor(2) ;
+  l4->Draw() ;
+
+
+  //now calculate bad map
+  //Exclude everething beyond the limits
+  AliPHOSEmcBadChannelsMap badMap;
+  for(Int_t i=1; i<=hSp->GetXaxis()->GetNbins(); i++){
+    for(Int_t j=1; j<=hSp->GetYaxis()->GetNbins(); j++){
+      Double_t soft = hSoft->GetBinContent(i,j) ;
+      Double_t hard= hHard->GetBinContent(i,j) ;
+      if(soft>=cutSoftMin   && soft <=cutSoftMax &&
+        hard>=cutHardMin && hard <=cutHardMax){
+        hBadMap->SetBinContent(i,j,1.) ;
+        hBadMap->SetBinContent(i,j,1.) ;
+      }
+      else{
+        hBadMap->SetBinContent(i,j,0.) ;
+        badMap.SetBadChannel(mod,i,j); //module, col,row
+      }
+    }
+  }
+  
+  TCanvas * cb = new TCanvas("BadMap","Bad map",15,630,500,300) ;
+  hBadMap->DrawClone("colz") ;
+
+  TFile * fout = new TFile("BadMap.root","recreate") ;
+  hBadMap->Write() ;
+  fout->Close() ;
+
+  //put now result to local CDB
+  AliCDBManager *CDB = AliCDBManager::Instance();
+  //  CDB->SetDefaultStorage("local://$ALICE_ROOT");
+  CDB->SetDefaultStorage("local://./");
+  //  CDB->SetSpecificStorage("local://./","PHOS");
+  
+  AliCDBMetaData *md= new AliCDBMetaData();
+  md->SetResponsible("Dmitri Peressounko");
+  md->SetComment("Dead channels for run 1234");
+  AliCDBId id("PHOS/Calib/EmcBadChannels",0,999999);
+  CDB->Put(&badMap,id, md);
+ 
+
+}

diff --git a/PHOS/macros/BadModules/scanESD.C b/PHOS/macros/BadModules/scanESD.C
new file mode 100644 (file)
index 0000000..a8435b0
--- /dev/null
+++ b/PHOS/macros/BadModules/scanESD.C
@@ -0,0 +1,66 @@
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include "TH3.h"
+#include "TChain.h"
+#include "TFile.h"
+#include "TTree.h"
+#include "AliCDBManager.h"
+#include "AliGeomManager.h"
+#include "AliRunLoader.h"
+#include "AliESDEvent.h"
+#include "AliESDCaloCluster.h"
+#include "AliPHOSGeometry.h"
+#endif
+void scanESD(Int_t irun=7727)
+{
+  //This macro scans ESD and fills 3D histogram:
+  //spectrum of clusters with center in each cristall
+  //this histogram can be used e.g. in bad modules selection
+  //Author: D.Peressounko Dmitri.Peressounko@cern.ch
+
+  //Uncomment the following if misalignement should be applied
+  AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT");
+  //   AliCDBManager::Instance()->SetDefaultStorage("local://./");
+  //   AliCDBManager::Instance()->SetSpecificStorage("PHOS/*","local:///data/prsnko/");
+  AliCDBManager::Instance()->SetRun(irun) ;
+
+  //The final histo
+  TH3D * hEsd = new TH3D("hEsd","Energy of all clusters",64,0.,64.,56,0.,56.,2000,0.,10.) ;
+
+  //Reading/creation of Geometry
+  AliGeomManager::LoadGeometry("geometry.root");
+  AliPHOSGeometry *phosgeom = AliPHOSGeometry::GetInstance("IHEP","") ;
+
+  //Create list of ESD files
+  TChain * chain = new TChain("esdTree") ;
+  char filename[255] ;
+  for(Int_t sec=1; sec<=17; sec++){
+    sprintf(filename,"Seq_%d0/AliESDs.root/esdTree",sec) ;
+    chain->AddFile(filename);
+  }
+  
+  AliESDEvent *event = new AliESDEvent();
+  event->ReadFromTree(chain);
+  
+  for(Int_t iEvent=1; iEvent<chain->GetEntries(); iEvent++){
+     if(iEvent%10000==0){
+       printf("Event %d \n",iEvent) ;
+     }
+    chain->GetEvent(iEvent);
+    Int_t multClu = event->GetNumberOfCaloClusters();
+    for (Int_t i0=0; i0<multClu; i0++) {
+      AliESDCaloCluster * clu1 = event->GetCaloCluster(i0);
+      Float_t xyz[3] = {0,0,0};
+      clu1->GetPosition(xyz);   //Gloabal position in ALICE system
+      Int_t absId ;
+      phosgeom->RelPosToAbsId(3,xyz[0],xyz[2],absId) ; //Here we assume that data taken with 3 module
+      Int_t relid[4] ;
+      phosgeom->AbsToRelNumbering(absId,relid) ;
+      hEsd->Fill(relid[2]*1.-0.5, relid[3]*1.-0.5,clu1->E()) ;
+    }   
+  }
+  
+  TFile *ff = new TFile("scan.root","recreate") ;
+  hEsd->Write() ;
+  ff->Close() ;
+
+}