StRoot: StEEmcPool/macros/runEEmcPreAnalysis.C Source File

00001 //
00002 // Runs StMuEEmcPreAnalysis maker, performs a (trivial) clustering
00003 //   around the high-tower in each event, and generates some histograms
00004 //   of tower, preshower and postshower response, and SMD response...
00005 //
00006 // Binning on the histograms may not be optimal...
00007 //
00008 
00009 
00010 // forward declarations
00011 class StChain;
00012 class StIOMaker;
00013 class StMuDstMaker;
00014 class StMuEEmcPreAnalysisMaker;
00015 class StMuDbReader;
00016 class StEEmcDb;
00017 class St_db_Maker;
00018 
00019 StChain                  *chain;
00020 StIOMaker                *ioMaker;
00021 StMuDstMaker             *muDstMaker;
00022 StMuEEmcPreAnalysisMaker *muEEmcAnal;
00023 StMuDbReader             *muDb;
00024 StEEmcDb                 *eemdDb;
00025 St_db_Maker              *starDb;
00026 
00027 #include <iomanip>
00028 
00029 void runEEmcPreAnalysis ( Int_t nevents = 100,  // Number of events
00030                       Int_t smdhist = 10,   // Prescale on SMD histograms
00031                       Char_t *muDst = "MuDst/mcPi0n_field_onpi0_10000_06TC05_5.MuDst.root",
00032                       Char_t *output = "output.root" // Output file
00033                       );
00034 
00036 
00037 void runEEmcPreAnalysis ( Int_t nevents, Int_t smdhist, Char_t *muDst, Char_t *output ) {
00038 
00039 
00040   if ( output != "" ) 
00041     TFile *f = new TFile( output, "RECREATE" );
00042 
00044   //
00045   // Preshower, postshower and tower 2D histos
00046   //
00047   TH2F *hPre1VsTow  = new TH2F("hPre1VsTow", "Preshower(1) energy deposit vs reco tower E", 100, 0., 20., 100, 0., 0.2);
00048   TH2F *hPre2VsTow  = new TH2F("hPre2VsTow", "Preshower(2) energy deposit vs reco tower E", 100, 0., 20., 100, 0., 0.2);
00049   TH2F *hPostVsTow  = new TH2F("hPostVsTow", "Postshower   energy deposit vs reco tower E", 100, 0., 20., 100, 0., 0.2);
00050   TH2F *hPre1VsPre2 = new TH2F("hPre1VsPre2","Preshower(1) energy deposit vs preshower(2) energy deposit", 100, 0., 0.2, 100, 0., 0.2);
00051 
00052 
00054   //
00055   // Load shared libraries
00056   //
00057   gROOT -> LoadMacro("$STAR/StRoot/StMuDSTMaker/COMMON/macros/loadSharedLibraries.C");
00058   loadSharedLibraries();
00059   gSystem->Load("StEEmcUtil.so");
00060   gSystem->Load("StMuEEmcPreAnalysisMaker.so");
00061   gSystem->Load("StDbLib");
00062   gSystem->Load("StDbBroker");
00063   gSystem->Load("St_db_Maker");
00064   gSystem->Load("StEEmcDbMaker");
00065 
00067   //
00068   // Create the chain and add the muDst maker
00069   //
00070   chain      = new StChain("StChain");
00071   muDstMaker = new StMuDstMaker(0,0,"",muDst,"MuDst.root",1);
00072   muDb       = StMuDbReader::instance();
00073 
00074   starDb = new St_db_Maker("StarDb", "MySQL:StarDb");
00075   new StEEmcDbMaker("eemcDb");
00076 
00077   muEEmcAnal = new StMuEEmcPreAnalysisMaker();
00078 
00080   //
00081   // Initialize the chain in preparation to loop over events
00082   //
00083   chain -> Init();
00084   chain -> ls(3);
00085   eemcDb = (StEEmcDb*)chain->GetDataSet("StEEmcDb");
00086   assert(eemcDb);
00087   //  starDb -> setTimeStampDay(20030516);
00088   eemcDb -> setTimeStampDay(20040101);
00089   eemcDb -> setPreferedFlavor( "set492", "eemcPMTcal" );
00090 
00092   // 
00093   // Event Loop
00094   //
00095   
00096   Int_t status = 0;
00097   Int_t ievent = 0;
00098   while ( !status && ievent < nevents ) {
00099 
00100     if ( !ievent % 100 ) std::cout << "Processing event number " << ievent << std::endl;
00101 
00102     // Clear all makers on the chain
00103     chain -> Clear();
00104 
00105     // Call all of the makers on the chain
00106     status = chain -> Make();
00107 
00108     // Get the high-tower from the EEMC analysis maker
00109     StMuEEmcTower *seed = muEEmcAnal -> getHighTower();
00110     //    StMuEEmcTower *neighbors[] = seed -> getNeighbors();
00111 
00112     // We take the 8 surrounding towers, plus the central
00113     //   "seed" tower, to be the cluster.  The following
00114     //   loop determines the total energy of the cluster,
00115     //   total energy deposited in pre and postshower
00116     //   layers, and determines the first and last U and
00117     //   V smd strips within the "fiducial volume" of the
00118     //   cluster.
00119 
00120     Float_t energy      = seed -> getEnergy(0);
00121     Float_t energy_pre1 = seed -> getEnergy(1);
00122     Float_t energy_pre2 = seed -> getEnergy(2);
00123     Float_t energy_post = seed -> getEnergy(3);
00124 
00125     StMuEEmcStrip *ufirst = seed -> getFirstStrip('U');
00126     StMuEEmcStrip *vfirst = seed -> getFirstStrip('V');
00127     StMuEEmcStrip *ulast  = seed -> getLastStrip('U');
00128     StMuEEmcStrip *vlast  = seed -> getLastStrip('V');
00129 
00130     std::cout << ievent << " seed:    " 
00131               << std::setw(7) << std::setprecision(3) << energy << " "
00132               << std::setw(7) << std::setprecision(3) << energy_pre1 << " "
00133               << std::setw(7) << std::setprecision(3) << energy_pre2 << " "
00134               << std::setw(7) << std::setprecision(3) << energy_post << " " 
00135               << std::setw(3) << ufirst -> getId() << " "
00136               << std::setw(3) << ulast -> getId() << " "
00137               << std::setw(3) << vfirst -> getId() << " "
00138               << std::setw(3) << vlast -> getId() << " "
00139               << std::endl;
00140 
00141     StMuEEmcTower *neighbor;
00142 
00143     for ( Int_t i = 0; i < 8; i++ ) {
00144 
00145       // Get a pointer to the neighboring tower, and 
00146       // punt if its NULL
00147       neighbor = seed -> getNeighbor(i);
00148       if ( !neighbor ) continue;
00149 
00150       energy      += neighbor -> getEnergy(0);
00151       energy_pre1 += neighbor -> getEnergy(1);
00152       energy_pre2 += neighbor -> getEnergy(2);
00153       energy_post += neighbor -> getEnergy(3);
00154 
00155       StMuEEmcStrip *mytmp = neighbor -> getFirstStrip('U');
00156       if ( mytmp -> getId() < ufirst -> getId() ) ufirst = mytmp;
00157       mytmp = neighbor -> getFirstStrip('V');
00158       if ( mytmp -> getId() < vfirst -> getId() ) vfirst = mytmp;
00159       mytmp = neighbor -> getLastStrip('U');
00160       if ( mytmp -> getId() > ulast -> getId() ) ulast = mytmp;
00161       mytmp = neighbor -> getLastStrip('V');
00162       if ( mytmp -> getId() > vlast -> getId() ) vlast = mytmp;
00163 
00164     }
00165 
00166     // Fill preshower/postshower histos
00167     hPre1VsTow -> Fill( energy, energy_pre1 );
00168     hPre2VsTow -> Fill( energy, energy_pre2 );
00169     hPostVsTow -> Fill( energy, energy_post );
00170     hPre1VsPre2-> Fill( energy_pre2, energy_pre1 );
00171 
00173     //
00174     // How often do we create a histogram for
00175     //   the SMD strips?
00176     //
00177     if ( !(ievent % smdhist) ) {
00178 
00179       f->cd();
00180 
00181       // Book and fill a histogram for SMD's (on event-by-event basis)
00182       TString myname = "uEvent"; myname += ievent;
00183       TString mytitle = "U SMD strip energy distribution";
00184       Int_t min = ufirst -> getId() - 1;
00185       Int_t max = ulast  -> getId() + 1;
00186       Int_t nbin = ( max - min + 1 );
00187       TH1F *uhist = new TH1F(myname,mytitle,nbin,min,max);
00188       
00189       StMuEEmcStrip *strip;
00190       for ( strip = ufirst; strip <= ulast; strip++ ) {
00191         //$$$uhist -> Fill( (Float_t)strip->getId(), strip -> getEnergy() );
00192 
00193         Int_t   id     = strip -> getId() - min;
00194         Float_t myenergy = strip -> getEnergy();
00195         Float_t nphoto = (myenergy>0.) ? strip -> getNumPhotoElectrons() : 1;
00196         
00197         uhist -> SetBinContent( id, myenergy );
00198         uhist -> SetBinError  ( id, myenergy / sqrt(nphoto) );
00199         
00200       }
00201       
00202       myname = "vEvent"; myname += ievent;
00203       mytitle = "V SMD strip energy distribution";
00204       min = vfirst -> getId() - 1;
00205       max = vlast  -> getId() + 1;
00206       nbin = ( max - min + 1 );
00207       TH1F *vhist = new TH1F(myname,mytitle,nbin,min,max);
00208       
00209       for ( strip = vfirst; strip <= vlast; strip++ ) {
00210         //$$$vhist -> Fill( (Float_t)strip->getId(), strip -> getEnergy() );
00211 
00212         Int_t   id     = strip -> getId() - min;
00213         Float_t myenergy = strip -> getEnergy();
00214         Float_t nphoto = (myenergy>0.) ? strip -> getNumPhotoElectrons() : 1;
00215         
00216         vhist -> SetBinContent( id, myenergy );
00217         vhist -> SetBinError  ( id, myenergy / sqrt(nphoto) );
00218 
00219       }
00220 
00221     }
00223 
00224 
00225 
00226     // Let the user know we're still alive
00227     std::cout << ievent << " cluster: " 
00228               << std::setw(7) << std::setprecision(3) << energy << " "
00229               << std::setw(7) << std::setprecision(3) << energy_pre1 << " "
00230               << std::setw(7) << std::setprecision(3) << energy_pre2 << " "
00231               << std::setw(7) << std::setprecision(3) << energy_post << " " 
00232               << std::setw(3) << ufirst -> getId() << " "
00233               << std::setw(3) << ulast -> getId() << " "
00234               << std::setw(3) << vfirst -> getId() << " "
00235               << std::setw(3) << vlast -> getId() << " "
00236               << std::endl
00237               << std::endl;
00238 
00239     // Increment the event counter
00240     ievent++; 
00241 
00242   }
00243 
00244   f -> Write();
00245 
00246 }