CMSSW Reference Manual

00001 /*
00002  * $Date: 2008/04/19 19:10:51 $
00003  * $Revision: 1.15 $
00004  *
00005  * \author: D. Giordano, domenico.giordano@cern.ch
00006  * Modified: M.De Mattia 2/3/2007 & R.Castello 5/4/2007 & Susy Borgia 15/11/07
00007  */
00008 
00009 #include "AnalysisExamples/SiStripDetectorPerformance/plugins/ClusterAnalysis.h"
00010 #include "DataFormats/TrackerRecHit2D/interface/ProjectedSiStripRecHit2D.h"
00011 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
00012 
00013 
00014 #include "AnalysisDataFormats/TrackInfo/src/TrackInfo.cc"
00015 
00016 #include "sstream"
00017 
00018 #include "TH3S.h"
00019 #include "TProfile.h"
00020 #include "TCanvas.h"
00021 #include "TPostScript.h"
00022 
00023 static const uint8_t _NUM_SISTRIP_SUBDET_ = 4;
00024 static TString SubDet[_NUM_SISTRIP_SUBDET_]={"_TIB","_TOB","_TID","_TEC"};
00025 static TString flags[3] = {"_onTrack","_offTrack","_All"};
00026 static TString width_flags[5] = {"","_width_1","_width_2","_width_3","_width_ge_4"};
00027 
00028 
00029 namespace cms{
00030   ClusterAnalysis::ClusterAnalysis(edm::ParameterSet const& conf): 
00031     conf_(conf),
00032     psfilename_(conf.getParameter<std::string>("psfileName")), 
00033     psfiletype_(conf.getParameter<int32_t>("psfiletype")),
00034     psfilemode_(conf.getUntrackedParameter<int32_t>("psfilemode",1)),
00035     Filter_src_( conf.getParameter<edm::InputTag>( "Filter_src" ) ),
00036     Track_src_( conf.getParameter<edm::InputTag>( "Track_src" ) ),
00037     Cluster_src_( conf.getParameter<edm::InputTag>( "Cluster_src" ) ),
00038     ModulesToBeExcluded_(conf.getParameter< std::vector<uint32_t> >("ModulesToBeExcluded")),
00039     EtaAlgo_(conf.getParameter<int32_t>("EtaAlgo")),
00040     NeighStrips_(conf.getParameter<int32_t>("NeighStrips")),
00041     not_the_first_event(false),
00042     tracksCollection_in_EventTree(true),
00043     trackAssociatorCollection_in_EventTree(true)
00044   {
00045   }
00046 
00047   ClusterAnalysis::~ClusterAnalysis(){
00048     std::cout << "Destructing object" << std::endl;
00049     //    delete Hlist;
00050   }
00051   
00052   void ClusterAnalysis::beginRun(const edm::Run& run, const edm::EventSetup& es ) {
00053 
00054     //get geom    
00055     es.get<TrackerDigiGeometryRecord>().get( tkgeom );
00056     edm::LogInfo("ClusterAnalysis") << "[ClusterAnalysis::beginJob] There are "<<tkgeom->detUnits().size() <<" detectors instantiated in the geometry" << std::endl;  
00057 
00058     es.get<SiStripDetCablingRcd>().get( SiStripDetCabling_ );
00059 
00060     es.get<SiStripQualityRcd>().get(SiStripQuality_);
00061 
00062     book();
00063   }
00064 
00065   void ClusterAnalysis::book() {
00066     
00067     TFileDirectory ClusterNoise = fFile->mkdir( "ClusterNoise" );
00068     TFileDirectory ClusterSignal = fFile->mkdir("ClusterSignal");
00069     TFileDirectory ClusterStoN = fFile->mkdir("ClusterStoN");
00070     TFileDirectory ClusterEta = fFile->mkdir("ClusterEta");
00071     TFileDirectory ClusterWidth = fFile->mkdir("ClusterWidth");
00072     TFileDirectory ClusterPos = fFile->mkdir("ClusterPos");
00073     TFileDirectory Tracks = fFile->mkdir("Tracks");
00074     TFileDirectory Layer = fFile->mkdir("Layer");
00075 
00076     const edm::ParameterSet mapSet = conf_.getParameter<edm::ParameterSet>("MapFlag");
00077     if( mapSet.getParameter<bool>("Map_ClusOccOn") ){
00078       tkMap_ClusOcc[0]=new TrackerMap( "ClusterOccupancy_onTrack" );
00079       tkMap_ClusOcc[1]=new TrackerMap( "ClusterOccupancy_offTrack" );
00080       tkMap_ClusOcc[2]=new TrackerMap( "ClusterOccupancy_All" );
00081     }  
00082     if( mapSet.getParameter<bool>("Map_InvHit") ) 
00083       tkInvHit=new TrackerMap("Invalid_Hit");
00084 
00085     //&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
00086     // get list of active detectors from SiStripDetCabling 
00087 
00088     std::vector<uint32_t> vdetId_;
00089     SiStripDetCabling_->addActiveDetectorsRawIds(vdetId_);
00090     //&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
00091 
00092     //Create histograms
00093     Hlist = new THashList();
00094 
00095     //Display 3D
00096     name = "ClusterGlobalPos";
00097     bookHlist("TH3","TH3ClusterGlobalPos",ClusterPos, name, "z (cm)","x (cm)","y (cm)"); 
00098 
00099     //&&&&&&&&&&&&&&&&&&&&&&&&
00100 
00101     // bookHlist(TObjArray, Name of the parameterset in the cfg, name of the TFileDirectory, name of the hitsogram, xbinnumber, xmin, xmax )
00102 
00103     name = "nTracks";
00104     bookHlist("TH1","TH1nTracks", Tracks, name, "N Tracks" );
00105     name = "nRecHits";
00106     bookHlist("TH1","TH1nRecHits",Tracks,  name, "N RecHits" );
00107 
00108     // Loop on onTrack, offTrack, All
00109     for (int j=0;j<3;j++){      
00110       //Number of Cluster 
00111       name="nClusters"+flags[j];
00112       bookHlist("TH1","TH1nClusters",Tracks, name, "N Clusters" );
00113 
00114       for (int i=0;i<_NUM_SISTRIP_SUBDET_;i++) {
00115         //Number of Cluster on each det
00116         name="nClusters"+SubDet[i]+flags[j];
00117         bookHlist("TH1","TH1nClusters",Tracks, name, "N Clusters" );
00118       }
00119     }
00120 
00121     // Loop on onTrack, offTrack, All
00122     for (int j=0;j<3;j++){
00123       //Histos for detector type
00124       for (int i=0;i<_NUM_SISTRIP_SUBDET_;i++){
00125         
00126         TString appString=SubDet[i]+flags[j];
00127 
00128         //Cluster Width
00129         name="cWidth"+appString;
00130         bookHlist("TH1","TH1ClusterWidth",ClusterWidth, name, "Nstrip" );
00131 
00132         //Loop for cluster width
00133         for (int iw=0;iw<5;iw++){
00134           
00135           appString=SubDet[i]+flags[j]+width_flags[iw];
00136         
00137           //Cluster Noise
00138           name="cNoise"+appString;
00139           bookHlist("TH1","TH1ClusterNoise",ClusterNoise, name, "ADC count" );
00140 
00141           //Cluster Signal
00142           name="cSignal"+appString;
00143           bookHlist("TH1","TH1ClusterSignal",ClusterSignal, name, "ADC count" );                 
00144 
00145           //Cluster Signal corrected
00146           if(j==0 && iw==0 ){
00147             name="cSignalCorr"+appString;
00148             bookHlist("TH1","TH1ClusterSignalCorr",ClusterSignal, name, "ADC count" );  
00149           }
00150           
00151           //Cluster StoN
00152           name="cStoN"+appString;
00153           bookHlist("TH1","TH1ClusterStoN",ClusterStoN, name );
00154 
00155           //Cluster SignaltoNoise corrected
00156           if(j==0 && iw==0 ){        
00157             name="cStoNCorr"+appString;
00158             bookHlist("TH1","TH1ClusterStoNCorr",ClusterStoN, name );  
00159           }
00160 
00161           //Cluster Position
00162           name="cPos"+appString;
00163           bookHlist("TH1","TH1ClusterPos",ClusterPos, name, "strip Num" );
00164 
00165           //Cluster StoN Vs Cluster Position
00166           name="cStoNVsPos"+appString;
00167           bookHlist("TH2","TH2ClusterStoNVsPos",ClusterPos, name, "strip Num");
00168 
00169           //Cluster Charge Division (only for study on Raw Data Runs)
00170           name="cEta"+appString;
00171           bookHlist("TH1","TH1ClusterEta",ClusterEta, name, "" );
00172 
00173           name="cEta_scatter"+appString;
00174           bookHlist("TH2","TH2ClusterEta",ClusterEta, name, "" , "");
00175         }//end loop on width 
00176 
00177         //cWidth Vs Angle
00178         name = "ClusterWidthVsAngle"+appString;
00179         bookHlist("TProfile","TProfileWidthAngle",ClusterWidth, name, "cos(angle_xz)", "clusWidth");
00180 
00181         //Residual Vs Angle
00182         name = "ResidualVsAngle"+appString;
00183         bookHlist("TProfile","TProfileResidualAngle",ClusterWidth, name, "Angle" , "Residual");
00184       
00185       } //end loop on det type 
00186     }//end loop on onTrack,offTrack,all
00187 
00188 
00189     //&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
00190     //Detector Detail Plots
00191 
00192     //Histos for each detector
00193     for (std::vector<uint32_t>::const_iterator detid_iter=vdetId_.begin();detid_iter!=vdetId_.end();detid_iter++){  
00194       uint32_t detid = *detid_iter;
00195       
00196       if (detid < 1){
00197         edm::LogError("ClusterAnalysis")<< "[" <<__PRETTY_FUNCTION__ << "] invalid detid " << detid<< std::endl;
00198         continue;
00199       }
00200       const StripGeomDetUnit* _StripGeomDetUnit = dynamic_cast<const StripGeomDetUnit*>(tkgeom->idToDetUnit(DetId(detid)));
00201       if (_StripGeomDetUnit==0){
00202         edm::LogError("SiStripCondObjDisplay")<< "[SiStripCondObjDisplay::beginJob] the detID " << detid << " doesn't seem to belong to Tracker" << std::endl; 
00203         continue;
00204       }
00205       
00206       
00207       //&&&&&&&&&&&&&&&&&
00208       // Insert here code to instantiate histos per detector
00209       
00210       unsigned int nstrips = _StripGeomDetUnit->specificTopology().nstrips();
00211       
00212       //      edm::LogError("ClusterAnalysis") << " Detid " << detid << " SubDet " << GetSubDetAndLayer(detid).first << " Layer " << GetSubDetAndLayer(detid).second << std::endl;   
00213       if (DetectedLayers.find(GetSubDetAndLayer(detid)) == DetectedLayers.end()){
00214 
00215         DetectedLayers[GetSubDetAndLayer(detid)]=true;
00216       }
00217  
00218       //&&&&&&&&&&&&&&
00219       // Retrieve information for the module
00220       //&&&&&&&&&&&&&&&&&&     
00221       char cdetid[128];
00222       sprintf(cdetid,"%d",detid);
00223       char aname[128];
00224       sprintf(aname,"%s_%d",_StripGeomDetUnit->type().name().c_str(),detid);
00225       char SubStr[128];
00226 
00227       SiStripDetId a(detid);
00228       if ( a.subdetId() == 3 ){
00229         TIBDetId b(detid);
00230         sprintf(SubStr,"_SingleDet_%d_TIB_%d_%d_%d_%d",detid,b.layer(),b.string()[0],b.string()[1],b.glued());
00231       } else if ( a.subdetId() == 4 ) {
00232         TIDDetId b(detid);
00233         sprintf(SubStr,"_SingleDet_%d_TID_%d_%d_%d_%d",detid,b.wheel(),b.ring(),b.side(),b.glued());
00234       } else if ( a.subdetId() == 5 ) {
00235         TOBDetId b(detid);
00236         sprintf(SubStr,"_SingleDet_%d_TOB_%d_%d_%d_%d",detid,b.layer(),b.rod()[0],b.rod()[1],b.glued());
00237       } else if ( a.subdetId() == 6 ) {
00238         TECDetId b(detid);
00239         sprintf(SubStr,"_SingleDet_%d_TEC_%d_%d_%d_%d_%d",detid,b.wheel(),b.ring(),b.side(),b.glued(),b.stereo());
00240       }
00241       
00242       TString appString=TString(SubStr);//+"_"+cdetid;
00243 
00244       TFileDirectory detid_dir = fFile->mkdir( cdetid );     
00245 
00246       //Cluster Noise
00247       name="cNoise"+appString;
00248       bookHlist("TH1","TH1ClusterNoise",detid_dir, name, "ADC count" );
00249 
00250       //Cluster Signal
00251       name="cSignal"+appString;
00252       bookHlist("TH1","TH1ClusterSignal",detid_dir, name, "ADC count" );
00253 
00254       //Cluster StoN
00255       name="cStoN"+appString;
00256       bookHlist("TH1","TH1ClusterStoN",detid_dir, name, "" );
00257 
00258       //Cluster Signal x Fiber
00259       name="cSignalxFiber"+appString+"_onTrack";
00260       bookHlist("TProfile","TProfileSignalxFiber",detid_dir, name, "ApvPair", "ADC count" );
00261 
00262       //Cluster Width
00263       name="cWidth"+appString;
00264       bookHlist("TH1","TH1ClusterWidth",detid_dir, name, "Nstrip" );
00265 
00266       //Cluster Position
00267       name="cPos"+appString;
00268       //bookHlist("TH1","TH1ClusterPos", name, "Nbinx", "xmin", "xmax" );
00269       Hlist->Add(new TH1F(name,name,nstrips,0,nstrips));
00270                 
00271       //Cluster StoN Vs Cluster Position
00272       name="cStoNVsPos"+appString;
00273       char labeln[128];
00274       sprintf(labeln,"strip Num (Ntot=%d)",nstrips);
00275       bookHlist("TH2","TH2ClusterStoNVsPos",detid_dir, name, labeln);
00276  
00277       //Cluster Charge Division (only for study on Raw Data Runs)
00278       name="cEta"+appString;
00279       bookHlist("TH1","TH1ClusterEta", detid_dir,name, "" );
00280 
00281       name="cEta_scatter"+appString;
00282       bookHlist("TH2","TH2ClusterEta", detid_dir,name, "" ,  "" );
00283 
00284       //      const StripGeomDetUnit* _StripGeomDetUnit = dynamic_cast<const StripGeomDetUnit*>(tkgeom->idToDetUnit(DetId(detid)));
00285       //      if (_StripGeomDetUnit==0){
00286       //        continue;
00287       //      }
00288           
00289       //      unsigned int nstrips = _StripGeomDetUnit->specificTopology().nstrips();
00290 
00291       name="DBPedestals"+appString;
00292       TH1F* pPed = detid_dir.make<TH1F>(name,name,nstrips,-0.5,nstrips-0.5);
00293       Hlist->Add(pPed);
00294     
00295       name="DBNoise"+appString;
00296       TH1F* pNoi = detid_dir.make<TH1F>(name,name,nstrips,-0.5,nstrips-0.5);
00297       Hlist->Add(pNoi);
00298     
00299       name="DBBadStrips"+appString;
00300       TH1F* pBad = detid_dir.make<TH1F>(name,name,nstrips,-0.5,nstrips-0.5);  
00301       Hlist->Add(pBad);
00302 
00303       //&&&&&&&&&&&&&&&&&
00304     }//end loop on detector     
00305 
00306 
00307     //&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
00308     //Layer Detail Plots
00309     
00310     for (std::map<std::pair<std::string,uint32_t>,bool>::const_iterator iter=DetectedLayers.begin(); iter!=DetectedLayers.end();iter++){
00311      
00312       char cApp[64];
00313       sprintf(cApp,"_Layer_%d",iter->first.second);
00314       TFileDirectory Layer = fFile->mkdir("Layer");
00315       // Loop on onTrack, offTrack, All
00316       for (int j=0;j<3;j++){
00317         
00318         TString appString="_"+TString(iter->first.first)+cApp+flags[j];
00319      
00320         //Cluster Noise
00321         name="cNoise"+appString;
00322         bookHlist("TH1","TH1ClusterNoise", Layer,name, "ADC count" );
00323 
00324         //Cluster Signal
00325         name="cSignal"+appString;
00326         bookHlist("TH1","TH1ClusterSignal",Layer, name, "ADC count" );
00327         
00328         //Cluster Signal corrected
00329         if(j==0){
00330           name="cSignalCorr"+appString;
00331           bookHlist("TH1","TH1ClusterSignalCorr",Layer, name, "ADC count" );
00332 
00333           //Cluster Signal Vs Angle
00334           name = "cSignalVsAngle"+appString;
00335           bookHlist("TProfile","TProfilecSignalVsAngle",Layer, name, "cos(angle_rz)" , "ADC count");
00336           name = "cSignalVsAngleH"+appString;
00337           bookHlist("TH2","TH2cSignalVsAngle",Layer, name, "cos(angle_rz)" , "ADC count");
00338         }
00339 
00340         //Cluster StoN
00341         name="cStoN"+appString;
00342         bookHlist("TH1","TH1ClusterStoN",Layer, name, "" );
00343         
00344         //Cluster SignaltoNoise corrected
00345         if(j==0){
00346           name="cStoNCorr"+appString;
00347           bookHlist("TH1","TH1ClusterStoNCorr",Layer, name, "" );
00348         }
00349 
00350         //Cluster Width
00351         name="cWidth"+appString;
00352         bookHlist("TH1","TH1ClusterWidth",Layer, name, "Nstrip" );
00353 
00354         //Cluster Position
00355         name="cPos"+appString;
00356         bookHlist("TH1","TH1ClusterPos",Layer, name, "strip Num" );
00357 
00358         //Cluster StoN Vs Cluster Position
00359         name="cStoNVsPos"+appString;
00360         bookHlist("TH2","TH2ClusterStoNVsPos",Layer, name, "strip Num");
00361 
00362         //residual
00363         name="res_x"+appString;
00364         bookHlist("TH1","TH1Residual_x",Layer, name, "" );
00365         
00366         //residual y
00367         name="res_y"+appString;
00368         bookHlist("TH1","TH1Residual_y",Layer, name, "" );
00369       
00370         //cWidth Vs Angle
00371         if(j==0){
00372           name = "ClusterWidthVsAngle"+appString;
00373           bookHlist("TProfile","TProfileWidthAngle",Layer, name, "angle_xz" , "clusWidth");
00374         
00375           //Residuals Vs Angle
00376           name = "ResidualVsAngle"+appString;
00377           bookHlist("TProfile","TProfileResidualAngle",Layer, name, "cos(angle_rz)" , "Residual");
00378 
00379           //Angle Vs phi
00380           name = "AngleVsPhi"+appString;
00381           bookHlist("TProfile","TProfileAngleVsPhi",Layer, name, "Phi (deg)" , "Impact angle angle_xz (deg)");
00382         }
00383       }
00384     }
00385   }
00386 
00387   //------------------------------------------------------------------------------------------
00388 
00389   void ClusterAnalysis::endJob() {  
00390     edm::LogInfo("ClusterAnalysis") << "[ClusterAnalysis::endJob] >>> saving histograms" << std::endl;
00391     
00392     if (psfilemode_>0){    
00393       
00394       edm::LogInfo("ClusterAnalysis")  << "... And now write on ps file " << psfiletype_ << std::endl;
00395       TPostScript ps(psfilename_.c_str(),psfiletype_);
00396       TCanvas Canvas("c","c");//("c","c",600,300);
00397       TIter lIter(Hlist);
00398       TObject *th = (TObject*)lIter.Next(); 
00399       while (th!=NULL){
00400         if (psfilemode_>1 || (strstr(th->GetName(),"SingleDet_")==NULL && strstr(th->GetName(),"_width_")==NULL)){  
00401           edm::LogInfo("ClusterAnalysis") << "Histos name " << th->GetName() << " title " <<  th->GetTitle() << std::endl;
00402           if (dynamic_cast<TH1F*>(th) !=NULL){
00403             if (dynamic_cast<TH1F*>(th)->GetEntries() != 0)
00404               th->Draw();
00405           }
00406           if (dynamic_cast<TH2F*>(th) !=NULL){
00407             if (dynamic_cast<TH2F*>(th)->GetEntries() != 0)
00408               th->Draw();
00409           }
00410           else if (dynamic_cast<TProfile*>(th) !=NULL){
00411             if (dynamic_cast<TProfile*>(th)->GetEntries() != 0)
00412               th->Draw();
00413           }   
00414           else if (dynamic_cast<TH3S*>(th) !=NULL){
00415             if (dynamic_cast<TH3S*>(th)->GetEntries() != 0)
00416               th->Draw();
00417           }
00418           Canvas.Update();
00419           ps.NewPage();
00420         }
00421         
00422       th = (TObject*)lIter.Next(); 
00423       }      
00424       ps.Close();      
00425     }    
00426 
00427     //TkMap->Save() and Print()
00428 
00429     const edm::ParameterSet mapSett = conf_.getParameter<edm::ParameterSet>("MapFlag");
00430     if( mapSett.getParameter<bool>("Map_ClusOccOn") ){
00431       tkMap_ClusOcc[0]->save(true,0,0,"ClusterOccupancyMap_onTrack.png");
00432       tkMap_ClusOcc[1]->save(true,0,0,"ClusterOccupancyMap_offTrack.png");
00433       tkMap_ClusOcc[2]->save(true,0,0,"ClusterOccupancyMap_All.png");
00434       
00435       tkMap_ClusOcc[0]->print(true,0,0,"ClusterOccupancyMap_onTrack");   
00436       tkMap_ClusOcc[1]->print(true,0,0,"ClusterOccupancyMap_offTrack");   
00437       tkMap_ClusOcc[2]->print(true,0,0,"ClusterOccupancyMap_All");
00438     }
00439 
00440     if( mapSett.getParameter<bool>("Map_InvHit")) {
00441       tkInvHit->save(true,0,0,"Inv_Hit_map.png");
00442       tkInvHit->print(true,0,0,"Inv_Hit_map");
00443     }
00444     
00445   }  
00446   //------------------------------------------------------------------------------------------
00447 
00448   void ClusterAnalysis::analyze(const edm::Event& e, const edm::EventSetup& es) {
00449 
00450     tracksCollection_in_EventTree=true;
00451     trackAssociatorCollection_in_EventTree=true;
00452     
00453     LogTrace("ClusterAnalysis") << "[ClusterAnalysis::analyse]  " << "Run " << e.id().run() << " Event " << e.id().event() << std::endl;
00454     runNb   = e.id().run();
00455     eventNb = e.id().event();
00456     edm::LogInfo("ClusterAnalysis") << "Processing run " << runNb << " event " << eventNb << std::endl;
00457 
00458     es.get<SiStripDetCablingRcd>().get( SiStripDetCabling_ );
00459 
00460     //get Pedestal and Noise  ES handle
00461     es.get<SiStripPedestalsRcd>().get(pedestalHandle);
00462     es.get<SiStripNoisesRcd>().get(noiseHandle);
00463 
00464     if (!not_the_first_event){
00465       fillPedNoiseFromDB();
00466       not_the_first_event=true;
00467     }
00468     
00469     //Get input
00470     e.getByLabel( Cluster_src_, dsv_SiStripCluster);    
00471     
00472     e.getByLabel(Track_src_, trackCollection);
00473     if(trackCollection.isValid()){
00474     }else{
00475       edm::LogError("ClusterAnalysis")<<"trackCollection not found "<<std::endl;
00476       tracksCollection_in_EventTree=false;
00477     }
00478     
00479     edm::InputTag TkiTag = conf_.getParameter<edm::InputTag>( "TrackInfo" );  
00480     
00481     e.getByLabel(TkiTag,tkiTkAssCollection); 
00482     if(tkiTkAssCollection.isValid()){
00483     }else{
00484       edm::LogError("ClusterAnalysis")<<"trackInfo not found "<<std::endl;
00485       trackAssociatorCollection_in_EventTree=false;
00486     }
00487     vPSiStripCluster.clear();
00488     countOn=0;
00489     countOff=0;
00490     countAll=0;
00491     // istart=oXZHitAngle.size();  
00492   
00493     // get geometry to evaluate local angles
00494     edm::ESHandle<TrackerGeometry> estracker;
00495     es.get<TrackerDigiGeometryRecord>().get(estracker);
00496     _tracker=&(* estracker);
00497 
00498     //Perform track study
00499     if (tracksCollection_in_EventTree || trackAssociatorCollection_in_EventTree) {
00500       LogTrace("ClusterAnalysis") << "\n["<<__PRETTY_FUNCTION__<<"] \n Executing trackStudy for Event = " << eventNb << std::endl; 
00501       trackStudy(es);
00502     }
00503    
00504     std::stringstream ss;
00505     ss << "\nList of SiStripClusterPointer\n";
00506     for (std::vector<const SiStripCluster*>::iterator iter=vPSiStripCluster.begin();iter!=vPSiStripCluster.end();iter++)
00507       ss << *iter << "\n";    
00508     LogTrace("ClusterAnalysis") << "\n["<<__PRETTY_FUNCTION__<<"] \n vPSiStripCluster.size()=" << vPSiStripCluster.size()<< ss.str() << std::endl;      
00509     
00510 
00511     //Perform Cluster Study (irrespectively to tracks)
00512     AllClusters(es);
00513 
00514     if (countAll != countOn+countOff)
00515       edm::LogWarning("ClusterAnalysis") << "\n["<<__PRETTY_FUNCTION__<<"] Counts (on, off, all) do not match" << countOn << " " << countOff << " " << countAll; 
00516 
00517     for (int j=0;j<3;j++){
00518       int nTot=0;
00519       for (int i=0;i<4;i++){
00520         ((TH1F*) Hlist->FindObject("nClusters"+SubDet[i]+flags[j]))
00521           ->Fill(NClus[i][j]);
00522         nTot+=NClus[i][j];
00523         NClus[i][j]=0;
00524       }
00525       ((TH1F*) Hlist->FindObject("nClusters"+flags[j]))->Fill(nTot);
00526     }
00527   }
00528   
00529   //------------------------------------------------------------------------
00530   
00531   void ClusterAnalysis::trackStudy( const edm::EventSetup& es){
00532     LogTrace("ClusterAnalysis") << "\n["<<__PRETTY_FUNCTION__<<"]" << std::endl;
00533     
00534     const reco::TrackCollection tC = *(trackCollection.product());
00535     
00536     int nTracks=tC.size();
00537     
00538     edm::LogInfo("ClusterAnalysis") << "Reconstructed "<< nTracks << " tracks" << std::endl ;
00539     
00540     int i=0;
00541     
00542     edm::LogInfo("ClusterAnalysis") <<"\t\t TrackCollection size "<< trackCollection->size() << std::endl;
00543     if(trackCollection->size()==0){
00544       ((TH1F*) Hlist->FindObject("nTracks"))->Fill(trackCollection->size());
00545       edm::LogInfo("ClusterAnalysis") <<"\t\t TrackCollection size zero!! Histo filled with a "<< trackCollection->size() << std::endl;
00546     }
00547     
00548     for (unsigned int track=0;track<trackCollection->size();++track){
00549       reco::TrackRef trackref = reco::TrackRef(trackCollection, track);
00550       
00551       //get the ref to the trackinfo
00552       reco::TrackInfoRef trackinforef=(*tkiTkAssCollection.product())[trackref];
00553       // loop on all the track hits
00554       edm:: LogInfo("ClusterAnalysis")
00555         << "Track number "<< i+1 
00556         << "\n\tmomentum: " << trackref->momentum()
00557         << "\n\tPT: " << trackref->pt()
00558         << "\n\tvertex: " << trackref->vertex()
00559         << "\n\timpact parameter: " << trackref->d0()
00560         << "\n\tcharge: " << trackref->charge()
00561         << "\n\tnormalizedChi2: " << trackref->normalizedChi2() 
00562         <<"\n\tFrom EXTRA : "
00563         <<"\n\t\touter PT "<< trackref->outerPt()<<std::endl;
00564       i++;
00565       
00566       int recHitsSize=trackref->recHitsSize();
00567       edm::LogInfo("ClusterAnalysis") <<"\t\tNumber of RecHits "<<recHitsSize<<std::endl;
00568       
00569       const edm::ParameterSet psett = conf_.getParameter<edm::ParameterSet>("TrackThresholds");
00570       if( trackref->normalizedChi2() < psett.getParameter<double>("maxChi2") && recHitsSize > psett.getParameter<double>("minRecHit") ){
00571         ((TH1F*) Hlist->FindObject("nRecHits"))->Fill(recHitsSize);
00572         ((TH1F*) Hlist->FindObject("nTracks"))->Fill(nTracks);
00573       }else{
00574         LogTrace("ClusterAnalysis") <<"\t\tSkipped track "<< i << std::endl;
00575         //continue;
00576       }
00577       
00578       //------------------------------RESIDUAL at the layer level ------------
00579       /*          
00580                  LocalPoint stateposition= 
00581                  LocalPoint rechitposition= 
00582                  fillTH1( stateposition.x() - rechitposition.x(),"res_x"+appString,0);
00583                  fillTH1( stateposition.y() - rechitposition.y(),"res_y"+appString,0);
00584                  ((TProfile*) Hlist->FindObject("ResidualVsAngle"))->Fill(angle,stateposition.x()- rechitposition.x(),1);
00585                  ((TProfile*) Hlist->FindObject("ResidualVsAngle"+appString+"_onTrack"))->Fill(angle,stateposition.x()- rechitposition.x(),1);
00586                  
00587       */
00588       
00589       TrackingRecHitRef rechitref=trackref->recHit(2);
00590       
00591       const uint32_t& ndetid = rechitref->geographicalId().rawId();
00592       if (find(ModulesToBeExcluded_.begin(),ModulesToBeExcluded_.end(),ndetid)!=ModulesToBeExcluded_.end()){
00593         LogTrace("ClusterAnalysis") << "Modules Excluded" << std::endl;
00594         return;
00595       } 
00596       
00597       const edm::ParameterSet mappSet = conf_.getParameter<edm::ParameterSet>("MapFlag");
00598       if( mappSet.getParameter<bool>("Map_InvHit") ){
00599         edm::LogInfo("TrackInfoAnalyzerExample") << "RecHit type: " << rechitref->getType() << std::endl;  
00600         if(!rechitref->isValid()){
00601           edm::LogInfo("TrackInfoAnalyzerExample") <<"The recHit is not valid"<< std::endl;
00602           //inactive??      
00603           if(rechitref->getType() == 2) {
00604             edm::LogInfo("TrackInfoAnalyzerExample") << "inactive rechit found on detid " << ndetid << std::endl;
00605             tkInvHit->fill(ndetid,1); 
00606             tkInvHit->fill(ndetid+1,1);
00607             tkInvHit->fill(ndetid+2,1);
00608           }  
00609         }else{
00610           edm::LogInfo("TrackInfoAnalyzerExample") <<"The recHit is valid"<< std::endl;
00611         }     
00612       }  
00613       reco::TrackInfo::TrajectoryInfo::const_iterator iter;
00614 
00615       for(iter=trackinforef->trajStateMap().begin();iter!=trackinforef->trajStateMap().end();iter++){
00616         
00617         //trajectory local direction and position on detector
00618         LocalVector statedirection=(trackinforef->stateOnDet(Updated,(*iter).first)->parameters()).momentum();
00619         LocalPoint  stateposition=(trackinforef->stateOnDet(Updated,(*iter).first)->parameters()).position();
00620         
00621         std::stringstream ss;
00622         ss <<"LocalMomentum: "<<statedirection
00623            <<"\nLocalPosition: "<<stateposition
00624            << "\nLocal x-z plane angle: "<<atan2(statedirection.x(),statedirection.z());
00625 
00626         if(trackinforef->type((*iter).first)==Matched){ // get the direction for the components
00627           
00628           const SiStripMatchedRecHit2D* matchedhit=dynamic_cast<const SiStripMatchedRecHit2D*>(&(*(iter)->first));
00629           if (matchedhit!=0){
00630             ss<<"\nMatched recHit found"<< std::endl;     
00631             //mono side
00632             statedirection= trackinforef->localTrackMomentumOnMono(Updated,(*iter).first);
00633             if(statedirection.mag() != 0)         RecHitInfo(es, matchedhit->monoHit(),statedirection,trackref);
00634             //stereo side
00635             statedirection= trackinforef->localTrackMomentumOnStereo(Updated,(*iter).first);
00636             if(statedirection.mag() != 0)         RecHitInfo(es, matchedhit->stereoHit(),statedirection,trackref);
00637             ss<<"\nLocalMomentum (stereo): "<<trackinforef->localTrackMomentumOnStereo(Updated,(*iter).first);
00638           }
00639         }
00640         else if (trackinforef->type((*iter).first)==Projected){//one should be 0
00641           ss<<"\nProjected recHit found"<< std::endl;
00642           const ProjectedSiStripRecHit2D* phit=dynamic_cast<const ProjectedSiStripRecHit2D*>(&(*(iter)->first));
00643           if(phit!=0){
00644             //mono side
00645             statedirection= trackinforef->localTrackMomentumOnMono(Updated,(*iter).first);
00646             if(statedirection.mag() != 0) RecHitInfo(es, &(phit->originalHit()),statedirection,trackref);
00647             //stereo side
00648             statedirection= trackinforef->localTrackMomentumOnStereo(Updated,(*iter).first);
00649             if(statedirection.mag() != 0)  RecHitInfo(es, &(phit->originalHit()),statedirection,trackref);
00650           }     
00651           
00652         }
00653         else {
00654           const SiStripRecHit2D* hit=dynamic_cast<const SiStripRecHit2D*>(&(*(iter)->first));
00655           if(hit!=0){
00656             ss<<"\nSingle recHit found"<< std::endl;      
00657             statedirection=(trackinforef->stateOnDet(Updated,(*iter).first)->parameters()).momentum();
00658             if(statedirection.mag() != 0) RecHitInfo(es, hit,statedirection,trackref);
00659 
00660           }
00661         }
00662         LogTrace("TrackInfoAnalyzerExample") <<ss.str() << std::endl;
00663       }
00664     }
00665   }
00666 
00667   void ClusterAnalysis::RecHitInfo( const edm::EventSetup& es, const SiStripRecHit2D* tkrecHit, LocalVector LV,reco::TrackRef track_ref ){
00668     
00669     if(!tkrecHit->isValid()){
00670       LogTrace("ClusterAnalysis") <<"\t\t Invalid Hit " << std::endl;
00671       return;  
00672     }
00673     
00674     const uint32_t& detid = tkrecHit->geographicalId().rawId();
00675     if (find(ModulesToBeExcluded_.begin(),ModulesToBeExcluded_.end(),detid)!=ModulesToBeExcluded_.end()){
00676       LogTrace("ClusterAnalysis") << "Modules Excluded" << std::endl;
00677       return;
00678     }
00679     
00680     LogTrace("ClusterAnalysis")
00681       <<"\n\t\tRecHit on det "<<tkrecHit->geographicalId().rawId()
00682       <<"\n\t\tRecHit in LP "<<tkrecHit->localPosition()
00683       <<"\n\t\tRecHit in GP "<<tkgeom->idToDet(tkrecHit->geographicalId())->surface().toGlobal(tkrecHit->localPosition()) 
00684       <<"\n\t\tRecHit trackLocal vector "<<LV.x() << " " << LV.y() << " " << LV.z() <<std::endl; 
00685     
00686     //Get SiStripCluster from SiStripRecHit
00687     if ( tkrecHit != NULL ){
00688       LogTrace("ClusterAnalysis") << "GOOD hit" << std::endl;
00689       const SiStripCluster* SiStripCluster_ = &*(tkrecHit->cluster());
00690       //const SiStripClusterInfo* SiStripClusterInfo_ = MatchClusterInfo(SiStripCluster_,detid);
00691       
00692       const edm::ParameterSet pset = conf_.getParameter<edm::ParameterSet>("TrackThresholds");
00693       if( track_ref->normalizedChi2() < pset.getParameter<double>("maxChi2") && track_ref->recHitsSize() > pset.getParameter<double>("minRecHit") ){                
00694         
00695         SiStripClusterInfo* clusterInfo = new SiStripClusterInfo( detid, *SiStripCluster_, es);
00696 
00697         if ( clusterInfos(clusterInfo,detid,"_onTrack", LV ) ) {
00698           vPSiStripCluster.push_back(SiStripCluster_);
00699           countOn++;
00700         }
00701         delete clusterInfo ; //CHECKME
00702       }
00703     }else{
00704       LogTrace("ClusterAnalysis") << "NULL hit" << std::endl;
00705     }     
00706   }
00707   
00708   //------------------------------------------------------------------------
00709   
00710   void ClusterAnalysis::AllClusters( const edm::EventSetup& es){
00711 
00712     LogTrace("ClusterAnalysis") << "Executing AllClusters" << std::endl;
00713     //Loop on Dets
00714     edm::DetSetVector<SiStripCluster>::const_iterator DSViter=dsv_SiStripCluster->begin();
00715     for (; DSViter!=dsv_SiStripCluster->end();DSViter++){
00716       uint32_t detid=DSViter->id;
00717 
00718       if (find(ModulesToBeExcluded_.begin(),ModulesToBeExcluded_.end(),detid)!=ModulesToBeExcluded_.end())
00719         continue;
00720       
00721       //Loop on Clusters
00722       LogTrace("ClusterAnalysis") << "\n["<<__PRETTY_FUNCTION__<<"] \n on detid "<< detid << " N Cluster= " << DSViter->data.size() <<std::endl;
00723       
00724       edm::DetSet<SiStripCluster>::const_iterator ClusIter = DSViter->data.begin();
00725       for(; ClusIter!=DSViter->data.end(); ClusIter++) {
00726 
00727         //const SiStripClusterInfo* SiStripClusterInfo_=MatchClusterInfo(&*ClusIter,detid);
00728         
00729         SiStripClusterInfo* clusterInfo = new SiStripClusterInfo( detid, *ClusIter, es);
00730 
00731         if ( clusterInfos(clusterInfo,detid,"_All", LV) ){ 
00732           countAll++;
00733 
00734           LogTrace("ClusterAnalysis") << "\n["<<__PRETTY_FUNCTION__<<"] ClusIter " << &*ClusIter << 
00735             "\t " << std::find(vPSiStripCluster.begin(),vPSiStripCluster.end(),&*ClusIter)-vPSiStripCluster.begin() << std::endl;
00736 
00737           if (std::find(vPSiStripCluster.begin(),vPSiStripCluster.end(),&*ClusIter) == vPSiStripCluster.end()){
00738             if ( clusterInfos(clusterInfo, detid,"_offTrack", LV) ) 
00739               countOff++;
00740           }
00741         }
00742         delete clusterInfo; // CHECKME
00743       }       
00744     }
00745   }
00746   
00747   //------------------------------------------------------------------------
00748 
00749  
00750   //------------------------------------------------------------------------
00751   
00752   bool ClusterAnalysis::clusterInfos( SiStripClusterInfo* clusterInfo, const uint32_t& detid,TString flag , const LocalVector LV ){
00753     LogTrace("ClusterAnalysis") << "\n["<<__PRETTY_FUNCTION__<<"]" << std::endl;
00754     
00755     if (clusterInfo==0) 
00756       return false;
00757    
00758 
00759     const  edm::ParameterSet ps_b = conf_.getParameter<edm::ParameterSet>("BadModuleStudies");
00760     if  ( ps_b.getParameter<bool>("Bad") ) {//it will perform Bad modules discrimination
00761       short n_Apv;
00762       switch((int)clusterInfo->getFirstStrip()/128){
00763       case 0:
00764         n_Apv=0;
00765         break;
00766       case 1:
00767         n_Apv=1;
00768         break;
00769       case 2:
00770         n_Apv=2;
00771         break;
00772       case 3:
00773         n_Apv=3;
00774         break;
00775       case 4:
00776         n_Apv=4;
00777         break;
00778       case 5:
00779         n_Apv=5;
00780         break;
00781       }
00782       if ( ps_b.getParameter<bool>("justGood") ){//it will analyze just good modules 
00783         LogTrace("SiStripQuality") << "Just good module selected " << std::endl;
00784         if(SiStripQuality_->IsModuleBad(detid)){
00785           LogTrace("SiStripQuality") << "\n Excluding bad module " << detid << std::endl;
00786           return false;
00787         }else if(SiStripQuality_->IsApvBad(detid, n_Apv)){
00788           LogTrace("SiStripQuality") << "\n Excluding bad module and APV " << detid << n_Apv << std::endl;
00789           return false;
00790         }
00791       }else{
00792         LogTrace("SiStripQuality") << "Just bad module selected " << std::endl;
00793         if(!SiStripQuality_->IsModuleBad(detid) || !SiStripQuality_->IsApvBad(detid, n_Apv)){
00794           LogTrace("SiStripQuality") << "\n Skipping good module " << detid << std::endl;
00795           return false;
00796         }
00797       }
00798     }
00799 
00800     const  edm::ParameterSet ps = conf_.getParameter<edm::ParameterSet>("clusterInfoConditions");
00801     if  ( ps.getParameter<bool>("On") 
00802           &&
00803           ( 
00804            //CHECKME
00805            clusterInfo->getCharge()/clusterInfo->getNoise() < ps.getParameter<double>("minStoN") 
00806            ||
00807            clusterInfo->getCharge()/clusterInfo->getNoise() > ps.getParameter<double>("maxStoN") 
00808            ||
00809            clusterInfo->getWidth() < ps.getParameter<double>("minWidth") 
00810            ||
00811            clusterInfo->getWidth() > ps.getParameter<double>("maxWidth") 
00812            )
00813           )
00814       return false;
00815 
00816     LogTrace("ClusterAnalysis") << "Executing clusterInfos for module " << detid << std::endl;
00817 
00818     const StripGeomDetUnit*_StripGeomDetUnit = dynamic_cast<const StripGeomDetUnit*>(tkgeom->idToDetUnit(DetId(detid)));
00819     //GeomDetEnumerators::SubDetector SubDet_enum=_StripGeomDetUnit->specificType().subDetector();
00820     int SubDet_enum=_StripGeomDetUnit->specificType().subDetector() -2;
00821     
00822     //&&&&&&&&&&&&&&&& GLOBAL POS &&&&&&&&&&&&&&&&&&&&&&&&
00823     const StripTopology &topol=(StripTopology&)_StripGeomDetUnit->topology();
00824     MeasurementPoint mp(clusterInfo->getPosition(),rnd.Uniform(-0.5,0.5));
00825     LocalPoint localPos = topol.localPosition(mp);
00826     GlobalPoint globalPos=(_StripGeomDetUnit->surface()).toGlobal(localPos);
00827     //&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
00828     
00829     //char cdetid[128];
00830     //sprintf(cdetid,"_%d",detid);
00831     int iflag;
00832     if (flag=="_onTrack")
00833       iflag=0;
00834     else if (flag=="_offTrack")
00835       iflag=1;
00836     else
00837       iflag=2;
00838     
00839     NClus[SubDet_enum][iflag]++;
00840 
00841     LogTrace("ClusterAnalysis") << "NClus on detid = " << detid << " " << flag << " is " << NClus[SubDet_enum][iflag] << std::endl;
00842     //    TrackerMap filling for each flag
00843     const edm::ParameterSet _mapSet = conf_.getParameter<edm::ParameterSet>("MapFlag");
00844     if( _mapSet.getParameter<bool>("Map_ClusOccOn") )
00845       tkMap_ClusOcc[iflag]->fill(detid,1);
00846  
00847  /* CHECKME
00848     std::stringstream ss;
00849     const_cast<SiStripClusterInfo*>(cluster)->print(ss);
00850     LogTrace("ClusterAnalysis") 
00851       << "\n["<<__PRETTY_FUNCTION__<<"]\n"
00852       << ss.str() 
00853       << "\n\t\tcluster LocalPos "     << localPos
00854       << "\n\t\tcluster GlobalPos "     << globalPos
00855       << std::endl;
00856 */
00857     long double tanXZ = -999;
00858     float cosRZ = -2;
00859     long double atanXZ = -200;
00860     LogTrace("ClusterAnalysis")<< "\n\tLV " << LV.x() << " " << LV.y() << " " << LV.z() << " " << LV.mag() << std::endl;
00861     if (LV.mag()!=0){
00862       double proj_yZ=_StripGeomDetUnit->surface().toGlobal(LocalVector(0.,1.,0.)).z();
00863       tanXZ= LV.x() /LV.z() * (-1) * proj_yZ/fabs(proj_yZ);
00864       cosRZ= fabs(LV.z())/LV.mag();
00865       atanXZ=atan(tanXZ)*189/Geom::pi();
00866 
00867       LogTrace("ClusterAnalysis")<< "\n\t tanXZ " << tanXZ << " cosRZ " << cosRZ << std::endl;
00868     }
00869 
00870     //Display
00871     ((TH3S*) Hlist->FindObject("ClusterGlobalPos"))
00872       ->Fill(globalPos.z(),globalPos.x(),globalPos.y());
00873 
00874     //&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
00875     //Cumulative Plots
00876 
00877     TString appString=SubDet[SubDet_enum]+flag;    
00878   
00879     fillTH1(clusterInfo->getCharge(),"cSignal"+appString,1,clusterInfo->getWidth());
00880 
00881     if (LV.mag()!=0){
00882       
00883       fillTH1(clusterInfo->getCharge()*cosRZ,"cSignalCorr"+appString,0); //Filled only for ontrack
00884       
00885       fillTProfile(atanXZ,clusterInfo->getWidth(),"ClusterWidthVsAngle"+appString,0); //Filled only for ontrack
00886 
00887       fillTH1((clusterInfo->getCharge()/clusterInfo->getNoise())*cosRZ,"cStoNCorr"+appString,0); //Filled only for ontrack
00888     } 
00889        
00890     fillTH1(clusterInfo->getNoise(),"cNoise"+appString,1,clusterInfo->getWidth());
00891     
00892     if (clusterInfo->getNoise()){
00893       fillTH1(clusterInfo->getCharge()/clusterInfo->getNoise(),"cStoN"+appString,1,clusterInfo->getWidth());
00894       
00895     }
00896       
00897     fillTH1(clusterInfo->getWidth(),"cWidth"+appString,0);
00898 
00899     fillTH1(clusterInfo->getPosition(),"cPos"+appString,1,clusterInfo->getWidth());
00900 
00901     fillTH2(clusterInfo->getPosition(),clusterInfo->getCharge()/clusterInfo->getNoise(),"cStoNVsPos"+appString,1,clusterInfo->getWidth());
00902    
00903    /*
00904     // to retrieve the (raw)digiamplitudesL and (raw)digiamplitudesR, one needs to provide
00905     // the rawdigicollection via e.g. a cfg-file. An example can be developped upon request
00906     // if needed by contacting the developper of SiStripClusterInfo(E.Delmeire)
00907      
00908     int neighbourStripNumber=3;
00909     std::string rawdigiProducer = "SiStripDigis"; 
00910     std::string rawdigiLabel = "VirginRaw"// "ProcessedRaw"
00911     edm::Handle< edm::DetSetVector<SiStripRawDigi> > rawDigiHandle;
00912     es.getByLabel(rawdigiProducer,rawdigiLabel,rawDigiHandle);
00913     
00914     vector<float> amplitudesL, amplitudesR;
00915     amplitudesL = clusterInfo->getRawDigiAmplitudesLR(neighbourStripNumber, 
00916                                                             *rawDigiHandle, 
00917                                                                dsv_SiStripCluster,
00918                                                               rawDigiLabel).first;
00919                                                               
00920     amplitudesR = clusterInfo->getRawDigiAmplitudesLR(neighbourStripNumber, 
00921                                                             *rawDigiHandle, 
00922                                                                 dsv_SiStripCluster,
00923                                                               rawDigiLabel).second;
00924                                               
00925     if (amplitudesL.size()!=0 ||  amplitudesR.size()!=0){
00926         
00927       float Ql=0;
00928       float Qr=0;
00929       float Qt=0;
00930 
00931       if (EtaAlgo_==1){
00932         Ql=clusterInfo->getChargeLR().first;
00933         Qr=clusterInfo->getChargeLR().second;
00934       
00935         for (std::vector<float>::const_iterator it=amplitudesL.begin(); it !=amplitudesL.end() && it-amplitudesL.begin()<NeighStrips_; it ++)
00936           { Ql += (*it);}
00937 
00938         
00939         for (std::vector<float>::const_iterator it=amplitudesR.begin(); it !=amplitudesR.end() && it-amplitudesR.begin()<NeighStrips_; it ++)
00940           { Qr += (*it);}
00941         
00942         Qt=Ql+Qr+clusterInfo->getMaxCharge();
00943       }
00944       else{
00945         
00946         int Nstrip=clusterInfo->getStripAmplitudes().size();
00947         float pos=clusterInfo->getPosition()-0.5;
00948         for(int is=0;is<Nstrip && clusterInfo->getFirstStrip()+is<=pos;is++)
00949           Ql+=clusterInfo->getStripAmplitudes()[is];
00950         
00951         Qr=clusterInfo->getCharge()-Ql;
00952 
00953         for (std::vector<float>::const_iterator it=amplitudesL.begin(); it !=amplitudesL.end() && it-amplitudesL.begin()<NeighStrips_; it ++)
00954           { Ql += (*it);}
00955         
00956         for (std::vector<float>::const_iterator it=amplitudesR.begin(); it !=amplitudesR.end() && it-amplitudesR.begin()<NeighStrips_; it ++)
00957           { Qr += (*it);}
00958 
00959         
00960         Qt=Ql+Qr;
00961       }
00962       
00963       LogTrace("ClusterAnalysis") << "\n["<<__PRETTY_FUNCTION__<<"] \n on detid "<< detid << " Ql=" << Ql << " Qr="<< Qr << " Qt="<<Qt<< " eta="<< Ql/Qt<< std::endl;
00964       
00965       fillTH1(Ql/Qt,"cEta"+appString,1,clusterInfo->getWidth());
00966     
00967       fillTH2((Ql-Qr)/Qt,(Ql+Qr)/Qt,"cEta_scatter"+appString,1,clusterInfo->getWidth());
00968       
00969     }
00970     */
00971     //&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
00972     //Detector Detail Plots    
00973     char aname[128];
00974     //sprintf(aname,"%s_%d",_StripGeomDetUnit->type().name().c_str(),detid);
00975     SiStripDetId a(detid);
00976     if ( a.subdetId() == 3 ){
00977       TIBDetId b(detid);
00978       sprintf(aname,"_SingleDet_%d_TIB_%d_%d_%d_%d",detid,b.layer(),b.string()[0],b.string()[1],b.glued());
00979     } else if ( a.subdetId() == 4 ) {
00980       TIDDetId b(detid);
00981       sprintf(aname,"_SingleDet_%d_TID_%d_%d_%d_%d",detid,b.wheel(),b.ring(),b.side(),b.glued());
00982     } else if ( a.subdetId() == 5 ) {
00983       TOBDetId b(detid);
00984       sprintf(aname,"_SingleDet_%d_TOB_%d_%d_%d_%d",detid,b.layer(),b.rod()[0],b.rod()[1],b.glued());
00985     } else if ( a.subdetId() == 6 ) {
00986       TECDetId b(detid);
00987       sprintf(aname,"_SingleDet_%d_TEC_%d_%d_%d_%d_%d",detid,b.wheel(),b.ring(),b.side(),b.glued(),b.stereo());
00988     }        
00989     appString=TString(aname);
00990         
00991     if(flag=="_All"){
00992 
00993       fillTH1(clusterInfo->getCharge(),"cSignal"+appString,0);
00994 
00995       fillTH1(clusterInfo->getNoise(),"cNoise"+appString,0);
00996 
00997       if (clusterInfo->getNoise()){
00998         fillTH1(clusterInfo->getCharge()/clusterInfo->getNoise(),"cStoN"+appString,0);
00999       }
01000       
01001       fillTH1(clusterInfo->getWidth(),"cWidth"+appString,0);
01002 
01003       fillTH1(clusterInfo->getPosition(),"cPos"+appString,0);
01004 
01005       fillTH2(clusterInfo->getPosition(),clusterInfo->getCharge()/clusterInfo->getNoise(),"cStoNVsPos"+appString,0);
01006     }
01007 
01008     if(flag=="_onTrack" && cosRZ>-2){
01009       fillTProfile((int)(clusterInfo->getPosition()-.5)/256,clusterInfo->getCharge()*cosRZ,"cSignalxFiber"+appString,0);
01010     }
01011 
01012     //&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
01013     // Layer Detail Plots
01014     char cApp[64];
01015     sprintf(cApp,"_Layer_%d",GetSubDetAndLayer(detid).second);
01016     appString="_"+TString(GetSubDetAndLayer(detid).first)+cApp+flag;
01017 
01018     fillTH1(clusterInfo->getCharge(),"cSignal"+appString,0);
01019   
01020     fillTH1(clusterInfo->getNoise(),"cNoise"+appString,0);
01021 
01022     if (clusterInfo->getNoise()){
01023       fillTH1(clusterInfo->getCharge()/clusterInfo->getNoise(),"cStoN"+appString,0);
01024     }
01025       
01026     fillTH1(clusterInfo->getWidth(),"cWidth"+appString,0);
01027     
01028     fillTH1(clusterInfo->getPosition(),"cPos"+appString,0);      
01029     
01030     if(flag=="_onTrack" && LV.mag()!=0){
01031       fillTProfile(atanXZ,clusterInfo->getWidth(),"ClusterWidthVsAngle"+appString+"_onTrack",0);
01032       
01033       fillTH1(clusterInfo->getCharge()*cosRZ,"cSignalCorr"+appString,0);
01034       
01035       fillTProfile(cosRZ,clusterInfo->getCharge(),"cSignalVsAngle"+appString,0);
01036       fillTH2(cosRZ,clusterInfo->getCharge(),"cSignalVsAngleH"+appString,0);
01037       
01038       if (clusterInfo->getNoise()){
01039         
01040         fillTH1((clusterInfo->getCharge()/clusterInfo->getNoise())*cosRZ,"cStoNCorr"+appString,0); 
01041       }
01042       
01043       
01044       //***Only for TIB and TOB rphi modules***//      
01045       if ( StripSubdetector(detid).stereo() == 0 ){
01046         
01047         fillTProfile(
01048                      tkgeom->idToDet(DetId(detid))->surface().toGlobal(LV).phi().degrees(),
01049                      atanXZ,
01050                      "AngleVsPhi"+appString,0
01051                      );
01052         
01053         //std::cout << " detid phi sinxz " <<  appString << " R " << tkgeom->idToDet(DetId(detid))->surface().toGlobal(_HitDir._TrackingRecHit->localPosition()).mag() << "\t |  phi " << tkgeom->idToDet(DetId(detid))->surface().toGlobal(_HitDir._TrackingRecHit->localPosition()).phi().degrees() << " anglexz " << atanXZ << " | tanxz " << tanXZ << " x " << _HitDir._LV.x() << " z " << _HitDir._LV.z() << " mag " << sqrt(_HitDir._LV.x()*_HitDir._LV.x()+_HitDir._LV.z()*_HitDir._LV.z()) << " \t | " << tkgeom->idToDet(DetId(detid))->surface().toGlobal(LocalVector(0,0,1)) << " " << tkgeom->idToDet(DetId(detid))->surface().toGlobal(LocalVector(1,0,0)) << std::endl;
01054         
01055         //LogTrace("ClusterAnalysis") << " det " << appString << " angle Ph " << tkgeom->idToDet(DetId(detid))->surface().toGlobal(_HitDir._TrackingRecHit->localPosition()).phi().degrees() << " " << _HitDir._TrackingRecHit->localPosition() << "   " <<  tkgeom->idToDet(DetId(detid))->surface().toGlobal(_HitDir._TrackingRecHit->localPosition()) << std::endl; 
01056       }      
01057     }
01058     return true;
01059   }
01060 
01061   //--------------------------------------------------------------------------------
01062   std::pair<std::string,uint32_t> ClusterAnalysis::GetSubDetAndLayer(const uint32_t& detid){
01063     
01064     std::string cSubDet;
01065     uint32_t layer=0;
01066     const StripGeomDetUnit* _StripGeomDetUnit = dynamic_cast<const StripGeomDetUnit*>(tkgeom->idToDetUnit(DetId(detid)));
01067     switch(_StripGeomDetUnit->specificType().subDetector())
01068       {
01069       case GeomDetEnumerators::TIB:
01070         cSubDet="TIB";
01071         layer=TIBDetId(detid).layer();
01072         break;
01073       case GeomDetEnumerators::TOB:
01074         cSubDet="TOB";
01075         layer=TOBDetId(detid).layer();
01076         break;
01077       case GeomDetEnumerators::TID:
01078         cSubDet="TID";
01079         layer=TIDDetId(detid).ring();
01080         // SUSY Modified me on 11/9/07
01081         LogTrace("ClusterAnalysis") << "SubDet found TID ring" << layer << std::endl;
01082         break;
01083       case GeomDetEnumerators::TEC:
01084         cSubDet="TEC";
01085         layer=TECDetId(detid).ring();
01086         LogTrace("ClusterAnalysis") << "SubDet found TEC ring" << layer << std::endl;
01087         break;
01088       default:
01089         edm::LogWarning("ClusterAnalysis") << "WARNING!!! this detid does not belong to tracker" << std::endl;
01090       }
01091     return std::make_pair(cSubDet,layer);
01092   }
01093 
01094 
01095   void ClusterAnalysis::fillTH1(float value,TString name,bool widthFlag,float cwidth){
01096 
01097     for (int iw=0;iw<5;iw++){
01098       if ( iw==0 || (iw==4 && cwidth>3) || ( iw>0 && iw<4 && cwidth==iw) ){     
01099         TH1F* hh = (TH1F*) Hlist->FindObject(name+width_flags[iw]);
01100         if (hh!=0)  
01101           hh->Fill(value);
01102       }
01103       if (!widthFlag)
01104         break;
01105     }
01106   }
01107 
01108   void ClusterAnalysis::fillTProfile(float xvalue,float yvalue,TString name,bool widthFlag,float cwidth){
01109 
01110     for (int iw=0;iw<5;iw++){
01111       if ( iw==0 || (iw==4 && cwidth>3) || ( iw>0 && iw<4 && cwidth==iw) ){     
01112         TProfile* hh = (TProfile*) Hlist->FindObject(name+width_flags[iw]);
01113         if (hh!=0)  
01114           hh->Fill(xvalue,yvalue);
01115       }
01116       if (!widthFlag)
01117         break;
01118     }
01119   }
01120 
01121   void ClusterAnalysis::fillTH2(float xvalue,float yvalue,TString name,bool widthFlag, float cwidth){
01122 
01123     for (int iw=0;iw<5;iw++){
01124       if ( iw==0 || (iw==4 && cwidth>3) || ( iw>0 && iw<4 && cwidth==iw) ){     
01125         TH2F* hh = (TH2F*) Hlist->FindObject(name+width_flags[iw]);
01126         if (hh!=0)  
01127           hh->Fill(xvalue,yvalue);
01128       }
01129       if (!widthFlag)
01130         break;
01131     }
01132   }
01133 
01134   void ClusterAnalysis::bookHlist(char* HistoType, char* ParameterSetLabel, TFileDirectory subDir, TString & HistoName, char* xTitle, char* yTitle, char* zTitle){
01135     if ( HistoType == "TH1" ) {
01136       Parameters =  conf_.getParameter<edm::ParameterSet>(ParameterSetLabel);
01137       TH1F* p = subDir.make<TH1F>(HistoName,HistoName,
01138                          Parameters.getParameter<int32_t>("Nbinx"),
01139                          Parameters.getParameter<double>("xmin"),
01140                          Parameters.getParameter<double>("xmax")
01141                          );
01142       if ( xTitle != "" )
01143         p->SetXTitle(xTitle);
01144       if ( yTitle != "" )
01145         p->SetYTitle(yTitle);
01146       Hlist->Add(p);
01147     }
01148     else if ( HistoType == "TH2" ) {
01149       Parameters =  conf_.getParameter<edm::ParameterSet>(ParameterSetLabel);
01150       TH2F* p = subDir.make<TH2F>(HistoName,HistoName,
01151                          Parameters.getParameter<int32_t>("Nbinx"),
01152                          Parameters.getParameter<double>("xmin"),
01153                          Parameters.getParameter<double>("xmax"),
01154                          Parameters.getParameter<int32_t>("Nbiny"),
01155                          Parameters.getParameter<double>("ymin"),
01156                          Parameters.getParameter<double>("ymax")
01157                          );
01158       if ( xTitle != "" )
01159         p->SetXTitle(xTitle);
01160       if ( yTitle != "" )
01161         p->SetYTitle(yTitle);
01162       if ( zTitle != "" )
01163         p->SetZTitle(zTitle);
01164       Hlist->Add(p);
01165     }
01166     else if ( HistoType == "TH3" ){
01167       Parameters =  conf_.getParameter<edm::ParameterSet>(ParameterSetLabel);
01168       TH3S* p = subDir.make<TH3S>(HistoName,HistoName,
01169                          Parameters.getParameter<int32_t>("Nbinx"),
01170                          Parameters.getParameter<double>("xmin"),
01171                          Parameters.getParameter<double>("xmax"),
01172                          Parameters.getParameter<int32_t>("Nbiny"),
01173                          Parameters.getParameter<double>("ymin"),
01174                          Parameters.getParameter<double>("ymax"),
01175                          Parameters.getParameter<int32_t>("Nbinz"),
01176                          Parameters.getParameter<double>("zmin"),
01177                          Parameters.getParameter<double>("zmax")
01178                          );
01179       if ( xTitle != "" )
01180         p->SetXTitle(xTitle);
01181       if ( yTitle != "" )
01182         p->SetYTitle(yTitle);
01183       if ( zTitle != "" )
01184         p->SetZTitle(zTitle);
01185       Hlist->Add(p);
01186     }
01187     else if ( HistoType == "TProfile" ){
01188       Parameters =  conf_.getParameter<edm::ParameterSet>(ParameterSetLabel);
01189       TProfile* p = subDir.make<TProfile>(HistoName,HistoName,
01190                                   Parameters.getParameter<int32_t>("Nbinx"),
01191                                   Parameters.getParameter<double>("xmin"),
01192                                   Parameters.getParameter<double>("xmax"),
01193                                   Parameters.getParameter<double>("ymin"),
01194                                   Parameters.getParameter<double>("ymax")
01195                                   );
01196       if ( xTitle != "" )
01197         p->SetXTitle(xTitle);
01198       if ( yTitle != "" )
01199         p->SetYTitle(yTitle);
01200       Hlist->Add(p);
01201     }
01202     else{
01203       edm::LogError("ClusterAnalysis")<< "[" <<__PRETTY_FUNCTION__ << "] invalid HistoType " << HistoType << std::endl;
01204     }
01205   }
01206 
01207 
01208 
01209   void ClusterAnalysis::fillPedNoiseFromDB(){
01210     std::vector<uint32_t> vdetId_;
01211     SiStripDetCabling_->addActiveDetectorsRawIds(vdetId_);
01212     
01213     for (std::vector<uint32_t>::const_iterator detid_iter=vdetId_.begin();detid_iter!=vdetId_.end();detid_iter++){
01214       
01215       uint32_t detid = *detid_iter;
01216       
01217       if (detid < 1){
01218         edm::LogError("ClusterAnalysis")<< "[" <<__PRETTY_FUNCTION__ << "] invalid detid " << detid<< std::endl;
01219         continue;
01220       }
01221       const StripGeomDetUnit* _StripGeomDetUnit = dynamic_cast<const StripGeomDetUnit*>(tkgeom->idToDetUnit(DetId(detid)));
01222       if (_StripGeomDetUnit==0){
01223         continue;
01224       }
01225       
01226       
01227       unsigned int nstrips = _StripGeomDetUnit->specificTopology().nstrips();
01228 
01229       //&&&&&&&&&&&&&&
01230       // Retrieve information for the module
01231       //&&&&&&&&&&&&&&&&&&     
01232       char cdetid[128];
01233       sprintf(cdetid,"%d",detid);
01234       char aname[128];
01235       sprintf(aname,"%s_%d",_StripGeomDetUnit->type().name().c_str(),detid);
01236       char SubStr[128];
01237       
01238       SiStripDetId a(detid);
01239       if ( a.subdetId() == 3 ){
01240         TIBDetId b(detid);
01241         sprintf(SubStr,"_SingleDet_%d_TIB_%d_%d_%d_%d",detid,b.layer(),b.string()[0],b.string()[1],b.glued());
01242       } else if ( a.subdetId() == 4 ) {
01243         TIDDetId b(detid);
01244         sprintf(SubStr,"_SingleDet_%d_TID_%d_%d_%d_%d",detid,b.wheel(),b.ring(),b.side(),b.glued());
01245       } else if ( a.subdetId() == 5 ) {
01246         TOBDetId b(detid);
01247         sprintf(SubStr,"_SingleDet_%d_TOB_%d_%d_%d_%d",detid,b.layer(),b.rod()[0],b.rod()[1],b.glued());
01248       } else if ( a.subdetId() == 6 ) {
01249         TECDetId b(detid);
01250         sprintf(SubStr,"_SingleDet_%d_TEC_%d_%d_%d_%d_%d",detid,b.wheel(),b.ring(),b.side(),b.glued(),b.stereo());
01251       }
01252       
01253       TString appString=TString(SubStr);
01254 
01255       SiStripNoises::Range noiseRange = noiseHandle->getRange(detid);
01256       SiStripPedestals::Range pedRange = pedestalHandle->getRange(detid);
01257       
01258       SiStripQuality::Range detQualityRange = SiStripQuality_->getRange(detid);
01259      
01260       for(size_t istrip=0;istrip<nstrips;istrip++){
01261         try{
01262           //Fill Pedestals
01263           TH1F* hh1 = (TH1F*) Hlist->FindObject("DBPedestals"+appString);
01264           if (hh1!=0) 
01265             hh1->Fill(istrip,pedestalHandle->getPed(istrip,pedRange)); 
01266           
01267           //Fill Noises
01268           TH1F* hh2 = (TH1F*) Hlist->FindObject("DBNoise"+appString);
01269           if (hh2!=0)   
01270             hh2->Fill(istrip,noiseHandle->getNoise(istrip,noiseRange));     
01271           
01272           //Fill BadStripsNoise
01273           TH1F* hh3 = (TH1F*) Hlist->FindObject("DBBadStrips"+appString);
01274           if (hh3!=0)  
01275             hh3->Fill(istrip,SiStripQuality_->IsStripBad(detQualityRange,istrip)?1.:0.);
01276           
01277         }catch(cms::Exception& e){
01278           edm::LogError("SiStripCondObjDisplay") << "[SiStripCondObjDisplay::endJob]  cms::Exception:  DetName " << name << " " << e.what() ;
01279         }
01280       }
01281       
01282     }
01283   }
01284 }
ClusterAnalysis.cc
