Using a batch script to simulate pileup


{

    //first, we define in the data base:
    //1.) The mean time difference between 2 events
    //    N.B. the time in Pluto is in mm/s, so we have to convert it
    //    let us start with 100000 events/s
    makeStaticData()->SetBatchValue("beam_particle_mean_time", 10e-5 * TMath::C() / 10e-3);
    
    //2.) The time of the first event (we choose 0)
    makeStaticData()->SetBatchValue("absolute_event_time", 0);
    
    //3.) The threshold for a pile-up event (just e.g. 50% of beam_particle_mean_time)
    makeStaticData()->SetBatchValue("pileup_time_theshold", 0.5 * 10e-5 * TMath::C() / 10e-3);

    //4.) Just a flag for pileup
    makeStaticData()->SetBatchValue("pileup_flag", 0);

    //this just constructs the temporary particles:
    makeDynamicData()->GetBatchParticle("p1");
    makeDynamicData()->GetBatchParticle("p2");
    makeDynamicData()->GetBatchParticle("eta1");

    //**************************************************************
    PReaction my_reaction("_T1=2.2", "p", "p", "p p eta", "pileup_eta");
    my_reaction.trackedParticles();  //this removes the "empty events"
    //**************************************************************
    
    //the following loop sets the absolute time for each particle
    my_reaction.Do("foreach(*); event_time = [*]->T(); [*]->SetT(event_time + absolute_event_time)");

    //set a event time difference based on a "flat" distribution
    my_reaction.Do("timediff = sampleFlat()*2*beam_particle_mean_time");
    
    //check if the last event was pileup and has been removed:
    my_reaction.Do("if !pileup_flag; goto check_pileup");
    my_reaction.Do("push(p1); push(p2); push(eta1);");
    my_reaction.Do("pileup_flag = 0");
    my_reaction.Do("goto event_ok");

    //check if we are below threshold
    my_reaction.Do("check_pileup:");
    my_reaction.Do("if timediff > pileup_time_theshold; goto event_ok");
    
    //if yes, store particles and kill event
    my_reaction.Do("p1 = [p,1]; p2 = [p,2]; eta1 = [eta];");
    my_reaction.Do("[p,1]->SetInActive(); [p,2]->SetInActive(); [eta]->SetInActive();");
    my_reaction.Do("pileup_flag = 1");
    

    my_reaction.Do("event_ok:");
    //set the new time for the next event
    my_reaction.Do("absolute_event_time = absolute_event_time + timediff");

    

    //this is just for debugging:
    my_reaction.Do("time = [p,1]->T(); echo $time");
    my_reaction.Do("foreach(*); [*]->Print();");
    
    my_reaction.Print();   //The "Print()" statement is optional
    my_reaction.Loop(100);
}
 eta_pileup.C:1
 eta_pileup.C:2
 eta_pileup.C:3
 eta_pileup.C:4
 eta_pileup.C:5
 eta_pileup.C:6
 eta_pileup.C:7
 eta_pileup.C:8
 eta_pileup.C:9
 eta_pileup.C:10
 eta_pileup.C:11
 eta_pileup.C:12
 eta_pileup.C:13
 eta_pileup.C:14
 eta_pileup.C:15
 eta_pileup.C:16
 eta_pileup.C:17
 eta_pileup.C:18
 eta_pileup.C:19
 eta_pileup.C:20
 eta_pileup.C:21
 eta_pileup.C:22
 eta_pileup.C:23
 eta_pileup.C:24
 eta_pileup.C:25
 eta_pileup.C:26
 eta_pileup.C:27
 eta_pileup.C:28
 eta_pileup.C:29
 eta_pileup.C:30
 eta_pileup.C:31
 eta_pileup.C:32
 eta_pileup.C:33
 eta_pileup.C:34
 eta_pileup.C:35
 eta_pileup.C:36
 eta_pileup.C:37
 eta_pileup.C:38
 eta_pileup.C:39
 eta_pileup.C:40
 eta_pileup.C:41
 eta_pileup.C:42
 eta_pileup.C:43
 eta_pileup.C:44
 eta_pileup.C:45
 eta_pileup.C:46
 eta_pileup.C:47
 eta_pileup.C:48
 eta_pileup.C:49
 eta_pileup.C:50
 eta_pileup.C:51
 eta_pileup.C:52
 eta_pileup.C:53
 eta_pileup.C:54
 eta_pileup.C:55
 eta_pileup.C:56
 eta_pileup.C:57
 eta_pileup.C:58
 eta_pileup.C:59
 eta_pileup.C:60
 eta_pileup.C:61
 eta_pileup.C:62
 eta_pileup.C:63
 eta_pileup.C:64