ROOT logo
//#include "UserTreeAnalysis.H"   // remove if copied to user working directory
#include <stdio.h>
#include <assert.h>
#include <math.h>
#include "MC4Vector.H"
#include "HEPParticle.H"
#include "TH1.h"
#include "Setup.H"
#include "TObjArray.h"
#include "TMath.h"
#define PI 3.141592653

inline bool ifSofty(int Id,int nparams, double *params){
  if (nparams<5 && Id==22) return true;   // to remove photons only
  for (int i=nparams-1; i>3; i--)
    if (Id==params[i]) return true;       // to remove all what is in params from nparams down to 4
  return false;
}

// very similar to  MC_FillUserHistogram from Generate.cxx
inline void fillUserHisto(char *name,double val, double weight=1.0, 
                          double min=Setup::bin_min[0][0], 
                          double max=Setup::bin_max[0][0]){

    TH1D *h=(TH1D*)(Setup::user_histograms->FindObject(name));
    if(!h){
      h=new TH1D(name,name,Setup::nbins[0][0],min,max);
      if(!h) return;
      Setup::user_histograms->Add(h);
      //      printf("user histogram created %s\n", name);
}
    h->Fill(val,weight);

}
double angle(double X, double Y){

 
    double  an=0.0; 
    double R=sqrt(X*X+Y*Y); 
    //    if(R<pow(10,-20)) printf(" angle this time X %f\n", 10000*X);
    if(R<pow(10,-20)) return an;

    if(TMath::Abs(X)/R<0.8) 
    { 
        an=acos(X/R);
	if(Y<0 && an>0) an=-an;
        if(Y>0 && an<0) an=-an; 
    }
    else 
    { 
      an=asin(Y/R);
      if(X<0 && an>=0.)  an=PI-an; 
      else if(X<0.)      an=-PI-an; 
       
    } 
  return an;
}

int counter=0;

int ZtautauAnalysis(HEPParticle *mother,HEPParticleList *stableDaughters, int nparams, double *params)
{
  // THIS IS EXAMPLE of userTreeAnalysis. It acts on list of particle X (under our MC-test) final
  // daughters. Of course in real life many options may need to be introduced by the user.
    // PARAMETERS:
    // params[0] threshold on Energy (or p_T) of particle expressed as a fraction of mother's 
    //		 Energy (or p_T) in mother's frame. If not specified - default is 0.05 
    //           Note that setting it to zero is possible.
    // params[1] maximum number of left soft-suspected particles. 
    //           0: (default) all listed particles are removed, even if hard
    //           1: 2:  one two etc  removable particles are kept (at most)
    //          -1: this option is off, all hard particles stay.
    // params[2] control tag on discriminating particles, 
    //           0: (default) Energy in decaying particle rest frame
    //           1: Energy in lab frame
    //           2: p_T with respect to beam direction in lab frame.
    // params[3] type of action to be applied on removed daughters 
    //           0: (default) nothing is done, they are lost
    //           1:  algorithm as in studies on PHOTOS is used, see papers P. Golonka and Z. Was.
    // params[4] from this paramter on PDG id-s of particles to be removed can be added. 
    //           Default is that only photons are removed.
    // To get to origin of our particle (mother) one need to go after UserEventAnalysis,
    // instructive example will be given later. 
    assert(mother!=0);
    assert(stableDaughters!=0);


    double threshold=0.05, leftmax=0.0, selector=0.0, actLost=0.0; 
    if (nparams>0 && params==0)  return -1;
    if (nparams>0) threshold=params[0]; 
    if (nparams>1) leftmax  =params[1]; 
    if (nparams>2) selector =params[2]; 
    if (nparams>3) actLost  =params[3]; 

    
    HEPParticleList *lostDaughters=new HEPParticleList;    
    double pt_limit=threshold*mother->GetM();

    HEPParticleListIterator daughters (*stableDaughters);
    int nphot=0;
    double ephot=pow(10,22);
    bool redo=false;
    for (HEPParticle *part=daughters.first(); part!=0; (redo)? part=part:part=daughters.next() ) {
      if(redo) redo=false;
	MC4Vector d4(part->GetE(),part->GetPx(),part->GetPy(),part->GetPz(),part->GetM());
	// boost to mother's frame:
	d4.Boost(mother->GetPx(),mother->GetPy(),mother->GetPz(),mother->GetE(),mother->GetM());



	double p_t=d4.GetX0();  // default
	switch ((int) selector)
        {
	case 1: p_t=part->GetE(); break;
	case 2: p_t=d4.Xt(); break;
	}


        if(ifSofty(part->GetPDGId(),nparams,params)) nphot++;
	if ( ifSofty(part->GetPDGId(),nparams,params) && p_t < pt_limit) {
	  //	  printf("Odrzucamy! %i\n",count++);
          nphot=0;
            lostDaughters->push_back(part);
	    stableDaughters->remove(part);
	    part=daughters.first();
            redo=true;
            continue;
	}
        if( ifSofty(part->GetPDGId(),nparams,params) && ephot>p_t) ephot=p_t;
    }
    while(nphot>(int) leftmax)
    {
      double ephot1=pow(10,22);
      redo=false;
      for (HEPParticle *part=daughters.first(); part!=0; (redo)? part=part:part=daughters.next() ) {
        if(redo) redo=false;
	MC4Vector d4(part->GetE(),part->GetPx(),part->GetPy(),part->GetPz(),part->GetM());
	// boost to mother's frame:
	d4.Boost(mother->GetPx(),mother->GetPy(),mother->GetPz(),mother->GetE(),mother->GetM());


	double p_t=d4.GetX0();  // default
	switch ((int) selector)
        {
	case 1: p_t=part->GetE(); break;
	case 2: p_t=d4.Xt(); break;
	}


	if ( ifSofty(part->GetPDGId(),nparams,params) && p_t == ephot) {

	  nphot--;
            lostDaughters->push_back(part);
	    stableDaughters->remove(part);
            ephot1=pow(10,22);
	     part=daughters.first();
             redo=true;
             continue;
	}
        if( ifSofty(part->GetPDGId(),nparams,params) && ephot1>p_t) ephot1=p_t;
      }
      ephot=ephot1;

    }

    // ##############################################################
    // here functionality of removig some daughters is finished
    // now we reconsider what to do with them
    // ##############################################################
    // delete lostDaughters;
    // return 1; 


    // ##############################################################
    // Now: What to do with lost daughters?
    // ##############################################################
    //    lostDaughters->ls();
    int version=(int) actLost;
        
    switch(version)
      {
      case 1:  // add lost to charged
	{
          HEPParticleListIterator lost (*lostDaughters);
	  for (HEPParticle *lostpart=lost.first(); lostpart!=0; lostpart=lost.next() ) {
            HEPParticle* Best=0;
	    double searchvirt=pow(10.0,22);
            MC4Vector VV;
	    for( HEPParticle *part=daughters.first(); part!=0; part=daughters.next() ){
	      if(part->GetCharge()==0) continue;
              VV=lostpart->GetP4()+part->GetP4();
	      VV.AdjustM();
	      if (VV.GetM()<searchvirt) {searchvirt=VV.GetM(); Best=part;}
	    }
	    if(Best) {
              Best->SetPx(Best->GetPx()+lostpart->GetPx());
              Best->SetPy(Best->GetPy()+lostpart->GetPy());
              Best->SetPz(Best->GetPz()+lostpart->GetPz());
              Best->SetE (Best->GetE ()+lostpart->GetE ());
	    }
	  }
	break;
	}
      default: break; // do nothing
      }

    delete lostDaughters;


	/*
		Create histogram
	*/
	double px=0,py=0,pz=0,e=0;
	double px1=0,py1=0,pz1=0,e1=0;

	//loop over all daughters
	for(HEPParticle *part=daughters.first();part!=0;part=daughters.next())
	{
		if(abs(part->GetPDGId())==211)
		{
			px1+=part->GetPx();
			py1+=part->GetPy();
			pz1+=part->GetPz();
			e1+=part->GetE();
		}
		if(abs(part->GetPDGId())==211) continue;
		//Sum e+ e- 4-vectors
		px+=part->GetPx();
		py+=part->GetPy();
		pz+=part->GetPz();
		e+=part->GetE();
	}
	double mass=e*e-px*px-py*py-pz*pz;
	double mass1=(e+e1)*(e+e1)-(px+px1)*(px+px1)-(py+py1)*(py+py1)-(pz+pz1)*(pz+pz1);
	mass=1-mass/mass1;
	char mm[] = "pi-energy";
	fillUserHisto(mm,mass,1.0,0.0,1.1);

    //For ZtautauAnalysis - other histograms are not needed
    bool activateUserHist=false;
    if(!activateUserHist) return 1;

    // segmet of user defined histograms

    double X=mother->GetPx(), Y=mother->GetPy(), Z=mother->GetPz(); 
    // double E=mother->GetE(),  MM=mother->GetM();


    double pt=sqrt(X*X+Y*Y);
    double eta=log((sqrt(pt*pt+Z*Z)+TMath::Abs(Z))/pt);
    if(Z<0 && eta>0) eta=-eta;
    if(Z>0 && eta<0) eta=-eta;
    double phi=angle(X,Y);
    char hist1[] = "mother-PT";
    char hist2[] = "mother-eta";
    char hist3[] = "mother-phi";
    fillUserHisto(hist1,pt,1.0,0.0,100.0);
    fillUserHisto(hist2,eta,1.0,-8.0,8.0);
    fillUserHisto(hist3,phi,1.0,-PI,PI);
    return 1;
}
 ZtautauAnalysis.C:1
 ZtautauAnalysis.C:2
 ZtautauAnalysis.C:3
 ZtautauAnalysis.C:4
 ZtautauAnalysis.C:5
 ZtautauAnalysis.C:6
 ZtautauAnalysis.C:7
 ZtautauAnalysis.C:8
 ZtautauAnalysis.C:9
 ZtautauAnalysis.C:10
 ZtautauAnalysis.C:11
 ZtautauAnalysis.C:12
 ZtautauAnalysis.C:13
 ZtautauAnalysis.C:14
 ZtautauAnalysis.C:15
 ZtautauAnalysis.C:16
 ZtautauAnalysis.C:17
 ZtautauAnalysis.C:18
 ZtautauAnalysis.C:19
 ZtautauAnalysis.C:20
 ZtautauAnalysis.C:21
 ZtautauAnalysis.C:22
 ZtautauAnalysis.C:23
 ZtautauAnalysis.C:24
 ZtautauAnalysis.C:25
 ZtautauAnalysis.C:26
 ZtautauAnalysis.C:27
 ZtautauAnalysis.C:28
 ZtautauAnalysis.C:29
 ZtautauAnalysis.C:30
 ZtautauAnalysis.C:31
 ZtautauAnalysis.C:32
 ZtautauAnalysis.C:33
 ZtautauAnalysis.C:34
 ZtautauAnalysis.C:35
 ZtautauAnalysis.C:36
 ZtautauAnalysis.C:37
 ZtautauAnalysis.C:38
 ZtautauAnalysis.C:39
 ZtautauAnalysis.C:40
 ZtautauAnalysis.C:41
 ZtautauAnalysis.C:42
 ZtautauAnalysis.C:43
 ZtautauAnalysis.C:44
 ZtautauAnalysis.C:45
 ZtautauAnalysis.C:46
 ZtautauAnalysis.C:47
 ZtautauAnalysis.C:48
 ZtautauAnalysis.C:49
 ZtautauAnalysis.C:50
 ZtautauAnalysis.C:51
 ZtautauAnalysis.C:52
 ZtautauAnalysis.C:53
 ZtautauAnalysis.C:54
 ZtautauAnalysis.C:55
 ZtautauAnalysis.C:56
 ZtautauAnalysis.C:57
 ZtautauAnalysis.C:58
 ZtautauAnalysis.C:59
 ZtautauAnalysis.C:60
 ZtautauAnalysis.C:61
 ZtautauAnalysis.C:62
 ZtautauAnalysis.C:63
 ZtautauAnalysis.C:64
 ZtautauAnalysis.C:65
 ZtautauAnalysis.C:66
 ZtautauAnalysis.C:67
 ZtautauAnalysis.C:68
 ZtautauAnalysis.C:69
 ZtautauAnalysis.C:70
 ZtautauAnalysis.C:71
 ZtautauAnalysis.C:72
 ZtautauAnalysis.C:73
 ZtautauAnalysis.C:74
 ZtautauAnalysis.C:75
 ZtautauAnalysis.C:76
 ZtautauAnalysis.C:77
 ZtautauAnalysis.C:78
 ZtautauAnalysis.C:79
 ZtautauAnalysis.C:80
 ZtautauAnalysis.C:81
 ZtautauAnalysis.C:82
 ZtautauAnalysis.C:83
 ZtautauAnalysis.C:84
 ZtautauAnalysis.C:85
 ZtautauAnalysis.C:86
 ZtautauAnalysis.C:87
 ZtautauAnalysis.C:88
 ZtautauAnalysis.C:89
 ZtautauAnalysis.C:90
 ZtautauAnalysis.C:91
 ZtautauAnalysis.C:92
 ZtautauAnalysis.C:93
 ZtautauAnalysis.C:94
 ZtautauAnalysis.C:95
 ZtautauAnalysis.C:96
 ZtautauAnalysis.C:97
 ZtautauAnalysis.C:98
 ZtautauAnalysis.C:99
 ZtautauAnalysis.C:100
 ZtautauAnalysis.C:101
 ZtautauAnalysis.C:102
 ZtautauAnalysis.C:103
 ZtautauAnalysis.C:104
 ZtautauAnalysis.C:105
 ZtautauAnalysis.C:106
 ZtautauAnalysis.C:107
 ZtautauAnalysis.C:108
 ZtautauAnalysis.C:109
 ZtautauAnalysis.C:110
 ZtautauAnalysis.C:111
 ZtautauAnalysis.C:112
 ZtautauAnalysis.C:113
 ZtautauAnalysis.C:114
 ZtautauAnalysis.C:115
 ZtautauAnalysis.C:116
 ZtautauAnalysis.C:117
 ZtautauAnalysis.C:118
 ZtautauAnalysis.C:119
 ZtautauAnalysis.C:120
 ZtautauAnalysis.C:121
 ZtautauAnalysis.C:122
 ZtautauAnalysis.C:123
 ZtautauAnalysis.C:124
 ZtautauAnalysis.C:125
 ZtautauAnalysis.C:126
 ZtautauAnalysis.C:127
 ZtautauAnalysis.C:128
 ZtautauAnalysis.C:129
 ZtautauAnalysis.C:130
 ZtautauAnalysis.C:131
 ZtautauAnalysis.C:132
 ZtautauAnalysis.C:133
 ZtautauAnalysis.C:134
 ZtautauAnalysis.C:135
 ZtautauAnalysis.C:136
 ZtautauAnalysis.C:137
 ZtautauAnalysis.C:138
 ZtautauAnalysis.C:139
 ZtautauAnalysis.C:140
 ZtautauAnalysis.C:141
 ZtautauAnalysis.C:142
 ZtautauAnalysis.C:143
 ZtautauAnalysis.C:144
 ZtautauAnalysis.C:145
 ZtautauAnalysis.C:146
 ZtautauAnalysis.C:147
 ZtautauAnalysis.C:148
 ZtautauAnalysis.C:149
 ZtautauAnalysis.C:150
 ZtautauAnalysis.C:151
 ZtautauAnalysis.C:152
 ZtautauAnalysis.C:153
 ZtautauAnalysis.C:154
 ZtautauAnalysis.C:155
 ZtautauAnalysis.C:156
 ZtautauAnalysis.C:157
 ZtautauAnalysis.C:158
 ZtautauAnalysis.C:159
 ZtautauAnalysis.C:160
 ZtautauAnalysis.C:161
 ZtautauAnalysis.C:162
 ZtautauAnalysis.C:163
 ZtautauAnalysis.C:164
 ZtautauAnalysis.C:165
 ZtautauAnalysis.C:166
 ZtautauAnalysis.C:167
 ZtautauAnalysis.C:168
 ZtautauAnalysis.C:169
 ZtautauAnalysis.C:170
 ZtautauAnalysis.C:171
 ZtautauAnalysis.C:172
 ZtautauAnalysis.C:173
 ZtautauAnalysis.C:174
 ZtautauAnalysis.C:175
 ZtautauAnalysis.C:176
 ZtautauAnalysis.C:177
 ZtautauAnalysis.C:178
 ZtautauAnalysis.C:179
 ZtautauAnalysis.C:180
 ZtautauAnalysis.C:181
 ZtautauAnalysis.C:182
 ZtautauAnalysis.C:183
 ZtautauAnalysis.C:184
 ZtautauAnalysis.C:185
 ZtautauAnalysis.C:186
 ZtautauAnalysis.C:187
 ZtautauAnalysis.C:188
 ZtautauAnalysis.C:189
 ZtautauAnalysis.C:190
 ZtautauAnalysis.C:191
 ZtautauAnalysis.C:192
 ZtautauAnalysis.C:193
 ZtautauAnalysis.C:194
 ZtautauAnalysis.C:195
 ZtautauAnalysis.C:196
 ZtautauAnalysis.C:197
 ZtautauAnalysis.C:198
 ZtautauAnalysis.C:199
 ZtautauAnalysis.C:200
 ZtautauAnalysis.C:201
 ZtautauAnalysis.C:202
 ZtautauAnalysis.C:203
 ZtautauAnalysis.C:204
 ZtautauAnalysis.C:205
 ZtautauAnalysis.C:206
 ZtautauAnalysis.C:207
 ZtautauAnalysis.C:208
 ZtautauAnalysis.C:209
 ZtautauAnalysis.C:210
 ZtautauAnalysis.C:211
 ZtautauAnalysis.C:212
 ZtautauAnalysis.C:213
 ZtautauAnalysis.C:214
 ZtautauAnalysis.C:215
 ZtautauAnalysis.C:216
 ZtautauAnalysis.C:217
 ZtautauAnalysis.C:218
 ZtautauAnalysis.C:219
 ZtautauAnalysis.C:220
 ZtautauAnalysis.C:221
 ZtautauAnalysis.C:222
 ZtautauAnalysis.C:223
 ZtautauAnalysis.C:224
 ZtautauAnalysis.C:225
 ZtautauAnalysis.C:226
 ZtautauAnalysis.C:227
 ZtautauAnalysis.C:228
 ZtautauAnalysis.C:229
 ZtautauAnalysis.C:230
 ZtautauAnalysis.C:231
 ZtautauAnalysis.C:232
 ZtautauAnalysis.C:233
 ZtautauAnalysis.C:234
 ZtautauAnalysis.C:235
 ZtautauAnalysis.C:236
 ZtautauAnalysis.C:237
 ZtautauAnalysis.C:238
 ZtautauAnalysis.C:239
 ZtautauAnalysis.C:240
 ZtautauAnalysis.C:241
 ZtautauAnalysis.C:242
 ZtautauAnalysis.C:243
 ZtautauAnalysis.C:244
 ZtautauAnalysis.C:245
 ZtautauAnalysis.C:246
 ZtautauAnalysis.C:247
 ZtautauAnalysis.C:248
 ZtautauAnalysis.C:249
 ZtautauAnalysis.C:250
 ZtautauAnalysis.C:251
 ZtautauAnalysis.C:252
 ZtautauAnalysis.C:253
 ZtautauAnalysis.C:254
 ZtautauAnalysis.C:255
 ZtautauAnalysis.C:256
 ZtautauAnalysis.C:257
 ZtautauAnalysis.C:258
 ZtautauAnalysis.C:259
 ZtautauAnalysis.C:260
 ZtautauAnalysis.C:261
 ZtautauAnalysis.C:262