StRoot: StTriggerUtilities/macros/rdMu2TrigSimu.C Source File

00001 // Before running macro create an output directory "outL2" if you want a copy of the L2 output histos
00002 // This macro runs on real data (flagMC=0) and MC files (flagMC=1)
00003 // Set flag ==1 for those detectors you want included in the trigger decisions 
00004 // Set configuration for BEMC (online, offline or custom).  EEMC is the same for offline and online
00005 // Choose the correct configuration year if you want to use L2
00006 
00007 int total=0;
00008 class StChain *chain=0; 
00009 char *eemcSetupPath="/afs/rhic.bnl.gov/star/users/kocolosk/public/StarTrigSimuSetup/";
00010 
00011 void rdMu2TrigSimu(char *file="/star/data47/reco/pp200/pythia6_410/9_11gev/cdf_a/y2006c/gheisha_on/p07ic/rcf1309_*_2000evts.MuDst.root"){
00012   
00013   int nevents = 200;
00014   int flagMC=1;  // 0/1 == Using Real/Simulation data files 
00015   int useEemc=1; // 0/1 == Exclude/Include EEMC in Trigger Decisions 
00016   int useBemc=1; // 0/1 == Exclude/Include BEMC in Trigger Decisions 
00017   int useL2=1;   // 0/1 == Exclude/Include L2 in Trigger Decisions 
00018   int L2ConfigYear=2006; // possible: 2006, 2008
00019   int bemcConfig=2; // Online==1, Offline==2, Expert==3
00020   int playConfig=0; // jan:100_199
00021   int emcEveDump=0; // extrating raw EMC data in a custom format
00022   int outputL2Histo=0;//output L2 histos to directory outL2
00023   TString outDir="./outL2/"; 
00024 
00025 
00026   gROOT->LoadMacro("$STAR/StRoot/StMuDSTMaker/COMMON/macros/loadSharedLibraries.C");
00027   loadSharedLibraries();
00028   assert( !gSystem->Load("StDetectorDbMaker"));
00029   assert( !gSystem->Load("StDbUtilities"));
00030   assert( !gSystem->Load("StDbBroker"));
00031   assert( !gSystem->Load("St_db_Maker"));
00032   assert( !gSystem->Load("StEEmcUtil")); // needed by eemcDb
00033   assert( !gSystem->Load("StEEmcDbMaker"));
00034   assert( !gSystem->Load("StDaqLib")); // needed by bemcDb
00035   assert( !gSystem->Load("StEmcRawMaker"));
00036   assert( !gSystem->Load("StEmcADCtoEMaker"));
00037   if (flagMC) {
00038     assert( !gSystem->Load("StMcEvent"));
00039     assert( !gSystem->Load("StMcEventMaker"));
00040     assert( !gSystem->Load("StEmcSimulatorMaker"));
00041     assert( !gSystem->Load("StEEmcSimulatorMaker"));
00042     assert( !gSystem->Load("StEpcMaker"));
00043   }
00044   assert( !gSystem->Load("StTriggerUtilities"));
00045 
00046   gROOT->Macro("LoadLogger.C");
00047   cout << " loading done " << endl;
00048   
00049   chain= new StChain("StChain"); 
00050 
00051   if (flagMC){
00052     TString geantFile;
00053     geantFile += file;
00054     geantFile.ReplaceAll("MuDst.root", "geant.root");
00055     printf("geantFile=%s\n", geantFile.Data());
00056     StIOMaker* ioMaker = new StIOMaker();
00057     ioMaker->SetFile(geantFile);
00058     ioMaker->SetIOMode("r");
00059     ioMaker->SetBranch("*",0,"0");             //deactivate all branches
00060     ioMaker->SetBranch("geantBranch",0,"r");   //activate geant Branch
00061     StMcEventMaker *evtMaker = new StMcEventMaker();
00062   }
00063 
00064   //Need MuDstMaker to get data
00065   printf(" Analyzing file=%s\n",file);  
00066   StMuDstMaker* muDstMaker =new StMuDstMaker(0,0,"",file,"",1000);
00067 
00068   //Database -- get a real calibration from the database
00069   St_db_Maker* dbMk =0;
00070   if(useEemc || useL2) // full DB access
00071     dbMk = new St_db_Maker("StarDb","MySQL:StarDb","MySQL:StarDb","$STAR/StarDb");
00072   else // only Barrel is uploaded, is faster 
00073     dbMk  = new St_db_Maker("Calibrations","MySQL:Calibrations_emc");
00074 
00075     
00076   //If MC then must set database time and date
00077   //If Endcap fast simu is used tower gains in DB do not matter,JB
00078   if (flagMC) dbMk->SetDateTime(20060522, 55000);//timestamp R7142018
00079   
00080   //Collect all output histograms 
00081   TObjArray* HList=new TObjArray; 
00082   
00083   //Endcap DB
00084   if(useEemc || useL2) new StEEmcDbMaker("eemcDb");
00085   
00086 
00087   //Get BEMC adc values
00088   if (flagMC && useBemc) {
00089     StEmcSimulatorMaker* emcSim = new StEmcSimulatorMaker(); //use this instead to "redo" converstion from geant->adc
00090     if (bemcConfig == 1) {
00091         emcSim->setCheckStatus(kBarrelEmcTowerId,false); //this returns hits regardless of offline tower status
00092     }
00093     emcSim->setCalibSpread(kBarrelEmcTowerId,0.15);//spread gains by 15%
00094   }
00095   if (flagMC==0 && useBemc){
00096     StEmcADCtoEMaker *bemcAdc = new StEmcADCtoEMaker();//for real data this sets calibration and status
00097     if (bemcConfig == 1) {
00098         bemcAdc->setCheckStatus(kBarrelEmcTowerId,false);
00099     }
00100   }
00101 
00102   //must use slow simulator to get pedestals correct for L2
00103   if (flagMC==1 && useEemc){
00104     StEEmcSlowMaker *slowSim = new StEEmcSlowMaker("slowSim");
00105     slowSim->setSamplingFraction(0.0384); // effectively scales all Tower energies with a factor of 1.3 (added by: Ilya Selyuzhenkov; April 11, 2008)
00106     slowSim->setAddPed(true);
00107     slowSim->setSmearPed(true);
00108   }
00109  
00110   //Get TriggerMaker
00111   StTriggerSimuMaker *simuTrig = new StTriggerSimuMaker("StarTrigSimu");
00112   simuTrig->setHList(HList);
00113   simuTrig->setMC(flagMC); // must be before individual detectors, to be passed
00114   simuTrig->useBbc();
00115   if(useEemc) {
00116     simuTrig->useEemc(0);//default=0:just process ADC, 1,2:comp w/trgData,see .
00117     simuTrig->eemc->setSetupPath(eemcSetupPath);
00118   } 
00119   if(useBemc){
00120     simuTrig->useBemc();
00121     simuTrig->bemc->setConfig(bemcConfig);
00122   }
00123 
00124   if(flagMC && useEemc==2){
00125     // pass one argument to M-C as generic switch    
00126     // Endcap specific params -- ok Jan you need to change this to a default "online" setup
00127     int eemcDsmSetup[20]; // see StEemcTriggerSimu::initRun() for definition
00128     memset(eemcDsmSetup, 0,sizeof(eemcDsmSetup));// clear all, may be a bad default
00129     eemcDsmSetup[0]=3;  // HTthr0
00130     eemcDsmSetup[1]=12; // HTthr1
00131     eemcDsmSetup[2]=22; // HTthr2
00132     eemcDsmSetup[3]=1;  // TPthr0
00133     eemcDsmSetup[4]=17; // TPthr1
00134     eemcDsmSetup[5]=31; // TPthr2
00135     eemcDsmSetup[10]=2; //HTTPthrSelc, 2=use_thres_#1
00136     simuTrig->eemc->setDsmSetup(eemcDsmSetup);    
00137   }
00138 
00139 
00140   if(useL2) {
00141     /* 
00142        reads all input/setup files from  L2setup-yyyymmdd/
00143        writes all output files to L2out-yyyymmdd 
00144        depending on the DB time stamp 
00145        both dierectiorie MUST exist, setup must be reasonable
00146     */
00147     StGenericL2Emulator* simL2Mk=0;
00148     if(L2ConfigYear==2006) simL2Mk= new StL2_2006EmulatorMaker;
00149     else if(L2ConfigYear==2008) simL2Mk= new StL2_2008EmulatorMaker;
00150     assert(simL2Mk);
00151     simL2Mk->setSetupPath(eemcSetupPath);
00152     simL2Mk->setOutPath(outDir.Data());
00153     if (flagMC) simL2Mk->setMC();
00154     simuTrig->useL2(simL2Mk);
00155   }
00156 
00157   
00158   //if(emcEveDump) new StJanEventMaker;
00159     
00160   StTriggerSimuPlayMaker *playMk= new StTriggerSimuPlayMaker; // to develope user  analysis of trigQA 
00161   playMk->setConfig(playConfig);
00162   playMk->setHList(HList);
00163   
00164   
00165   chain->ls(3);
00166   chain->Init();
00167  
00168   for (Int_t iev=0;iev<nevents; iev++) {
00169     cout << "\n****************************************** " << endl;
00170     cout << "Working on eventNumber:\t" << iev <<"\tof:\t"<<nevents<<endl;
00171     cout << "****************************************** " << endl;
00172     chain->Clear();
00173     int iret = chain->Make(iev);
00174     total++;   
00175     if (iret % 10 == kStEOF || iret % 10 == kStFatal)  {
00176       cout << "Bad return code!" << endl;
00177       break;
00178     }
00179 
00180     int trigID[3]={127213,117211,137611};
00181     StMuDst *muDst = muDstMaker->muDst();
00182     StMuEvent *muEvent = muDst->event();    
00183     StMuTriggerIdCollection trig = muEvent -> triggerIdCollection();
00184     StTriggerId l1trig = trig.nominal();
00185     if( l1trig.isTrigger(trigID[0])) {
00186       cout<<" SimuTrigger 127213 ="<<simuTrig->isTrigger(trigID[0])<<" BEMC="<<simuTrig->bemc->triggerDecision(trigID[0])<<" L2="<<simuTrig->lTwo->triggerDecision(trigID[0])<<endl;
00187      }
00188     if( l1trig.isTrigger(trigID[1])) {
00189       cout<<" SimuTrigger 117211 ="<<simuTrig->isTrigger(trigID[1])<<" BEMC="<<simuTrig->bemc->triggerDecision(trigID[1])<<" L2="<<simuTrig->lTwo->triggerDecision(trigID[1])<<endl;
00190      }
00191     if( l1trig.isTrigger(trigID[2])) {
00192       cout<<" SimuTrigger 137611 ="<<simuTrig->isTrigger(trigID[2])<<" BEMC="<<simuTrig->bemc->triggerDecision(trigID[2])<<" L2="<<simuTrig->lTwo->triggerDecision(trigID[2])<<endl;
00193     }
00194 
00195     
00196     StTriggerSimuResult trigResult = simuTrig->detailedResult(trigID[2]);
00197     if (trigResult.bemcDecision()==1){
00198        vector<short> towerId = trigResult.highTowerIds();
00199       for (unsigned i=0; i<towerId.size(); i++) {
00200         cout<<" LO Trigger BEMC Tower="<<towerId[i]<<" adc="<<trigResult.highTowerAdc(towerId[i])<<endl;
00201       }
00202     }
00203 
00204 
00205    if (trigResult.l2Decision()==1){
00206        vector<short> towerId = trigResult.highTowerIds();
00207       for (unsigned i=0; i<towerId.size(); i++) {
00208         cout<<" L2 Trigger BEMC Tower="<<towerId[i]<<" adc="<<trigResult.highTowerAdc(towerId[i])<<endl;
00209       }
00210     }
00211  
00212   }
00213 
00214  
00215   chain->Finish();
00216   cout << "****************************************** " << endl;
00217   cout << "total number of events  " << total << endl;
00218   cout << "****************************************** " << endl;
00219 
00220 
00221   if (outputL2Histo==1) {
00222   
00223     TString fileMu=file;
00224     printf("=%s=\n",fileMu.Data());
00225     if(fileMu.Contains(".lis")) fileMu.ReplaceAll(".lis",".trgSim");
00226     if(fileMu.Contains(".MuDst.root")) fileMu.ReplaceAll(".MuDst.root",".trgSim");
00227     TString outF=outDir+fileMu;
00228     outF+=".hist.root";
00229     printf("=%s=\n",outF.Data());
00230     hf=new TFile(outF,"recreate");
00231     if(hf->IsOpen()) {
00232       HList->Write();
00233       printf("\n Histo saved -->%s<\n",outF.Data());
00234     } 
00235     else {
00236       printf("\n Failed to open Histo-file -->%s<, continue\n",outF.Data());
00237     }
00238   }
00239   
00240   //cout <<Form("sorting done %d of   nEve=%d, CPU rate=%.1f Hz, total time %.1f minute(s) \n\n",total,nEntries,rate,tMnt)<<endl;
00241  
00242 }
00243 
00244 
00245 
00246 
00247 //========================================
00248 //========================================
00249 #if 0
00250 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
00251 clean it up, Jan
00252 
00253   const int mxTr=13;
00254   int myTrgList[mxTr]    ={127580,127551,127271,127571,127821,127831,127575,127622,127221,127611,117705,127501,127652};
00255   // prescale as for run 7101015
00256   char* myTrgListN[mxTr]={"eemc-http*116","eemc-jp0-mb*221", "eemc-jp1-mb*1", "bemc-jp1*433", 
00257                           "bemc-http-mb-fast*1", "eemc-http-mb-fast*1", "bemc-jp0-etot*979" , 
00258                           "bemc-jp0-etot-mb-L2jet*2", "bemc-jp1-mb*1", "bemc-http-mb-l2gamma*1",
00259                           "jpsi-mb*20", "bemc-jp0-mb*909","eemc-jp0-etot-mb-L2jet*2"};
00260   int nR[mxTr],nS[mxTr],nRS[mxTr];
00261   memset(nR,0, sizeof(nR));
00262   memset(nS,0, sizeof(nR));
00263   memset(nRS,0, sizeof(nR));
00264 
00265   int BL1_ADC[6];
00266     cout<<Form(" Simu trgSize=%d, firedID:", simuTrig->mTriggerList.size())<<endl;
00267     int j;
00268     for(j=0; j<simuTrig->mTriggerList.size();j++) 
00269       cout<<Form("s%d, ", simuTrig->mTriggerList[j]);
00270     cout<<"\n";
00271       
00272 
00273      // Access to muDst .......................
00274     StMuEvent* muEve = muDstMaker->muDst()->event();
00275 
00276     StEventInfo &info=muEve->eventInfo();
00277     StMuTriggerIdCollection &tic=muEve->triggerIdCollection();
00278     vector<unsigned int> trgL=tic.nominal().triggerIds();
00279     cout<<Form(" real trgSize=%d, firedID:",trgL.size());
00280     for(j=0; j<trgL.size();j++) 
00281       cout<<Form("r%d, ", trgL[j]);
00282     cout<<"\n";
00283 
00284     //.. compare Simu vs.  real
00285     for(j=0;j<mxTr;j++) {
00286       int trgId=myTrgList[j];
00287       bool realT=tic.nominal().isTrigger(trgId);
00288       bool simT=simuTrig->isTrigger(trgId); // endcap only, tmp
00289       if(realT) nR[j]++;
00290       if(realT && simT) nRS[j]++;
00291       if(simT) nS[j]++;
00292         cout <<Form("C:j=%d  trg=%d  R=%d S=%d  RS=%d",j, myTrgList[j],realT,simT,realT && simT )<<endl;
00293      }
00294 
00295    for(j=0;j<mxTr;j++) {
00296      cout <<Form("SUM  trg=%d  nR=%d nS=%d  nRS=%d  %s", myTrgList[j], nR[j],nS[j],nRS[j],myTrgListN[j])<<endl;
00297    }
00298  
00299 
00300 #endif