StRoot: StEEmcSimulatorMaker/macros/runEEmcSlowMaker.C Source File

00001 class StChain;
00002 class StMuEmcCollection;
00003 
00004 class StEEmcDb;
00005 class StMuDstMaker;
00006 class TChain;
00007 class TObjArray;
00008 
00009 class StEEmcA2EMaker;
00010 
00011 
00012 StEEmcDb *myDb;
00013 StMuDstMaker* muMk;
00014 StChain *chain=0;
00015 TObjArray *HList;
00016 
00017 StEEmcA2EMaker *mEEbefore;
00018 StEEmcA2EMaker *mEEafter;
00019 
00020 // Comment out the line below to read EEMC data from StEvent (default is MuDst)
00021 #define useMuDst
00022 
00023 int runEEmcSlowMaker( int nEve=2000 ){
00024 
00025   int firstSec=1;
00026   int lastSec=12;
00027 
00028   Int_t nFiles=10; 
00029 
00030 #ifdef useMuDst
00031   char* file="rcf1308_01_2000evts.MuDst.root";
00032 #else
00033   char* file="rcf1308_01_2000evts.event.root";
00034 #endif
00035   char* inDir="/star/data47/reco/pp200/pythia6_410/11_15gev/cdf_a/y2006c/gheisha_on/p07ic/";
00036 
00037 #ifdef useMuDst
00038   char *outname="mudst";
00039 #else
00040   char *outname="stevent";
00041 #endif
00042   TString fullName=file;
00043 
00044   gROOT->LoadMacro("$STAR/StRoot/StMuDSTMaker/COMMON/macros/loadSharedLibraries.C");
00045   loadSharedLibraries();
00046   gROOT->Macro("LoadLogger.C");
00047   cout << " loading done " << endl;
00048   //gSystem->Load("StDbLib");
00049   gSystem->Load("StDbBroker");
00050   gSystem->Load("St_db_Maker");
00051   gSystem->Load("StEEmcUtil");
00052   gSystem->Load("StEEmcDbMaker");
00053   gSystem->Load("StEEmcSimulatorMaker");
00054   gSystem->Load("StEEmcA2EMaker");
00055 
00056   // create chain    
00057   chain = new StChain("StChain"); 
00058 
00059   printf("adding MuDst/StEvent from run '%s' ....\n",fullName.Data());
00060 
00061   // Now we add Makers to the chain...
00062 #ifdef useMuDst
00063   StMuDstMaker* muMk = new StMuDstMaker(0,0,inDir,fullName,".MuDst.root",nFiles);
00064   TChain* tree=muMk->chain(); assert(tree); 
00065   int nEntries=tree->GetEntries();
00066   printf("total eve in chain =%d\n",nEntries);
00067 #else
00068   StIOMaker* ioMaker = new StIOMaker;
00069   ioMaker->SetFile(gSystem->ConcatFileName(inDir, file));
00070   ioMaker->SetBranch("*", 0, "0");
00071   ioMaker->SetBranch("eventBranch", 0, "r");
00072   ioMaker->SetBranch("geantBranch", 0, "r");
00073   Int_t nEntries=nEve;
00074 #endif
00075 
00076   St_db_Maker *stDb = new St_db_Maker("StarDb", "MySQL:StarDb");
00077   stDb->SetDateTime(20031120,0);
00078   stDb->SetFlavor("sim","eemcPMTcal");
00079   stDb->SetFlavor("sim","eemcPIXcal");
00080   stDb->SetFlavor("sim","eemcPMTped");
00081   stDb->SetFlavor("sim","eemcPMTstat");
00082   stDb->SetFlavor("sim","eemcPMTname");
00083   stDb->SetFlavor("sim","eemcADCconf");
00084 
00085   new StEEmcDbMaker("eemcDb");
00086   HList=new TObjArray;
00087 
00088   // ADC 2 energy before slow simulator
00089   mEEbefore = new StEEmcA2EMaker("before");
00090   mEEbefore->database("eemcDb");
00091 #ifdef useMuDst
00092   mEEbefore->source("MuDst",1);
00093 #else
00094   mEEbefore->source("StEvent",2);
00095 #endif
00096 
00097 
00098 
00099   /*  
00100    * The slow simulator is setup by default with the best
00101    * estimates we have for the single p.e. resolution, energy
00102    * deposited by a single mip in an smd/pre/postshower layer,
00103    * and number of photoelectrons produced per mip.  If one
00104    * really wants to play with fire, these are documented
00105    * in the header file.
00106    *
00107    */
00108   StEEmcSlowMaker *slowSim=new StEEmcSlowMaker();
00109   slowSim->setDropBad(0);   // 0=no action, 1=drop chn marked bad in db
00110   slowSim->setAddPed(0);    // 0=no action, 1=ped offset from db
00111   slowSim->setSmearPed(0);  // 0=no action, 1=gaussian ped, width from db
00112   slowSim->setOverwrite(1); // 0=no action, 1=overwrite muDst/StEvent
00113 
00114 
00115 #ifdef useMuDst
00116   slowSim->setSource("MuDst");
00117 #else
00118   slowSim->setSource("StEvent");
00119 #endif
00120 
00122   slowSim->SetHList(HList);
00123 
00124 
00125   // ADC 2 energy after slow simulator
00126   mEEafter = new StEEmcA2EMaker("after");
00127   mEEafter->database("eemcDb");
00128 #ifdef useMuDst
00129   mEEafter->source("MuDst",1);
00130 #else
00131   mEEafter->source("StEvent",2);
00132 #endif
00133 
00134 
00135 
00136   // some test histograms
00137   TH2F *h_compare=new TH2F("h_compare","ADC after slow maker vs before; uncorrected ADC; corrected ADC",512,0.,1024.,512,0.,1024. );
00138   TH2F *h_compare_pre1=new TH2F("h_compare_pre1","ADC after slow maker vs before",256,0.,1024.,256,0.,256. );
00139   TH2F *h_compare_pre2=new TH2F("h_compare_pre2","ADC after slow maker vs before",256,0.,1024.,256,0.,256. );
00140   TH2F *h_testit=new TH2F("h_testit","ADC after-before vs ADC pre1+pre2; pre1+pre2 AD; corrected-uncorrected",256,0.,512.,100,-50.,50.);
00141   TH2F *h_perinc=new TH2F("h_perinc","percent increase vs before; ADC [before simu]; % increase",128,0.,1024.,120,-10.,50.);
00142   TH2F *h_perincE=new TH2F("h_perincE","percent increase vs original E; Energy [before simu]; % increase",150,0.,30.,120,-10.,50.);
00143   TH2F *h_perincEZ=new TH2F("h_perincEZ","percent increase vs original E (zeros supressed); Energy [before simu]; % increase",150,0.,30.,120,-10.,50.);
00144 
00145 
00146 
00147 
00148   
00149   // Put your analysis maker here
00150   // StMyAnaMaker *bhla =new  StMyAnaMaker;
00151 
00152   chain->ls(3);
00153   chain->Init();
00154 
00155   StEEmcDb* myDb = (StEEmcDb*)chain->GetDataSet("StEEmcDb");
00156   myDb->setSectors(firstSec, lastSec);
00157 
00158   int eventCounter=0;
00159   int stat=0;
00160   int t1=time(0);
00161   //---------------------------------------------------
00162   while ( stat==0 ) {// loop over events
00163     if(eventCounter>=nEve) break;
00164     eventCounter++;
00165     chain->Clear();
00166     stat = chain->Make();
00167     if(eventCounter%100==0)
00168       printf("\n====================%d  processing  ============\n", eventCounter);
00169 
00170     for ( Int_t ii=0;ii<mEEafter->numberOfHitTowers(0);ii++ )
00171       {
00172         StEEmcTower tafter=mEEafter->hittower(ii,0);
00173         StEEmcTower tbefore=mEEbefore->tower(tafter.index(),0);
00174         StEEmcTower pbefore=mEEbefore->tower(tafter.index(),1);
00175         StEEmcTower qbefore=mEEbefore->tower(tafter.index(),1);
00176         h_compare->Fill( tbefore.adc(), tafter.adc());
00177 
00178         h_testit->Fill( pbefore.adc()+pbefore.adc(), tafter.adc()-tbefore.adc() );
00179 
00180         h_perinc->Fill( tbefore.adc(), 100.0* (tafter.adc()-tbefore.adc() )/tafter.adc() );
00181         h_perincE->Fill( tbefore.energy(), 100.0* (tafter.adc()-tbefore.adc() )/tafter.adc() );
00182         if ( TMath::Abs(tbefore.adc() - tafter.adc()) >1. )
00183           h_perincEZ->Fill( tbefore.energy(), 100.0* (tafter.adc()-tbefore.adc() )/tafter.adc() );
00184 
00185       }
00186 
00187     
00188   }
00189   printf("sorting done, nEve=%d of %d\n",nEve, nEntries);
00190   //  return;
00191   int t2=time(0);
00192   if(t2==t1) t2=t1+1;
00193   float tMnt=(t2-t1)/60.;
00194   float rate=1.*nEve/(t2-t1);
00195   printf("sorting done %d of   nEve=%d, CPU rate=%.1f Hz, total time %.1f minute(s) \n",eventCounter,nEntries,rate,tMnt);
00196 
00197   chain->Finish();
00198 
00199   //  return;
00200   // save output histograms 
00201   // HList.ls();
00202 
00203 
00204 
00206 #if 1  
00207   TString out="./";
00208   out+=outname;
00209   out+=".hist.root";
00210   TFile f( out,"recreate");
00211   assert(f.IsOpen());
00212   printf("%d histos are written  to '%s' ...\n",HList->GetEntries(),out.Data());
00213   HList->Write();
00214   h_compare->Write();
00215   h_testit->Write();
00216   h_perinc->Write();h_perincE->Write();h_perincEZ->Write();
00217   f.Close();
00218 #endif  
00219 
00220 }