loon: AlgAtNuRecoMCTruth.cxx Source File

00001 // AlgAtNuRecoMCTruth.cxx
00002 
00003 #include "AlgAtNuRecoMCTruth.h"
00004 
00005 #include "Algorithm/AlgConfig.h"
00006 #include "Algorithm/AlgFactory.h"
00007 #include "Algorithm/AlgHandle.h"
00008 
00009 #include "RerootExodus/RerootExodus.h"
00010 #include "REROOT_Classes/REROOT_NeuKin.h"
00011 #include "REROOT_Classes/REROOT_FLSHit.h"
00012 
00013 #include "Record/SimSnarlRecord.h"
00014 #include "Record/SimSnarlHeader.h"
00015 #include "Digitization/DigiSignal.h"
00016 
00017 #include "CandData/CandRecord.h"
00018 #include "Candidate/CandContext.h"
00019 #include "JobControl/JobCModuleRegistry.h"
00020 #include "JobControl/JobCommand.h"
00021 #include "MessageService/MsgService.h"
00022 #include "MinosObjectMap/MomNavigator.h"
00023 #include "UgliGeometry/UgliGeomHandle.h"
00024 #include "Validity/VldContext.h"
00025 #include "Plex/PlexPlaneId.h"
00026 
00027 #include "RecoBase/CandSliceListHandle.h"
00028 #include "RecoBase/CandSliceHandle.h"
00029 #include "RecoBase/CandStripListHandle.h"
00030 #include "RecoBase/CandStripHandle.h"
00031 #include "CandDigit/CandDeMuxDigitHandle.h"
00032 #include "CandDigit/CandDeMuxDigitListHandle.h"
00033 
00034 #include "CandShowerAtNuHandle.h"
00035 #include "CandShowerAtNuListHandle.h"
00036 #include "CandTrackAtNuHandle.h"
00037 #include "CandTrackAtNuListHandle.h"
00038 
00039 #include "HitAtNu.h"
00040 #include "TrackAtNu.h"
00041 #include "ShowerAtNu.h"
00042 
00043 #include "CandAtNuRecoHandle.h"
00044 
00045 #include "TClonesArray.h"
00046 
00047 #include <sys/time.h>
00048 
00049 //
00050 // $Log $
00051 //
00052 
00053 ClassImp(AlgAtNuRecoMCTruth)
00054 
00055 CVSID("$Id: AlgAtNuRecoMCTruth.cxx,v 1.7 2006/03/28 13:54:11 blake Exp $");
00056 
00057 AlgAtNuRecoMCTruth::AlgAtNuRecoMCTruth()
00058 {
00059 
00060 }
00061   
00062 AlgAtNuRecoMCTruth::~AlgAtNuRecoMCTruth()
00063 {
00064 
00065 }
00066 
00067 void AlgAtNuRecoMCTruth::RunAlg(AlgConfig& ac, CandHandle &ch, CandContext &cx)
00068 {
00069 
00070   MSG("AlgAtNuRecoMCTruth", Msg::kDebug) << "AlgAtNuRecoMCTruth::RunAlg(...)" << endl;
00071 
00072   // Get CandAtNuRecoHandle
00073   CandAtNuRecoHandle& atnu_handle = dynamic_cast<CandAtNuRecoHandle&>(ch);
00074 
00075   // Unpack AlgConfig
00076   TString fListOutTrk,fListOutShw;
00077   fListOutTrk = ac.GetCharString("ListOutTrk");
00078   fListOutShw = ac.GetCharString("ListOutShw");
00079  
00080   // Unpack CandContext
00081   CandSliceListHandle* slice_list = (CandSliceListHandle*)(cx.GetDataIn()); 
00082 
00083   // clear timebanks
00084   Double_t dtime;
00085   Double_t fTime_TrackAtNu,fTime_ShowerAtNu;
00086 
00087   struct timeval tpbefore;
00088   gettimeofday(&tpbefore,0);
00089 
00090   Int_t i,j,k;
00091   Int_t pln,bpln,epln;
00092   Int_t trkplns,shwplns;
00093   TObjArray tmptrk,tmpshw,tmpmu;
00094   Int_t itr,idnu,ipdg,muflag,trkflag,shwflag;
00095   Int_t trueview,packedPEC,trueplane,truecell,trueext,truestrip;
00096 
00097   // Get Geometry
00098   CandRecord* candrectmp = (CandRecord*)(cx.GetCandRecord()); 
00099   VldContext *vldc = (VldContext*)(candrectmp->GetVldContext());
00100   UgliGeomHandle ugh(*vldc);
00101 
00102   Int_t StrpExtr[8]={ 0, 28, 56, 76, 96, 116, 136, 164 };
00103   Int_t StrpCell[28]={ 0, 1, 2, 3, 4, 5, 6, 
00104                        7, 8, 9, 10, 11, 12, 13, 
00105                        14, 15, 16, 17, 18, 19, 20, 
00106                        21, 22, 23, 24, 25, 26, 27 };
00107 
00108   // Get FLSHits
00109   const TClonesArray* FLShits = (TClonesArray*)(RerootExodus::GetFLSHitList());
00110   const TClonesArray* NeuKins = (TClonesArray*)(RerootExodus::GetNeuKinList());
00111   REROOT_NeuKin* neukin = (REROOT_NeuKin*)(NeuKins->At(0));
00112   idnu = neukin->INu();
00113   for(i=0;i<1+FLShits->GetLast();i++){
00114     REROOT_FLSHit* FLS_hit = (REROOT_FLSHit*)FLShits->At(i);
00115     ipdg = FLS_hit->IPDG();                   
00116     if( (idnu==0&&(ipdg==13||ipdg==-13))
00117       ||(idnu==14&&ipdg==13)||(idnu==-14&&ipdg==-13) ){
00118       tmpmu.Add(FLS_hit);
00119     }
00120   }
00121 
00122   // Iterate over slices
00123   TIter sliceitr(slice_list->GetDaughterIterator());
00124   while(CandSliceHandle* slice = dynamic_cast<CandSliceHandle*>(sliceitr())){
00125     atnu_handle.AddDaughterLink(*slice);
00126 
00127     TIter stripitr(slice->GetDaughterIterator());
00128     while(CandStripHandle* strip = (CandStripHandle*)(stripitr())){
00129       HitAtNu* hit = new HitAtNu(strip); pln=hit->GetPlane(); 
00130       hitplnbnk[pln].Add(hit); hitbnk.Add(hit);
00131     }
00132 
00133     trkplns=0; shwplns=0;
00134     for(i=0;i<500;i++){
00135       trkflag=0; shwflag=0;
00136       for(j=0;j<1+hitplnbnk[i].GetLast();j++){
00137         HitAtNu* hit = (HitAtNu*)(hitplnbnk[i].At(j));
00138         if(hit->GetCharge()>1.0){
00139           muflag=0;
00140           for(k=0;k<1+tmpmu.GetLast();k++){
00141             if(!muflag){
00142               REROOT_FLSHit* FLS_hit = (REROOT_FLSHit*)(tmpmu.At(k));
00143               packedPEC=FLS_hit->IPackedPEC();
00144               trueplane=(Int_t)(packedPEC/65536.0);
00145               if(trueplane==hit->GetPlane()){
00146                 trueext=(Int_t)((packedPEC-trueplane*65536)/256.0);
00147                 truecell=packedPEC-trueplane*65536-trueext*256;
00148                 truestrip=StrpExtr[trueext-1]+StrpCell[truecell-1];
00149                 PlexPlaneId trueplnid(vldc->GetDetector(),trueplane,0);
00150                 UgliScintPlnHandle trueplnhandle = ugh.GetScintPlnHandle(trueplnid);
00151                 if(trueplnhandle.GetPlaneView()==PlaneView::kU) trueview = 0;
00152                 if(trueplnhandle.GetPlaneView()==PlaneView::kV) trueview = 1;   
00153                 if(truestrip==hit->GetStrip()) muflag=1;
00154               }
00155             }
00156           }
00157           if(muflag){
00158             trkflag=1; tmptrk.Add(hit); hit->SetTrkFlag(2);
00159           }
00160           else{
00161             shwflag=1; tmpshw.Add(hit); hit->SetShwFlag(2);
00162           }
00163         }
00164       }
00165       if(trkflag) trkplns++; if(shwflag) shwplns++;
00166     }
00167 
00168 
00169     if(1+tmptrk.GetLast()>0){
00170       bpln=-1; epln=-1;
00171       for(j=0;j<1+tmptrk.GetLast();j++){
00172         HitAtNu* hit = (HitAtNu*)(tmptrk.At(j));
00173         if(bpln<0||hit->GetPlane()<bpln) bpln=hit->GetPlane();
00174         if(epln<0||hit->GetPlane()>epln) epln=hit->GetPlane();
00175       }
00176 
00177       if(trkplns>5 && 1+epln-bpln>5){
00178         TrackAtNu* trku = new TrackAtNu(slice);
00179         TrackAtNu* trkv = new TrackAtNu(slice);
00180         trkbnk[0].Add(trku); trkbnk[1].Add(trkv);
00181 
00182         for(j=0;j<1+tmptrk.GetLast();j++){
00183           HitAtNu* hit = (HitAtNu*)(tmptrk.At(j));
00184           if(hit->GetPlaneView()==0) trku->AddHit(hit);
00185           if(hit->GetPlaneView()==1) trkv->AddHit(hit);
00186         }
00187         
00188         trku->SetPartner(trkv); trkv->SetPartner(trku);
00189       }
00190 
00191     }
00192     
00193 
00194     if(1+tmpshw.GetLast()>0){
00195       bpln=-1; epln=-1;
00196       for(j=0;j<1+tmpshw.GetLast();j++){
00197         HitAtNu* hit = (HitAtNu*)(tmpshw.At(j));
00198         if(bpln<0||hit->GetPlane()<bpln) bpln=hit->GetPlane();
00199         if(epln<0||hit->GetPlane()>epln) epln=hit->GetPlane();
00200       }
00201 
00202       if(shwplns>1 && 1+epln-bpln>1){
00203         ShowerAtNu* shwu = new ShowerAtNu(slice);
00204         ShowerAtNu* shwv = new ShowerAtNu(slice);
00205         shwbnk[0].Add(shwu); shwbnk[1].Add(shwv);
00206 
00207         for(j=0;j<1+tmpshw.GetLast();j++){
00208           HitAtNu* hit = (HitAtNu*)(tmpshw.At(j));
00209           if(hit->GetPlaneView()==0) shwu->AddHit(hit);
00210           if(hit->GetPlaneView()==1) shwv->AddHit(hit);
00211         }
00212         
00213         shwu->SetPartner(shwv); shwv->SetPartner(shwu);
00214       }
00215 
00216     }
00217 
00218 
00219     for(i=0;i<500;i++){
00220       hitplnbnk[i].Clear();
00221     }
00222 
00223   }
00224 
00225   
00226   
00227   // *********************************************************
00228   // * F O R M A T I O N   O F   C A N D X X X H A N D L E S *
00229   // *********************************************************
00230   MSG("AlgAtNuRecoMCTruth", Msg::kDebug) << " Formation of CandXXXHandles " << endl;
00231 
00232   // create Singleton Instance of AlgFactory
00233   AlgFactory &af = AlgFactory::GetInstance();
00234 
00235   // create CandTrackAtNuListHandle
00236   TObjArray* objtrk = 0;
00237   for(itr=0;itr<1+trkbnk[0].GetLast();itr++){
00238     TrackAtNu* trku = (TrackAtNu*)(trkbnk[0].At(itr)); 
00239     TrackAtNu* trkv = (TrackAtNu*)(trkbnk[1].At(itr));
00240     objtrk = new TObjArray(); 
00241     objtrk->Add(trku);
00242     objtrk->Add(trkv); 
00243     trkcxt.Add(objtrk);
00244   }
00245 
00246   MSG("AlgAtNuRecoMCTruth", Msg::kDebug) << "*** CREATE CandTrackAtNuListHandle *** " << endl;
00247   AlgHandle ah_trk = af.GetAlgHandle("AlgTrackAtNuList", "default");
00248   CandContext cx_trk(this, cx.GetMom());
00249   cx_trk.SetCandRecord(cx.GetCandRecord());
00250   cx_trk.SetDataIn(&trkcxt);
00251   CandTrackAtNuListHandle trk_list = CandTrackAtNuList::MakeCandidate(ah_trk, cx_trk);
00252   trk_list.SetName(fListOutTrk.Data());
00253   trk_list.SetTitle(TString("Created by AtNuFindModule from ").Append(slice_list->GetName()));
00254   
00255   TIter trkitr(trk_list.GetDaughterIterator());
00256   while(CandTrackAtNuHandle* trk = dynamic_cast<CandTrackAtNuHandle*>(trkitr())){
00257     atnu_handle.AddTrack(trk);
00258   }
00259 
00260   //  create CandShowerAtNuListHandle
00261   TObjArray* objshw = 0;
00262   for(itr=0;itr<1+shwbnk[0].GetLast();itr++){
00263     ShowerAtNu* shwu = (ShowerAtNu*)(shwbnk[0].At(itr));
00264     ShowerAtNu* shwv = (ShowerAtNu*)(shwbnk[1].At(itr));
00265     objshw = new TObjArray();
00266     objshw->Add(shwu);
00267     objshw->Add(shwv);
00268     shwcxt.Add(objshw);
00269   }
00270 
00271   MSG("AlgAtNuRecoMCTruth", Msg::kDebug) << "*** CREATE CandShowerAtNuListHandle *** " << endl;
00272   AlgHandle ah_shw = af.GetAlgHandle("AlgShowerAtNuList", "default");
00273   CandContext cx_shw(this, cx.GetMom());
00274   cx_shw.SetCandRecord(cx.GetCandRecord());
00275   cx_shw.SetDataIn(&shwcxt);
00276   CandShowerAtNuListHandle shw_list = CandShowerAtNuList::MakeCandidate(ah_shw, cx_shw);
00277   shw_list.SetName(fListOutShw.Data());
00278   shw_list.SetTitle(TString("Created by AtNuFindModule from ").Append(slice_list->GetName())); 
00279   
00280   TIter shwitr(shw_list.GetDaughterIterator());
00281   while(CandShowerAtNuHandle* shw = dynamic_cast<CandShowerAtNuHandle*>(shwitr())){
00282     atnu_handle.AddShower(shw);
00283   }
00284 
00285 
00286   // ***************************************************
00287   // * D E L E T E   R E M A I N I N G   O B J E C T S *
00288   // ***************************************************
00289 
00290   hitbnk.Delete(); // hit bank
00291 
00292   for(itr=0;itr<2;itr++){ 
00293     trkbnk[itr].Delete(); // track bank
00294   }
00295   trkcxt.Delete();      // candcontext for tracks
00296 
00297   for(itr=0;itr<2;itr++){ 
00298     shwbnk[itr].Delete(); // shower bank
00299   }
00300   shwcxt.Delete();      // candcontext for showers
00301 
00302   // *********************************
00303   // * C A L C U L A T E   T I M E S *
00304   // *********************************
00305 
00306   struct timeval tpafter;
00307   gettimeofday(&tpafter,0);
00308   dtime=1000.0*(tpafter.tv_sec-tpbefore.tv_sec)+0.001*(tpafter.tv_usec-tpbefore.tv_usec);
00309   fTime_TrackAtNu=0.5*dtime; fTime_ShowerAtNu=0.5*dtime; 
00310 
00311   trk_list.SetCPUTime(fTime_TrackAtNu);
00312   MSG("AlgAtNuRecoMCTruth", Msg::kDebug) << " *** TRACK RECO TIME = " << fTime_TrackAtNu << endl; 
00313 
00314   shw_list.SetCPUTime(fTime_ShowerAtNu);
00315   MSG("AlgAtNuRecoMCTruth", Msg::kDebug) << " *** SHOWER RECO TIME = " << fTime_ShowerAtNu << endl;
00316 
00317   // ***************************************************************
00318   // * G I V E   C A N D H A N D L E S   T O   C A N D R E C O R D *
00319   // ***************************************************************
00320 
00321   CandRecord* candrec = (CandRecord*)(cx.GetCandRecord());
00322   candrec->SecureCandHandle(trk_list);
00323   candrec->SecureCandHandle(shw_list);
00324 
00325   return;
00326 
00327 }
00328 
00329 void AlgAtNuRecoMCTruth::Trace(const char * /* c */) const
00330 {
00331 
00332 }
00333 
00334 
00335 
00336