! ! ---------------------------- MC for dalitz-decays ------------------------------- ! ! Genesis program for pair analysis ! ! input defined in runcard file which needs to be set as link befor ! execution of program. ( ln -s 'file' runcard ) ! ! existing runcard 'files': runcard-pBe ! runcard-Pb ! ! output file, hbook file, and resource file are defined in the runcards ! ! ! 3-10-97 Axel ! ! -------------------------------variable and parameter declarations ------------ ! use genesis_m use mom_m use runcard_m implicit none ! ! parameters and variablers for histograming ! integer, parameter :: nwpawc = 4000000 integer :: istat, icycle, hmemor, iquest integer, parameter :: lendata = 12 integer, parameter :: lunhb = 10, lunout = 20 integer, parameter :: idp = 2 real :: data(lendata) character(len=8) :: ptags(lendata) common /pawc/ hmemor(nwpawc) common /quest/ iquest(100) ! ! other variables ! integer :: ievent, ipart, idec, ichild ! loop index integer :: pairs = 0 ! counter type(tree), pointer :: ptree, pchild, pele, ppos ! pointer to particles in tree real :: yrap,ppt,the,thp,wt,eta real, external :: atg ! ! for control of cpu-time ! integer :: clock_at_start, clock_at_end integer :: clock_now, clock_rate ! ! ------------------------------- end of declarations ----------------------------------- ! ! read runcards ! call GeGetRuncard ! ! ------------------------------- initialization ---------------------------------------- ! ! Init Hbook ! if (write_hb) then print *,'booking histos...' call hlimit(nwpawc) iquest(10) = 65000 print "(' open histo file: ',a)", hbfile call hropen(lunhb,'HISTOS',hbfile,'NQ',1024,istat) if (istat.ne.0) then print *,'pairs: error opening hbook file. stop.' stop endif ptags(1) = 'mass' ! invariant mass of pair ptags(2) = 'pe' ! p of electron ptags(3) = 'the' ! theta of electron ptags(4) = 'phe' ! phi of electron ptags(5) = 'pp' ! p of positron ptags(6) = 'thp' ! theta of positron ptags(7) = 'php' ! phi of positron ptags(8) = 'p' ! p of parent particle ptags(9) = 'th' ! theta of parent particle ptags(10) = 'ph' ! phi of parent particle ptags(11) = 'dec' ! decay id ptags(12) = 'wt' ! weight factor for each decay call hbookn(idp,'pairs',lendata,'HISTOS',1024,ptags) call hbook1(8121,' pion0 pt ',100,0.,4.,0.) call hbook1(8122,' eta pt ',50,0.,2.,0.) call hbook1(8123,' omega pt ',50,0.,2.,0.) call hbook1(8124,' rho pt ',50,0.,2.,0.) call hbook1(8125,' etap pt ',50,0.,2.,0.) call hbook1(8126,' phi pt ',50,0.,2.,0.) call hbook1(9121,' pion0 y ',100,-2.,8.,0.) call hbook1(9122,' eta y ',100,-2.,8.,0.) call hbook1(9123,' omega y ',100,-2.,8.,0.) call hbook1(9124,' rho y ',100,-2.,8.,0.) call hbook1(9125,' etap y ',100,-2.,8.,0.) call hbook1(9126,' phi y ',100,-2.,8.,0.) call hbook1(9221,' pion0 eta ',100,-2.,8.,0.) call hbook1(9222,' eta eta ',100,-2.,8.,0.) call hbook1(9223,' omega eta ',100,-2.,8.,0.) call hbook1(9224,' rho eta ',100,-2.,8.,0.) call hbook1(9225,' etap eta ',100,-2.,8.,0.) call hbook1(9226,' phi eta ',100,-2.,8.,0.) endif ! ! Open data file ! if (write_data) then print "(' open output file: ',a)", outfile open(lunout,file=outfile,form='unformatted',status='unknown',iostat=istat) if (istat.ne.0) then print *,'pairs: Error opening output file. stop.' stop endif endif ! ! Init Genesis ! call GeInit(resFile, .true.) ! ! Print out runing conditions ! print *,'Run started: ',eurodate() print *,'----------------------------------------------------' print *,'Events................ ',nevents print *,'Multiplicity.......... ',nch print *,'Hbook file ........... ',trim(hbfile) print *,'Acceptance (theta).... ',theta_lower,theta_upper print *,'Pt-cut................ ',pt_min,' GeV/c' print *,'----------------------------------------------------' ! ! Timing ! call system_clock(count=clock_at_start) ! ! ------------------------------- end of initialization --------------------------------- ! ! ------------------------------- start of event loop ---------------------------------- ! event_loop: do ievent = 1, nevents call GeTriggerEvent(nch,bpar) ! trigger one event if (mod(ievent,100) == 0) then ! print-out call system_clock(count=clock_now,count_rate=clock_rate) print '("event: ",i5,5x,"still ",i4," min")', & ievent, & nint((float(nevents)/ievent-1.)* & float(clock_now-clock_at_start)/clock_rate/60.) endif particle_loop: do ipart = 1, lasttree ! loop all particles ptree => treeList(ipart) ! pointer to tree ! if (ptree%parent == 0) then ! take only primary produced particles ! if (ptree%parent /= 0 .and. associated(ptree%decaytype)) then ! take only secondary decays ! ! Assign decay id ! idec = 0 if (ptree%decaytype%name == 'eta->electron+positron+photon') idec = 1 if (ptree%decaytype%name == 'rho->electron+positron') idec = 2 if (ptree%decaytype%name == 'omega->pion0+electron+positron') idec = 3 if (ptree%decaytype%name == 'omega->electron+positron') idec = 4 if (ptree%decaytype%name == 'phi->electron+positron') idec = 5 if (ptree%decaytype%name == 'phi->eta+electron+positron') idec = 6 if (ptree%decaytype%name == 'pion0->photon+electron+positron') idec = 7 if (ptree%decaytype%name == "eta'->electron+positron+photon") idec = 8 if (idec.eq.0) cycle particle_loop ! cycle loop if not dilepton decay ! ! Loop children to get pointer to electron and positron ! child_loop: do ichild = 1, ptree%decaytype%nBody pchild => treeList(ptree%child(ichild)) if (pchild%parttype%name == "electron") then pele => treeList(ptree%child(ichild)) else if (pchild%parttype%name == "positron") then ppos => treeList(ptree%child(ichild)) endif enddo child_loop ! ! histogram parent particles ! yrap = yof(ptree%p4) ! particle rapidity ppt = ptof(ptree%p4%p) ! and pt eta = etaof(ptree%p4%p) ! and pseudo rapidity wt = ptree%mass_weight if (ptree%parttype%name == 'pion0') then if (3.1= theta_upper ) cycle particle_loop if (thp <= theta_lower .or. thp >= theta_upper ) cycle particle_loop ppt = ptof(pele%p4%p) if ( ppt <= pt_min) cycle particle_loop ! ppt = sqrt(ppos%p4%p%px**2 + ppos%p4%p%py**2) ppt = ptof(ppos%p4%p) if ( ppt <= pt_min) cycle particle_loop ! ! store kinematics of pair ! data(1) = abs(pele%p4 + ppos%p4) ! invariant mass data(2) = abs(pele%p4%p) ! electron momentum data(3) = the ! electron theta data(4) = atg(sngl(pele%p4%p%py),sngl(pele%p4%p%px)) ! electron phi data(5) = abs(ppos%p4%p) ! positron momentum data(6) = thp ! positron theta data(7) = atg(sngl(ppos%p4%p%py),sngl(ppos%p4%p%px)) ! positron phi data(8) = abs(ptree%p4%p) ! parent p data(9) = thetaof(ptree%p4%p) ! parent theta data(10) = atg(sngl(ptree%p4%p%py),sngl(ptree%p4%p%px)) ! parent phi data(11) = idec ! decay data(12) = pele%mass_weight ! weight if (write_hb) call hfn (idp, data) ! store in ntuple if (write_data) write(lunout) data ! write data to output file if requested pairs = pairs + 1 ! count number of pairs endif enddo particle_loop enddo event_loop ! ! ------------------------------- end of event loop ---------------------------------- ! ! ------------------------------- store output and clean up --------------------------- ! ! Close files ! if (write_hb) then call hrout (0,icycle,' ') call hrend ('HISTOS') close (lunhb) endif if (write_data) close (lunout) ! ! Timing ! call system_clock(count=clock_at_end,count_rate=clock_rate) ! ! Say goodbye ! print *,'Run ended: ',eurodate() print *,'----------------------------------------------------' print *,'Events processed...... ',ievent-1 print *,'Pairs generated...... ',pairs print *,'CPU total (min)....... ', & float(clock_at_end-clock_at_start)/clock_rate/60. print *,'CPU/event (s)......... ', & float(clock_at_end-clock_at_start)/clock_rate/ievent print *,'CPU/pair (ms)......... ', & 1000.*float(clock_at_end-clock_at_start)/clock_rate/pairs print *,'----------------------------------------------------' contains function eurodate() character(len=19) :: eurodate character(len=8) :: cdate character(len=10) :: ctime call date_and_time(date=cdate,time=ctime) eurodate = cdate(7:8)//'/'//cdate(5:6)//'/'//cdate(1:4)//' '// & ctime(1:2)//':'//ctime(3:4)//':'//ctime(5:6) end function eurodate end program