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