loon: UberDST.cxx Source File

00001 //#include <strstream>
00002 #include <cmath>
00003 #include <sstream>
00004 #include <set>
00005 #include "TTree.h"
00006 #include "TFile.h"
00007 #include "MINF_Classes/MINFast.h"
00008 #include "CandData/CandRecord.h"
00009 #include "CandData/CandHeader.h"
00010 #include "CalDetSI/CandCalDetSIHandle.h"
00011 #include "CandDigit/CandDigitListHandle.h"
00012 #include "CandDigit/CandDigitHandle.h"
00013 #include "REROOT_Classes/REROOT_Event.h"
00014 #include "REROOT_Classes/REROOT_FLSDigit.h"
00015 #include "REROOT_Classes/REROOT_StdHep.h"
00016 #include "RerootExodus/RerootExodus.h"
00017 #include "Plex/PlexStripEndId.h"
00018 #include "Plex/PlexSEIdAltLItem.h"
00019 #include "Plex/PlexHandle.h"
00020 #include "Plex/PlexPixelSpotId.h"
00021 #include "Navigation/NavSet.h"
00022 #include "Calibrator/Calibrator.h"
00023 #include "Conventions/CalDigitType.h"
00024 #include "Conventions/CalTimeType.h"
00025 #include "Conventions/Munits.h"
00026 #include "Conventions/ReadoutType.h"
00027 #include "MinosObjectMap/MomNavigator.h"
00028 #include "JobControl/JobCResult.h"
00029 #include "JobControl/JobCModuleRegistry.h" // JOBMODULE registration macro
00030 #include "JobControl/JobCommand.h"
00031 #include "MessageService/MsgService.h"
00032 #include "CalDetSI/Helpers.h"
00033 #include "CalDetDST/UberHit.h"
00034 #include "CalDetDST/UberEvent.h"
00035 #include "CalDetDST/UberHeader.h"
00036 #include "CalDetDST/UberMC.h"
00037 #include "CalDetDST/UberDST.h"
00038 
00039 CVSID("$Id: UberDST.cxx,v 1.11 2006/06/14 04:43:55 rhatcher Exp $");
00040 JOBMODULE(UberDST, "UberDST", "A CalDet ntuple");
00041 
00042 using namespace std;
00043 
00044 //const Float_t UberDST::fMC_PE_Conversion = 1.0/100.0; // 1pe == 100 mc adc counts
00045 //const Float_t UberDST::fMC_MIP_Conversion = 1.0/(100.0*3.567) ; // 3.567 pes/mip 
00046 
00047 const Float_t UberDST::fMC_PE_Conversion = 1.0; // taken care of by calibrator
00048 const Float_t UberDST::fMC_MIP_Conversion = 1.0/(60.0*3.567) ; // 3.567 pes/mip 
00049 
00050 
00051 //________________________________________________________________________________
00052 UberDST::UberDST():
00053   fMadeTrees(kFALSE),
00054   fFilledHeader(kFALSE),
00055   fIsMC(kFALSE),
00056   fFirstEvent(kTRUE),
00057   fMCDecision(kFALSE),
00058   snarltree(0),
00059   headertree(0),
00060   mctree(0),
00061   fRootFile(0),
00062   fCDLH(0),
00063   fCDSI(0)
00064 {
00065 
00066   MSG("UberDST",Msg::kDebug)<<"In UberDST creator"<<endl;
00067   snarldata = new UberEvent();
00068   header = new UberHeader();
00069   mc = new UberMC();
00070 }//end UberDST()
00071 //________________________________________________________________________________
00072 UberDST::~UberDST()
00073 {
00074   delete snarldata;
00075   delete header;
00076   delete mc;
00077 }//end ~UberDST()
00078 //________________________________________________________________________________
00079 JobCResult UberDST::Get(MomNavigator* mom)
00080 {
00081   MSG("UberDST",Msg::kDebug)<<"In UberDST get"<<endl;
00082 
00083   //ask mom  for a candidate record
00084   CandRecord* candrec=dynamic_cast<CandRecord*>
00085     (mom->GetFragment("CandRecord","PrimaryCandidateRecord"));
00086   if(candrec==0){
00087     MSG("UberDST",Msg::kError)<<"candrec==0"<<endl;
00088     return JobCResult(JobCResult::kFailed);
00089   }
00090   if(fRunNumber==0){
00091     const CandHeader *ch = candrec->GetCandHeader();
00092     if(ch!=0){
00093       fRunNumber=ch->GetRun();
00094       header->SetRunNo(fRunNumber);
00095       MSG("UberDST",Msg::kDebug)<<"RUN NUMBER "<<fRunNumber<<endl;
00096     }
00097   }
00098 
00099   fCDLH = dynamic_cast<CandDigitListHandle *>
00100     (candrec->FindCandHandle("CandDigitListHandle","canddigitlist"));
00101   if(fCDLH == 0) {//failure
00102     MSG("UberDST",Msg::kError)<<"fCDLH==0"<<endl;
00103     return JobCResult(JobCResult::kFailed);
00104   }
00105 
00106   //now ask candrecord for a CandCalDetSIHandle
00107   fCDSI = dynamic_cast<CandCalDetSIHandle *>
00108     (candrec->FindCandHandle("CandCalDetSIHandle"));
00109   if(fCDSI == 0) {//failure
00110     MSG("UberDST",Msg::kError)<<"fCDSI==0"<<endl;
00111     return JobCResult(JobCResult::kFailed);
00112   }
00113 
00114   MSG("UberDST",Msg::kDebug)<<"Returning kPassed"<<endl;
00115   
00116   return JobCResult(JobCResult::kPassed); //success
00117 }//end Get
00118 //________________________________________________________________________________
00119 JobCResult UberDST::Ana(const MomNavigator * /* mom */)
00120 {
00121   MSG("UberDST",Msg::kDebug)<<"In UberDST Ana"<<endl;
00122   //reset event data members
00123   snarldata->ResetEvent();
00124   mc->ResetEvent();
00125 
00126   //get the time since beginning of timeframe from from candigitlisthandle
00127   snarldata->triggertime = fCDLH->GetAbsTime();
00128 
00129   //get pid info from the CandCalDetSIHandle
00130   snarldata->snarlno = fCDSI->GetSnarl();
00131   snarldata->triggerword = fCDSI->GetTrigSource();
00132   snarldata->ceradc[0]=fCDSI->GetKovADC2();
00133   snarldata->ceradc[1]=fCDSI->GetKovADC1();
00134   snarldata->ceradc[2]=fCDSI->GetKovADC3();
00135 
00136   //get time relative to trigger time of cerenkov hits
00137   //note:  CalDetConstants::FARTIMECONVERT()--magic number to convert ticks to 
00138   //nanoseconds for the FAR DETECTOR electronics
00139   //we will need to fix this if the cerenkov is ever read out
00140   //using Near detector electronics (ie 2003)
00141   if(fCDSI->GetKovTimeStamp2()!=0){
00142     snarldata->certime[0]=fCDSI->GetKovTimeStamp2()
00143       *CalDetConstants::FARTIMECONVERT-snarldata->triggertime/Munits::nanosecond;
00144     MSG("UberDST",Msg::kDebug)<<"Kov 2 timestamp "<<fCDSI->GetKovTimeStamp2()
00145                               <<" time convert "<<CalDetConstants::FARTIMECONVERT
00146                               <<" trigger time "<<snarldata->triggertime
00147                               <<" in ns "
00148                               <<snarldata->triggertime/Munits::nanosecond
00149                               <<endl;
00150   }
00151   if(fCDSI->GetKovTimeStamp1()!=0){
00152     snarldata->certime[1]=fCDSI->GetKovTimeStamp1()
00153       *CalDetConstants::FARTIMECONVERT-snarldata->triggertime/Munits::nanosecond;
00154   }
00155   if(fCDSI->GetKovTimeStamp3()!=0){
00156     snarldata->certime[2]=fCDSI->GetKovTimeStamp3()
00157       *CalDetConstants::FARTIMECONVERT-snarldata->triggertime/Munits::nanosecond;
00158   }
00159   snarldata->toftdc[0]=fCDSI->GetTofTDC0();
00160   snarldata->toftdc[1]=fCDSI->GetTofTDC1();
00161   snarldata->toftdc[2]=fCDSI->GetTofTDC2();
00162   snarldata->tofadc[0]=fCDSI->GetTofADC0();
00163   snarldata->tofadc[1]=fCDSI->GetTofADC1();
00164   snarldata->tofadc[2]=fCDSI->GetTofADC2();
00165   if(fCDSI->GetTofTimeStamp()!=0){
00166     snarldata->toftime=fCDSI->GetTofTimeStamp()
00167       *CalDetConstants::FARTIMECONVERT-snarldata->triggertime/Munits::nanosecond;
00168         MSG("UberDST",Msg::kDebug)<<"tof timestamp "<<fCDSI->GetTofTimeStamp()
00169                               <<" time convert "<<CalDetConstants::FARTIMECONVERT
00170                               <<" trigger time "<<snarldata->triggertime
00171                               <<" in ns "
00172                               <<snarldata->triggertime/Munits::nanosecond
00173                                   <<" tof time "<<snarldata->toftime<<endl;
00174 
00175   }
00176   if(fCDSI->GetTofADCTimeStamp0()!=0){
00177     snarldata->tofhittime[0]=fCDSI->GetTofADCTimeStamp0()
00178       *CalDetConstants::FARTIMECONVERT-snarldata->triggertime/Munits::nanosecond;
00179   }
00180   if(fCDSI->GetTofADCTimeStamp1()!=0){
00181     snarldata->tofhittime[1]=fCDSI->GetTofADCTimeStamp1()
00182       *CalDetConstants::FARTIMECONVERT-snarldata->triggertime/Munits::nanosecond;
00183   }
00184   if(fCDSI->GetTofADCTimeStamp2()!=0){
00185     snarldata->tofhittime[2]=fCDSI->GetTofADCTimeStamp2()
00186       *CalDetConstants::FARTIMECONVERT-snarldata->triggertime/Munits::nanosecond;
00187   }
00188   //get an iterator over the digits
00189   CandDigitHandleItr hi(fCDLH->GetDaughterIterator());
00190   if(!hi.IsValid()){
00191     MSG("UberDST",Msg::kError)<<"DaugherIter not valid"<<endl;
00192     return JobCResult(JobCResult::kFailed);
00193   }
00194 
00195   //Get a validity context and a plex handle
00196   const VldContext *vc = (*hi)->GetVldContext();
00197   //add the number of seconds since the beginning of the run to the
00198   //trigger time to make trigger time absolute since beginning of run
00199   if(fFirstEvent){
00200     header->SetStartTime(vc->GetTimeStamp().GetSec());
00201     fFirstEvent=kFALSE;
00202   }
00203 
00204   snarldata->triggertime+=vc->GetTimeStamp().GetSec()-header->GetStartTime();
00205   MSG("UberDST",Msg::kDebug)<<"Trigger time "<<snarldata->triggertime<<endl
00206                             <<"vld context time "
00207                             <<vc->GetTimeStamp().GetSec()<<endl;
00208 
00209   PlexHandle ph(*vc, kTRUE);
00210   if(!fMCDecision){
00211     if(vc->GetSimFlag()==SimFlag::kReroot||vc->GetSimFlag()==SimFlag::kMC){
00212       fIsMC = kTRUE;
00213     }
00214     fMCDecision = kTRUE;
00215   }
00216   //ReInit the mip calibrator
00217   Calibrator::Instance().Reset(*vc);
00218   Bool_t gotsigcorrconv = kFALSE;
00219 
00220 
00221   //check dead chip list in the CandCalDetSIHandle to find mindeadplane
00222   //and and number of dead planes
00223   set<UShort_t> deadplanemap;
00224   vector<RawChannelId>::const_iterator dcit(fCDSI->GetDeadChips().begin());
00225   MSG("UberDST", Msg::kDebug)<<"Starting to get dead chip vector"<<endl
00226                                <<"Size of Dead Chip vector "
00227                                <<fCDSI->GetDeadChips().size()<<endl;
00228   while(dcit!=fCDSI->GetDeadChips().end()){
00229     //make a rawchannel id out of the dead chip (use va channel 3)
00230     RawChannelId newid((*dcit).GetDetector(), 
00231                        (*dcit).GetElecType(), 
00232                        (*dcit).GetCrate(), 
00233                        (*dcit).GetVarcId(), 
00234                        (*dcit).GetVmm(), 
00235                        (*dcit).GetVaAdcSel(), 
00236                        (*dcit).GetVaChip(),
00237                        3);
00238     if(ph.GetSEIdAltL(newid).GetPlane()>=0){
00239       deadplanemap.insert(ph.GetSEIdAltL(newid).GetPlane());
00240     }
00241     if(newid.IsSameVAChip(fCDSI->GetCerenkovChannel3())){
00242       deadplanemap.insert(65);
00243     }
00244     if(newid.IsSameVAChip(fCDSI->GetCerenkovChannel1())){
00245       deadplanemap.insert(66);
00246     }
00247     if(newid.IsSameVAChip(fCDSI->GetCerenkovChannel2())){
00248       deadplanemap.insert(67);
00249     }
00250     //make another possible rawchannel id out of the dead chip(use va channel 11)
00251     //which could possibly be on the next plane
00252     RawChannelId newid2((*dcit).GetDetector(), 
00253                         (*dcit).GetElecType(), 
00254                         (*dcit).GetCrate(), 
00255                         (*dcit).GetVarcId(), 
00256                         (*dcit).GetVmm(), 
00257                         (*dcit).GetVaAdcSel(), 
00258                         (*dcit).GetVaChip(), 
00259                         11);
00260     if(ph.GetSEIdAltL(newid2).GetPlane()>=0){
00261       deadplanemap.insert(ph.GetSEIdAltL(newid2).GetPlane());
00262     }
00263     dcit++;
00264   }
00265   if(deadplanemap.size()!=0){
00266     snarldata->mindeadplaneno = *(deadplanemap.begin());
00267     snarldata->ndeadplanes = deadplanemap.size();
00268   }
00269   MSG("UberDST",Msg::kDebug)<<"First dead plane "<<snarldata->mindeadplaneno
00270                             <<" No. dead planes "<<snarldata->ndeadplanes<<endl;
00271   
00272 
00273 
00274   //sort hits in event by plane
00275   CandDigitHandleKeyFunc* kf = hi.CreateKeyFunc();
00276   kf->SetFun(CalHelpers::KeyFromPlane);
00277   hi.GetSet()->AdoptSortKeyFunc(kf,kTRUE,kFALSE);
00278   //further sort hits by strip
00279   kf=hi.CreateKeyFunc();
00280   kf->SetFun(CalHelpers::KeyFromStrip);
00281   hi.GetSet()->AdoptSortKeyFunc(kf,kFALSE,kFALSE);
00282   //and finally, sort by strip end
00283   kf=hi.CreateKeyFunc();
00284   kf->SetFun(CalHelpers::KeyFromEnd);
00285   hi.GetSet()->AdoptSortKeyFunc(kf,kFALSE,kTRUE);
00286   kf = 0;
00287   //MSG("UberDST",Msg::kDebug)<<"Done sorting"<<endl;
00288   
00289   //declare variables to find maxes
00290   Float_t maxplanemip = 0.;
00291   Float_t maxp0stripmip=0.;
00292   Float_t maxp1stripmip=0.;
00293   Float_t maxp0hitmip=0.;
00294   Float_t maxp1hitmip=0.;
00295   Float_t totmipeven=0.;
00296   Float_t totmipodd=0.;
00297 
00298   map<PlexStripEndId, REROOT_FLSDigit> FLSDigitMap;
00299   int nmchitstrips=0;
00300   if(fIsMC){
00301     const TClonesArray *hep = RerootExodus::GetStdHepList();
00302     mc->SetStdHep(hep);
00303     MSG("UberDST",Msg::kDebug)<<"size of hep is "<<hep->GetEntries()<<endl;
00304     REROOT_Event *revt = gMINFast->GetREROOTEvent();
00305     mc->nhep = revt->n_stdheps();
00306     MSG("UberDST",Msg::kDebug)<<"and mc nhep is  "<<mc->nhep<<endl;
00307     REROOT_StdHep *stdhep = static_cast<REROOT_StdHep *>(hep->At(0));
00308     if(stdhep!=0){
00309       mc->mctype = stdhep->ID();
00310       mc->mcpx = stdhep->Px();
00311       mc->mcpy = stdhep->Py();
00312       mc->mcpz = stdhep->Pz();
00313       mc->mcenergy = stdhep->E();
00314       mc->mcvx = stdhep->Xmm();
00315       mc->mcvy = stdhep->Ymm();
00316       mc->mcvz = stdhep->Zmm();
00317     }
00318     const TClonesArray *flsd = RerootExodus::GetFLSDigitList();
00319     int i=0;
00320     TObject *obj;
00321     while((obj=flsd->At(i))){
00322       REROOT_FLSDigit *mcdig = static_cast<REROOT_FLSDigit*>(obj);
00323       PlexStripEndId mcseid = RerootExodus::PECAB2SEId(mcdig->IPln(),
00324                                                        mcdig->IExtr(),
00325                                                        mcdig->ICell(),0);
00326       FLSDigitMap[mcseid]=*mcdig;
00327       i++;
00328     }
00329   }
00330 
00331     
00332 
00333   
00334   //begin looping over planes
00335   for(hi.Reset();hi.IsValid();hi.NextKey()){ 
00336     UShort_t plane = (*hi)->GetPlexSEIdAltL().GetPlane();
00337     Bool_t iscosmic = kFALSE;
00338     if(plane>=CalDetConstants::PLANECONST){ //if not a detector plane, continue
00339       iscosmic=kTRUE;
00340     }
00341     if(!iscosmic){
00342       snarldata->nhitplanes++; //increment nhitplanes
00343     }
00344     Float_t planemip=0.;
00345 
00346     //begin looping over strips
00347     for(CandDigitHandleItr siter(hi, kTRUE);
00348         siter.IsValid();siter.NextKey()){ 
00349       UShort_t strip = (*siter)->GetPlexSEIdAltL().GetBestSEId().GetStrip();
00350 
00351       if(!iscosmic){
00352         snarldata->nhitstrips++;  //increment nhitstrips
00353         snarldata->AddNextHit(plane,strip);//make a hit in UberEvent's TClonesArray
00354         if(fIsMC){
00355           mc->AddNextHit(plane,strip);
00356         }
00357       }
00358       if(iscosmic){
00359         snarldata->AddNextCosmicHit(plane,strip);
00360       }
00361       Float_t p0stripmip=0.;
00362       Float_t p1stripmip=0.;
00363 
00364       //begin looping over ends
00365       for(CandDigitHandleItr eiter(siter,kTRUE);
00366           eiter.IsValid();eiter.Next()){ //loop over ends
00367 
00368         /*
00369         MSG("UberDST",Msg::kDebug)<<"RAW CHANNEL ID FOR THIS TIME ROUND:"<<endl;
00370         (*eiter)->GetChannelId().Print();
00371         MSG("UberDST",Msg::kDebug)<<endl;
00372         (*eiter)->GetPlexSEIdAltL().Print();
00373         MSG("UberDST",Msg::kDebug)<<" readouttype "
00374                                   <<ph.GetReadoutType((*eiter)->GetChannelId())<<endl;
00375         */
00376 
00377         if((*eiter)->GetPlexSEIdAltL().size()<=0){
00378           MSG("UberDST",Msg::kDebug)<<"plexseidaltl is EMPTY"<<endl;
00379           //(*eiter)->GetChannelId().Print();
00380           //MSG("UberDST",Msg::kDebug)<<endl;
00381           //(*eiter)->GetPlexSEIdAltL().Print();
00382           //MSG("UberDST",Msg::kDebug)<<endl;
00383           continue;
00384         }
00385         /*
00386         //for debugging the reading of the near detector stuff:
00387         if((*eiter)->GetChannelId().GetElecType()==ElecType::kQIE){
00388           MSG("UberDST",Msg::kDebug)<<" FOUND A QIE DIGIT "<<endl;
00389         }
00390         */
00391         PlexStripEndId pse = (*eiter)->GetPlexSEIdAltL().GetBestSEId();
00392         StripEnd::StripEnd_t se = pse.GetEnd();
00393         PlexSEIdAltLItem pseitem = (*eiter)->GetPlexSEIdAltL().GetBestItem();
00394 
00395         //get hit info
00396         Int_t adc = (int)(*eiter)->GetCharge();
00397         Float_t siglin = (*eiter)->GetCharge(CalDigitType::kSigLin);
00398         Float_t npe = pseitem.GetPE();
00399         Float_t sigcorr = (*eiter)->GetCharge(CalDigitType::kSigCorr);
00400         Float_t mip = Calibrator::Instance().GetMIP((*eiter)->GetCharge(CalDigitType::kSigCorr),
00401                                                     pse);
00402         Float_t time = (*eiter)->GetTime(CalTimeType::kT0)/Munits::nanosecond;
00403         PlexPixelSpotId pps=ph.GetPixelSpotId(pse);
00404         Int_t agg = pps.GetEncoded();
00405         // added by mak March 25, 2003
00406         if(fIsMC){
00407           // make mips right for the mc
00408           npe*=fMC_PE_Conversion;
00409           mip*=fMC_MIP_Conversion;
00410         }
00411 
00412         if(iscosmic){
00413           snarldata->AddNextCosmicHitValues(se, adc, time, agg);
00414           continue;
00415         }
00416 
00417         snarldata->nhits++;  //increment nhits
00418         if(!gotsigcorrconv){
00419           snarldata->sigcorrconv = sigcorr/mip;
00420           gotsigcorrconv=kTRUE;
00421         }
00422 
00423         //add the hit info to UberEvent's TClonesArray
00424 //CandDigit now removes the pedestal offset!!!!!, no need to do anything special
00425 //for qie digits anymore
00426 /*
00427         if((*eiter)->GetChannelId().GetElecType()==ElecType::kQIE){
00428           float Nd = 1.*adc/siglin;
00429           float Nc = 1.*siglin/mip;
00430           snarldata->AddNextHitValues(se, adc-50, siglin, 
00431                                       (adc-50)*1./(1.*adc/npe), (adc-50)/(Nd*Nc),
00432                                       time, agg);
00433         }
00434         else{
00435           snarldata->AddNextHitValues(se, adc, siglin, npe, mip, time, agg);
00436         }
00437 */
00438         snarldata->AddNextHitValues(se, adc, siglin, npe, mip, time, agg);
00439         MSG("UberDST",Msg::kDebug)<<"STRIP END "<<se<<" adc "<<adc<<endl;
00440         planemip+=mip;  //increment the mips deposited in current plane
00441 
00442         if(pse.GetPlane()%2==0){  //calculate center in even planes
00443           snarldata->mipweighcentereven+=mip*pse.GetStrip();
00444           totmipeven+=mip;
00445         }
00446         else{ //calculate center in odd planes
00447           snarldata->mipweighcenterodd+=mip*pse.GetStrip();
00448           totmipodd+=mip;
00449         }
00450 
00451         //look for max hit adc, npe, mip, time
00452         if(adc>snarldata->maxadc){
00453           snarldata->maxadc = adc;
00454         }
00455         if(npe>snarldata->maxnpe){
00456           snarldata->maxnpe = npe;
00457         }
00458         if(mip>snarldata->maxmip){
00459           snarldata->maxmip = mip;
00460         }
00461         if(time>snarldata->maxtime){
00462           snarldata->maxtime = time;
00463         }
00464 
00465         //add mip in all planes but plane 0
00466         if(pse.GetPlane()!=0){
00467           snarldata->totmip+=mip;
00468         }
00469 
00470         //find mips in plane 0
00471         if(pse.GetPlane()==0){
00472           snarldata->p0totmip+=mip;
00473           p0stripmip+=mip;
00474           if(mip>maxp0hitmip){
00475             maxp0hitmip = mip;
00476             snarldata->p0maxmiptstamp = time;
00477           } 
00478         }
00479 
00480         //find mips in plane 1
00481         if(pse.GetPlane()==1){
00482           snarldata->p1totmip+=mip;
00483           p1stripmip+=mip;
00484           if(mip>maxp1hitmip){
00485             maxp1hitmip = mip;
00486             snarldata->p1maxmiptstamp = time;
00487           } 
00488         }
00489       }//end loop over ends
00490       
00491       if(fIsMC&&!iscosmic){
00492         PlexStripEndId mcpse = (*siter)->GetPlexSEIdAltL().GetBestSEId();
00493         mcpse.SetEnd(StripEnd::kNegative);
00494         REROOT_FLSDigit mcd = FLSDigitMap.find(mcpse)->second;
00495         MSG("UberDST",Msg::kDebug)<<"pse.GetPlane() "
00496                                   <<(*siter)
00497           ->GetPlexSEIdAltL().GetBestSEId().GetPlane()
00498                                   <<" mc plane "<<mcd.IPln()<<endl;
00499         mc->AddNextHitValues(mcd.HitBits(),mcd.SumETrue(),mcd.TPos(),mcd.RawB(),
00500                              mcd.RawA(), mcd.CorrB(), mcd.CorrA(), mcd.TDCB(),
00501                              mcd.TDCA(), mcd.AveDistTrueB(),mcd.AveDistTrueA());
00502         mc->mcenergydep += mcd.SumETrue();
00503         nmchitstrips++;
00504       }
00505       
00506       //find max mip in an strip in plane 0
00507       if(p0stripmip>maxp0stripmip){ 
00508         maxp0stripmip = p0stripmip;
00509         snarldata->p0stripmaxmip = strip;
00510       }
00511       //find max mip in a strip in plane 1
00512       if(p1stripmip>maxp1stripmip){
00513         maxp1stripmip = p1stripmip;
00514         snarldata->p1stripmaxmip = strip;
00515       }
00516     }//end loop over strips
00517 
00518     snarldata->mipweighaveplane+=planemip*plane;  //calculate mip weigh ave plane
00519     if(planemip>maxplanemip){  //find showermax
00520       maxplanemip = planemip;
00521       snarldata->mipshowermax = maxplanemip;
00522       snarldata->showermax = plane;
00523     }
00524   }//end loop over planes
00525 
00526   snarldata->mipweighaveplane/=(snarldata->totmip+snarldata->p0totmip);
00527   if(totmipeven!=0){
00528      snarldata->mipweighcentereven/=totmipeven;
00529   }
00530   if(totmipodd!=0){
00531      snarldata->mipweighcenterodd/=totmipodd;
00532   }
00533 
00534   //loop again to calculate radius of event:
00535   for(hi.Reset();hi.IsValid();hi.NextKey()){
00536     UShort_t plane = ((*hi)->GetPlexSEIdAltL()).GetPlane();
00537     if(plane>=CalDetConstants::PLANECONST){
00538       continue;
00539     }
00540     float center = 0.;
00541     if(plane%2==0){
00542       center=snarldata->mipweighcentereven;
00543     }
00544     else{
00545       center = snarldata->mipweighcenterodd;
00546     }
00547     for(CandDigitHandleItr siter(hi, kTRUE);
00548         siter.IsValid();siter.NextKey()){ //loop over strips
00549       for(CandDigitHandleItr eiter(siter,kTRUE);
00550           eiter.IsValid();eiter.Next()){ //loop over ends
00551         PlexStripEndId pse = (*eiter)->GetPlexSEIdAltL().GetBestSEId();
00552         if(pse.GetStrip()>=CalDetConstants::STRIPCONST){
00553           continue;
00554         }
00555         Float_t mip = Calibrator::Instance().GetMIP((*eiter)->GetCharge(CalDigitType::kSigCorr),
00556                                                     pse);
00557         if(fIsMC){
00558           mip*=fMC_MIP_Conversion;
00559         }
00560         snarldata->mipweighrad += fabs(center-pse.GetStrip())*mip;
00561       }//end end loop
00562     }//end strip loop
00563   }//end planeloop
00564   snarldata->mipweighrad/=(snarldata->totmip+snarldata->p0totmip);
00565 
00566   if(!fFilledHeader){
00567     header->SetBeamMomentum(header->ReadBeamMomentum());
00568     header->SetComments(header->ReadComments());
00569     header->SetTemperature(header->ReadTemperature());
00570     header->SetRunType(header->ReadRunType());
00571   }
00572 
00573   if(!fMadeTrees){
00574     MakeTrees();
00575   }
00576 
00577   snarltree->Fill();
00578   if(!fFilledHeader){
00579     headertree->Fill();
00580     fFilledHeader=kTRUE;
00581   }
00582   if(fIsMC){
00583     mctree->Fill();
00584   }
00585   MSG("UberDST",Msg::kDebug)<<"Size of mc map: "<<FLSDigitMap.size()
00586                             <<" Number of hitstrips "<<snarldata->nhitstrips
00587                             <<" Nuber of mchitstrips "<<nmchitstrips<<endl;
00588   return JobCResult(JobCResult::kPassed);
00589 }//end Ana()
00590 //________________________________________________________________________________
00591 void UberDST::EndJob()
00592 {
00593 
00594   MSG("UberDST",Msg::kDebug)<<"In end job "
00595                             <<"comments "<<header->GetComments()<<endl;
00596   fRootFile->cd();
00597   headertree->Write();
00598   snarltree->Write();
00599   if(fIsMC){
00600     mctree->Write();
00601   }
00602   fRootFile->Close();
00603 
00604 }//end EndJob()
00605 //________________________________________________________________________________
00606 JobCResult UberDST::MakeTrees()
00607 {
00608   MSG("UberDST",Msg::kDebug)<<"In UberDST MakeTrees"<<endl;
00609   //try to open a file with name UberDST-runnumber.root
00610   string basedir;
00611   char *bd = getenv("UBERDST_OUT_DIR");
00612   if(bd==NULL){
00613   basedir ="./";
00614   }
00615   else{
00616     basedir=bd;
00617   }
00618   MSG("UberDST",Msg::kDebug)<<"UBERDSTPATH "<<basedir<<" bd "<<bd<<endl;
00619   stringstream rf;
00620   rf<<basedir<<"/UberDST-"<<fRunNumber<<".root"<<'\0';
00621   string srf = rf.str();
00622   //  rf.freeze(false);
00623   fRootFile = new TFile(srf.c_str(),"CREATE");
00624   //if a file with that name already exists, try appending a cycle number
00625   //if there are more than 10 files with this name, give up
00626   int i=1;
00627   while(!fRootFile->IsOpen()&&i<10){
00628     fRootFile = 0;
00629     stringstream rf2;
00630     rf2<<basedir<<"UberDST-"<<fRunNumber<<"-"<<i<<".root"<<'\0';
00631     string srf2 = rf2.str();
00632     //    rf2.freeze(false);
00633     fRootFile = new TFile(srf2.c_str(),"CREATE");
00634     i++;
00635   }
00636 
00637   //if we still don't have a file open, return kError
00638   if(!fRootFile->IsOpen()){
00639     MSG("UberDST",Msg::kWarning)<<"Could not open "<<srf<<" for writing"<<endl;
00640     return JobCResult(JobCResult::kFatal);
00641   }
00642   MSG("UberDST",Msg::kDebug)<<"About to make trees"<<endl;
00643   snarltree = new TTree("snarltree","snarltree");
00644   headertree = new TTree("headertree","headertree");
00645 
00646   MSG("UberDST",Msg::kDebug)<<"About to make branches"<<endl;
00647   snarltree->Branch("snarldata","UberEvent",&snarldata);
00648   headertree->Branch("header","UberHeader",&header);
00649 
00650   MSG("UberDST",Msg::kDebug)<<"Setting Directories"<<endl;
00651   snarltree->SetDirectory(fRootFile);
00652   headertree->SetDirectory(fRootFile);
00653 
00654   if(fIsMC){
00655     mctree = new TTree("mctree","mctree");
00656     mctree->Branch("mctree","UberMC",&mc);
00657     mctree->SetDirectory(fRootFile);
00658   }
00659 
00660   fMadeTrees = kTRUE;
00661   MSG("UberDST",Msg::kDebug)<<"Returning"<<endl;
00662   return JobCResult(JobCResult::kAOK);
00663 }//end MakeTrees()
00664 //_______________________________________________________________________________
00665 
00666 
00667 
00668 
00669 
00670 
00671 
00672 
00673 
00674 
00675 
00676