StRoot: StPmdDiscriminatorMaker/StPmdDiscriminatorMaker.cxx Source File

00001 /***************************************************************
00002  *
00003  * $Id: StPmdDiscriminatorMaker.cxx,v 1.10 2007/04/26 04:13:19 perev Exp $
00004  * Author: Subhasis Chattopadhyay
00005  ***************************************************************
00006  *
00007  * Description: StPmdDiscriminationMaker is class for photon/hadron
00008  * discrimination 
00009  *
00010  ****************************************************************
00011  * $Log: StPmdDiscriminatorMaker.cxx,v $
00012  * Revision 1.10  2007/04/26 04:13:19  perev
00013  * Remove StBFChain dependency
00014  *
00015  * Revision 1.9  2005/10/21 00:54:06  subhasis
00016  * fillStEvent crash due to mismatch in length of pmdclust and phmdclust fixed
00017  *
00018  * Revision 1.8  2005/02/23 05:17:18  subhasis
00019  * DiscbyNN option fixed
00020  *
00021  * Revision 1.7  2004/10/30 00:08:00  subhasis
00022  * TranFlag added and set to 0 in ctor
00023  *
00024  * Revision 1.6  2004/07/16 14:29:32  subhasis
00025  * more checks on edep Discriminate
00026  *
00027  * Revision 1.5  2004/07/13 12:59:11  subhasis
00028  *  Ncluster() is put with check in PrepareInput()
00029  *
00030  * Revision 1.4  2003/09/02 17:58:48  perev
00031  * gcc 3.2 updates + WarnOff
00032  *
00033  * Revision 1.3  2003/08/04 18:53:44  perev
00034  * warnOff
00035  *
00036  * Revision 1.2  2003/05/29 13:12:51  subhasis
00037  * several changes to include NN
00038  *
00039  ****************************************************************/
00040 #include<Stiostream.h>
00041 #include<assert.h>
00042 #include<math.h>
00043 #include"TROOT.h"
00044 #include<TRandom.h>
00045 #include<TBrowser.h>
00046 #include<TPad.h>
00047 #include<StMessMgr.h>
00048 #include<TFile.h>
00049 
00050 #include "StPmdDiscriminator.h"
00051 #include "StPmdDiscriminatorMaker.h"
00052 #include "StPmdDiscriminatorNN.h"
00053 #include "StPmdUtil/StPmdCollection.h"
00054 #include "StPmdUtil/StPmdDetector.h"
00055 #include "StPmdUtil/StPmdClusterCollection.h"
00056 #include "StPmdUtil/StPmdCluster.h"
00057 #include "StNNCluster.h"
00058 #include "StEventTypes.h"
00059 
00060 StPmdCl NNContainer;
00061 Int_t EventProcessed=0;
00062 Int_t EtaPhiArr_CPV[400][1400];
00063 
00064 ClassImp(StPmdDiscriminatorMaker)
00065   TDataSet *clusterIn;
00066 StPmdDiscriminatorMaker *discm; 
00067 StPmdCollection *cluster_hit;
00069 
00070 bool discbyNN=false;
00071 StPmdDiscriminatorMaker::StPmdDiscriminatorMaker(const char *name):StMaker(name){
00072   mApplyFlagNN=0;
00073   mTrainFlag=0;
00074 }
00075 
00076 StPmdDiscriminatorMaker::~StPmdDiscriminatorMaker()
00077 {
00078   // !destructor
00079 }
00080 
00081 Int_t StPmdDiscriminatorMaker::Init()
00082 {
00083   bookHistograms();
00084   return StMaker::Init();
00085 }
00086 
00087 void StPmdDiscriminatorMaker::bookHistograms()
00088 {
00089   mEdepPmd = new TH1F("EdepPmd","Energy deposited",100,0.,0.02);
00090   mEtaPmd  = new TH1F("EtaPmd","Eta Distribution",100,-6., 6.);
00091   mPhiPmd  = new TH1F("PhiPmd","Phi Distribution",100,-3.14,3.14);
00092   mMCPID  = new TH1F("McPID","Mc PID",11,-0.5,10.5);
00093   mcpvmatched  = new TH1F("cpvmatched","CPV Matched",100,0,100.);
00094   mEtadiff  = new TH1F("DeltaEta","Eta Difference",200,0.,200.);
00095   mPhidiff  = new TH1F("DeltaPhi","Phi Difference",4000,0.,3000.);
00096   mEtaPhi = new TH2F("DeltaEtaDetaPhi","(Delta) Eta Phi Distribution",100,0,200.,4000,0.,3000.);
00097   mClusterPID  = new TH1F("ClusterPID","ClusterPID based on matching",10,0.5,10.5);
00098   mClusterEdepPID  = new TH1F("ClusterEdepPID","ClusterPID based on Edep",10,0.5,10.5);
00099   mEdepVsPID  = new TH2F("EdepVsPID","PID based on Edep",100,0.,200.,100,0.,10.);
00100   mNNoutput  = new TH1F("NNTrainOut","Training output from NN",100,0.,2.);
00101 }
00102 
00103 Int_t StPmdDiscriminatorMaker::Make()
00104 {
00105   EventProcessed++;
00106   
00108   StEvent *currevent = (StEvent*)GetInputDS("StEvent");
00109   if(!currevent){
00110     cout << "PMDDISC ***** Can not get StEvent pointer . sorry \n";                     return kStErr;
00111   }
00112   
00113   StPhmdCollection* phmdcl = (StPhmdCollection*)currevent->phmdCollection();        
00114   clusterIn = GetDataSet("PmdSimulator");
00115   cluster_hit = (StPmdCollection*)clusterIn->Find("PmdCollection");
00116   StPmdDetector* pmd_det;
00117   StPmdDetector* cpv_det;
00118   if(cluster_hit){
00119     pmd_det = cluster_hit->detector(Int_t(0));
00120     cpv_det = cluster_hit->detector(Int_t(1));
00121   }
00122   
00123   if(phmdcl){
00124     StPhmdClusterCollection* cluscoll;
00125     StPhmdClusterCollection* cpvcluscoll;
00126     for(Int_t d=0;d<2;d++){
00127       StDetectorId pdet=static_cast<StDetectorId>(kPhmdCpvId+d);
00128       StPhmdDetector* detector=(StPhmdDetector*)phmdcl->detector(StDetectorId(pdet));
00129       if(detector){
00130         if(d==1)cluscoll = detector->cluster();
00131         if(d==0)cpvcluscoll = detector->cluster();
00132         if(d==1){
00133           if(discbyNN){
00134             Int_t ret=PrepareInputforNN(pmd_det,cpv_det,cluscoll,cpvcluscoll);
00135             if(ret==1)goto DISC;
00136           }
00137         }
00138         if(d==1&&cluscoll)cout<<"ncls "<<cluscoll->numberOfclusters()<<endl;
00139       }
00140     }
00141   }
00142   cout<<"DISCMAKER*****"<<mApplyFlagNN<<endl;
00143   if(mApplyFlagNN==1){
00144     
00145     StPmdDiscriminatorNN *discNN = new StPmdDiscriminatorNN(NNContainer);
00146     discNN->setApplyFlag(1);
00147     discNN->setDisMaker(this);
00148 
00149     if(discNN)
00150       {
00151         discNN->Discriminate();
00152       }
00153     NNContainer.clear();
00154     
00155   }
00156   
00157  DISC:
00158   
00159   if(cluster_hit){
00160     StPmdDetector * pmd_det = cluster_hit->detector(Int_t(0));
00161     StPmdDetector * cpv_det = cluster_hit->detector(Int_t(1));
00162     StPmdDiscriminator *disc = new StPmdDiscriminator(mEdepThreshold, pmd_det, cpv_det);
00163     
00164     if(disc){
00165       disc->SetEdepcut(0.0000063);  
00166       cout<<"calling edep_disc "<<endl;
00167 
00168       disc->Discriminate();   
00169     //  disc->Print();  //! Print photon like hits
00170     }
00171     fillHistograms(pmd_det,cpv_det);
00172   }     
00173   
00175   
00176   fillStEvent(cluster_hit,phmdcl);
00177   
00179   
00180   return kStOK;
00181 }
00182 
00183 
00184 void StPmdDiscriminatorMaker::Browse(TBrowser *b)
00185 {
00186   TDataSet::Browse(b);
00187 }
00188 
00189 void StPmdDiscriminatorMaker::fillHistograms(StPmdDetector* pmd_det, StPmdDetector* cpv_det)
00190 {
00191   StPmdClusterCollection* clusters = (StPmdClusterCollection*)pmd_det->cluster();
00192   Int_t nclust=0;
00193   if(clusters){nclust = clusters->Nclusters();}
00194 
00195   if(nclust>0){
00196   TIter next(clusters->Clusters());
00197           StPmdCluster *spmcl1;
00198   for(Int_t i=0; i<nclust ; i++)
00199     {
00200       spmcl1 = (StPmdCluster*)next();
00201       //      Float_t eta=spmcl1->CluEta();
00202       //      Float_t phi=spmcl1->CluPhi();
00203       Float_t edep=spmcl1->CluEdep();
00204       Int_t PID=spmcl1->CluPID();
00205       Int_t EdepPID=spmcl1->CluEdepPID();
00206       
00207          mEdepVsPID->Fill(edep*1000000,Float_t(EdepPID));
00208          mClusterPID->Fill(Float_t(PID));
00209          mClusterEdepPID->Fill(Float_t(EdepPID));
00210     }
00211   }
00212 }
00213 
00214 
00215 Int_t StPmdDiscriminatorMaker::PrepareInputforNN(StPmdDetector* pmd_det,StPmdDetector* cpv_det,StPhmdClusterCollection* cluscoll,StPhmdClusterCollection* cpvcluscoll)
00216   
00217 {
00218   cout<<"prepareInput***"<<endl;
00219   Double_t etapmd,phipmd,edeppmd;
00220   Double_t etacpv,phicpv,edepcpv;
00221 
00222   Bool_t mcpvdone;
00223   
00224   for(Int_t ibineta=0;ibineta<400;ibineta++){
00225     for(Int_t ibinphi=0;ibinphi<1400;ibinphi++){
00226       EtaPhiArr_CPV[ibineta][ibinphi]=-999;
00227     }
00228   }
00229   
00230   StPmdClusterCollection* clusters = (StPmdClusterCollection*)pmd_det->cluster();
00231   StPmdClusterCollection* clusters1 = (StPmdClusterCollection*)cpv_det->cluster();
00232   //  Int_t nclust = clusters->Nclusters(); 
00233 
00234   if(cluscoll){
00235     Int_t Ncluster0=cluscoll->numberOfclusters();
00236     if(Ncluster0>0)
00237       {
00238         const StSPtrVecPhmdCluster& pmdclusters= cluscoll->clusters();
00239         const StSPtrVecPhmdCluster& cpvclusters= cpvcluscoll->clusters();
00240         
00241         cout<<"PMD cluster size "<<pmdclusters.size()<<endl;
00242         cout<<"CPV cluster size "<<cpvclusters.size()<<endl;
00243         
00244 if(clusters1){
00245           Int_t nclust1=0;
00246           nclust1 = clusters1->Nclusters(); 
00247         TIter nextCPV(clusters1->Clusters());
00248         StPmdCluster *spmcl2;
00249         for(Int_t icpv=0; icpv<nclust1 ; icpv++)
00250           {
00251             Int_t etabin,phibin;
00252             spmcl2 = (StPmdCluster*)nextCPV();
00253             etacpv=spmcl2->CluEta();
00254             phicpv=spmcl2->CluPhi();
00255             edepcpv=spmcl2->CluEdep();
00256             
00257             if(fabs(etacpv)>2. && fabs(etacpv)<4.){etabin=Int_t((fabs(etacpv)-2)/mDeltaEta);}
00258             if(fabs(phicpv)<3.14){phibin=Int_t((phicpv+3.14)/mDeltaPhi);}
00259             EtaPhiArr_CPV[etabin][phibin]=icpv;
00260           }
00261   }     
00262 
00263 
00264         Int_t cpvmatched=0;
00265         Int_t hadmatched=0;
00266   StPmdCluster *spmcl1;
00267   if(clusters){
00268           TIter next(clusters->Clusters());
00269         for(UInt_t i=0;i<pmdclusters.size();i++)
00270           {
00271             Int_t etabin,phibin;
00272             mcpvdone=true;
00273             spmcl1 = (StPmdCluster*)next();
00274             StPhmdCluster *cl1=(StPhmdCluster*)pmdclusters[i];
00275             cl1->setMcPid(spmcl1->McCluPID());
00276             
00277             StNNCluster* nncls = new StNNCluster();
00278             nncls->setPmdCluster(cl1);
00279             
00280             etapmd=spmcl1->CluEta();
00281             phipmd=spmcl1->CluPhi();
00282             edeppmd=spmcl1->CluEdep();
00283             mMCPID->Fill(Float_t(spmcl1->McCluPID()));
00284             
00285             if(fabs(etapmd)>2. && fabs(etapmd)<4.)etabin=Int_t((fabs(etapmd)-2)/mDeltaEta);
00286             if(fabs(phipmd)<3.14)phibin=Int_t((phipmd+3.14)/mDeltaPhi);
00287             
00288             for(Int_t ibineta=0;ibineta<400;ibineta++){
00289               for(Int_t ibinphi=0;ibinphi<1400;ibinphi++){
00290                 if(EtaPhiArr_CPV[ibineta][ibinphi]!=-999 && EtaPhiArr_CPV[ibineta][ibinphi]!=0){
00291                   Float_t etadiff=abs(ibineta-etabin);
00292                   Float_t phidiff=abs(ibinphi-phibin);
00293                   mEtaPhi->Fill(etadiff,phidiff);
00294                   mEtadiff->Fill(etadiff);
00295                   mPhidiff->Fill(phidiff);
00296                 }
00297               }
00298             }
00299             
00300             if(EtaPhiArr_CPV[etabin][phibin]!=-999){
00301               Int_t index=EtaPhiArr_CPV[etabin][phibin];
00302               nncls->setCpvCluster((StPhmdCluster*)cpvclusters[index]);
00303               if(spmcl1->McCluPID()==8){
00304                 hadmatched++;
00305                 spmcl1->setCluPID(8); //PID 8 set for photon
00306               }
00307               if(spmcl1->McCluPID()==1){
00308                 hadmatched++;
00309                 spmcl1->setCluPID(1); //PID 8 set for photon
00310               }
00311               cpvmatched++;
00312             }
00313             
00314             
00315             NNContainer.push_back(nncls); 
00316           }
00317           }
00318         if(cpvmatched>0)mcpvmatched->Fill(cpvmatched);
00319         
00320       }
00321   }
00322   return 0;
00323 }
00324 
00325 Int_t StPmdDiscriminatorMaker::Finish()
00326 {
00327   if(mApplyFlagNN!=1 && mTrainFlag==1){
00328     StPmdDiscriminatorNN *discNN = new StPmdDiscriminatorNN(NNContainer);
00329     discNN->setApplyFlag(0);
00330     discNN->setDisMaker(this);
00331     if(discNN)
00332       {
00333         discNN->Discriminate();
00334       }         
00335   }
00336   return StMaker::Finish();
00337 }
00338 
00339 
00340 
00341 void StPmdDiscriminatorMaker::fillStEvent(StPmdCollection* cluster_hit,StPhmdCollection* phmdcl)
00342 {
00343   
00344   if(cluster_hit){
00345     StPmdDetector * pmd_det = cluster_hit->detector(Int_t(0));
00346     StPmdClusterCollection* clusters = (StPmdClusterCollection*)pmd_det->cluster();
00347     if(phmdcl){
00348 
00349       StDetectorId pdet=static_cast<StDetectorId>(kPhmdId);
00350       StPhmdDetector* detector=(StPhmdDetector*)phmdcl->detector(StDetectorId(pdet));
00351       if(detector){
00352         StPhmdClusterCollection* cluscoll = detector->cluster();
00353         if(cluscoll){
00354           const StSPtrVecPhmdCluster& pmdclusters= cluscoll->clusters();
00355           Int_t nclust = clusters->Nclusters();
00356           TIter next(clusters->Clusters());
00357           StPmdCluster *spmcl1;
00358           int nclust_phmd=pmdclusters.size();
00359 
00360           int loopcls=0;
00361           if(nclust==nclust_phmd)loopcls=nclust;
00362           if(nclust<nclust_phmd)loopcls=nclust;
00363           if(nclust>nclust_phmd)loopcls=nclust_phmd;
00364 
00365           for(Int_t i=0; i<loopcls ; i++)
00366             {
00367               spmcl1 = (StPmdCluster*)next();
00368               int edeppid=spmcl1->CluEdepPID();
00369               StPhmdCluster *cl1=(StPhmdCluster*)pmdclusters[i];
00370                 if(!cl1)break;
00371               if(cl1)cl1->setEnergyPid(edeppid);
00372             }
00373         }
00374       }
00375     }
00376   }
00377 }
00378