Reference Manual

Go to the documentation of this file.
00001 // Original Author:  Loic QUERTENMONT
00002 //         Created:  Mon Nov  16 08:55:18 CET 2009
00003 
00004 #include <memory>
00005 
00006 #include "FWCore/Framework/interface/Frameworkfwd.h"
00007 #include "FWCore/Framework/interface/EDAnalyzer.h"
00008 #include "FWCore/Framework/interface/Event.h"
00009 #include "FWCore/Framework/interface/MakerMacros.h"
00010 #include "FWCore/Framework/interface/ESHandle.h"
00011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00012 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00013 #include "FWCore/Utilities/interface/Exception.h"
00014 
00015 #include "Geometry/CommonDetUnit/interface/GeomDetUnit.h"
00016 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
00017 #include "Geometry/CommonTopologies/interface/StripTopology.h"
00018 #include "DataFormats/GeometrySurface/interface/TrapezoidalPlaneBounds.h"
00019 #include "DataFormats/GeometrySurface/interface/RectangularPlaneBounds.h"
00020 
00021 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
00022 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
00023 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
00024 #include "Geometry/TrackerGeometryBuilder/interface/PixelGeomDetUnit.h"
00025 #include "Geometry/TrackerNumberingBuilder/interface/GeometricDet.h"
00026 #include "Geometry/CommonDetUnit/interface/TrackingGeometry.h"
00027 
00028 #include "DataFormats/SiStripCluster/interface/SiStripClusterCollection.h"
00029 
00030 #include "CalibFormats/SiStripObjects/interface/SiStripDetCabling.h"
00031 #include "CalibTracker/Records/interface/SiStripDetCablingRcd.h"
00032 
00033 #include "DataFormats/FEDRawData/interface/FEDNumbering.h"
00034 #include "DataFormats/TrackReco/interface/Track.h"
00035 #include "DataFormats/TrackReco/interface/TrackFwd.h"
00036 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h"
00037 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
00038 #include "DataFormats/SiStripDetId/interface/SiStripSubStructure.h"
00039 #include "DataFormats/DetId/interface/DetId.h"
00040 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
00041 #include "DataFormats/SiStripDetId/interface/TECDetId.h"
00042 #include "DataFormats/SiStripDetId/interface/TIBDetId.h"
00043 #include "DataFormats/SiStripDetId/interface/TIDDetId.h"
00044 #include "DataFormats/SiStripDetId/interface/TOBDetId.h"
00045 #include "DataFormats/TrackReco/interface/DeDxHit.h"
00046 #include "DataFormats/TrackReco/interface/TrackDeDxHits.h"
00047 
00048 #include "CommonTools/ConditionDBWriter/interface/ConditionDBWriter.h"
00049 #include "CondFormats/SiStripObjects/interface/SiStripApvGain.h"
00050 
00051 #include "TrackingTools/PatternTools/interface/Trajectory.h"
00052 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
00053 
00054 #include "DQMServices/Core/interface/DQMStore.h"
00055 #include "DQMServices/Core/interface/MonitorElement.h"
00056 #include "FWCore/ServiceRegistry/interface/Service.h"
00057 #include "CommonTools/UtilAlgos/interface/TFileService.h"
00058 
00059 #include "CalibFormats/SiStripObjects/interface/SiStripGain.h"
00060 #include "CalibTracker/Records/interface/SiStripGainRcd.h"
00061 
00062 #include "CalibFormats/SiStripObjects/interface/SiStripQuality.h"
00063 #include "CalibTracker/Records/interface/SiStripQualityRcd.h"
00064 
00065 
00066 #include "TFile.h"
00067 #include "TObjString.h"
00068 #include "TString.h"
00069 #include "TH1F.h"
00070 #include "TH2F.h"
00071 #include "TProfile.h"
00072 #include "TF1.h"
00073 #include "TROOT.h"
00074 #include "TTree.h"
00075 #include "TChain.h"
00076 
00077 #include <ext/hash_map>
00078 
00079 
00080 
00081 using namespace edm;
00082 using namespace reco;
00083 using namespace std;
00084 using namespace __gnu_cxx;
00085 using __gnu_cxx::hash_map;
00086 using __gnu_cxx::hash;
00087 
00088 struct stAPVGain{
00089    unsigned int Index; 
00090    unsigned int Bin;
00091    unsigned int DetId;
00092    unsigned int APVId;
00093    unsigned int SubDet;
00094    float        x;
00095    float        y;
00096    float        z;
00097    float        Eta;
00098    float        R;
00099    float        Phi;
00100    float        Thickness;
00101    double       FitMPV;
00102    double       FitMPVErr;
00103    double       FitWidth;
00104    double       FitWidthErr;
00105    double       FitChi2;
00106    double       Gain;
00107    double       CalibGain;
00108    double       PreviousGain;
00109    double       NEntries;
00110    TH1F*        HCharge;
00111    TH1F*        HChargeN;
00112    bool         isMasked;
00113 };
00114 
00115 class SiStripGainFromCalibTree : public ConditionDBWriter<SiStripApvGain> {
00116    public:
00117       explicit SiStripGainFromCalibTree(const edm::ParameterSet&);
00118       ~SiStripGainFromCalibTree();
00119 
00120 
00121    private:
00122       virtual void algoBeginJob      (const edm::EventSetup&);
00123       virtual void algoEndJob        ();
00124       virtual void algoAnalyze       (const edm::Event &, const edm::EventSetup &);
00125 
00126               void algoAnalyzeTheTree();
00127               void algoComputeMPVandGain();
00128 
00129               void getPeakOfLandau(TH1* InputHisto, double* FitResults, double LowRange=50, double HighRange=5400);
00130               bool IsGoodLandauFit(double* FitResults); 
00131               void storeOnTree();
00132 
00133               void MakeCalibrationMap();
00134 
00135 
00136       SiStripApvGain* getNewObject();
00137       edm::Service<TFileService> tfs;
00138 
00139       double       MinNrEntries;
00140       double       MaxMPVError;
00141       double       MaxChi2OverNDF;
00142       double       MinTrackMomentum;
00143       double       MaxTrackMomentum;
00144       double       MinTrackEta;
00145       double       MaxTrackEta;
00146       unsigned int MaxNrStrips;
00147       unsigned int MinTrackHits;
00148       double       MaxTrackChiOverNdf;
00149       bool         AllowSaturation;
00150       bool         FirstSetOfConstants;
00151       bool         Validation;
00152       bool         OldGainRemoving;
00153       int          CalibrationLevel;
00154 
00155       bool         useCalibration;
00156       string       m_calibrationPath;
00157 
00158       std::string  OutputGains;
00159       vector<string> VInputFiles;
00160 
00161       TH2F*        Charge_Vs_Index;
00162       TH2F*        Charge_Vs_Index_Absolute;
00163       TH2F*        Charge_Vs_PathlengthTIB;
00164       TH2F*        Charge_Vs_PathlengthTOB;
00165       TH2F*        Charge_Vs_PathlengthTIDP;
00166       TH2F*        Charge_Vs_PathlengthTIDM;
00167       TH2F*        Charge_Vs_PathlengthTECP1;
00168       TH2F*        Charge_Vs_PathlengthTECP2;
00169       TH2F*        Charge_Vs_PathlengthTECM1;
00170       TH2F*        Charge_Vs_PathlengthTECM2;
00171 
00172       unsigned int NEvent;    
00173       unsigned int NTrack;
00174       unsigned int NCluster;
00175       unsigned int SRun;
00176       unsigned int ERun;
00177       unsigned int GOOD;
00178       unsigned int BAD;
00179 
00180    private :
00181       class isEqual{
00182          public:
00183                  template <class T> bool operator () (const T& PseudoDetId1, const T& PseudoDetId2) { return PseudoDetId1==PseudoDetId2; }
00184       };
00185 
00186       std::vector<stAPVGain*> APVsCollOrdered;
00187       __gnu_cxx::hash_map<unsigned int, stAPVGain*,  __gnu_cxx::hash<unsigned int>, isEqual > APVsColl;
00188 };
00189 
00190 SiStripGainFromCalibTree::SiStripGainFromCalibTree(const edm::ParameterSet& iConfig) : ConditionDBWriter<SiStripApvGain>(iConfig)
00191 {
00192    OutputGains         = iConfig.getParameter<std::string>("OutputGains");
00193 
00194    MinNrEntries        = iConfig.getUntrackedParameter<double>  ("minNrEntries"       ,  20);
00195    MaxMPVError         = iConfig.getUntrackedParameter<double>  ("maxMPVError"        ,  500.0);
00196    MaxChi2OverNDF      = iConfig.getUntrackedParameter<double>  ("maxChi2OverNDF"     ,  5.0);
00197    MinTrackMomentum    = iConfig.getUntrackedParameter<double>  ("minTrackMomentum"   ,  3.0);
00198    MaxTrackMomentum    = iConfig.getUntrackedParameter<double>  ("maxTrackMomentum"   ,  99999.0);
00199    MinTrackEta         = iConfig.getUntrackedParameter<double>  ("minTrackEta"        , -5.0);
00200    MaxTrackEta         = iConfig.getUntrackedParameter<double>  ("maxTrackEta"        ,  5.0);
00201    MaxNrStrips         = iConfig.getUntrackedParameter<unsigned>("maxNrStrips"        ,  2);
00202    MinTrackHits        = iConfig.getUntrackedParameter<unsigned>("MinTrackHits"       ,  8);
00203    MaxTrackChiOverNdf  = iConfig.getUntrackedParameter<double>  ("MaxTrackChiOverNdf" ,  3);
00204    AllowSaturation     = iConfig.getUntrackedParameter<bool>    ("AllowSaturation"    ,  false);
00205    FirstSetOfConstants = iConfig.getUntrackedParameter<bool>    ("FirstSetOfConstants",  true);
00206    Validation          = iConfig.getUntrackedParameter<bool>    ("Validation"         ,  false);
00207    OldGainRemoving     = iConfig.getUntrackedParameter<bool>    ("OldGainRemoving"    ,  false);
00208 
00209    CalibrationLevel    = iConfig.getUntrackedParameter<int>     ("CalibrationLevel"   ,  0);
00210    VInputFiles         = iConfig.getParameter<vector<string> >  ("InputFiles");
00211 
00212 
00213    useCalibration      = iConfig.getUntrackedParameter<bool>("UseCalibration", false);
00214    m_calibrationPath   = iConfig.getUntrackedParameter<string>("calibrationPath");
00215 
00216 
00217 }
00218 
00219 
00220 
00221 
00222 void
00223 SiStripGainFromCalibTree::algoBeginJob(const edm::EventSetup& iSetup)
00224 {
00225    Charge_Vs_Index           = tfs->make<TH2F>("Charge_Vs_Index"          , "Charge_Vs_Index"          , 72785, 0   , 72784,1000,0,2000);
00226    Charge_Vs_Index_Absolute  = tfs->make<TH2F>("Charge_Vs_Index_Absolute" , "Charge_Vs_Index_Absolute" , 72785, 0   , 72784, 500,0,2000);
00227    Charge_Vs_PathlengthTIB   = tfs->make<TH2F>("Charge_Vs_PathlengthTIB"  , "Charge_Vs_PathlengthTIB"  , 20   , 0.3 , 1.3  , 250,0,2000);
00228    Charge_Vs_PathlengthTOB   = tfs->make<TH2F>("Charge_Vs_PathlengthTOB"  , "Charge_Vs_PathlengthTOB"  , 20   , 0.3 , 1.3  , 250,0,2000);
00229    Charge_Vs_PathlengthTIDP  = tfs->make<TH2F>("Charge_Vs_PathlengthTIDP" , "Charge_Vs_PathlengthTIDP" , 20   , 0.3 , 1.3  , 250,0,2000);
00230    Charge_Vs_PathlengthTIDM  = tfs->make<TH2F>("Charge_Vs_PathlengthTIDM" , "Charge_Vs_PathlengthTIDM" , 20   , 0.3 , 1.3  , 250,0,2000);
00231    Charge_Vs_PathlengthTECP1 = tfs->make<TH2F>("Charge_Vs_PathlengthTECP1", "Charge_Vs_PathlengthTECP1", 20   , 0.3 , 1.3  , 250,0,2000);
00232    Charge_Vs_PathlengthTECP2 = tfs->make<TH2F>("Charge_Vs_PathlengthTECP2", "Charge_Vs_PathlengthTECP2", 20   , 0.3 , 1.3  , 250,0,2000);
00233    Charge_Vs_PathlengthTECM1 = tfs->make<TH2F>("Charge_Vs_PathlengthTECM1", "Charge_Vs_PathlengthTECM1", 20   , 0.3 , 1.3  , 250,0,2000);
00234    Charge_Vs_PathlengthTECM2 = tfs->make<TH2F>("Charge_Vs_PathlengthTECM2", "Charge_Vs_PathlengthTECM2", 20   , 0.3 , 1.3  , 250,0,2000);
00235 
00236    edm::ESHandle<TrackerGeometry> tkGeom;
00237    iSetup.get<TrackerDigiGeometryRecord>().get( tkGeom );
00238    vector<GeomDet*> Det = tkGeom->dets();
00239 
00240    edm::ESHandle<SiStripGain> gainHandle;
00241    iSetup.get<SiStripGainRcd>().get(gainHandle);
00242    if(!gainHandle.isValid()){printf("\n#####################\n\nERROR --> gainHandle is not valid\n\n#####################\n\n");exit(0);}
00243  
00244 //   size_t numberOfTag = gainHandle->getNumberOfTags();
00245    for(unsigned int i=0;i<gainHandle->getNumberOfTags();i++){
00246       printf("Reccord %i --> Rcd Name = %s    Label Name = %s\n",i,gainHandle->getRcdName(i).c_str(), gainHandle->getLabelName(i).c_str());
00247    }
00248  
00249 
00250    edm::ESHandle<SiStripQuality> SiStripQuality_;
00251    iSetup.get<SiStripQualityRcd>().get(SiStripQuality_);
00252 
00253 
00254    unsigned int Index=0;
00255    for(unsigned int i=0;i<Det.size();i++){
00256       DetId  Detid  = Det[i]->geographicalId(); 
00257       int    SubDet = Detid.subdetId();
00258 
00259       if( SubDet == StripSubdetector::TIB ||  SubDet == StripSubdetector::TID ||
00260           SubDet == StripSubdetector::TOB ||  SubDet == StripSubdetector::TEC  ){
00261 
00262           StripGeomDetUnit* DetUnit     = dynamic_cast<StripGeomDetUnit*> (Det[i]);
00263           if(!DetUnit)continue;
00264 
00265           const StripTopology& Topo     = DetUnit->specificTopology();  
00266           unsigned int         NAPV     = Topo.nstrips()/128;
00267           for(unsigned int j=0;j<NAPV;j++){
00268                 stAPVGain* APV = new stAPVGain;
00269                 APV->Index         = Index;
00270                 APV->Bin           = Charge_Vs_Index->GetXaxis()->FindBin(APV->Index);
00271                 APV->DetId         = Detid.rawId();
00272                 APV->APVId         = j;
00273                 APV->SubDet        = SubDet;
00274                 APV->FitMPV        = -1;
00275                 APV->FitMPVErr     = -1;
00276                 APV->FitWidth      = -1;
00277                 APV->FitWidthErr   = -1;
00278                 APV->FitChi2       = -1;
00279                 APV->Gain          = -1;
00280                 APV->PreviousGain  = 1;
00281                 APV->x             = DetUnit->position().basicVector().x();
00282                 APV->y             = DetUnit->position().basicVector().y();
00283                 APV->z             = DetUnit->position().basicVector().z();
00284                 APV->Eta           = DetUnit->position().basicVector().eta();
00285                 APV->Phi           = DetUnit->position().basicVector().phi();
00286                 APV->R             = DetUnit->position().basicVector().transverse();
00287                 APV->Thickness     = DetUnit->surface().bounds().thickness();
00288                 APV->NEntries      = 0;
00289                 APV->isMasked      = SiStripQuality_->IsApvBad(Detid.rawId(),j);
00290 
00291                 if(!FirstSetOfConstants){
00292                    if(gainHandle->getNumberOfTags()!=2){printf("ERROR: NUMBER OF GAIN TAG IS EXPECTED TO BE 2\n");fflush(stdout);exit(0);};                
00293                    APV->PreviousGain  = gainHandle->getApvGain(APV->APVId,gainHandle->getRange(APV->DetId, 1),1);
00294                    printf("DETID = %7i APVID=%1i Previous Gain=%8.4f\n",APV->DetId,APV->APVId,APV->PreviousGain);
00295                 }
00296 
00297                 APVsCollOrdered.push_back(APV);
00298                 APVsColl[(APV->DetId<<3) | APV->APVId] = APV;
00299                 Index++;
00300           }
00301       }
00302    }
00303 
00304    MakeCalibrationMap();
00305 
00306 
00307    NEvent     = 0;
00308    NTrack     = 0;
00309    NCluster   = 0;
00310    SRun       = 1<<31;
00311    ERun       = 0;
00312    GOOD       = 0;
00313    BAD        = 0;
00314 
00315    algoAnalyzeTheTree();
00316    algoComputeMPVandGain();
00317 }
00318 
00319 
00320 void
00321 SiStripGainFromCalibTree::algoAnalyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
00322 {
00323 }
00324 
00325 void 
00326 SiStripGainFromCalibTree::algoEndJob() {
00327 }
00328 
00329 
00330 void SiStripGainFromCalibTree::getPeakOfLandau(TH1* InputHisto, double* FitResults, double LowRange, double HighRange)
00331 { 
00332    FitResults[0]         = -0.5;  //MPV
00333    FitResults[1]         =  0;    //MPV error
00334    FitResults[2]         = -0.5;  //Width
00335    FitResults[3]         =  0;    //Width error
00336    FitResults[4]         = -0.5;  //Fit Chi2/NDF
00337 
00338    if( InputHisto->GetEntries() < MinNrEntries)return;
00339 
00340    // perform fit with standard landau
00341    TF1* MyLandau = new TF1("MyLandau","landau",LowRange, HighRange);
00342    MyLandau->SetParameter(1,300);
00343    InputHisto->Fit("MyLandau","0QR WW");
00344 
00345    // MPV is parameter 1 (0=constant, 1=MPV, 2=Sigma)
00346    FitResults[0]         = MyLandau->GetParameter(1);  //MPV
00347    FitResults[1]         = MyLandau->GetParError(1);   //MPV error
00348    FitResults[2]         = MyLandau->GetParameter(2);  //Width
00349    FitResults[3]         = MyLandau->GetParError(2);   //Width error
00350    FitResults[4]         = MyLandau->GetChisquare() / MyLandau->GetNDF();  //Fit Chi2/NDF
00351 
00352    delete MyLandau;
00353 }
00354 
00355 bool SiStripGainFromCalibTree::IsGoodLandauFit(double* FitResults){
00356    if(FitResults[0] < 2             )return false;
00357    if(FitResults[1] > MaxMPVError   )return false;
00358    if(FitResults[4] > MaxChi2OverNDF)return false;
00359    return true;   
00360 }
00361 
00362 
00363 void SiStripGainFromCalibTree::algoAnalyzeTheTree()
00364 {
00365    for(unsigned int i=0;i<VInputFiles.size();i++){
00366       printf("Openning file %3i/%3i --> %s\n",i+1, (int)VInputFiles.size(), (char*)(VInputFiles[i].c_str())); fflush(stdout);
00367       TChain* tree = new TChain("gainCalibrationTree/tree");
00368       tree->Add(VInputFiles[i].c_str());
00369 
00370       TString EventPrefix("");
00371       TString EventSuffix("");
00372 
00373       TString TrackPrefix("track");
00374       TString TrackSuffix("");
00375 
00376       TString CalibPrefix("GainCalibration");
00377       TString CalibSuffix("");
00378 
00379       unsigned int                 eventnumber    = 0;    tree->SetBranchAddress(EventPrefix + "event"          + EventSuffix, &eventnumber   , NULL);
00380       unsigned int                 runnumber      = 0;    tree->SetBranchAddress(EventPrefix + "run"            + EventSuffix, &runnumber     , NULL);
00381       std::vector<bool>*           TrigTech    = 0;    tree->SetBranchAddress(EventPrefix + "TrigTech"    + EventSuffix, &TrigTech   , NULL);
00382 
00383       std::vector<double>*         trackchi2ndof  = 0;    tree->SetBranchAddress(TrackPrefix + "chi2ndof"       + TrackSuffix, &trackchi2ndof , NULL);
00384       std::vector<float>*          trackp         = 0;    tree->SetBranchAddress(TrackPrefix + "momentum"       + TrackSuffix, &trackp        , NULL);
00385       std::vector<float>*          trackpt        = 0;    tree->SetBranchAddress(TrackPrefix + "pt"             + TrackSuffix, &trackpt       , NULL);
00386       std::vector<double>*         tracketa       = 0;    tree->SetBranchAddress(TrackPrefix + "eta"            + TrackSuffix, &tracketa      , NULL);
00387       std::vector<double>*         trackphi       = 0;    tree->SetBranchAddress(TrackPrefix + "phi"            + TrackSuffix, &trackphi      , NULL);
00388       std::vector<unsigned int>*   trackhitsvalid = 0;    tree->SetBranchAddress(TrackPrefix + "hitsvalid"      + TrackSuffix, &trackhitsvalid, NULL);
00389 
00390       std::vector<int>*            trackindex     = 0;    tree->SetBranchAddress(CalibPrefix + "trackindex"     + CalibSuffix, &trackindex    , NULL);
00391       std::vector<unsigned int>*   rawid          = 0;    tree->SetBranchAddress(CalibPrefix + "rawid"          + CalibSuffix, &rawid         , NULL);
00392       std::vector<float>*          localdirx      = 0;    tree->SetBranchAddress(CalibPrefix + "localdirx"      + CalibSuffix, &localdirx     , NULL);
00393       std::vector<float>*          localdiry      = 0;    tree->SetBranchAddress(CalibPrefix + "localdiry"      + CalibSuffix, &localdiry     , NULL);
00394       std::vector<float>*          localdirz      = 0;    tree->SetBranchAddress(CalibPrefix + "localdirz"      + CalibSuffix, &localdirz     , NULL);
00395       std::vector<unsigned short>* firststrip     = 0;    tree->SetBranchAddress(CalibPrefix + "firststrip"     + CalibSuffix, &firststrip    , NULL);
00396       std::vector<unsigned short>* nstrips        = 0;    tree->SetBranchAddress(CalibPrefix + "nstrips"        + CalibSuffix, &nstrips       , NULL);
00397       std::vector<bool>*           saturation     = 0;    tree->SetBranchAddress(CalibPrefix + "saturation"     + CalibSuffix, &saturation    , NULL);
00398       std::vector<bool>*           overlapping    = 0;    tree->SetBranchAddress(CalibPrefix + "overlapping"    + CalibSuffix, &overlapping   , NULL);
00399       std::vector<bool>*           farfromedge    = 0;    tree->SetBranchAddress(CalibPrefix + "farfromedge"    + CalibSuffix, &farfromedge   , NULL);
00400       std::vector<unsigned int>*   charge         = 0;    tree->SetBranchAddress(CalibPrefix + "charge"         + CalibSuffix, &charge        , NULL);
00401       std::vector<float>*          path           = 0;    tree->SetBranchAddress(CalibPrefix + "path"           + CalibSuffix, &path          , NULL);
00402       std::vector<float>*          chargeoverpath = 0;    tree->SetBranchAddress(CalibPrefix + "chargeoverpath" + CalibSuffix, &chargeoverpath, NULL);
00403       std::vector<unsigned char>*  amplitude      = 0;    tree->SetBranchAddress(CalibPrefix + "amplitude"      + CalibSuffix, &amplitude     , NULL);
00404       std::vector<double>*         gainused       = 0;    tree->SetBranchAddress(CalibPrefix + "gainused"       + CalibSuffix, &gainused      , NULL);
00405 
00406 
00407       printf("Number of Events = %i + %i = %i\n",NEvent,(unsigned int)tree->GetEntries(),(unsigned int)(NEvent+tree->GetEntries()));
00408       printf("Progressing Bar              :0%%       20%%       40%%       60%%       80%%       100%%\n");
00409       printf("Looping on the Tree          :");
00410       int TreeStep = tree->GetEntries()/50;if(TreeStep<=1)TreeStep=1;
00411       for (unsigned int ientry = 0; ientry < tree->GetEntries(); ientry++) {
00412 //      for (unsigned int ientry = 0; ientry < tree->GetEntries() && ientry<50000; ientry++) {
00413       if(ientry%TreeStep==0){printf(".");fflush(stdout);}
00414          tree->GetEntry(ientry);
00415 
00416          if(runnumber<SRun)SRun=runnumber;
00417          if(runnumber>ERun)ERun=runnumber;
00418 
00419          NEvent++;
00420          NTrack+=(*trackp).size();
00421 
00422          unsigned int FirstAmplitude=0;
00423          for(unsigned int i=0;i<(*chargeoverpath).size();i++){
00424             FirstAmplitude+=(*nstrips)[i];
00425                  int    TI = (*trackindex)[i];
00426             if((*tracketa      )[TI]  < MinTrackEta        )continue;
00427             if((*tracketa      )[TI]  > MaxTrackEta        )continue;
00428             if((*trackp        )[TI]  < MinTrackMomentum   )continue;
00429             if((*trackp        )[TI]  > MaxTrackMomentum   )continue;
00430             if((*trackhitsvalid)[TI]  < MinTrackHits       )continue;
00431             if((*trackchi2ndof )[TI]  > MaxTrackChiOverNdf )continue;
00432             if((*farfromedge)[i]        == false           )continue;
00433             if((*overlapping)[i]        == true            )continue;
00434             if((*saturation )[i]        && !AllowSaturation)continue;
00435             if((*nstrips    )[i]      > MaxNrStrips        )continue;
00436             NCluster++;
00437 
00438             stAPVGain* APV = APVsColl[((*rawid)[i]<<3) | ((*firststrip)[i]/128)];
00439 
00440 //            printf("detId=%7i run=%7i event=%9i charge=%5i cs=%3i\n",(*rawid)[i],runnumber,eventnumber,(*charge)[i],(*nstrips)[i]);
00441 
00442 
00443 
00444             //double trans = atan2((*localdiry)[i],(*localdirx)[i])*(180/3.14159265);
00445             //double alpha = acos ((*localdirx)[i] / sqrt( pow((*localdirx)[i],2) +  pow((*localdirz)[i],2) ) ) * (180/3.14159265);
00446             //double beta  = acos ((*localdiry)[i] / sqrt( pow((*localdirx)[i],2) +  pow((*localdirz)[i],2) ) ) * (180/3.14159265);
00447           
00448             //printf("NStrip = %i : Charge = %i --> Path = %f  --> ChargeOverPath=%f\n",(*nstrips)[i],(*charge)[i],(*path)[i],(*chargeoverpath)[i]);
00449             //printf("Amplitudes: ");
00450             //for(unsigned int a=0;a<(*nstrips)[i];a++){printf("%i ",(*amplitude)[FirstAmplitude+a]);}
00451             //printf("\n");
00452 
00453             
00454             int Charge = 0;
00455             if(useCalibration || !FirstSetOfConstants){
00456                bool Saturation = false;
00457                for(unsigned int s=0;s<(*nstrips)[i];s++){
00458                   int StripCharge =  (*amplitude)[FirstAmplitude-(*nstrips)[i]+s];
00459                   if(useCalibration && !FirstSetOfConstants){ StripCharge=(int)(StripCharge*(APV->PreviousGain/APV->CalibGain));
00460                   }else if(useCalibration){                   StripCharge=(int)(StripCharge/APV->CalibGain);
00461                   }else if(!FirstSetOfConstants){             StripCharge=(int)(StripCharge*APV->PreviousGain);}
00462                   if(StripCharge>1024){
00463                      StripCharge = 255;
00464                      Saturation = true;
00465                   }else if(StripCharge>254){
00466                      StripCharge = 254;
00467                      Saturation = true;
00468                   }
00469                   Charge += StripCharge;
00470                }
00471                if(Saturation && !AllowSaturation)continue;
00472             }else{
00473                Charge = (*charge)[i];
00474             }
00475 
00476             //printf("ChargeDifference = %i Vs %i with Gain = %f\n",(*charge)[i],Charge,APV->CalibGain);
00477 
00478             double ClusterChargeOverPath   =  ( (double) Charge )/(*path)[i] ;
00479             if(Validation)     {ClusterChargeOverPath/=(*gainused)[i];}
00480             if(OldGainRemoving){ClusterChargeOverPath*=(*gainused)[i];}
00481             Charge_Vs_Index_Absolute->Fill(APV->Index,Charge);   
00482             Charge_Vs_Index         ->Fill(APV->Index,ClusterChargeOverPath);
00483 
00484 
00485                   if(APV->SubDet==StripSubdetector::TIB){ Charge_Vs_PathlengthTIB  ->Fill((*path)[i],Charge); 
00486             }else if(APV->SubDet==StripSubdetector::TOB){ Charge_Vs_PathlengthTOB  ->Fill((*path)[i],Charge);
00487             }else if(APV->SubDet==StripSubdetector::TID){
00488                      if(APV->Eta<0){                      Charge_Vs_PathlengthTIDM ->Fill((*path)[i],Charge);
00489                }else if(APV->Eta>0){                      Charge_Vs_PathlengthTIDP ->Fill((*path)[i],Charge);
00490                }
00491             }else if(APV->SubDet==StripSubdetector::TEC){
00492                      if(APV->Eta<0){
00493                         if(APV->Thickness<0.04){          Charge_Vs_PathlengthTECM1->Fill((*path)[i],Charge);
00494                   }else if(APV->Thickness>0.04){          Charge_Vs_PathlengthTECM2->Fill((*path)[i],Charge);
00495                   }
00496                }else if(APV->Eta>0){
00497                         if(APV->Thickness<0.04){          Charge_Vs_PathlengthTECP1->Fill((*path)[i],Charge);
00498                   }else if(APV->Thickness>0.04){          Charge_Vs_PathlengthTECP2->Fill((*path)[i],Charge);
00499                   }
00500                }
00501             }
00502 
00503          }// END OF ON-CLUSTER LOOP
00504       }printf("\n");// END OF EVENT LOOP
00505 
00506    }
00507 }
00508 
00509 
00510 
00511 void SiStripGainFromCalibTree::algoComputeMPVandGain() {
00512    unsigned int I=0;
00513    TH1D* Proj = NULL;
00514    double FitResults[5];
00515    double MPVmean = 300;
00516 
00517    printf("Progressing Bar              :0%%       20%%       40%%       60%%       80%%       100%%\n");
00518    printf("Fitting Charge Distribution  :");
00519    int TreeStep = APVsColl.size()/50;
00520    for(__gnu_cxx::hash_map<unsigned int, stAPVGain*,  __gnu_cxx::hash<unsigned int>, isEqual >::iterator it = APVsColl.begin();it!=APVsColl.end();it++,I++){
00521    if(I%TreeStep==0){printf(".");fflush(stdout);}
00522    //if(I>1000)break;
00523       stAPVGain* APV = it->second;
00524 
00525       Proj = (TH1D*)(Charge_Vs_Index->ProjectionY("",APV->Bin,APV->Bin,"e"));
00526       if(!Proj)continue;
00527 
00528       if(CalibrationLevel==0){
00529       }else if(CalibrationLevel==1){
00530          int SecondAPVId = APV->APVId;
00531          if(SecondAPVId%2==0){    SecondAPVId = SecondAPVId+1; }else{ SecondAPVId = SecondAPVId-1; }
00532          stAPVGain* APV2 = APVsColl[(APV->DetId<<3) | SecondAPVId];
00533          TH1D* Proj2 = (TH1D*)(Charge_Vs_Index->ProjectionY("",APV2->Bin,APV2->Bin,"e"));
00534          if(Proj2){Proj->Add(Proj2,1);delete Proj2;}
00535       }else if(CalibrationLevel==2){
00536           for(unsigned int i=0;i<6;i++){
00537             __gnu_cxx::hash_map<unsigned int, stAPVGain*,  __gnu_cxx::hash<unsigned int>, isEqual >::iterator tmpit;
00538             tmpit = APVsColl.find((APV->DetId<<3) | i);
00539             if(tmpit==APVsColl.end())continue;
00540             stAPVGain* APV2 = tmpit->second;
00541             if(APV2->DetId != APV->DetId || APV2->APVId == APV->APVId)continue;            
00542             TH1D* Proj2 = (TH1D*)(Charge_Vs_Index->ProjectionY("",APV2->Bin,APV2->Bin,"e"));
00543 //            if(Proj2 && APV->DetId==369171124)printf("B) DetId %6i APVId %1i --> NEntries = %f\n",APV2->DetId, APV2->APVId, Proj2->GetEntries());
00544             if(Proj2){Proj->Add(Proj2,1);delete Proj2;}
00545           }          
00546       }else{
00547          CalibrationLevel = 0;
00548          printf("Unknown Calibration Level, will assume %i\n",CalibrationLevel);
00549       }
00550 
00551       getPeakOfLandau(Proj,FitResults);
00552       APV->FitMPV      = FitResults[0];
00553       APV->FitMPVErr   = FitResults[1];
00554       APV->FitWidth    = FitResults[2];
00555       APV->FitWidthErr = FitResults[3];
00556       APV->FitChi2     = FitResults[4];
00557       APV->NEntries    = Proj->GetEntries();
00558 
00559       if(APV->FitMPV>0){
00560           APV->Gain = APV->FitMPV / MPVmean;
00561           GOOD++;
00562       }else{
00563           APV->Gain = 1;
00564           BAD++;
00565       }
00566       if(APV->Gain<=0)           APV->Gain  = 1;
00567 //      if(!FirstSetOfConstants)   APV->Gain *= APV->PreviousGain;
00568 
00569 
00570       //printf("%5i/%5i:  %6i - %1i  %5E Entries --> MPV = %f +- %f\n",I,APVsColl.size(),APV->DetId, APV->APVId, Proj->GetEntries(), FitResults[0], FitResults[1]);fflush(stdout);
00571       delete Proj;
00572    }printf("\n");
00573 
00574    storeOnTree();
00575 }
00576 
00577 
00578 void SiStripGainFromCalibTree::storeOnTree()
00579 {
00580    unsigned int  tree_Index;
00581    unsigned int  tree_Bin;
00582    unsigned int  tree_DetId;
00583    unsigned char tree_APVId;
00584    unsigned char tree_SubDet;
00585    float         tree_x;
00586    float         tree_y;
00587    float         tree_z;
00588    float         tree_Eta;
00589    float         tree_R;
00590    float         tree_Phi;
00591    float         tree_Thickness;
00592    float         tree_FitMPV;
00593    float         tree_FitMPVErr;
00594    float         tree_FitWidth;
00595    float         tree_FitWidthErr;
00596    float         tree_FitChi2NDF;
00597    double        tree_Gain;
00598    double        tree_PrevGain;
00599    double        tree_NEntries;
00600    bool          tree_isMasked;
00601 
00602    TTree*         MyTree;
00603    MyTree = tfs->make<TTree> ("APVGain","APVGain");
00604    MyTree->Branch("Index"             ,&tree_Index      ,"Index/i");
00605    MyTree->Branch("Bin"               ,&tree_Bin        ,"Bin/i");
00606    MyTree->Branch("DetId"             ,&tree_DetId      ,"DetId/i");
00607    MyTree->Branch("APVId"             ,&tree_APVId      ,"APVId/b");
00608    MyTree->Branch("SubDet"            ,&tree_SubDet     ,"SubDet/b");
00609    MyTree->Branch("x"                 ,&tree_x          ,"x/F"); 
00610    MyTree->Branch("y"                 ,&tree_y          ,"y/F");   
00611    MyTree->Branch("z"                 ,&tree_z          ,"z/F");   
00612    MyTree->Branch("Eta"               ,&tree_Eta        ,"Eta/F");
00613    MyTree->Branch("R"                 ,&tree_R          ,"R/F");
00614    MyTree->Branch("Phi"               ,&tree_Phi        ,"Phi/F");
00615    MyTree->Branch("Thickness"         ,&tree_Thickness  ,"Thickness/F");
00616    MyTree->Branch("FitMPV"            ,&tree_FitMPV     ,"FitMPV/F");
00617    MyTree->Branch("FitMPVErr"         ,&tree_FitMPVErr  ,"FitMPVErr/F");
00618    MyTree->Branch("FitWidth"          ,&tree_FitWidth   ,"FitWidth/F");
00619    MyTree->Branch("FitWidthErr"       ,&tree_FitWidthErr,"FitWidthErr/F");
00620    MyTree->Branch("FitChi2NDF"        ,&tree_FitChi2NDF ,"FitChi2NDF/F");
00621    MyTree->Branch("Gain"              ,&tree_Gain       ,"Gain/D");
00622    MyTree->Branch("PrevGain"          ,&tree_PrevGain   ,"PrevGain/D");
00623    MyTree->Branch("NEntries"          ,&tree_NEntries   ,"NEntries/D");
00624    MyTree->Branch("isMasked"          ,&tree_isMasked   ,"isMasked/O");
00625 
00626 
00627    FILE* Gains = fopen(OutputGains.c_str(),"w");
00628    fprintf(Gains,"NEvents   = %i\n",NEvent);
00629    fprintf(Gains,"NTracks   = %i\n",NTrack);
00630    fprintf(Gains,"NClusters = %i\n",NCluster);
00631    fprintf(Gains,"Number of APVs = %lu\n",static_cast<unsigned long>(APVsColl.size()));
00632    fprintf(Gains,"GoodFits = %i BadFits = %i ratio = %f\n",GOOD,BAD,(100.0*GOOD)/(GOOD+BAD));
00633 
00634    for(unsigned int a=0;a<APVsCollOrdered.size();a++){
00635       stAPVGain* APV = APVsCollOrdered[a];
00636       if(APV==NULL)continue;
00637 //     printf(      "%i | %i | PreviousGain = %7.5f NewGain = %7.5f\n", APV->DetId,APV->APVId,APV->PreviousGain,APV->Gain);
00638       fprintf(Gains,"%i | %i | PreviousGain = %7.5f NewGain = %7.5f\n", APV->DetId,APV->APVId,APV->PreviousGain,APV->Gain);
00639 
00640       tree_Index      = APV->Index;
00641       tree_Bin        = APV->Bin;
00642       tree_DetId      = APV->DetId;
00643       tree_APVId      = APV->APVId;
00644       tree_SubDet     = APV->SubDet;
00645       tree_x          = APV->x;
00646       tree_y          = APV->y;
00647       tree_z          = APV->z;
00648       tree_Eta        = APV->Eta;
00649       tree_R          = APV->R;
00650       tree_Phi        = APV->Phi;
00651       tree_Thickness  = APV->Thickness;
00652       tree_FitMPV     = APV->FitMPV;
00653       tree_FitMPVErr  = APV->FitMPVErr;
00654       tree_FitWidth   = APV->FitWidth;
00655       tree_FitWidthErr= APV->FitWidthErr;
00656       tree_FitChi2NDF = APV->FitChi2;
00657       tree_Gain       = APV->Gain;
00658       tree_PrevGain   = APV->PreviousGain;
00659       tree_NEntries   = APV->NEntries;
00660       tree_isMasked   = APV->isMasked;
00661 
00662       MyTree->Fill();
00663    }
00664    fclose(Gains);
00665 
00666 
00667 }
00668 
00669 
00670 
00671 SiStripApvGain* SiStripGainFromCalibTree::getNewObject() 
00672 {
00673    SiStripApvGain* obj = new SiStripApvGain();
00674    std::vector<float>* theSiStripVector = NULL;
00675    unsigned int PreviousDetId = 0; 
00676    for(unsigned int a=0;a<APVsCollOrdered.size();a++)
00677    {
00678       stAPVGain* APV = APVsCollOrdered[a];
00679       if(APV==NULL){ printf("Bug\n"); continue; }
00680       if(APV->DetId != PreviousDetId){
00681          if(theSiStripVector!=NULL){
00682             SiStripApvGain::Range range(theSiStripVector->begin(),theSiStripVector->end());
00683             if ( !obj->put(PreviousDetId,range) )  printf("Bug to put detId = %i\n",PreviousDetId);
00684          }
00685          theSiStripVector = new std::vector<float>;
00686          PreviousDetId = APV->DetId;
00687       }
00688       theSiStripVector->push_back(APV->Gain);
00689    }
00690    if(theSiStripVector!=NULL){
00691       SiStripApvGain::Range range(theSiStripVector->begin(),theSiStripVector->end());
00692       if ( !obj->put(PreviousDetId,range) )  printf("Bug to put detId = %i\n",PreviousDetId);
00693    }
00694 
00695    return obj;
00696 }
00697 
00698 
00699 SiStripGainFromCalibTree::~SiStripGainFromCalibTree()
00700 { 
00701 }
00702 
00703 void SiStripGainFromCalibTree::MakeCalibrationMap(){
00704    if(!useCalibration)return;
00705 
00706   
00707    TChain* t1 = new TChain("SiStripCalib/APVGain");
00708    t1->Add(m_calibrationPath.c_str());
00709 
00710    unsigned int  tree_DetId;
00711    unsigned char tree_APVId;
00712    double        tree_Gain;
00713 
00714    t1->SetBranchAddress("DetId"             ,&tree_DetId      );
00715    t1->SetBranchAddress("APVId"             ,&tree_APVId      );
00716    t1->SetBranchAddress("Gain"              ,&tree_Gain       );
00717 
00718    for (unsigned int ientry = 0; ientry < t1->GetEntries(); ientry++) {
00719       t1->GetEntry(ientry);
00720        stAPVGain* APV = APVsColl[(tree_DetId<<3) | (unsigned int)tree_APVId];
00721        APV->CalibGain = tree_Gain;
00722    }
00723 
00724 }
00725 
00726 
00727 
00728 
00729 DEFINE_FWK_MODULE(SiStripGainFromCalibTree);
CMSSW_4_4_3_patch1/src/CalibTracker/SiStripChannelGain/plugins/SiStripGainFromCalibTree.cc
