end-of-line normalization
[u/mrichter/AliRoot.git] / ITS / GetGainModuleLevel.C
index 1421fcd..d9e0fa7 100644 (file)
-#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("PWGPPdEdxSSDQA/output");\r
-       if(!listin)     \r
-               listin=(TList*)file_data->Get("PWGPPdEdxSSDQA/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
+#include "TMath.h"
+#include "TList.h"
+#include "TH2F.h"
+#include "TH1F.h"
+#include "TF1.h"
+#include "TH1D.h"
+#include "TProfile.h"
+#include "TFile.h"
+#include "TCanvas.h"
+#include <Riostream.h>
+#include <TROOT.h>
+#include <TString.h>
+#include "TStyle.h"
+#include "TGrid.h"
+class TList;
+class TH2F;
+class TH1F;
+class TH1D;
+class TProfile;
+class TFile;
+class TCanvas;
+class TStyle;
+class TF1;
+class TGrid;
+
+
+
+
+
+Double_t convolution(Double_t* x,Double_t *par)
+{
+         //Fit parameters:
+   //par[0]=Width (scale) parameter of Landau density
+   //par[1]=Most Probable (MP, location) parameter of Landau density
+   //par[2]=Total area (integral -inf to inf, normalization constant)
+   //par[3]=Width (sigma) of convoluted Gaussian function
+   //
+   //In the Landau distribution (represented by the CERNLIB approximation),
+   //the maximum is located at x=-0.22278298 with the location parameter=0.
+   //This shift is corrected within this function, so that the actual
+   //maximum is identical to the MP parameter.
+
+      // Numeric constants
+      Double_t invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
+      Double_t mpshift  = -0.22278298;       // Landau maximum location
+
+      // Control constants
+      Double_t np = 200.0;      // number of convolution steps
+      Double_t sc =   5.0;      // convolution extends to +-sc Gaussian sigmas
+
+      // Variables
+      Double_t xx;
+      Double_t mpc;
+      Double_t fland;
+      Double_t sum = 0.0;
+      Double_t xlow,xupp;
+      Double_t step;
+      Double_t i;
+
+
+      // MP shift correction
+      mpc = par[1] - mpshift * par[0];
+
+      // Range of convolution integral
+      xlow = x[0] - sc * par[3];
+      xupp = x[0] + sc * par[3];
+
+      step = (xupp-xlow) / np;
+
+      // Convolution integral of Landau and Gaussian by sum
+      for(i=1.0; i<=np/2; i++) 
+      {
+                xx = xlow + (i-.5) * step;
+               fland = TMath::Landau(xx,mpc,par[0]) / par[0];
+               sum += fland * TMath::Gaus(x[0],xx,par[3]);
+
+               xx = xupp - (i-.5) * step;
+               fland = TMath::Landau(xx,mpc,par[0]) / par[0];
+               sum += fland * TMath::Gaus(x[0],xx,par[3]);
+      }
+
+      return (par[2] * step * sum * invsq2pi / par[3]);
+}
+
+
+
+
+void GetGainModuleLevel(TString filename,Bool_t normal=1,Int_t ntofit=500, Bool_t grid=0)
+{
+       gROOT->SetStyle("Plain");
+       gStyle->SetPalette(1,0);
+       gStyle->SetOptStat(111111);
+       gStyle->SetOptFit(1);   
+       if(grid)
+               TGrid::Connect("alien://");
+       
+       TFile* file_data=TFile::Open(filename);
+       if(!file_data)
+               return;
+       TList *listin=0x0;
+       listin=(TList*)file_data->Get("output");
+       if(!listin)
+               listin=(TList*)file_data->Get("PWGPPdEdxSSDQA/output");
+       if(!listin)     
+               listin=(TList*)file_data->Get("PWGPPdEdxSSDQA/SSDdEdxQA");
+       if(!listin)     
+               listin=(TList*)file_data->Get("SSDdEdxQA");
+       if(!listin)     
+               return;
+       TH2F* fHistQ=0x0;
+       if (normal)             
+               fHistQ=(TH2F*)listin ->FindObject("QACharge");
+       else
+               fHistQ=(TH2F*)listin ->FindObject("QAChargeCorrected");
+       if(!fHistQ)
+               return;
+       TH2F* fHistCR=(TH2F*)listin ->FindObject("QAChargeRatio");
+       if(!fHistCR)
+               return;
+       TList *listout1=new TList();
+       
+       TList *listout2=new TList();
+       
+       TH1F* fHistMPVs=new TH1F("HistMPVS","HistMPVs;MPV;N",75,70,95);
+       fHistMPVs->SetDirectory(0);
+       listout2->Add(fHistMPVs);
+       
+       TH1F* fHistSL=new TH1F("HistSL","#sigma_{Landau};#sigma_{Landau};N",40,0,16);
+       fHistSL->SetDirectory(0);
+       listout2->Add(fHistSL);
+       
+       TH1F* fHistSG=new TH1F("HistSG","#sigma_{Gaus};#sigma_{Gaus};N",40,0,16);
+       fHistSG->SetDirectory(0);
+       listout2->Add(fHistSG);
+       
+       TH1F* fHistCRmean=new TH1F("HistCRmean","HistCRmean;CRmean;N",200,-1,1);
+       fHistCRmean->SetDirectory(0);
+       listout2->Add(fHistCRmean);
+       
+       TH1F* fHistCRRMS=new TH1F("HistCRRMS","HistCRRMS;CRRMS;N",100,0,1);
+       fHistCRRMS->SetDirectory(0);
+       listout2->Add(fHistCRRMS);
+       
+       TH1F* fHistGainP=new TH1F("HistGainP","HistGainP;CRGainPcorr;N",120,0.5,2.0);
+       fHistGainP->SetDirectory(0);
+       listout2->Add(fHistGainP);
+       
+       TH1F* fHistGainN=new TH1F("HistGainN","HistGainN;CRGainNcorr;N",120,0.5,2.0);
+       fHistGainN->SetDirectory(0);
+       listout2->Add(fHistGainN);
+       
+       TH1F *fMPVGraph = new TH1F("MPVgraph","MPVgraph;Module number;MPV",1698,-0.5,1697.5);
+       fMPVGraph->SetMarkerColor(kRed);
+       fMPVGraph->SetMarkerStyle(22);
+       listout2->Add(fMPVGraph);
+       
+       TH1F *fCRmeanGraph = new TH1F("CRmeangraph","CRmeangraph;Module number;MPV",1698,-0.5,1697.5);
+       fCRmeanGraph->SetMarkerColor(kBlue);
+       fCRmeanGraph->SetMarkerStyle(23);
+       listout2->Add(fCRmeanGraph);
+       
+       Float_t gainP[1698];
+       Float_t gainN[1698];
+       Float_t mpv[1698];
+       Int_t flag[1698];       
+       
+       ofstream outfiletxt;
+       outfiletxt.open("gain.txt");
+        outfiletxt.width(10) ;
+       outfiletxt.setf(outfiletxt.left);
+       outfiletxt<<"MODULE"<<"\t";
+       outfiletxt.width(10);
+       outfiletxt.setf(outfiletxt.left);
+       outfiletxt<<"FLAG"<<"\t";
+       outfiletxt.width(10) ;
+       outfiletxt.setf(outfiletxt.left);
+       outfiletxt<<"GainPcorr"<<"\t";
+       outfiletxt.width(10) ;
+       outfiletxt.setf(outfiletxt.left);
+       outfiletxt<<"GainNcorr"<<"\t";
+       outfiletxt.width(10) ;
+       outfiletxt.setf(outfiletxt.left);
+       outfiletxt<<"MPV"<<endl;
+       
+       
+       ofstream outfiletxtbad;
+       outfiletxtbad.open("badModules.txt");
+       
+       
+       
+       
+       for (int i =0;i<1698;i++)
+       {
+               cout<<i<<endl;
+               TString tmpQ("Q");
+               tmpQ+=i;
+               TString tmpCR("CR");
+               tmpCR+=i;
+               TH1D* fHist1DCR= fHistCR->ProjectionY(tmpCR,i+1,i+1);
+               Double_t mean=fHist1DCR->GetMean();
+               if(!(TMath::Abs(mean)<1.0)||fHist1DCR->GetEntries()<10)
+               {               
+                       flag[i]=-2;
+                       gainN[i]=1.0;
+                       gainP[i]=1.0;
+                       mpv[i]=1.0;
+                       continue;
+               }
+               fHistCRmean->Fill(mean);
+               fHistCRRMS->Fill(fHist1DCR->GetRMS());
+               gainN[i]=1.0/(1.0+mean);
+               gainP[i]=1.0/(1.0-mean);
+               fHistGainP->Fill(gainP[i]);
+               fHistGainN->Fill(gainN[i]);
+               fCRmeanGraph->SetBinContent(i+1,mean);
+               fCRmeanGraph->SetBinError(i+1,fHist1DCR->GetRMS());
+               
+               
+               TH1D* fHist1DQ=fHistQ->ProjectionY(tmpQ,i+1,i+1);
+                fHist1DQ->SetDirectory(0);
+                listout1->Add(fHist1DQ);
+               if(fHist1DQ->GetEntries()<ntofit)
+               {
+                       flag[i]=-1;
+                       mpv[i]=1.0;
+                        outfiletxtbad<<"Low statistic \t module= "<<i<<" netries="<<fHist1DQ->GetEntries()<<endl;
+                       continue;
+               }
+               else
+               {
+                       tmpQ+="fit";
+                       Float_t range=fHist1DQ->GetBinCenter(fHist1DQ->GetMaximumBin());
+                       TF1 *f1 = new TF1(tmpQ,convolution,range*0.45,range*3.0,4);
+                       f1->SetParameters(7.0,range,1.0,5.5);
+                       Float_t normalization=fHist1DQ->GetEntries()*fHist1DQ->GetXaxis()->GetBinWidth(2)/f1->Integral(range*0.45,range*3.0);
+                       f1->SetParameters(7.0,range,normalization,5.5);
+                       //f1->SetParameters(7.0,range,fHist1DQ->GetMaximum(),5.5);
+                       f1->SetParNames("sigma Landau","MPV","N","sigma Gaus");
+                       f1->SetParLimits(0,2.0,100.0);
+                       f1->SetParLimits(3,0.0,100.0);
+                       if(fHist1DQ->Fit(tmpQ,"BRQ")==0)
+                       {
+                               mpv[i]=f1->GetParameter(1);
+                               fHistMPVs->Fill(mpv[i]);
+                               fHistSL->Fill(f1->GetParameter(0));
+                               fHistSG->Fill(f1->GetParameter(3));
+                               flag[i]=1;              
+                               fMPVGraph->SetBinContent(i+1,f1->GetParameter(1));
+                               fMPVGraph->SetBinError(i+1,f1->GetParError(1));
+                               if(mpv[i]<75.0)
+                               {
+                                       outfiletxtbad<<"MPV lower than 75 \t module="<<i<<endl;
+                                       flag[i]=0;
+                               }       
+                               if(mpv[i]>100.0)
+                               {
+                                       outfiletxtbad<<"MPV higher than 100 \t module="<<i<<endl;
+                                       flag[i]=0;
+                                       
+                               }
+                               if(f1->GetParError(1)>1.0)
+                               {
+                                       outfiletxtbad<<"MPV high error on MPV  \t module="<<i<<endl;                            
+                                       flag[i]=0;
+                               }
+                       }
+                       else
+                       {
+                               mpv[i]=1;
+                               flag[i]=0;
+                                outfiletxtbad<<"BAD FIT \t module="<<i<<endl;
+                               //continue;
+                       }       
+               }       
+       }       
+       
+       for (int i=0;i<1698;i++)
+       {       
+               outfiletxt.setf(outfiletxt.scientific);
+               outfiletxt.precision(2);        
+                outfiletxt.width(10) ;
+               outfiletxt.setf(outfiletxt.left);
+               outfiletxt<<i<<"\t";
+                outfiletxt.width(10) ;
+               outfiletxt.setf(outfiletxt.left);
+               outfiletxt<<flag[i]<<"\t";
+                outfiletxt.width(10) ;
+               outfiletxt.setf(outfiletxt.left);
+               outfiletxt<<gainP[i]<<"\t";
+                outfiletxt.width(10) ;
+               outfiletxt.setf(outfiletxt.left);
+               outfiletxt<<gainN[i]<<"\t";
+                outfiletxt.width(10) ;
+               outfiletxt.setf(outfiletxt.left);
+               outfiletxt<<mpv[i]<<endl;
+       }
+                
+       TCanvas *c1 = new TCanvas("1","1",1200,800);
+       c1->Divide(2,1);
+       c1->cd(1);
+       fHistQ->Draw("colz");
+       c1->cd(2);
+       fHistCR->Draw("colz");
+       
+       
+       
+       TFile* fout1=TFile::Open("gain_all_fits.root","recreate");
+       listout1->Write("output",TObject::kSingleKey);  
+       fout1->Close();
+       
+       TFile* fout2=TFile::Open("gain_control_plots.root","recreate");
+       listout2->Write("output",TObject::kSingleKey);  
+       fout2->Close();
+       
+       
+       
+}