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