macros to extract SSD gains from QA train output and to update the OCDB (M. Chojnacki)
authormasera <masera@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sun, 27 Feb 2011 09:38:26 +0000 (09:38 +0000)
committermasera <masera@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sun, 27 Feb 2011 09:38:26 +0000 (09:38 +0000)
ITS/GetGainModuleLevel.C [new file with mode: 0644]
ITS/UpdateSSDGainOCDB.C [new file with mode: 0644]

diff --git a/ITS/GetGainModuleLevel.C b/ITS/GetGainModuleLevel.C
new file mode 100644 (file)
index 0000000..8b341f4
--- /dev/null
@@ -0,0 +1,319 @@
+#include "TMath.h"\r
+#include "TList.h"\r
+#include "TH2F.h"\r
+#include "TH1F.h"\r
+#include "TF1.h"\r
+#include "TH1D.h"\r
+#include "TProfile.h"\r
+#include "TFile.h"\r
+#include "TCanvas.h"\r
+#include <Riostream.h>\r
+#include <TROOT.h>\r
+#include <TString.h>\r
+#include "TStyle.h"\r
+#include "TGrid.h"\r
\r
\r
+class TList;\r
+class TH2F;\r
+class TH1F;\r
+class TH1D;\r
+class TProfile;\r
+class TFile;\r
+class TCanvas;\r
+class TStyle;\r
+class TF1;\r
+class TGrid;\r
+\r
+\r
+\r
+\r
+\r
+Double_t convolution(Double_t* x,Double_t *par)\r
+{\r
+         //Fit parameters:\r
+   //par[0]=Width (scale) parameter of Landau density\r
+   //par[1]=Most Probable (MP, location) parameter of Landau density\r
+   //par[2]=Total area (integral -inf to inf, normalization constant)\r
+   //par[3]=Width (sigma) of convoluted Gaussian function\r
+   //\r
+   //In the Landau distribution (represented by the CERNLIB approximation),\r
+   //the maximum is located at x=-0.22278298 with the location parameter=0.\r
+   //This shift is corrected within this function, so that the actual\r
+   //maximum is identical to the MP parameter.\r
+\r
+      // Numeric constants\r
+      Double_t invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)\r
+      Double_t mpshift  = -0.22278298;       // Landau maximum location\r
+\r
+      // Control constants\r
+      Double_t np = 200.0;      // number of convolution steps\r
+      Double_t sc =   5.0;      // convolution extends to +-sc Gaussian sigmas\r
+\r
+      // Variables\r
+      Double_t xx;\r
+      Double_t mpc;\r
+      Double_t fland;\r
+      Double_t sum = 0.0;\r
+      Double_t xlow,xupp;\r
+      Double_t step;\r
+      Double_t i;\r
+\r
+\r
+      // MP shift correction\r
+      mpc = par[1] - mpshift * par[0];\r
+\r
+      // Range of convolution integral\r
+      xlow = x[0] - sc * par[3];\r
+      xupp = x[0] + sc * par[3];\r
+\r
+      step = (xupp-xlow) / np;\r
+\r
+      // Convolution integral of Landau and Gaussian by sum\r
+      for(i=1.0; i<=np/2; i++) \r
+      {\r
+                xx = xlow + (i-.5) * step;\r
+               fland = TMath::Landau(xx,mpc,par[0]) / par[0];\r
+               sum += fland * TMath::Gaus(x[0],xx,par[3]);\r
+\r
+               xx = xupp - (i-.5) * step;\r
+               fland = TMath::Landau(xx,mpc,par[0]) / par[0];\r
+               sum += fland * TMath::Gaus(x[0],xx,par[3]);\r
+      }\r
+\r
+      return (par[2] * step * sum * invsq2pi / par[3]);\r
+}\r
+\r
+\r
+\r
+\r
+void GetGainModuleLevel(TString filename,Bool_t normal=1,Int_t ntofit=500, Bool_t grid=0)\r
+{\r
+       gROOT->SetStyle("Plain");\r
+       gStyle->SetPalette(1,0);\r
+       gStyle->SetOptStat(111111);\r
+       gStyle->SetOptFit(1);   \r
+       if(grid)\r
+               TGrid::Connect("alien://");\r
+       \r
+       TFile* file_data=TFile::Open(filename);\r
+       if(!file_data)\r
+               return;\r
+       TList *listin=0x0;\r
+       listin=(TList*)file_data->Get("output");\r
+       if(!listin)\r
+               listin=(TList*)file_data->Get("PWG1dEdxSSDQA/output");\r
+       if(!listin)     \r
+               listin=(TList*)file_data->Get("PWG1dEdxSSDQA/SSDdEdxQA");\r
+       if(!listin)     \r
+               listin=(TList*)file_data->Get("SSDdEdxQA");\r
+       if(!listin)     \r
+               return;\r
+       TH2F* fHistQ=0x0;\r
+       if (normal)             \r
+               fHistQ=(TH2F*)listin ->FindObject("QACharge");\r
+       else\r
+               fHistQ=(TH2F*)listin ->FindObject("QAChargeCorrected");\r
+       if(!fHistQ)\r
+               return;\r
+       TH2F* fHistCR=(TH2F*)listin ->FindObject("QAChargeRatio");\r
+       if(!fHistCR)\r
+               return;\r
+       TList *listout1=new TList();\r
+       \r
+       TList *listout2=new TList();\r
+       \r
+       TH1F* fHistMPVs=new TH1F("HistMPVS","HistMPVs;MPV;N",75,70,95);\r
+       fHistMPVs->SetDirectory(0);\r
+       listout2->Add(fHistMPVs);\r
+       \r
+       TH1F* fHistSL=new TH1F("HistSL","#sigma_{Landau};#sigma_{Landau};N",40,0,16);\r
+       fHistSL->SetDirectory(0);\r
+       listout2->Add(fHistSL);\r
+       \r
+       TH1F* fHistSG=new TH1F("HistSG","#sigma_{Gaus};#sigma_{Gaus};N",40,0,16);\r
+       fHistSG->SetDirectory(0);\r
+       listout2->Add(fHistSG);\r
+       \r
+       TH1F* fHistCRmean=new TH1F("HistCRmean","HistCRmean;CRmean;N",200,-1,1);\r
+       fHistCRmean->SetDirectory(0);\r
+       listout2->Add(fHistCRmean);\r
+       \r
+       TH1F* fHistCRRMS=new TH1F("HistCRRMS","HistCRRMS;CRRMS;N",100,0,1);\r
+       fHistCRRMS->SetDirectory(0);\r
+       listout2->Add(fHistCRRMS);\r
+       \r
+       TH1F* fHistGainP=new TH1F("HistGainP","HistGainP;CRGainPcorr;N",120,0.5,2.0);\r
+       fHistGainP->SetDirectory(0);\r
+       listout2->Add(fHistGainP);\r
+       \r
+       TH1F* fHistGainN=new TH1F("HistGainN","HistGainN;CRGainNcorr;N",120,0.5,2.0);\r
+       fHistGainN->SetDirectory(0);\r
+       listout2->Add(fHistGainN);\r
+       \r
+       TH1F *fMPVGraph = new TH1F("MPVgraph","MPVgraph;Module number;MPV",1698,-0.5,1697.5);\r
+       fMPVGraph->SetMarkerColor(kRed);\r
+       fMPVGraph->SetMarkerStyle(22);\r
+       listout2->Add(fMPVGraph);\r
+       \r
+       TH1F *fCRmeanGraph = new TH1F("CRmeangraph","CRmeangraph;Module number;MPV",1698,-0.5,1697.5);\r
+       fCRmeanGraph->SetMarkerColor(kBlue);\r
+       fCRmeanGraph->SetMarkerStyle(23);\r
+       listout2->Add(fCRmeanGraph);\r
+       \r
+       Float_t gainP[1698];\r
+       Float_t gainN[1698];\r
+       Float_t mpv[1698];\r
+       Int_t flag[1698];       \r
+       \r
+       ofstream outfiletxt;\r
+       outfiletxt.open("gain.txt");\r
+        outfiletxt.width(10) ;\r
+       outfiletxt.setf(outfiletxt.left);\r
+       outfiletxt<<"MODULE"<<"\t";\r
+       outfiletxt.width(10);\r
+       outfiletxt.setf(outfiletxt.left);\r
+       outfiletxt<<"FLAG"<<"\t";\r
+       outfiletxt.width(10) ;\r
+       outfiletxt.setf(outfiletxt.left);\r
+       outfiletxt<<"GainPcorr"<<"\t";\r
+       outfiletxt.width(10) ;\r
+       outfiletxt.setf(outfiletxt.left);\r
+       outfiletxt<<"GainNcorr"<<"\t";\r
+       outfiletxt.width(10) ;\r
+       outfiletxt.setf(outfiletxt.left);\r
+       outfiletxt<<"MPV"<<endl;\r
+       \r
+       \r
+       ofstream outfiletxtbad;\r
+       outfiletxtbad.open("badModules.txt");\r
+       \r
+       \r
+       \r
+       \r
+       for (int i =0;i<1698;i++)\r
+       {\r
+               cout<<i<<endl;\r
+               TString tmpQ("Q");\r
+               tmpQ+=i;\r
+               TString tmpCR("CR");\r
+               tmpCR+=i;\r
+               TH1D* fHist1DCR= fHistCR->ProjectionY(tmpCR,i+1,i+1);\r
+               Double_t mean=fHist1DCR->GetMean();\r
+               if(!(TMath::Abs(mean)<1.0)||fHist1DCR->GetEntries()<10)\r
+               {               \r
+                       flag[i]=-2;\r
+                       gainN[i]=1.0;\r
+                       gainP[i]=1.0;\r
+                       mpv[i]=1.0;\r
+                       continue;\r
+               }\r
+               fHistCRmean->Fill(mean);\r
+               fHistCRRMS->Fill(fHist1DCR->GetRMS());\r
+               gainN[i]=1.0/(1.0+mean);\r
+               gainP[i]=1.0/(1.0-mean);\r
+               fHistGainP->Fill(gainP[i]);\r
+               fHistGainN->Fill(gainN[i]);\r
+               fCRmeanGraph->SetBinContent(i+1,mean);\r
+               fCRmeanGraph->SetBinError(i+1,fHist1DCR->GetRMS());\r
+               \r
+               \r
+               TH1D* fHist1DQ=fHistQ->ProjectionY(tmpQ,i+1,i+1);\r
+                fHist1DQ->SetDirectory(0);\r
+                listout1->Add(fHist1DQ);\r
+               if(fHist1DQ->GetEntries()<ntofit)\r
+               {\r
+                       flag[i]=-1;\r
+                       mpv[i]=1.0;\r
+                        outfiletxtbad<<"Low statistic \t module= "<<i<<" netries="<<fHist1DQ->GetEntries()<<endl;\r
+                       continue;\r
+               }\r
+               else\r
+               {\r
+                       tmpQ+="fit";\r
+                       Float_t range=fHist1DQ->GetBinCenter(fHist1DQ->GetMaximumBin());\r
+                       TF1 *f1 = new TF1(tmpQ,convolution,range*0.45,range*3.0,4);\r
+                       f1->SetParameters(7.0,range,1.0,5.5);\r
+                       Float_t normalization=fHist1DQ->GetEntries()*fHist1DQ->GetXaxis()->GetBinWidth(2)/f1->Integral(range*0.45,range*3.0);\r
+                       f1->SetParameters(7.0,range,normalization,5.5);\r
+                       //f1->SetParameters(7.0,range,fHist1DQ->GetMaximum(),5.5);\r
+                       f1->SetParNames("sigma Landau","MPV","N","sigma Gaus");\r
+                       f1->SetParLimits(0,2.0,100.0);\r
+                       f1->SetParLimits(3,0.0,100.0);\r
+                       if(fHist1DQ->Fit(tmpQ,"BRQ")==0)\r
+                       {\r
+                               mpv[i]=f1->GetParameter(1);\r
+                               fHistMPVs->Fill(mpv[i]);\r
+                               fHistSL->Fill(f1->GetParameter(0));\r
+                               fHistSG->Fill(f1->GetParameter(3));\r
+                               flag[i]=1;              \r
+                               fMPVGraph->SetBinContent(i+1,f1->GetParameter(1));\r
+                               fMPVGraph->SetBinError(i+1,f1->GetParError(1));\r
+                               if(mpv[i]<75.0)\r
+                               {\r
+                                       outfiletxtbad<<"MPV lower than 75 \t module="<<i<<endl;\r
+                                       flag[i]=0;\r
+                               }       \r
+                               if(mpv[i]>100.0)\r
+                               {\r
+                                       outfiletxtbad<<"MPV higher than 100 \t module="<<i<<endl;\r
+                                       flag[i]=0;\r
+                                       \r
+                               }\r
+                               if(f1->GetParError(1)>1.0)\r
+                               {\r
+                                       outfiletxtbad<<"MPV high error on MPV  \t module="<<i<<endl;                            \r
+                                       flag[i]=0;\r
+                               }\r
+                       }\r
+                       else\r
+                       {\r
+                               mpv[i]=1;\r
+                               flag[i]=0;\r
+                                outfiletxtbad<<"BAD FIT \t module="<<i<<endl;\r
+                               //continue;\r
+                       }       \r
+               }       \r
+       }       \r
+       \r
+       for (int i=0;i<1698;i++)\r
+       {       \r
+               outfiletxt.setf(outfiletxt.scientific);\r
+               outfiletxt.precision(2);        \r
+                outfiletxt.width(10) ;\r
+               outfiletxt.setf(outfiletxt.left);\r
+               outfiletxt<<i<<"\t";\r
+                outfiletxt.width(10) ;\r
+               outfiletxt.setf(outfiletxt.left);\r
+               outfiletxt<<flag[i]<<"\t";\r
+                outfiletxt.width(10) ;\r
+               outfiletxt.setf(outfiletxt.left);\r
+               outfiletxt<<gainP[i]<<"\t";\r
+                outfiletxt.width(10) ;\r
+               outfiletxt.setf(outfiletxt.left);\r
+               outfiletxt<<gainN[i]<<"\t";\r
+                outfiletxt.width(10) ;\r
+               outfiletxt.setf(outfiletxt.left);\r
+               outfiletxt<<mpv[i]<<endl;\r
+       }\r
+                \r
+       TCanvas *c1 = new TCanvas("1","1",1200,800);\r
+       c1->Divide(2,1);\r
+       c1->cd(1);\r
+       fHistQ->Draw("colz");\r
+       c1->cd(2);\r
+       fHistCR->Draw("colz");\r
+       \r
+       \r
+       \r
+       TFile* fout1=TFile::Open("gain_all_fits.root","recreate");\r
+       listout1->Write("output",TObject::kSingleKey);  \r
+       fout1->Close();\r
+       \r
+       TFile* fout2=TFile::Open("gain_control_plots.root","recreate");\r
+       listout2->Write("output",TObject::kSingleKey);  \r
+       fout2->Close();\r
+       \r
+       \r
+       \r
+}\r
diff --git a/ITS/UpdateSSDGainOCDB.C b/ITS/UpdateSSDGainOCDB.C
new file mode 100644 (file)
index 0000000..734bdba
--- /dev/null
@@ -0,0 +1,143 @@
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include <Riostream.h>
+#include <TFile.h>
+#include <TString.H>
+#include "AliITSGainSSDv2.h"
+#include "AliCDBEntry.h"
+#include "AliCDBManager.h"
+#include "AliCDBMetaData.h"
+#endif
+Int_t GainCalibration(AliCDBManager * man, Float_t* gainP, Float_t* gainN) 
+{
+       const Int_t fgkSSDMODULES = 1698;
+       static const Int_t fgkDefaultNStripsSSD = 768;
+       AliITSGainSSDv2 *gainSSD = new AliITSGainSSDv2();
+       AliCDBEntry *entryGainSSD = man->Get("ITS/Calib/GainSSD");
+       TObject *empty = (TObject *)entryGainSSD->GetObject();
+       TString objectname = empty->GetName();
+       if(objectname=="AliITSGainSSDv2") 
+       {
+               cout<<"Reading the new format of the calibration file..."<<endl;
+               gainSSD = (AliITSGainSSDv2 *)entryGainSSD->GetObject();
+       }
+       else
+       {
+               cout<<"Now OCDB object"<<endl;
+               return 0 ;
+       }
+       for (Int_t i = 0; i < fgkSSDMODULES; i++) 
+       {
+                       //for (Int_t i = 0; i < 1; i++) {
+                       //AliITSgeomTGeo::GetModuleId(i+500,layer,ladder,module);
+               for(Int_t j = 0; j < fgkDefaultNStripsSSD; j++) 
+               {
+                       gainP[i*fgkDefaultNStripsSSD+j] = gainSSD->GetGainP(i,j);
+                       gainN[i*fgkDefaultNStripsSSD+j]= gainSSD->GetGainN(i,j);
+       
+               }//strip loop
+       }//module loop
+       return 1;
+}
+
+
+void UpdateSSDGainOCDB(TString gainfile,Float_t globalMPV=84.0,const char* type = "alien", Int_t runNumber = 0, Int_t runRangemin=0 ,Int_t runRangemax= AliCDBRunRange::Infinity()) 
+{
+       if((runRangemax<runRangemin))
+       {
+               cout<<"wrong run range "<<endl;
+               return;
+       }       
+       const Int_t fgkPNStripsPerModule = 768;
+       const Int_t fgkNumberOfSSDModules = 1698;
+       Float_t gainPold[fgkPNStripsPerModule* fgkNumberOfSSDModules];
+       Float_t gainNold[fgkPNStripsPerModule* fgkNumberOfSSDModules];
+       
+       for(int i=0;i<fgkPNStripsPerModule* fgkNumberOfSSDModules;i++)
+       {
+                       gainPold[i]=0.0;
+                       gainNold[i]=0.0;
+       }
+       
+       TString gType = type;
+       AliCDBManager * man = AliCDBManager::Instance();
+       if(gType == "alien") 
+                man->SetDefaultStorage("alien://folder=/alice/data/2010/OCDB/");
+       else if(gType == "local") 
+               man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
+       else 
+       {
+               cout<<"Allowed types: local or alien!!!"<<endl;
+               return;
+       }
+         man->SetRun(runNumber);
+       if(!GainCalibration(man,gainPold,gainNold))
+               return;
+               
+       char buffer[256];
+       ifstream ingainfile(gainfile);
+       if(!ingainfile.good())
+               return;
+       
+       
+        ingainfile.getline(buffer,256);
+       cout<<buffer<<endl;
+       
+       
+       Int_t module=0;
+       Int_t gainflag=0;
+       Float_t gainP=0;
+       Float_t gainN=0;
+       Float_t fMPV=0;
+       Float_t corrP[1698];
+       Float_t corrN[1698];
+       
+       for(int jj=0;jj<1698;jj++)
+       {
+               ingainfile>>module;
+               ingainfile>>gainflag;
+               ingainfile>>gainP;
+               ingainfile>>gainN;
+               ingainfile>>fMPV;
+               if(gainflag==1)
+               {
+                       corrP[module]=gainP*globalMPV/fMPV;
+                       corrN[module]=gainN*globalMPV/fMPV;
+                       
+               }
+               else
+               {
+                       corrP[module]=1.0;
+                       corrN[module]=1.0;
+               }
+                       
+       }
+       
+       
+       TObjArray *array = new TObjArray();
+       array->SetName("Gain");
+       AliITSGainSSDv2  *mc;
+       mc = new AliITSGainSSDv2();
+       for(Int_t mind = 0; mind < fgkNumberOfSSDModules; mind++) 
+       {
+               for (Int_t i = 0; i < fgkPNStripsPerModule; i++) 
+               {
+                       mc->AddGainP(mind,i, gainPold[mind* fgkPNStripsPerModule+i]*corrP[mind]);
+                       mc->AddGainN(mind,i,gainNold[mind* fgkPNStripsPerModule+i]*corrN[mind]);                        
+               }
+       }
+
+       man = AliCDBManager::Instance();
+       man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
+       man->SetRun(0);
+       //AliCDBId id("ITS/Calib/GainSSD",0,AliCDBRunRange::Infinity());
+       AliCDBId id("ITS/Calib/GainSSD",runRangemin,runRangemax);
+       AliCDBMetaData *metadata= new AliCDBMetaData();
+       metadata->SetResponsible("Marek.Chojnacki@cern.ch");
+       metadata->SetComment("Default values for the SSD gain calibration");
+       //man->Put(array,id,metadata);
+       man->Put(mc,id,metadata);
+        TFile* fout=TFile::Open("AliITSGainSSD_new.root","recreate");
+       mc ->Write();   
+               fout->Close();
+       
+}