StRoot: StPmdDiscriminatorMaker/StPmdDiscriminatorNN.cxx Source File

00001 /****************************************************
00002  *
00003  * $Id: StPmdDiscriminatorNN.cxx,v 1.4 2007/04/26 04:13:46 perev Exp $
00004  *
00005  * Author: Subhasis Chattopadhyay
00006  *
00007  ******************************************************
00008  *
00009  * Description: Class for discrimination through energy
00010  * cut is defined.
00011  *
00012  ******************************************************
00013  *
00014  * $Log: StPmdDiscriminatorNN.cxx,v $
00015  * Revision 1.4  2007/04/26 04:13:46  perev
00016  * Remove StBFChain dependency
00017  *
00018  * Revision 1.3  2003/09/02 17:58:48  perev
00019  * gcc 3.2 updates + WarnOff
00020  *
00021  * Revision 1.2  2003/08/04 18:53:44  perev
00022  * warnOff
00023  *
00024  * Revision 1.1  2003/05/29 13:21:05  subhasis
00025  * NN discriminator
00026  *
00027  *
00028  ******************************************************/
00029 
00030 #include<Stiostream.h>
00031 #include"Stiostream.h"
00032 #include<assert.h>
00033 #include<math.h>
00034 #include"TROOT.h"
00035 #include<TRandom.h>
00036 #include<TBrowser.h>
00037 #include<TPad.h>
00038 #include<StMessMgr.h>
00039 #include<TFile.h>
00040 
00041 #include "StPmdUtil/StPmdGeom.h"
00042 #include "StPmdUtil/StPmdDetector.h"
00043 #include "StPmdDiscriminatorMaker.h"
00044 #include "StPmdDiscriminatorNN.h"
00045 #include "StPmdUtil/StPmdClusterCollection.h"
00046 #include "StPmdUtil/StPmdCluster.h"
00047 #include "StPmdNeuNet.h" 
00048 #include "StEventTypes.h" 
00049 #include "StNNCluster.h" 
00050 
00051 ClassImp(StPmdDiscriminatorNN)
00052   ofstream fileo("test_NN",ios::app);   
00053 
00054 Int_t npmdvalue=0;
00055 Int_t ncpvvalue=0;
00056 Int_t Trained=0;
00057 Int_t NTrain=0;
00058 
00059 StPmdDiscriminatorNN::StPmdDiscriminatorNN(StPmdCl cl_con)
00060 {
00061   mClContainer=cl_con;
00062   mApplyFlagNN=0;
00063   cout<<"inside discNN, size**"<<cl_con.size()<<" "<<mClContainer.size()<<endl;
00064 }
00065 
00066 
00067 StPmdDiscriminatorNN::StPmdDiscriminatorNN(StPmdDetector* pmd_det, StPmdDetector* cpv_det)
00068 {
00069   m_PmdDet=pmd_det;
00070   m_CpvDet=cpv_det;
00071 }
00072 
00073 StPmdDiscriminatorNN::~StPmdDiscriminatorNN()
00074 {
00075   //destructor
00076 }
00077 
00078 void StPmdDiscriminatorNN::Discriminate()
00079 {
00080   StPmdNeuNet *sneu=new StPmdNeuNet("for PMD",4,"4",1);
00081   sneu->setDiscMaker(m_DiscMaker);
00082   sneu->SetLearnParam(0.2); // the learning parameter (<1)
00083   sneu->SetInitParam(-2,2); // bounds for the initialisation of weights
00084   sneu ->SetUseBiases();
00085   sneu->Init(); // initialisation of the kernel
00086   sneu->PrintS(); // printing of network structure
00087   Int_t NNSize=mClContainer.size();
00088   sneu->SetNTrainEvents(NNSize);
00089   sneu->SetArraySize(NNSize);
00090   cout<<"nTrainevts  "<<sneu->GetNTrainEvents()<<endl;
00091   Input(sneu);
00092   cout<<"No of INputs **"<<NTrain<<"NNFlag  "<<mApplyFlagNN<<endl;
00093   
00094   if(mApplyFlagNN!=1){
00095     if(Trained!=1){
00096       sneu->TrainNCycles(100);
00097       Trained=1;
00098       sneu->Export("NNTrain.out"); // printing of network structure
00099     }
00100   }
00101   
00102   cout<<" NNTrain to be imported "<<endl;
00103   sneu->Import("NNTrain.out"); // printing of network structure
00104   if(mApplyFlagNN==1)
00105     {
00106       Float_t Teach[20000];
00107       Float_t Value[20000];
00108       for(Int_t i=0;i<20000;i++){Teach[i]=999; Value[i]=0;}
00109       
00110       cout<<"NN Valid() Called "<<endl;
00111       sneu->ApplyWeights(Teach,Value);
00112       for(Int_t i=0;i<20000;i++){
00113         if(Teach[i]!=999)fileo<<Teach[i]<<" "<<Value[i]<<endl;
00114       }
00115     }
00116 }
00117 
00118 void StPmdDiscriminatorNN::Input(StPmdNeuNet* sneu)
00119 {
00120   TFile * file=new TFile("nninput.root","RECREATE");
00121   m_NNedep_ph=new TH1F("nnedp_ph","(ph) PMD edep",100,0.,1000.);
00122   m_NNncell_ph=new TH1F("nn_ncell_ph","(ph) PMD ncell",100,0.,20.);
00123   m_NNsigma_ph=new TH1F("nn_sigma_ph","(ph) PMD sigma",100,0.,20.);
00124   m_NNedep_cpv_ph=new TH1F("nnedp_cpv_ph","(ph) CPV edep",100,0.,100.);
00125   m_NNedep_had=new TH1F("nnedp_had","(had) PMD edep",100,0.,1000.);
00126   m_NNncell_had=new TH1F("nn_ncell_had","(had) PMD ncell",100,0.,20.);
00127   m_NNsigma_had=new TH1F("nn_sigma_had","(had) PMD sigma",100,0.,20.);
00128   m_NNedep_cpv_had=new TH1F("nnedp_cpv_had","(had) CPV edep",100,0.,100.);
00129   
00130   
00131   Int_t totno = 0,totcpvno = 0;
00132   Float_t aveEnergy = 0., aveNcell = 0., aveSigma = 0.,aveCpvEnergy = 0;
00133   Float_t totEnergy = 0., totNcell = 0., totSigma = 0.,totCpvEnergy = 0;
00134   //  cout<<"CONTAINER SIZE ********* "<<mClContainer.size()<<endl;
00135   for(UInt_t i=0;i<mClContainer.size();i++)
00136     {
00137       StPhmdCluster *cl1=(StPhmdCluster*)(mClContainer[i]->PmdCluster());
00138       StPhmdCluster *cl2=(StPhmdCluster*)(mClContainer[i]->CpvCluster());
00140       
00141       totEnergy = totEnergy + cl1->energy();
00142       totNcell = totNcell +  cl1->numberOfCells();
00143       totSigma = totSigma + cl1->sigma();
00144       totno++;
00145       if(cl2){
00146         totCpvEnergy = totCpvEnergy + cl2->energy();
00147         totcpvno++;
00148       }
00149     }
00150   aveEnergy = totEnergy/totno;
00151   aveNcell =  totNcell/totno;
00152   aveSigma = totSigma/totno;
00153   aveCpvEnergy = totCpvEnergy/totcpvno;
00154   
00155   if(mApplyFlagNN==1){
00156     // Put avergaes by hand
00157     
00158     /*Single particle
00159       aveEnergy =0.000253643 ;
00160       aveNcell =2.54512;
00161       aveSigma =0.467905;
00162       aveCpvEnergy =1.8775e-05;
00163     */
00164     // AuAu
00165     aveEnergy =0.000157085 ;
00166     aveNcell =2.20981;
00167     aveSigma =0.476074;
00168     aveCpvEnergy =1.26994e-05;
00169     
00170   }
00171   
00172   fileo<<" AveEnergy "<<aveEnergy<<" AveNcell "<<aveNcell<<" aveSigma "<<aveSigma<<" aveCpvEne "<<aveCpvEnergy<<" TotNo "<<totno<<" TotCPvNo "<<totcpvno<<endl; 
00173   
00174   Float_t outEnergy,outNcell,outSigma,outCpvEnergy;
00175   
00176   for(UInt_t i=0;i<mClContainer.size();i++)
00177     {
00178       StPhmdCluster *cl1=(StPhmdCluster*)(mClContainer[i]->PmdCluster());
00179       StPhmdCluster *cl2=(StPhmdCluster*)(mClContainer[i]->CpvCluster());
00180       Float_t target;
00181       if(cl1->mcPid()==1)target=1.;
00182       if(cl1->mcPid()==8)target=0.;
00183       sneu->fillArrayOut(target,i,0);
00184       
00185 //VP      Int_t sm=cl1->module();
00186       Float_t energy=cl1->energy();
00187       InputRange(energy,aveEnergy,outEnergy);
00188       sneu->FillArray(i,0,outEnergy);
00189       
00190       Int_t ncell=cl1->numberOfCells();
00191       InputRange(ncell,aveNcell,outNcell);
00192 
00193       sneu->FillArray(i,1,outNcell);
00194       
00195       Float_t sigma=cl1->sigma();
00196       InputRange(sigma,aveSigma,outSigma);
00197 
00198       sneu->FillArray(i,2,outSigma);
00199       
00200       Float_t cpv_energy=0.;
00201       if(cl2){
00202         cpv_energy=cl2->energy();
00203         InputRange(cpv_energy,aveCpvEnergy,outCpvEnergy);
00204 
00205         sneu->FillArray(i,3,outCpvEnergy);
00206       }
00207 
00208       if(target==1)m_NNncell_ph->Fill(Float_t(ncell));
00209       if(target==0)m_NNncell_had->Fill(Float_t(ncell));
00210       if(target==1)m_NNsigma_ph->Fill(sigma);
00211       if(target==0)m_NNsigma_had->Fill(sigma);
00212       if(target==1&&cpv_energy>0)m_NNedep_cpv_ph->Fill(cpv_energy);
00213       if(target==0&&cpv_energy>0)m_NNedep_cpv_had->Fill(cpv_energy);
00214       if(target==1)m_NNedep_ph->Fill(energy);
00215       if(target==0)m_NNedep_had->Fill(energy);
00216       NTrain++;
00217     }
00218   cout<<"In Input **"<<NTrain<<endl;
00219   m_NNedep_ph->Write();
00220   m_NNncell_ph->Write();
00221   m_NNsigma_ph->Write();
00222   m_NNedep_cpv_ph->Write();
00223   
00224   m_NNedep_had->Write();
00225   m_NNncell_had->Write();
00226   m_NNsigma_had->Write();
00227   m_NNedep_cpv_had->Write();
00228   file->Close();
00229 }
00230 
00231 void StPmdDiscriminatorNN::setFormula()
00232 {
00233   cout<<"In the Neural Network"<<npmdvalue<<ncpvvalue<<endl;    
00234 }
00235 
00236 //function used for scaling the input variables to (-1,1)
00237 
00238 Float_t StPmdDiscriminatorNN::InputRange(Float_t Input,Float_t aveInput, Float_t& Output)
00239 {
00240   Float_t fx;
00241   Float_t ax = 1.;
00242   if(aveInput !=0.){
00243     if((Input/aveInput)<10.){
00244       fx = (2./(1. + exp(-ax*Input/aveInput))) - 1;
00245     }
00246     else
00247       {
00248         fx = (2./(1. + exp(-10.))) - 1.;
00249       }
00250     Output = 1 - 2. * fx;
00251   }
00252   return Output;
00253 }
00254         
00255