PHOS/macros/BadModules/makeBadModESD.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2 #include "TH3.h"
   3 #include "TF1.h"
   4 #include "TFile.h"
   5 #include "TLine.h"
   6 #include "TCanvas.h"
   7 #include "AliPHOSEmcBadChannelsMap.h"
   8 #include "AliCDBId.h"
   9 #include "AliCDBMetaData.h"
  10 #include "AliCDBEntry.h"
  11 #include "AliCDBManager.h"
  12 #endif
  13 void makeBadModESD(Int_t mod=3){
  14   //Here we use 3D distributions of pedestals and spectra
  15   //produced by the macro scanESD.C
  16   //Having these distributions we calculate distributions
  17   //over number of clusters in "Soft" and "Hard" regions and
  18   //exlude cells, havig considerable deviations from avrage
  19   //Result is stored in TH2S in root file and as CDB root file
  20   //
  21   //Currect version is valid for one PHOS module (e.g. cosmic run)
  22   //Author D.Peressounko
  23
  24   TFile * f = new TFile("scan.root") ;
  25   TH3D * hSp = (TH3D*)f->Get("hEsd") ;
  26
  27   //Boundary of "Soft" and "Hard" regions
  28   //Soft region used to exclude dead or almost dead cells
  29   //Hard region used to exclude "singers"
  30   //All energies in GeV
  31   Double_t aSoftMin=0.22 ;
  32   Double_t aSoftMax=0.5 ;
  33   Double_t aHardMin=1. ;
  34   Double_t aHardMax=2. ;
  35
  36   //Final bad map
  37   TH2D * hBadMap = new TH2D("hBadMap","Bad Modules map",64,0.,64.,56,0.,56.) ;
  38
  39   TH2D * hSoft = new TH2D("hSoft","Number of soft clusters",64,0.,64.,56,0.,56.) ;
  40   TH2D * hHard = new TH2D("hHard","Number of hard clusters",64,0.,64.,56,0.,56.) ;
  41
  42   //The range of these histos depends on statistics, one probaly will have to increase it
  43   TH1D * hAllSoft = new TH1D("hallSoft","Summary of Soft",200,0.,200.) ;
  44   TH1D * hAllHard = new TH1D("hallHard","Summary of Hard",100,0.,100.) ;
  45
  46   TAxis * ax = hSp->GetZaxis() ;
  47   Int_t iSoftMin=ax->FindBin(aSoftMin) ;
  48   Int_t iSoftMax=ax->FindBin(aSoftMax) ;
  49   Int_t iHardMin=ax->FindBin(aHardMin) ;
  50   Int_t iHardMax=ax->FindBin(aHardMax) ;
  51
  52   //Fill histograms with number of counts per cell
  53   for(Int_t i=1; i<=hSp->GetXaxis()->GetNbins(); i++){
  54     for(Int_t j=1; j<=hSp->GetYaxis()->GetNbins(); j++){
  55       Double_t soft=0 ;
  56       for(Int_t k=iSoftMin; k<=iSoftMax;k++){
  57          soft+=hSp->GetBinContent(i,j,k) ;
  58       }
  59       Double_t hard=0. ;
  60       for(Int_t k=iHardMin; k<=iHardMax;k++){
  61          hard+=hSp->GetBinContent(i,j,k) ;
  62       }
  63       hSoft->SetBinContent(i,j,soft) ;
  64       hHard->SetBinContent(i,j,hard) ;
  65       hAllSoft->Fill(soft) ;
  66       hAllHard->Fill(hard) ;
  67     }
  68   }
  69
  70
  71   //=====================================================================================
  72   //
  73   //Usually the part below needs some manual intervention to properly select good modules region
  74   //
  75   //Fit overall distribution with Gauss and calculate limits of "goodness"
  76   //usually +-3-4 sigma for soft and ~4-5 sigma for hard (less statistics)
  77   TF1 * gaus = new TF1("gaus","[0]*exp(-(x-[1])*(x-[1])/2./[2]/[2])",0.,1000.) ;
  78   gaus->SetParameters(100.,10.,20.) ;
  79   gaus->FixParameter(1,0.) ;
  80   hAllSoft->Fit(gaus,"","",10,80) ;
  81   Double_t m = gaus->GetParameter(1) ;
  82   Double_t s = TMath::Abs(gaus->GetParameter(2)) ;
  83   //one can use either sigma or hardwired parameterization
  84   //of limits of goodness
  85   Double_t cutSoftMin= 8. ; //TMath::Max(0.,m-4.*s) ;
  86   Double_t cutSoftMax= 80. ; //m+4.*s ;
  87   Double_t softH= gaus->GetParameter(0) ;
  88
  89   //repeat same procedure for hard
  90   gaus->SetParameters(100.,2.,10.) ;
  91   gaus->FixParameter(1,0.) ;
  92   hAllHard->Fit(gaus,"","",1.,20.) ;
  93   m = gaus->GetParameter(1) ;
  94   s = TMath::Abs(gaus->GetParameter(2)) ;
  95   Double_t cutHardMin= TMath::Max(0.,m-5.*s) ;
  96   Double_t cutHardMax= m+5.*s ;
  97   Double_t hardH= gaus->GetParameter(0) ;
  98
  99   //=====================================================================================
 100
 101   //Draw  results
 102   TCanvas * cm = new TCanvas("Soft","Soft",15,22,500,300) ;
 103   hSoft->Draw("colz") ;
 104
 105   TCanvas * cr = new TCanvas("Hard","Hard",520,22,500,300) ;
 106   hHard->Draw("colz") ;
 107
 108   TCanvas * cma = new TCanvas("AllSoft","Soft Multiplicity",15,330,500,300) ;
 109   cma->SetLogy() ;
 110   hAllSoft->Draw() ;
 111   TLine * l1 = new TLine(cutSoftMin,0.1,cutSoftMin,2.*softH) ;
 112   l1->SetLineColor(2) ;
 113   l1->Draw() ;
 114   TLine * l2 = new TLine(cutSoftMax,0.1,cutSoftMax,2.*softH) ;
 115   l2->SetLineColor(2) ;
 116   l2->Draw() ;
 117
 118   TCanvas * cra = new TCanvas("AllHard","Hard Multiplicity",520,330,500,300) ;
 119   cra->SetLogy() ;
 120   hAllHard->Draw() ;
 121   printf("Min=%f, max=%f \n",cutHardMin,cutHardMax) ;
 122   TLine * l3 = new TLine(cutHardMin,0.1,cutHardMin,2.*hardH) ;
 123   l3->SetLineColor(2) ;
 124   l3->Draw() ;
 125   TLine * l4 = new TLine(cutHardMax,0.1,cutHardMax,2.*hardH) ;
 126   l4->SetLineColor(2) ;
 127   l4->Draw() ;
 128
 129
 130   //now calculate bad map
 131   //Exclude everething beyond the limits
 132   AliPHOSEmcBadChannelsMap badMap;
 133   for(Int_t i=1; i<=hSp->GetXaxis()->GetNbins(); i++){
 134     for(Int_t j=1; j<=hSp->GetYaxis()->GetNbins(); j++){
 135       Double_t soft = hSoft->GetBinContent(i,j) ;
 136       Double_t hard= hHard->GetBinContent(i,j) ;
 137       if(soft>=cutSoftMin   && soft <=cutSoftMax &&
 138         hard>=cutHardMin && hard <=cutHardMax){
 139         hBadMap->SetBinContent(i,j,1.) ;
 140         hBadMap->SetBinContent(i,j,1.) ;
 141       }
 142       else{
 143         hBadMap->SetBinContent(i,j,0.) ;
 144         badMap.SetBadChannel(mod,i,j); //module, col,row
 145       }
 146     }
 147   }
 148
 149   TCanvas * cb = new TCanvas("BadMap","Bad map",15,630,500,300) ;
 150   hBadMap->DrawClone("colz") ;
 151
 152   TFile * fout = new TFile("BadMap.root","recreate") ;
 153   hBadMap->Write() ;
 154   fout->Close() ;
 155
 156   //put now result to local CDB
 157   AliCDBManager *CDB = AliCDBManager::Instance();
 158   //  CDB->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
 159   CDB->SetDefaultStorage("local://./");
 160   //  CDB->SetSpecificStorage("local://./","PHOS");
 161
 162   AliCDBMetaData *md= new AliCDBMetaData();
 163   md->SetResponsible("Dmitri Peressounko");
 164   md->SetComment("Dead channels for run 1234");
 165   AliCDBId id("PHOS/Calib/EmcBadChannels",0,999999);
 166   CDB->Put(&badMap,id, md);
 167
 168
 169 }