loon: ReweightHelpers.cxx Source File

00001 #include "TClonesArray.h"
00002 #include "MCNtuple/NtpMCRecord.h"
00003 #include "MCNtuple/NtpMCTruth.h"
00004 #include "MCNtuple/NtpMCStdHep.h"
00005 #include "MessageService/MsgService.h"
00006 #include "StandardNtuple/NtpStRecord.h"
00007 #include "AnalysisNtuples/ANtpTruthInfoBeam.h"
00008 #include "MCReweight/ReweightHelpers.h"
00009 #include "MCReweight/MCEventInfo.h"
00010 
00011 CVSID("$Id: ReweightHelpers.cxx,v 1.3 2007/03/01 17:13:02 rhatcher Exp $");
00012 //......................................................................
00013 
00014 
00015 void ReweightHelpers::EventRegistryFilla(NtpMCRecord *mc, int evtno, Registry &fReg)
00016 {
00017    if(evtno>=mc->mc->GetEntries()){
00018       MSG("ReweightHelper",Msg::kError)<<"EventRegistryFilla Error, event no: "<<evtno
00019                                        <<" size "<<mc->mc->GetEntries()<<std::endl;
00020       return;
00021    }
00022 
00023    NtpMCTruth *mcth = static_cast<NtpMCTruth *>((*mc->mc)[evtno]);
00024 
00025    Int_t initial_state=0;  
00026    Int_t had_fs=0;
00027    bool foundfs=false;
00028    TClonesArray& heparray = *(mc->stdhep);
00029    Int_t nhep = heparray.GetEntries();
00030 
00031    int protneut = -1;  // 0 = neutron, 1 = proton, 2 = N, 3 = A
00032    int nubarnu = 0;    // +1 = neutrino, -1 = antineutrino
00033    
00034    for(int i=0;i<nhep;i++){    
00035       NtpMCStdHep *sh = static_cast<NtpMCStdHep *>(heparray[i]);
00036       if(sh->mc==evtno){
00037          if(sh->IstHEP==0){  //initial state particle
00038             if(abs(sh->IdHEP)==12 || 
00039                abs(sh->IdHEP)==14 || 
00040                abs(sh->IdHEP)==16){   //(anti)neutrino    
00041                nubarnu = sh->IdHEP/abs(sh->IdHEP);  //get sign
00042             }
00043          }
00044          if(sh->IstHEP==11){    //target nucleon
00045             if(sh->IdHEP==2212) protneut = 1;  //proton
00046             else if(sh->IdHEP==2112) protneut = 0;  //neutron
00047             else if(abs(sh->IdHEP)>1000000000) protneut = 2;  //nucleus
00048             else protneut = 3; //atom - probably never get here since IdHEP A>N?
00049          }
00050          if(sh->IstHEP==3&&!foundfs){
00051            if(mcth->iresonance!=1002 ||   //if not RES event
00052               TMath::Abs(sh->IdHEP)!=15){ //or if it is RES, IdHEP!=tau lepton
00053              had_fs = (sh->IdHEP%1000);
00054              foundfs=true;
00055            }
00056          }
00057       }
00058    }
00059    
00060    if(protneut==1 && nubarnu==1)  initial_state=1;  //p + v
00061    if(protneut==0 && nubarnu==1)  initial_state=2;  //n + v
00062    if(protneut==1 && nubarnu==-1) initial_state=3;  //p + vbar
00063    if(protneut==0 && nubarnu==-1) initial_state=4;  //n + vbar
00064    if(protneut==2 && nubarnu==1)  initial_state=5;  //N + v
00065    if(protneut==3 && nubarnu==1)  initial_state=6;  //A + v
00066    if(protneut==2 && nubarnu==-1) initial_state=7;  //N + vbar
00067    if(protneut==3 && nubarnu==-1) initial_state=8;  //A + vbar
00068 
00069    Int_t nucleus=FindNucleusNumber((int)mcth->z,(int)mcth->a);
00070 
00071    fReg.UnLockKeys();
00072    fReg.UnLockValues();
00073    
00074    fReg.Set("event:nu_en",1.*mcth->p4neu[3]);
00075    fReg.Set("event:nu_px",1.*mcth->p4neu[0]);
00076    fReg.Set("event:nu_py",1.*mcth->p4neu[1]);
00077    fReg.Set("event:nu_pz",1.*mcth->p4neu[2]);
00078 
00079    fReg.Set("event:tar_en",1.*mcth->p4tgt[3]);
00080    fReg.Set("event:tar_px",1.*mcth->p4tgt[0]);
00081    fReg.Set("event:tar_py",1.*mcth->p4tgt[1]);
00082    fReg.Set("event:tar_pz",1.*mcth->p4tgt[2]);
00083 
00084    fReg.Set("event:y",mcth->y);
00085    fReg.Set("event:x",mcth->x);
00086    fReg.Set("event:q2",mcth->q2);
00087    fReg.Set("event:w2",mcth->w2);
00088 
00089    fReg.Set("event:iaction",mcth->iaction);
00090    fReg.Set("event:inu",mcth->inu);
00091    fReg.Set("event:iresonance",mcth->iresonance);
00092    fReg.Set("event:initial_state",initial_state);
00093 
00094    fReg.Set("event:nucleus",nucleus);
00095 
00096    fReg.Set("event:had_fs",had_fs);
00097 
00098    fReg.LockValues();
00099    fReg.LockKeys();
00100 }
00101 
00102 
00103 void ReweightHelpers::EventRegistryFilla(ANtpTruthInfoBeam *ir, Registry &fReg)
00104 {
00105    fReg.UnLockKeys();
00106    fReg.UnLockValues();
00107    
00108    fReg.Set("event:nu_en",ir->nuEnergy);
00109    fReg.Set("event:nu_px",ir->nuDCosX*sqrt(ir->nuEnergy*ir->nuEnergy));
00110    fReg.Set("event:nu_py",ir->nuDCosY*sqrt(ir->nuEnergy*ir->nuEnergy));
00111    fReg.Set("event:nu_pz",ir->nuDCosZ*sqrt(ir->nuEnergy*ir->nuEnergy));
00112 
00113    fReg.Set("event:tar_en",ir->targetEnergy);
00114    fReg.Set("event:tar_px",ir->targetPX);
00115    fReg.Set("event:tar_py",ir->targetPY);
00116    fReg.Set("event:tar_pz",ir->targetPZ);
00117 
00118    fReg.Set("event:y",ir->hadronicY);
00119    fReg.Set("event:x",ir->bjorkenX);
00120    fReg.Set("event:q2",ir->q2);
00121    fReg.Set("event:w2",ir->w2);
00122 
00123    fReg.Set("event:iaction",ir->interactionType);
00124    fReg.Set("event:inu",ir->nuFlavor);
00125    fReg.Set("event:iresonance",ir->resonanceCode);
00126    fReg.Set("event:initial_state",ir->initialState);
00127 
00128    Int_t nucleus=FindNucleusNumber((int)ir->atomicNumber,
00129                                    (int)ir->atomicWeight);
00130    fReg.Set("event:nucleus",nucleus);
00131    fReg.Set("event:had_fs",ir->hadronicFinalState);
00132 
00133 
00134    fReg.LockValues();
00135    fReg.LockKeys();
00136 
00137 
00138 }
00139 
00140 void ReweightHelpers::BeamRegistryFilla(ANtpTruthInfoBeam *ir, Registry &fReg)
00141 {
00142   //note: this function assumes the parent hadron information 
00143   //in ANtpTruthInfoBeam
00144   //has already been filled 
00145   //however, this information does not get filled
00146   //automatically by calling the filler object for 
00147   //ANtpTruthInfoBeam, instead, it gets filled
00148   //by a call to GnumiInterface::FillANtpTruth
00149   fReg.UnLockKeys();
00150   fReg.UnLockValues();
00151 
00152   fReg.Set("event:nuparent_x",ir->parentX);
00153   fReg.Set("event:nuparent_y",ir->parentY);
00154   fReg.Set("event:nuparent_z",ir->parentZ);
00155   fReg.Set("event:nuparent_px",ir->parentPX);
00156   fReg.Set("event:nuparent_py",ir->parentPY);
00157   fReg.Set("event:nuparent_pz",ir->parentPZ);
00158   fReg.Set("event:nuparent_pid",ir->parentPID);
00159   fReg.Set("event:nuparent_gen",ir->parentGen);
00160   
00161   fReg.LockValues();
00162   fReg.LockKeys();
00163   
00164   return;
00165 }
00166 
00167 
00168 void ReweightHelpers::MCEventInfoFilla(MCEventInfo *ei, NtpStRecord *st, 
00169                                        int evtno)
00170 {
00171   if(evtno>=st->mc->GetEntries()){
00172     MSG("ReweightHelper",Msg::kError)<<"EventRegistryFilla Error, "
00173                                      <<"event no: "<<evtno
00174                                      <<" size "
00175                                      <<st->mc->GetEntries()<<std::endl;
00176       return;
00177    }
00178 
00179    NtpMCTruth *mcth = static_cast<NtpMCTruth *>((*st->mc)[evtno]);
00180 
00181    Int_t initial_state=0;  
00182    Int_t had_fs=0;
00183    bool foundfs=false;
00184    TClonesArray& heparray = *(st->stdhep);
00185    Int_t nhep = heparray.GetEntries();
00186 
00187    int protneut = -1;  // 0 = neutron, 1 = proton, 2 = N, 3 = A
00188    int nubarnu = 0;    // +1 = neutrino, -1 = antineutrino
00189    
00190    for(int i=0;i<nhep;i++){    
00191       NtpMCStdHep *sh = static_cast<NtpMCStdHep *>(heparray[i]);
00192       if(sh->mc==evtno){
00193          if(sh->IstHEP==0){  //initial state particle
00194             if(abs(sh->IdHEP)==12 || 
00195                abs(sh->IdHEP)==14 || 
00196                abs(sh->IdHEP)==16){   //(anti)neutrino    
00197                nubarnu = sh->IdHEP/abs(sh->IdHEP);  //get sign
00198             }
00199          }
00200          if(sh->IstHEP==11){    //target nucleon
00201             if(sh->IdHEP==2212) protneut = 1;  //proton
00202             else if(sh->IdHEP==2112) protneut = 0;  //neutron
00203             else if(abs(sh->IdHEP)>1000000000) protneut = 2;  //nucleus
00204             else protneut = 3; //atom - probably never get here since IdHEP A>N?
00205          }
00206          if(sh->IstHEP==3&&!foundfs){
00207            if(mcth->iresonance!=1002 ||   //if not RES event
00208               TMath::Abs(sh->IdHEP)!=15){ //or if it is RES, IdHEP!=tau lepton
00209              had_fs = (sh->IdHEP%1000);
00210              foundfs=true;
00211            }
00212          }
00213       }
00214    }
00215    
00216    if(protneut==1 && nubarnu==1)  initial_state=1;  //p + v
00217    if(protneut==0 && nubarnu==1)  initial_state=2;  //n + v
00218    if(protneut==1 && nubarnu==-1) initial_state=3;  //p + vbar
00219    if(protneut==0 && nubarnu==-1) initial_state=4;  //n + vbar
00220    if(protneut==2 && nubarnu==1)  initial_state=5;  //N + v
00221    if(protneut==3 && nubarnu==1)  initial_state=6;  //A + v
00222    if(protneut==2 && nubarnu==-1) initial_state=7;  //N + vbar
00223    if(protneut==3 && nubarnu==-1) initial_state=8;  //A + vbar
00224 
00225    Int_t nucleus=FindNucleusNumber((int)mcth->z,(int)mcth->a);
00226 
00227    ei->nuE=1.*mcth->p4neu[3];
00228    ei->nuPx=1.*mcth->p4neu[0];
00229    ei->nuPy=1.*mcth->p4neu[1];
00230    ei->nuPz=1.*mcth->p4neu[2];
00231 
00232    ei->tarE=1.*mcth->p4tgt[3];
00233    ei->tarPx=1.*mcth->p4tgt[0];
00234    ei->tarPy=1.*mcth->p4tgt[1];
00235    ei->tarPz=1.*mcth->p4tgt[2];
00236 
00237    ei->y=mcth->y;
00238    ei->x=mcth->x;
00239    ei->q2=mcth->q2;
00240    ei->w2=mcth->w2;
00241 
00242    ei->iaction=mcth->iaction;
00243    ei->inu=mcth->inu;
00244    ei->iresonance=mcth->iresonance;
00245    ei->initial_state=initial_state;
00246 
00247    ei->nucleus=nucleus;
00248 
00249    ei->had_fs=had_fs;
00250 }
00251 
00252 int ReweightHelpers::FindNucleusNumber(int z, int a)
00253 {
00254   Int_t nucleus=1;
00255 
00256   switch (z) {
00257     //  case 1:
00258     //nucleus=0;   // free nucleon
00259     //break;
00260   case 1:
00261     switch (a) {
00262     case 1:
00263       nucleus=256;   // hydrogen1
00264       break;
00265     case 2:
00266       nucleus=257;   // hydrogen2
00267       break;
00268     case 3:
00269       nucleus=258;   // hydrogen2
00270       break;
00271     default:
00272       nucleus=256;   // hydrogen1
00273       break;
00274     }
00275     break;
00276   case 6:
00277     switch (a) {
00278     case 11:
00279       nucleus=274; // carbon11   
00280       break;
00281     case 12:
00282       nucleus=275; // carbon12
00283       break;
00284     case 13:
00285       nucleus=276; // carbon13
00286       break;
00287     case 14:
00288       nucleus=277; // carbon14
00289       break;
00290     default:
00291       nucleus=275; // carbon12
00292       break;
00293     }
00294     break;
00295   case 7:
00296     switch (a) {
00297     case 13:
00298       nucleus=278; // nitrogen13   
00299       break;
00300     case 14:
00301       nucleus=279; // nitrogen14
00302       break;
00303     case 15:
00304       nucleus=280; // nitrogen15
00305       break;
00306     case 16:
00307       nucleus=281; // nitrogen16
00308       break;
00309     case 17:
00310       nucleus=282; // nitrogen17
00311       break;
00312     default:
00313       nucleus=279; // nitrogen14
00314       break;
00315     }
00316     break;
00317   case 8:
00318     switch (a) {
00319     case 15:
00320       nucleus=283;   // oxygen15
00321       break;
00322     case 16:
00323       nucleus=284;   // oxygen16
00324       break;
00325     case 17:
00326       nucleus=285;   // oxygen17
00327       break;
00328     case 18:
00329       nucleus=286;   // oxygen18
00330       break;
00331     default:
00332       nucleus=284;   // oxygen16
00333       break;
00334     }
00335     break;
00336   case 13:
00337     switch (a) {
00338     case 26:
00339       nucleus=303;   // aluminium26
00340       break;
00341     case 27:
00342       nucleus=304;   // aluminium27
00343       break;
00344     case 28:
00345       nucleus=305;   // aluminium28
00346       break;
00347     case 29:
00348       nucleus=306;   // aluminium29
00349       break;
00350     default:
00351       nucleus=304;   // aluminium27
00352       break;
00353     }
00354     break;
00355   case 14:
00356     switch (a) {
00357     case 27:
00358       nucleus=307;   // silicon27
00359       break;
00360     case 28:
00361       nucleus=308;   // silicon28
00362       break;
00363     case 29:
00364       nucleus=309;   // silicon29
00365       break;
00366     case 30:
00367       nucleus=310;   // silicon30
00368       break;
00369     default:         
00370       nucleus=308;   // silicon28
00371       break;
00372     }
00373     break;
00374   case 15:
00375     switch (a) {
00376     case 30:
00377       nucleus=311;   //phosphorus30
00378       break;
00379     case 31:
00380       nucleus=312;   //phosphorus31
00381       break;
00382     case 32:
00383       nucleus=313;   //phosphorus32
00384       break;
00385     case 33:
00386       nucleus=314;   //phosphorus33
00387       break;
00388     default:
00389       nucleus=312;   //phosphorus31
00390       break;
00391     }
00392     break;
00393   case 16:
00394     switch (a) {
00395     case 31:
00396       nucleus=315;   //sulphur31
00397       break;
00398     case 32:
00399       nucleus=316;   //sulphur32
00400       break;
00401     case 33:
00402       nucleus=317;   //sulphur33
00403       break;
00404     case 34:
00405       nucleus=318;   //sulphur34
00406       break;
00407     case 35:
00408       nucleus=319;   //sulphur35
00409       break;
00410     case 36:
00411       nucleus=320;   //sulphur36
00412       break;
00413     default:
00414       nucleus=316;   //sulphur32
00415       break;
00416     }
00417     break;
00418   case 22:
00419     switch (a) {
00420     case 45:
00421       nucleus=347;   //titanium45
00422       break;
00423     case 46:
00424       nucleus=348;   //titanium46
00425       break;
00426     case 47:
00427       nucleus=349;   //titanium47
00428       break;
00429     case 48:
00430       nucleus=350;   //titanium48
00431       break;
00432     case 49:
00433       nucleus=351;   //titanium49
00434       break;
00435     case 50:
00436       nucleus=352;   //titanium50
00437       break;
00438     default:
00439       nucleus=350;   //titanium48
00440       break;
00441     }
00442     break;
00443   case 23:
00444     switch (a) {
00445     case 49:
00446       nucleus=353;   //vanadium49
00447       break;
00448     case 50:
00449       nucleus=354;   //vanadium50
00450       break;
00451     case 51:
00452       nucleus=355;   //vanadium51
00453       break;
00454     case 52:
00455       nucleus=356;   //vanadium52
00456       break;
00457     case 53:
00458       nucleus=357;   //vanadium53
00459       break;
00460     default:
00461       nucleus=355;   //vanadium51
00462       break;
00463     }
00464     break;
00465   case 24:
00466     switch (a) {
00467     case 49:
00468       nucleus=358;   //chromium49
00469       break;
00470     case 50:
00471       nucleus=359;   //chromium50
00472       break;
00473     case 51:
00474       nucleus=360;   //chromium51
00475       break;
00476     case 52:
00477       nucleus=361;   //chromium52
00478       break;
00479     case 53:
00480       nucleus=362;   //chromium53
00481       break;
00482     case 54:
00483       nucleus=363;   //chromium54
00484       break;
00485     default:
00486       nucleus=361;   //chromium52
00487       break;
00488     }
00489     break;
00490   case 25:
00491     switch (a) {
00492     case 53:
00493       nucleus=364;   //manganese53
00494       break;
00495     case 54:
00496       nucleus=365;   //manganese54
00497       break;    
00498     case 55:
00499       nucleus=366;   //manganese55
00500       break;    
00501     case 56:
00502       nucleus=367;   //manganese56
00503       break;    
00504     case 57:
00505       nucleus=368;   //manganese57
00506       break;    
00507     default:
00508       nucleus=366;   //manganese55
00509       break;
00510     }
00511     break;
00512   case 26:
00513     switch (a) {
00514     case 53:
00515       nucleus=369;   //iron53
00516       break;
00517     case 54:
00518       nucleus=370;   //iron54
00519       break;
00520     case 55:
00521       nucleus=371;   //iron55
00522       break;
00523     case 56:
00524       nucleus=372;   //iron56
00525       break;
00526     case 57:
00527       nucleus=373;   //iron57
00528       break;
00529     case 58:
00530       nucleus=374;   //iron58
00531       break;
00532     default:
00533       nucleus=372;   //iron56
00534       break;
00535     }
00536     break;
00537   case 28:
00538     switch (a) {
00539     case 57:
00540       nucleus=382;   //nickel57
00541       break;
00542     case 58:
00543       nucleus=383;   //nickel58
00544       break;
00545     case 59:
00546       nucleus=384;   //nickel59
00547       break;
00548     case 60:
00549       nucleus=385;   //nickel60
00550       break;
00551     case 61:
00552       nucleus=386;   //nickel61
00553       break;
00554     case 62:
00555       nucleus=387;   //nickel62
00556       break;
00557     case 63:
00558       nucleus=388;   //nickel63
00559       break;
00560     case 64:
00561       nucleus=389;   //nickel64
00562       break;
00563     default:
00564       nucleus=383;   //nickel58
00565       break;
00566     }
00567     break;
00568   case 29:
00569     switch (a) {
00570     case 62:
00571       nucleus=390;   //copper62
00572       break;
00573     case 63:
00574       nucleus=391;   //copper63
00575       break;
00576     case 64:
00577       nucleus=392;   //copper64
00578       break;
00579     case 65:
00580       nucleus=393;   //copper65
00581       break;
00582     case 66:
00583       nucleus=394;   //copper66
00584       break;
00585     case 67:
00586       nucleus=395;   //copper67
00587       break;
00588     default:
00589       nucleus=392;   //copper64
00590       break;
00591     }
00592     break;
00593   default:
00594     nucleus=1;   // unknown
00595     break;
00596   }
00597 
00598   return nucleus;
00599 }