#define musrAnalysis_cxx
#include "musrAnalysis.hh"
#include "musrTH.hh"
//#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <algorithm>
#include <TF1.h>
//#include <math.h>

void musrAnalysis::Loop(char* runChar, char* v1190FileName, Int_t nrEvents)
{
//   In a ROOT session, you can do:
//      Root > .L musrAnalysis.C
//      Root > musrAnalysis t
//      Root > t.GetEntry(12); // Fill t data members with entry number 12
//      Root > t.Show();       // Show values of entry 12
//      Root > t.Show(16);     // Read and show values of entry 16
//      Root > t.Loop();       // Loop on all entries
//

//     This is the loop skeleton where:
//    jentry is the global entry number in the chain
//    ientry is the entry number in the current Tree
//  Note that the argument to GetEntry must be:
//    jentry for TChain::GetEntry
//    ientry for TTree::GetEntry and TBranch::GetEntry
//
//       To read only selected branches, Insert statements like:
// METHOD1:
//    fChain->SetBranchStatus("*",0);  // disable all branches
//    fChain->SetBranchStatus("branchname",1);  // activate branchname
// METHOD2: replace line
//    fChain->GetEntry(jentry);       //read all branches
//by  b_branchname->GetEntry(ientry); //read only this branch
   if (fChain == 0) return;

   hInfo = new TH1D("hInfo","Different values that need to be passed to ploting program",10000,1.,10001.);
   ReadInInputParameters(v1190FileName);
   CreateHistograms();

   nentries = fChain->GetEntriesFast();
   if (nrEvents!=0)  nentries = nrEvents;

   Long64_t nbytes = 0, nb = 0;
   for (Long64_t jentry=0; jentry<nentries;jentry++) {
      Long64_t ientry = LoadTree(jentry);
      if (ientry < 0) break;
      nb = fChain->GetEntry(jentry);  InitialiseEvent();  nbytes += nb;
      // if (Cut(ientry) < 0) continue;
      if (jentry%100000==0) std::cout << "Analysing event nr. jentry=" << jentry << std::endl;
      AnalyseEvent(jentry);
   }
   SaveHistograms(runChar,v1190FileName);
}



//================================================================

void musrAnalysis::ReadInInputParameters(char* charV1190FileName) {
  std::cout<<"musrAnalysis::ReadInInputParameters()"<<std::endl;
 
  //  std::string V1190_FileName = "2.v1190";
  //  FILE *fSteeringFile=fopen(V1190_FileName.c_str(),"r");
  FILE *fSteeringFile=fopen(charV1190FileName,"r");
  if (fSteeringFile==NULL) {
    if (strcmp(charV1190FileName,"Unset")==0) {
      std::cout<<"musrDetectorConstruction:  No input macro file specified"<<std::endl;
    }
    else {
      std::cout << "E R R O R :    Failed to open detector configuration file " << charV1190FileName << std::endl;
    }
    std::cout << "S T O P    F O R C E D" << std::endl;
    exit(1);
  }
  std::cout << "Configuration file \"" << charV1190FileName << "\" was opened."<<std::endl;
  char  line[1001];

  // Write out the configuration file into the output file:
  std::cout << "\n\n....oooOO0OOooo........oooOO0OOooo......Configuration file used for this run....oooOO0OOooo........oooOO0OOooo......"<<std::endl; 
  while (!feof(fSteeringFile)) {
    fgets(line,1000,fSteeringFile);
    //    if ((line[0]!='#')&&(line[0]!='\n')&&(line[0]!='\r'))  // TS: Do not print comments
    std::cout << line;
  }
  std::cout <<"....oooOO0OOooo........oooOO0OOooo......End of the configuration file.....oooOO0OOooo........oooOO0OOooo......\n\n"<<std::endl;

  //  Loop
  rewind (fSteeringFile);
  instrument = "Unset";
  description = "Unset";
  tdchwtype = "Unset";
  tdcresolution = 195.3125*picosecond;
  mdelay    = 0;
  pdelay    = 800;
  mcoincwin = 40;
  pcoincwin = 40;
  vcoincwin = 80;
  muonRateFactor = 1.;
  rewindTimeBins  =  1000000000;
  dataWindowMin   = -0.2;
  dataWindowMax   = 10.0;
  pileupWindowMin = -10.5;
  pileupWindowMax = 10.5;
  musrMode = 'D';
  overallBinDelay = 0;
  boolInfinitelyLowMuonRate = false;
  //  globTime  = 0.;
  
  int    tmpivar;
  float  tmpfvar;
  //  int    tmpCoincidenceMultimap[maxChannels][maxChannels]; 
  std::multimap<int, int> tmpCoincidenceMultimap;
  //  for (int i=0; i<maxChannels; i++) {
  //    for (int j=0; j<maxChannels; j++) {
  //      //  tmpCoincidenceMultimap[i][j]=0;
  //      tmpCoincidenceMultimap.insert(pair<int, int>("a", 1));
  //
  //    }
  //  }
  while (!feof(fSteeringFile)) {
    fgets(line,1000,fSteeringFile);
    if ((line[0]!='#')&&(line[0]!='\n')&&(line[0]!='\r')&&(line[0]!='$')&&(line[0]!='!')) {
      char tmpString0[200]="Unset";
      char tmpString1[200]="Unset";
      sscanf(&line[0],"%s",tmpString0);
      if (strncmp(tmpString0,"INSTRUMENT=",strlen("INSTRUMENT="))==0) {
	instrument=&line[strlen("INSTRUMENT=")];
      }
      else if (strncmp(tmpString0,"DESCRIPTION=",strlen("DESCRIPTION="))==0) {
	description=&line[strlen("DESCRIPTION=")];
      }
      else if (strncmp(tmpString0,"TYPE=",strlen("TYPE="))==0) {
	tdchwtype=&line[strlen("TYPE=")];
      }
      else if (strncmp(tmpString0,"RESOLUTION=",strlen("RESOLUTION="))==0) {
	sscanf(&line[strlen("RESOLUTION=")],"%g",&tmpfvar);
	tdcresolution = tmpfvar * picosecond;
	if (tdcresolution<1*picosecond) {
	  std::cout<<"\n !!! ERROR !!! Time resolution of TDC must be larger than 1 ps (you requested "<<tdcresolution/picosecond<<" ps)"<<std::endl;
	  std::cout<<"   ===> too low time resolution might lead to numerical problems in musrCounter ==> S T O P"<<std::endl;
	  exit(1);
	}
      }
      else if (strncmp(tmpString0,"MDELAY=",strlen("MDELAY="))==0) {
	sscanf(&line[strlen("MDELAY=")],"%d",&tmpivar);
	mdelay = tmpivar;
      }
      else if (strncmp(tmpString0,"PDELAY=",strlen("PDELAY="))==0) {
	sscanf(&line[strlen("PDELAY=")],"%d",&tmpivar);
	pdelay = tmpivar;
      }
      else if (strncmp(tmpString0,"MCOINCIDENCEW=",strlen("MCOINCIDENCEW="))==0) {
	sscanf(&line[strlen("MCOINCIDENCEW=")],"%d",&tmpivar);
	mcoincwin = tmpivar;
      }
      else if (strncmp(tmpString0,"PCOINCIDENCEW=",strlen("PCOINCIDENCEW="))==0) {
	sscanf(&line[strlen("PCOINCIDENCEW=")],"%d",&tmpivar);
	pcoincwin = tmpivar;
      }
      else if (strncmp(tmpString0,"VCOINCIDENCEW=",strlen("VCOINCIDENCEW="))==0) {
	sscanf(&line[strlen("VCOINCIDENCEW=")],"%d",&tmpivar);
	vcoincwin = tmpivar;
      }
      else if (strncmp(tmpString0,"MUONRATEFACTOR=",strlen("MUONRATEFACTOR="))==0) {
	sscanf(&line[strlen("MUONRATEFACTOR=")],"%g",&tmpfvar);
	muonRateFactor = tmpfvar;
      }
      else if (strncmp(tmpString0,"INFINITELYLOWMUONRATE=",strlen("INFINITELYLOWMUONRATE"))==0) {
	boolInfinitelyLowMuonRate=true;
      }
      else if (strncmp(tmpString0,"REWINDTIMEBINS=",strlen("REWINDTIMEBINS="))==0) {
	sscanf(&line[strlen("REWINDTIMEBINS=")],"%d",&tmpivar);
	rewindTimeBins = tmpivar; std::cout<<" JUHELAK: rewindTimeBins="<<rewindTimeBins<<std::endl;
      }
      else if (strncmp(tmpString0,"DATAWINDOWMIN=",strlen("DATAWINDOWMIN="))==0) {
	sscanf(&line[strlen("DATAWINDOWMIN=")],"%g",&tmpfvar);
	dataWindowMin = tmpfvar;
      }
      else if (strncmp(tmpString0,"DATAWINDOWMAX=",strlen("DATAWINDOWMAX="))==0) {
	sscanf(&line[strlen("DATAWINDOWMAX=")],"%g",&tmpfvar);
	dataWindowMax = tmpfvar;
      }
      else if (strncmp(tmpString0,"PILEUPWINDOWMIN=",strlen("PILEUPWINDOWMIN="))==0) {
	sscanf(&line[strlen("PILEUPWINDOWMIN=")],"%g",&tmpfvar);
	pileupWindowMin = tmpfvar;
      }
      else if (strncmp(tmpString0,"PILEUPWINDOWMAX=",strlen("PILEUPWINDOWMAX="))==0) {
	sscanf(&line[strlen("PILEUPWINDOWMAX=")],"%g",&tmpfvar);
	pileupWindowMax = tmpfvar;
      }
      else if (strncmp(tmpString0,"MUSRMODE=",strlen("MUSRMODE="))==0) {
	sscanf(&line[strlen("MUSRMODE=")],"%s",&tmpString1);
	musrMode = tmpString1[0];
      }
      else if (strncmp(tmpString0,"musrTH",strlen("musrTH"))==0) {
	// Definition of the histograms - either musrTH1D or musrTH2D
	int beginningOfHistoTitle=0, endOfHistoTitle =0;
	for (Int_t i=0; i<(int)strlen(line); i++) {
	  if (line[i]=='"') {
	    if (beginningOfHistoTitle==0) beginningOfHistoTitle=i+1;
	    else                          endOfHistoTitle=i-1;
	  }
	}
	char histoName[200];
	sscanf(&line[0],"%*s %s",histoName);
	char histoTitle[501];
	CopySubstring(line,beginningOfHistoTitle,endOfHistoTitle,histoTitle);
	int nrBinsX, nrBinsY;
	float minX, maxX, minY, maxY;
	char varXName[500];  char varYName[500];
	char furtherOption[200]; 
	float rot_ref_frequency, rot_ref_phase;
	if (strcmp(tmpString0,"musrTH1D")==0) {
	  sscanf(&line[endOfHistoTitle+2],"%d %g %g %s %s %g %g",&nrBinsX,&minX,&maxX,varXName,furtherOption,&rot_ref_frequency,&rot_ref_phase);
	}
	else if (strcmp(tmpString0,"musrTH2D")==0) {
	  sscanf(&line[endOfHistoTitle+2],"%d %g %g %d %g %g %s %s",&nrBinsX,&minX,&maxX,&nrBinsY,&minY,&maxY,varXName,varYName);
	}
	else {
	  std::cout<<" !!! ERROR:  keyword \""<<tmpString0<<"\" in the steering file not understood ==> S T O P !!!"<<std::endl;
	  exit(1);
	}
	std::string stringVariableName1 = varXName;
	std::string stringVariableName2 = varYName;
	variableMapType::iterator iter1 = variableMap.find(stringVariableName1);
	variableMapType::iterator iter2;
	if (iter1==variableMap.end()) {
	  std::cout<<" !!!  ERROR: Variable with the name \""<<stringVariableName1<<"\" not defined in variableMap ==> ignored, histogram not defined!   S T O P !!!"<<std::endl;
	  exit(1);
	}
	if (strcmp(tmpString0,"musrTH2D")==0) {
	  iter2 = variableMap.find(stringVariableName2);
	  if (iter2==variableMap.end()) {
	    std::cout<<" !!!  ERROR: Variable with the name \""<<stringVariableName2<<"\" not defined in variableMap ==> ignored, histogram not defined!   S T O P !!!"<<std::endl;
	    exit(1);
	  }
	}
	
	Double_t* pointerToVariable1 = iter1->second;
	if (strcmp(tmpString0,"musrTH1D")==0) {
	  //	  std::cout<<"pointerToVariable1="<<pointerToVariable1<<"  pointer to muDecayPosZ="<<&muDecayPosZ<<std::endl;
	  musrTH* myTH = new musrTH("1D",histoName,histoTitle,nrBinsX,minX,maxX,nrBinsY,minY,maxY,pointerToVariable1,NULL);
	  listOfAllHistograms1D.push_back(myTH);
	  if (strcmp(furtherOption,"rotreference")==0) {
	    myTH->SetRotatingReferenceFrame(rot_ref_frequency, rot_ref_phase);
	    hInfo->Fill(21, rot_ref_frequency);  // value may be overwritten - just the last rot. ref. frame will be saved
	    hInfo->Fill(22, rot_ref_phase);      // value may be overwritten - just the last rot. ref. frame will be saved
	  }
	}
	else {
	  Double_t* pointerToVariable2 = iter2->second;
	  musrTH* myTH = new musrTH("2D",histoName,histoTitle,nrBinsX,minX,maxX,nrBinsY,minY,maxY,pointerToVariable1,pointerToVariable2);
	  listOfAllHistograms2D.push_back(myTH);
	}
      } 
      else if  (strncmp(tmpString0,"humanDecayHistograms",strlen("humanDecayHistograms"))==0) {
	char motherHistoName[200];
	char motherHistoPileupName[200];
	sscanf(&line[0],"%*s %s %s",motherHistoName,motherHistoPileupName);
	std::cout<<" Defining humanDecayHistograms (motherHistoName="<<motherHistoName<<")"<<std::endl;
	for(std::list<musrTH*>::const_iterator it = listOfAllHistograms1D.begin(); it != listOfAllHistograms1D.end(); ++it) {
	  if ((*it)->IsThisHistoNamed(motherHistoName)) {
	    motherOfHumanDecayHistograms = *it;
	  }
	  if ((*it)->IsThisHistoNamed(motherHistoPileupName)) {
	    motherOfHumanDecayPileupHistograms = *it;
	  }
	}
	if (motherOfHumanDecayHistograms==NULL) {
	  std::cout<<"\n\nmusrAnalysis::ReadInInputParameters  Mother histogram (motherOfHumanDecayHistograms)"
		   <<"of the \"Human decay histogram\" not found!"<<std::endl;
	  std::cout<<"                 "<<line<<std::endl;
	  std::cout<<"  ==> S T O P "<<std::endl;
	  exit(1);
	}
	if (motherOfHumanDecayPileupHistograms==NULL) {
	  std::cout<<"\n\nmusrAnalysis::ReadInInputParameters  Mother histogram (motherOfHumanDecayPileupHistograms)"
		   <<"of the \"Pileup human decay histogram\" not found!"<<std::endl;
	  std::cout<<"                 "<<line<<std::endl;
	  std::cout<<"  ==> S T O P "<<std::endl;
	  exit(1);
	}

	char* pch0 = strstr(line,motherHistoName)+strlen(motherHistoName);
	char* pch  = strstr(pch0,motherHistoPileupName)+strlen(motherHistoPileupName);

	int nscan; int N1; char NAME1[100];
	nscan = sscanf(pch,"%d %s",&N1,NAME1);
	//	indexHuman=0;
	Int_t iBinHuman=0;
	Int_t nBinHuman=0;
	do {
	  Bool_t nameDefined=false;
	  for (humanDecayMapType::const_iterator it = humanDecayMap.begin(); it!=humanDecayMap.end(); ++it) {
	    if ((it->second)==(std::string) NAME1) {
	      nameDefined = true;
	      iBinHuman   = it->first;
	    }
	  }
	  if (!nameDefined) {
	    nBinHuman++;
	    iBinHuman=nBinHuman;
	    humanDecayMap.insert(std::pair<Int_t,std::string>(iBinHuman,NAME1));
	  }
	  humanDecayMultimap.insert(std::pair<Int_t,Int_t>(iBinHuman,N1));
	  char* pch2 = strstr(pch,NAME1)+strlen(NAME1);
	  pch=pch2;
	  nscan = sscanf(pch,"%d %s",&N1,NAME1);
	  //	  indexHuman++;
	} while (nscan==2);
	
	Int_t nrBinsHumanHist=humanDecayMap.size();
	//	Int_t nrBinsHumanHist=0;
	//	std::string oldString="blaBlaoldStringUndefiNED";
	//	for (humanDecayMultimapType::const_iterator it = humanDecayMultimap.begin(); it!=humanDecayMultimap.end(); ++it) {
	//	  if (oldString!=(it->first)) nrBinsHumanHist++;
	//	}
	std::cout<<"\nnrBinsHumanHist="<<nrBinsHumanHist<<std::endl;
	//	humanDecayHistograms = new musrTH("1D","humanDecayHistograms","Where the muons stop and decay",indexHuman,0,float(indexHuman),0,0,0,NULL,NULL);
	humanDecayHistograms = new musrTH("1D","humanDecayHistograms","Where the muons stop and decay",nrBinsHumanHist,0,float(nrBinsHumanHist),0,0,0,NULL,NULL);
	humanDecayPileupHistograms = new musrTH("1D","humanDecayPileupHistograms","Where the muons contr. to pileup stop and decay",nrBinsHumanHist,0,float(nrBinsHumanHist),0,0,0,NULL,NULL);
	//	for (Int_t k=0; k<indexHuman; k++) {
	//	  humanDecayHistograms->SetBinLabel1D(k+1,sHumanDecayArray[k]);
	//	}
	//	oldString="blaBlaoldStringUndefiNED";
	for (humanDecayMapType::const_iterator it = humanDecayMap.begin(); it!=humanDecayMap.end(); ++it) {
	  humanDecayHistograms        -> SetBinLabel1D(it->first,it->second);
	  humanDecayPileupHistograms -> SetBinLabel1D(it->first,it->second);
	}
      }
      else if  (strncmp(tmpString0,"condition",strlen("condition"))==0) {
	// Selection of the predefined conditions to be applied
	int  iConditionTMP;
	char conditionNameTMP[200];
	sscanf(&line[0],"%*s %d %s",&iConditionTMP,conditionNameTMP);
	if (iConditionTMP<0) {
	  std::cout<<" !!!  ERROR: Condition number must be >= 0  (not "<<iConditionTMP<<")"<< std::endl; 
	  exit(1);
	}
	if (iConditionTMP>=nrConditions) {
	  std::cout<<" !!!  ERROR: Condition number must be < "<<nrConditions<<"  (not "<<iConditionTMP<<")"<< std::endl; 
	  exit(1);
	}
	if      (strcmp(conditionNameTMP,"alwaysTrue")==0)              conditionMap[iConditionTMP]=&alwaysTrue;
	else if (strcmp(conditionNameTMP,"oncePerEvent")==0)            conditionMap[iConditionTMP]=&oncePerEvent;
	else if (strcmp(conditionNameTMP,"muonDecayedInSample_gen")==0) conditionMap[iConditionTMP]=&muonDecayedInSample_gen;
	else if (strcmp(conditionNameTMP,"muonTriggered_gen")==0)       conditionMap[iConditionTMP]=&muonTriggered_gen;
	else if (strcmp(conditionNameTMP,"muonTriggered_det")==0)       conditionMap[iConditionTMP]=&muonTriggered_det;
	else if (strcmp(conditionNameTMP,"positronHit_det")==0)         conditionMap[iConditionTMP]=&positronHit_det;
	else if (strcmp(conditionNameTMP,"goodEvent_det")==0)           conditionMap[iConditionTMP]=&goodEvent_det;
	else if (strcmp(conditionNameTMP,"goodEvent_gen")==0)           conditionMap[iConditionTMP]=&goodEvent_gen;
	else if (strcmp(conditionNameTMP,"goodEvent_det_AND_goodEvent_gen")==0)  conditionMap[iConditionTMP]=&goodEvent_det_AND_goodEvent_gen;
	else if (strcmp(conditionNameTMP,"pileupEventCandidate")==0)    conditionMap[iConditionTMP]=&pileupEventCandidate;
	else if (strcmp(conditionNameTMP,"pileupEvent")==0)             conditionMap[iConditionTMP]=&pileupEvent;
	else if (strcmp(conditionNameTMP,"goodEvent_det_AND_muonDecayedInSample_gen")==0) conditionMap[iConditionTMP]=&goodEvent_det_AND_muonDecayedInSample_gen;
	else {
	  std::cout<<" !!!  ERROR: Condition of the name \""<<conditionNameTMP<<"\" not predefined ==> Add it in the musrAnalysis.cxx   S T O P !!!"<<std::endl;
	  exit(1);
	}
      }
      else if  (strcmp(tmpString0,"draw")==0) {
	// Which histograms should be drawn
	int  iConditionTMP;
	char histoNameTMP[200];
	sscanf(&line[0],"%*s %s %d",histoNameTMP,&iConditionTMP);
	if (iConditionTMP<-1) {
	  std::cout<<" !!!  ERROR: draw: condition number must be >= -1  (not "<<iConditionTMP<<")"<< std::endl; 
	  exit(1);
	}
	if (iConditionTMP>=nrConditions) {
	  std::cout<<" !!!  ERROR: draw: condition number must be < "<<nrConditions<<"  (not "<<iConditionTMP<<")"<< std::endl; 
	  exit(1);
	}
	for(std::list<musrTH*>::const_iterator it = listOfAllHistograms1D.begin(); it != listOfAllHistograms1D.end(); ++it) {
	  if ((*it)->IsThisHistoNamed(histoNameTMP)) {
	    (*it)->SetDrawListOfHistograms(iConditionTMP);
	  }
	}
	for(std::list<musrTH*>::const_iterator it = listOfAllHistograms2D.begin(); it != listOfAllHistograms2D.end(); ++it) {
	  if ((*it)->IsThisHistoNamed(histoNameTMP)) {
	    (*it)->SetDrawListOfHistograms(iConditionTMP);
	  }
	}
	if (humanDecayHistograms->IsThisHistoNamed(histoNameTMP)) {
	  humanDecayHistograms->SetDrawListOfHistograms(iConditionTMP);
	}
	if (humanDecayPileupHistograms->IsThisHistoNamed(histoNameTMP)) {
	  humanDecayPileupHistograms->SetDrawListOfHistograms(iConditionTMP);
	}
      }
      else if  (strcmp(tmpString0,"counterPhaseShifts")==0) {
	int nscan; int N1; float PHASE_SHIFT; char NAME1[100]; char NAME2[100];
	char *pch = line + strlen("counterPhaseShifts");
	do {
	  nscan = sscanf(pch,"%d %g",&N1,&PHASE_SHIFT);
	  phaseShiftMap.insert(std::pair<int,Double_t>(N1,PHASE_SHIFT/180*pi));
	  std::cout<<"N1="<<N1<<", PHASE_SHIFT="<<PHASE_SHIFT/180<<" pi"<<std::endl;
	  nscan = sscanf(pch,"%s %s",NAME1,NAME2);
	  char* pch1 = strstr(pch ,NAME1)+strlen(NAME1);
	  char* pch2 = strstr(pch1,NAME2)+strlen(NAME2);
	  pch=pch2;
	  nscan = sscanf(pch,"%d %g",&N1,&PHASE_SHIFT);
	} while (nscan==2);
      }
      else if  (strcmp(tmpString0,"fit")==0) {
	char histoName[100]; char functionName[100]; float xMin; float xMax; float p0, p1, p2, p3, p4, p5, p6;
	sscanf(&line[0],"%*s %s %s %g %g %g %g %g %g %g",histoName,functionName,&xMin,&xMax,&p0,&p1,&p2,&p3,&p4,&p5,&p6);

	TF1 *funct;
	
	if (strcmp(functionName,"funct1")==0) {
	  funct = new TF1("funct1","[3]*exp(([4]-x)/2.19703)*(1+[2]*cos(x*[0]+[1]))");
	  funct -> SetParameter(0,p0);
	  funct -> SetParameter(1,p1);
	  funct -> SetParameter(2,p2);
	  funct -> SetParameter(3,p3);
	  funct -> FixParameter(4,xMin);
	}
	else if (strcmp(functionName,"funct2")==0) {
	  funct = new TF1("funct2","[3]*exp(([4]-x)/2.19703)*(1+[2]*cos(x*[0]+[1]))+[5]");
	  funct -> SetParameter(0,p0);
	  funct -> SetParameter(1,p1);
	  funct -> SetParameter(2,p2);
	  funct -> SetParameter(3,p3);
	  funct -> FixParameter(4,xMin);
	  funct -> SetParameter(5,p4);
	}
	else if (strcmp(functionName,"funct3")==0) {
	  funct = new TF1("funct3","[3]*exp(([4]-x)/2.19703)*(1+[2]*cos(x*[0]+[1]))+[5]");
	  funct -> SetParameter(0,p0);
	  funct -> SetParameter(1,p1);
	  funct -> SetParameter(2,p2);
	  funct -> SetParameter(3,p3);
	  funct -> FixParameter(4,xMin);
	  funct -> FixParameter(5,p4);
	}
	else if (strcmp(functionName,"funct4")==0) {
	  funct = new TF1("funct4","[3]*exp(-x/2.19703)*(1+[2]*cos(x*[0]+[1]))+[4]");
	  funct -> SetParameter(0,p0);
	  funct -> SetParameter(1,p1);
	  funct -> SetParameter(2,p2);
	  funct -> SetParameter(3,p3);
	  funct -> SetParameter(4,p4);
	}
	else if (strcmp(functionName,"pol0")==0) {
	  funct = new TF1("pol0","pol0");
	}
	else if (strcmp(functionName,"simpleExpoPLUSconst")==0) {
	  funct = new TF1("simpleExpoPLUSconst","[0]*exp(-x/2.19703)+[1]");
	  funct -> SetParameter(0,p0);
	  funct -> SetParameter(1,p1);
	}
	else if (strcmp(functionName,"rotFrameTime20")==0) {
	  funct = new TF1("rotFrameTime20","[2]*exp(-x/2.19703)*cos(x*[0]+[1]) ");
	  funct -> SetParameter(0,p0);
	  funct -> SetParameter(1,p1);
	  funct -> SetParameter(2,p2);
	}
	else {
	  std::cout<<"musrAnalysis::ReadInInputParameters:  function \""<<functionName<<"\" not defined! ==> S T O P"<<std::endl;
	  exit(1);
	}
	
	musrTH* tmpHistograms;
	for(std::list<musrTH*>::const_iterator it = listOfAllHistograms1D.begin(); it != listOfAllHistograms1D.end(); ++it) {
	  if ((*it)->IsThisHistoNamed(histoName)) {
	    tmpHistograms = *it;
	  }
	}
	if (tmpHistograms==NULL) {
	  std::cout<<"\n\nmusrAnalysis::ReadInInputParameters: Fit required, but the corresponding histogram "
		   <<"not found!"<<std::endl;
	  std::cout<<"                 "<<line<<std::endl;
	  std::cout<<"  ==> S T O P "<<std::endl;
	  exit(1);
	}

	tmpHistograms->AssignFunction(funct, xMin, xMax);
      }
      
      else { 
	//------------------------------------------------------------------
	// Second part of the configuration file
	//
	// Replace spaces between quotation marks by '_' 
	// and remove quatation marks and semicolon by ' '. 
	// Remember position of the beginning and end of the defintion of the coincidence counters
	int nrOfSemicolons=0, beginningOfCoincidenceArray=0, endOfCoincidenceArray =0;
	for (Int_t i=0; i<(int)strlen(line); i++) {
	  if (line[i]=='"') {
	    line[i] = ' ';
	    Int_t ii=i;
	    do {
	      ii++;
	      if (isspace(line[ii])) line[ii] = '_';
	    } while (line[ii]!='"');
	    line[ii]= ' ';
	  }
	  else if (line[i] == ';')   {
	    line[i] = ' ';
	    nrOfSemicolons++; 
	    if (nrOfSemicolons==5) beginningOfCoincidenceArray=i;
	    if (nrOfSemicolons==6)       endOfCoincidenceArray=i;
	  }
	}

	// Read in the values from the configuration file 
	int tmpChannel=-1, tmpTimeShift=0; 
	char tmpName[200], tmpType[200];
	float tmpThreshold =0.;
	sscanf(&line[0],"%d %s %s %g %d",&tmpChannel,tmpName,tmpType,&tmpThreshold,&tmpTimeShift);
	musrCounter* counter = new musrCounter(tmpChannel,tmpName,tmpType[0],tmpThreshold,tmpTimeShift);
	if      (tmpType[0]=='M') {mCounter               =counter; allCounterMap[tmpChannel]=counter; overallBinDelay=tmpTimeShift;}
	else if (tmpType[0]=='P') {pCounterMap[tmpChannel]=counter; allCounterMap[tmpChannel]=counter;}
	else if (tmpType[0]=='K') {kCounterMap[tmpChannel]=counter; allCounterMap[tmpChannel]=counter;}
	else if (tmpType[0]=='V') {vCounterMap[tmpChannel]=counter; allCounterMap[tmpChannel]=counter;}

	// Read in the values for the coincidence detectors
	// (Tricky while the number of detectors in coincidence is not known)
	if ((endOfCoincidenceArray>0)&&(endOfCoincidenceArray!=beginningOfCoincidenceArray)) {
	  char tmpCoinc[100]; strcpy(tmpCoinc,"987654321");
	  int jj=beginningOfCoincidenceArray;
	  do {
	    do {jj++;} while (isspace(line[jj]));  // skip space characters;
	    if (jj>endOfCoincidenceArray) {break;}
	    sscanf(&line[jj],"%s",tmpCoinc);
	    jj+=strlen(tmpCoinc);
	    int iCoinc=atoi(tmpCoinc);
	    if (iCoinc==987654321) continue;
//ckdel 20.9.2010	    if ((maxChannels<=tmpChannel)||(maxChannels<=iCoinc)) {
//ckdel 20.9.2010	      std::cout<<" Channel number ("<<tmpChannel<<") or channel coincidence required ("<<iCoinc
//ckdel 20.9.2010		       <<") larger than maxChannels ("<<maxChannels<<") ==> S T O P"<<std::endl;
//ckdel 20.9.2010	      exit(1);
//ckdel 20.9.2010	    }
	    tmpCoincidenceMultimap.insert(std::pair<int, int>(tmpChannel,iCoinc));
	    //	    tmpCoincidenceMultimap[tmpChannel][abs(iCoinc)]=iCoinc;
	  } while (jj<endOfCoincidenceArray);
	  // Read in the last part of the configuration line (after the coincidence definition)
	  if (tmpType[0]=='P') {
	    char tmpHistoName[500];     strcpy(tmpHistoName,   "Unset");
	    char tmpHistoNameAdd[500];  strcpy(tmpHistoNameAdd,"Unset");
	    int tmpT0=0, tmpT1=0, tmpT2=0, tmpHistoaddNr=0;
	    sscanf(&line[endOfCoincidenceArray],"%s %d %d %d %d %s",tmpHistoName,&tmpT0,&tmpT1,&tmpT2,&tmpHistoaddNr,tmpHistoNameAdd);
	    counter->SetTDChistogram(tmpHistoName,tmpT0,tmpT1,tmpT2,tmpHistoaddNr,tmpHistoNameAdd);
	  }
	}
      }
      //      std::cout << line;
    }
  }
  std::cout << "INSTRUMENT=" << instrument << std::endl;
  std::cout << "DESCRIPTION=" << description << std::endl;
  std::cout << "TYPE=" << tdchwtype << std::endl;
  std::cout << "RESOLUTION="<<tdcresolution<<" ps"<<std::endl;
  std::cout << "MDELAY="<<mdelay<<std::endl;
  std::cout << "PDELAY="<<pdelay<<std::endl;
  std::cout << "MCOINCIDENCEW="<<mcoincwin<<std::endl;
  std::cout << "PCOINCIDENCEW="<<pcoincwin<<std::endl;
  std::cout << "VCOINCIDENCEW="<<vcoincwin<<std::endl;
  std::cout << "DATAWINDOWMIN="<<dataWindowMin<<std::endl;
  std::cout << "DATAWINDOWMAX="<<dataWindowMax<<std::endl;
  std::cout << "PILEUPWINDOWMIN="<<pileupWindowMin<<std::endl;
  std::cout << "PILEUPWINDOWMAX="<<pileupWindowMax<<std::endl;
  std::cout << "REWINDTIMEBINS="<<rewindTimeBins<<std::endl;


  // Loop over all counters and set maximum coincidence window (maxCoincidenceTimeWindow):
  // The "K" counters are tricky - they have to be set according to the detector, in which they
  // are in coincidence.
 
  pileupWindowBinMin = Long64_t(pileupWindowMin/tdcresolution);
  pileupWindowBinMax = Long64_t(pileupWindowMax/tdcresolution);
  dataWindowBinMin   = Long64_t(  dataWindowMin/tdcresolution);
  dataWindowBinMax   = Long64_t(  dataWindowMax/tdcresolution);
  for (counterMapType::const_iterator it = allCounterMap.begin(); it!=allCounterMap.end(); ++it) {
    char DetectorType = (it->second)->GetCounterType();
    if      (DetectorType=='M') {
      (it->second)->SetMaxCoincidenceTimeWindow(pileupWindowBinMin);
      (it->second)->SetCoincidenceTimeWindow(pileupWindowBinMin,pileupWindowBinMax);
      (it->second)->SetCoincidenceTimeWindowOfAllCoincidenceDetectors(-mcoincwin,-mcoincwin,mcoincwin);
    }
    else if (DetectorType=='P') {
      (it->second)->SetMaxCoincidenceTimeWindow(dataWindowBinMin);
      (it->second)->SetCoincidenceTimeWindow(dataWindowBinMin,dataWindowBinMax);
      (it->second)->SetCoincidenceTimeWindowOfAllCoincidenceDetectors(-pcoincwin,-pcoincwin,pcoincwin);
    }
    else if (DetectorType=='V') {
      (it->second)->SetMaxCoincidenceTimeWindow(-vcoincwin);
      (it->second)->SetCoincidenceTimeWindow(-vcoincwin,vcoincwin);
    }
    int iChanNr = (it->second)->GetCounterNr();
    for (std::multimap<int,int>::iterator itCoinc = tmpCoincidenceMultimap.begin(); itCoinc != tmpCoincidenceMultimap.end();  ++itCoinc)  {
      if ((*itCoinc).first == iChanNr) {
	int iChn2 = itCoinc->second;
	counterMapType::const_iterator itTwo = allCounterMap.find(abs(iChn2));
        if (itTwo==allCounterMap.end()) { 
	  std::cout<<" Pointer to coincidence counter ("<<iChn2<<") not found!  This should never happen!!! ==> S T O P"<<std::endl;
	  exit(1);
	}
	else {
	  (it->second)->SetCoincidenceCounter(itTwo->second,iChn2);
	}
      }
    }
  }

    //    for (int j=0; j<maxChannels; j++) {
    //    for (multimap<int,int>::iterator itCoinc = tmpCoincidenceMultimap.begin(); itCoinc != tmpCoincidenceMultimap.end();  ++itCoinc)  {
    //      std::cout << "   houby  [" << (*itCoinc).first << ", " << (*itCoinc).second << "]" << std::endl;
 
    //
    //      if (tmpCoincidenceMultimap[iChanNr][j]!=0) {
//      int j = (*itCoinc).first;
//        counterMapType::const_iterator itTwo = allCounterMap.find(j);
//        if (itTwo==allCounterMap.end()) { 
//	  std::cout<<" Pointer to coincidence counter ("<<j<<") not found!  This should never happen!!! ==> S T O P"<<std::endl;
//	  exit(1);
//	}
//	else {
//	  //	  int iChn = tmpCoincidenceMultimap[iChanNr][j];
//	  int iChn = itCoinc->second;
//	  (it->second)->SetCoincidenceCounter(itTwo->second,iChn);
//	}
//	//      }
//    }
//  }


  
 //      
  //
  //      char tmpString1[100]="Unset",tmpString2[100]="Unset",tmpString3[100]="Unset";
  //      sscanf(&line[0],"%*s %s %s",tmpString1,tmpString2);

}

//================================================================

void musrAnalysis::CreateHistograms() {
  std::cout<<"musrAnalysis::CreateHistograms()"<<std::endl;
  // safeTimeWindow=2*(std::max(datawindow,pileupwindow));
  safeTimeWindow = (dataWindowMax>pileupWindowMax)? 3*dataWindowMax: 3*pileupWindowMax;
  std::cout<<"musrAnalysis::CreateHistograms():  safeTimeWindow = "<<safeTimeWindow<<" (bins="<<Long64_t(safeTimeWindow/tdcresolution)<<")"<<std::endl;
  currentTime=0.;
  //  nextEventTime=0.;
  nextUnfilledEventTime=0.;
  numberOfRewinds=0;
  numberOfGoodMuons=0;
  lastPreprocessedEntry=-1;
  //  currentEventID =0;
  fChain->GetEntry(1); Double_t fieldValue = fieldNomVal[0]; if (nFieldNomVal>1) fieldValue += fieldNomVal[1];
  //  omega = 851.372*fabs(fieldValue);
  //  omega = 851.372*fieldValue;     // value set originally in this program
  //  omega = 850.62*fieldValue;  // value used in Geant ?
  omega = 851.610*fieldValue;       // value from the fits of the data


  hInfo->Fill(1, fieldValue);
  hInfo->Fill(6, runID);
  hInfo->Fill(7, hGeantParameters->GetBinContent(7));
  std::cout<<"musrAnalysis::CreateHistograms():   fieldValue="<<fieldValue<<"T,  omega="<<omega<<std::endl;  
}
//================================================================

void musrAnalysis::SaveHistograms(char* runChar, char* v1190FileName) {

  for(std::list<musrTH*>::const_iterator it = listOfAllHistograms1D.begin(); it != listOfAllHistograms1D.end(); ++it) {
    //    (*it)->ListHistograms();
    (*it)->FitHistogramsIfRequired(omega);
    (*it)->DrawTH1DdrawList("");
    //    (*it)->DrawTH1D("hist",0);
    //    (*it)->DrawTH1D("hist",1);
  }
  for(std::list<musrTH*>::const_iterator it = listOfAllHistograms2D.begin(); it != listOfAllHistograms2D.end(); ++it) {
    (*it)->DrawTH2DdrawList("hist");
  }
  for (counterMapType::const_iterator it = pCounterMap.begin(); it!=pCounterMap.end(); ++it) {
    (it->second)->DrawTDChistogram();
  }
  if (humanDecayHistograms) {
    //    humanDecayHistograms->FillHumanDecayArray(motherOfHumanDecayHistograms,indexHuman,nHumanDecayArray);
    humanDecayHistograms->FillHumanDecayArray(motherOfHumanDecayHistograms,humanDecayMap,humanDecayMultimap);
    humanDecayHistograms->DrawTH1DdrawList("text");
    humanDecayPileupHistograms->FillHumanDecayArray(motherOfHumanDecayPileupHistograms,humanDecayMap,humanDecayMultimap);
    humanDecayPileupHistograms->DrawTH1DdrawList("text");
  }

  Long64_t numberOfMuonCandidates = mCounter->GetNumberOfMuonCandidates();
  hInfo->Fill(5,(Double_t) numberOfMuonCandidates);  // number of "triggers"

  //==============================
  //  Write out the histograms
  std::cout<<"musrAnalysis::SaveHistograms()"<<std::endl;
  TDirectory *dir = gDirectory;
  char HistogramOutputFileName[500];
  sprintf(HistogramOutputFileName,"data/his_%s_%s.root",runChar,v1190FileName);
  TFile *fHistograms = new TFile(HistogramOutputFileName,"recreate");
  TIter next(dir->GetList());
  while ( TObject *obj = next() ) {
    if (obj->InheritsFrom(TH1::Class()) ) {obj->Write();}
  }
  fHistograms->Close();
  delete fHistograms;


  //==============================

  std::cout<<"==================================================================="<<std::endl;
  std::cout<<"Number of raw trigger counts (ignoring pileup gate, vetos and coincidences): "<<numberOfMuonCandidates<<std::endl;
  std::cout<<"Number of triggered events (i.e. only good \"muons\"): "<<numberOfGoodMuons<<std::endl;
  Double_t durationOfExperiment = (numberOfRewinds*rewindTimeBins*tdcresolution+currentTime)/1000000;
  std::cout<<"Duration of the \"experiment\": "<<durationOfExperiment<<" second"<<std::endl;
  //  std::cout<<"                        (numberOfRewinds="<<numberOfRewinds<<")"<<std::endl;
  std::cout<<"In this run, MUONRATEFACTOR was set to "<<muonRateFactor<<", and the raw trigger rate was "
	   <<numberOfMuonCandidates/durationOfExperiment<<" muons/second"<<std::endl;
  std::cout<<"               To get event rate of 30000 events/second, set MUONRATEFACTOR to "
	   << muonRateFactor*( (numberOfMuonCandidates/durationOfExperiment)/30000 )<<std::endl;
  std::cout<<"========================== E N D =================================="<<std::endl;
}

//================================================================

void musrAnalysis::AnalyseEvent(Long64_t iiiEntry) {
  //  std::cout<<"\n______________________________________________________________________________________________\n";
  //  std::cout<<"musrAnalysis::AnalyseEvent("<<iiiEntry<<")   currentTime="<<currentTime<<" (bin="<<Long64_t(currentTime/tdcresolution)<<")"<<std::endl;
  

  // Reset some initial variables

  // Loop over several next event and preprocess them (i.e. fill
  // them into the lists/maps of the class musrCounter).
  while (((iiiEntry>lastPreprocessedEntry)||(((nextUnfilledEventTime-currentTime)<safeTimeWindow))&&(!boolInfinitelyLowMuonRate)) && (lastPreprocessedEntry+1<nentries)) {
    Double_t deltaT = PreprocessEvent(lastPreprocessedEntry+1);
    nextUnfilledEventTime+=deltaT;
  };
  
  fChain->GetEntry(iiiEntry);  InitialiseEvent();
  
  // Loop over all interesting "moments", which are:
  //    1) any hit in the muon counter
  //    2) any new event (even if nothing was detected in the muon counter)

  Long64_t currentEventID = eventID;
  Double_t timeOfM0Hit;
  Long64_t evID;
  Bool_t thisEventAlreadyFilled = false;

  //  std::cout<<"___________________\n";
  //  for (counterMapType::const_iterator it = allCounterMap.begin(); it!=allCounterMap.end(); ++it) {
  //    (*it).second->myPrintThisCounter(eventID);
  //  }
  FillHistograms(iiiEntry);


  currentTime+=timeToNextEvent*muonRateFactor;
  Long64_t currentTimeBin = Long64_t(currentTime/tdcresolution);
  if (!boolInfinitelyLowMuonRate) RemoveOldHitsFromCounters(currentTimeBin-1);              // Remove all hits that can not influence the next event
  else                            RemoveOldHitsFromCounters(std::numeric_limits<Long64_t>::max()/2);
  if (currentTimeBin>rewindTimeBins) {  // When the time variables start to be too large, rewind the time info everywhere; 
    RewindAllTimeInfo(rewindTimeBins);
  }
}

//================================================================

//void musrAnalysis::RemoveOldHitsFromCounters(Double_t timeLimit) {
void musrAnalysis::RemoveOldHitsFromCounters(Long64_t timeBinLimit) {
  // Loop over all Counters and remove hits that happen at or before the "timeBinLimit".
  // std::cout<<"musrAnalysis::RemoveOldHitsFromCounters: timeBinLimit="<<timeBinLimit<<std::endl;
  for (counterMapType::const_iterator it = allCounterMap.begin(); it!=allCounterMap.end(); ++it) {
    (*it).second->RemoveHitsInCounter(timeBinLimit);
  }
}

//================================================================

//void musrAnalysis::RewindAllTimeInfo(Double_t timeToRewind) {
void musrAnalysis::RewindAllTimeInfo(Long64_t timeBinsToRewind) {
  currentTime           -= timeBinsToRewind*tdcresolution;
  nextUnfilledEventTime -= timeBinsToRewind*tdcresolution;
  numberOfRewinds++;
  // Loop over all timing information and rewind it by "timeToRewind";
  //  std::cout<<"musrAnalysis::RewindAllTimeInfo()"<<std::endl;
  for (counterMapType::const_iterator it = allCounterMap.begin(); it!=allCounterMap.end(); ++it) {
    (*it).second->RewindHitsInCounter(timeBinsToRewind);
  }
}

//================================================================

void musrAnalysis::FillHistograms(Int_t iiiEntry) {
  //  std::cout<<"musrAnalysis::FillHistograms() event="<<eventID<<" , bool="<<generatedInfo<<","<<detectorInfo<<std::endl;

  oncePerEvent = true;
  Bool_t mCounterHitExistsForThisEventID = false;
  Bool_t pCounterHitExistsForThisEventID = false;
  Long64_t timeBinOfThePreviouslyProcessedHit = -100000000;
  Long64_t timeBin0                    = -100000000;
  Long64_t timeBin1                    = -100000000;
  Long64_t timeBin2                    = -100000000;
  Int_t    kEntry                      = 0;
  Int_t    posEntry                      = 0;
  Int_t    idetM                       = 0;
  Int_t    idetM_ID                    = 0;
  Double_t idetM_edep                  = 0.;
  Int_t    idetP                       = 0;
  Int_t    idetP_ID                    = 0;
  Double_t idetP_edep                  = 0.;
  Bool_t   doubleHitM                  = false;
  Bool_t   doubleHitP                  = false;

  //  std::cout<<" FillHistograms: timeBinOfThePreviouslyProcessedHit = "<<timeBinOfThePreviouslyProcessedHit<<std::endl;
  //  mCounterHitExistsForThisEventID = mCounter->GetNextGoodHitInThisEvent(eventID,timeBinOfThePreviouslyProcessedHit,0,'M',timeBin0,kEntry,idetM,idetM_ID,idetM_edep,doubleHitM);  
  mCounterHitExistsForThisEventID = mCounter->GetNextGoodMuon(eventID,timeBinOfThePreviouslyProcessedHit,timeBin0,kEntry,idetM,idetM_ID,idetM_edep,doubleHitM);
  timeBinOfThePreviouslyProcessedHit = timeBin0;

  //___________________________________________________________ 

  do {
    //    std::cout<<"                                                   timeBin0 ="<<timeBin0<<std::endl;
    // Check whether there was good hit in the Positron counter
    //    Long64_t dataBinMin = (mCounterHitExistsForThisEventID) ? timeBin0+dataWindowBinMin : timeBinOfThePreviouslyProcessedHit-100000000;
    //    Long64_t dataBinMax = (mCounterHitExistsForThisEventID) ? timeBin0+dataWindowBinMax : timeBinOfThePreviouslyProcessedHit+100000000;
    pileup_muDecayDetID_double = -1001;
    if (mCounterHitExistsForThisEventID) {
      numberOfGoodMuons++;
      Long64_t dataBinMin = timeBin0+dataWindowBinMin;
      Long64_t dataBinMax = timeBin0+dataWindowBinMax;
      pCounterHitExistsForThisEventID = PositronCounterHit(eventID,dataBinMin,dataBinMax,timeBin1,timeBin2,posEntry,idetP,idetP_ID,idetP_edep,doubleHitP);
      //cDEL      if (pCounterHitExistsForThisEventID) std::cout<<"        timeBin1-timeBin2 = "<<timeBin1<<"-"<<timeBin2<<"="<<timeBin1-timeBin2 <<std::endl;
      if (pCounterHitExistsForThisEventID && (kEntry>0)&&(posEntry>0)&&(kEntry!=posEntry)) {
	// This must be a pileup event (positron counter hit comes from the different event than the muon counter hit)
	fChain->GetEntry(posEntry);
	pileup_muDecayDetID_double = muDecayDetID;
	//	if (pileup_muDecayDetID_double==-1000) {
	//	  std::cout<<"DEBUG:  pileup_muDecayDetID_double==-1000, posEntry="<<posEntry<<", eventID="<<eventID<<", idetP_edep="<<idetP_edep<<", idetP="<<idetP<<",  idetP_ID="<<idetP_ID<<std::endl;
	//	}
	fChain->GetEntry(iiiEntry);
      }
    }
    //    if (doubleHitP) {
    //      std::cout<<"yyyyyyyyyyyyyyyyy Total double hit nr="<<testIVar1++<<std::endl;
    //    }

    // CALCULATE VARIABLES
    // Initial positron
    pos_Trans_Momentum = sqrt(posIniMomX*posIniMomX+posIniMomY*posIniMomY);
    pos_Momentum = sqrt(pos_Trans_Momentum*pos_Trans_Momentum+posIniMomZ*posIniMomZ);
    pos_Radius = pos_Trans_Momentum/(-BFieldAtDecay_Bz)/0.3;
    pos_Theta = myAtan2(pos_Trans_Momentum,posIniMomZ);
    pos_Phi   = myAtan2(posIniMomY,posIniMomX);
    // Detector info
    det_m0edep = (mCounterHitExistsForThisEventID) ? idetM_edep : -1000.;
    //   det_time0 = timeBin0*tdcresolution;
    //    det_time1 = timeBin1*tdcresolution;
    det_time10 = ((timeBin0==-100000000)|| (timeBin1==-100000000)) ? -100000000 : (timeBin1-timeBin0)*tdcresolution;
    det_time20 = ((timeBin0==-100000000)|| (timeBin2==-100000000)) ? -100000000 : (timeBin2-timeBin0)*tdcresolution;
    //cks    if ((det_time10>-1)&&(det_time10<0.01)) {std::cout<<" eventID="<<eventID<<"   det_time10="<< det_time10<<std::endl;}
    gen_time10 = muDecayTime-muM0Time;
    det_time10_MINUS_gen_time10 = (det_time10 - gen_time10)/picosecond;
    det_posEdep = (pCounterHitExistsForThisEventID) ?  idetP_edep : -1000.;
    //    std::cout<<eventID<<"  det_time10="<<det_time10<<" t1="<<timeBin1*tdcresolution<<" t0="<<timeBin0*tdcresolution<<"  gen_time10="<<gen_time10<<std::endl;
    
    // CALCULATE CONDITIONS
    alwaysTrue = true;
    muonDecayedInSample_gen = muDecayDetID==201;
    muonTriggered_gen = muM0Time > -1000;
    muonTriggered_det = mCounterHitExistsForThisEventID;
    positronHit_det   = pCounterHitExistsForThisEventID;
    goodEvent_det     = muonTriggered_det && positronHit_det;
    goodEvent_gen     = (muDecayTime>-999)&&(muM0Time>-999);
    goodEvent_det_AND_goodEvent_gen = goodEvent_det && goodEvent_gen;
    pileupEventCandidate = ((kEntry>0)&&(posEntry>0)&&(kEntry!=posEntry)) ? true:false;
    pileupEvent       = pileupEventCandidate&&goodEvent_det;
    goodEvent_det_AND_muonDecayedInSample_gen = goodEvent_det && muonDecayedInSample_gen;
    //    if (pileupEvent) std::cout<<"pileupEvent (eventID="<<eventID<<std::endl;

    // Fill pileup-variables, but only if positron comes from different muon than the trigger signal
    //      if (mCounterHitExistsForThisEventID&&pCounterHitExistsForThisEventID)
    //      std::cout<<eventID<<" "<<mCounterHitExistsForThisEventID<<":  iiiEntry (trigger) = "<<iiiEntry<<"     kEntry (muon hit)="<<kEntry<<" "<<pCounterHitExistsForThisEventID<<"    posEntry (positron hit)="<<posEntry<<std::endl;    
    //    if (goodEvent_det) {
    //      std::cout<<"wwwwwwwwwwwwww \nEvent Nr: "<<eventID<<" muDecayTime="<< muDecayTime <<"  muM0Time="<<muM0Time<<"  det_time_start[idetP]="<<det_time_start[idetP]<<"   t1-t0="<<(timeBin1-timeBin0)*tdcresolution<<std::endl;
    //    }
    for (Int_t i=0; i<nrConditions; i++) { condition[i]=false; }
    //    MyPrintConditions();
    for (conditionMapType::const_iterator it = conditionMap.begin(); it!=conditionMap.end(); ++it) {
      int i = it->first;
      condition[i]=*(it->second);
    }
    
    // A special condition can me added here:
    // condition[29] = muIniMomZ>27;
    
    // Fill all "Detector Histograms"
    for(std::list<musrTH*>::const_iterator it = listOfAllHistograms1D.begin(); it != listOfAllHistograms1D.end(); ++it) {
      (*it)->FillTH1D(wght,&condition[0]);
    }
    for(std::list<musrTH*>::const_iterator it = listOfAllHistograms2D.begin(); it != listOfAllHistograms2D.end(); ++it) {
      (*it)->FillTH2D(wght,&condition[0]);
    }
    if ((goodEvent_det)&&(pCounterHitExistsForThisEventID)) {
      //	musrCounter* pCounterHitInThisEvent 
      counterMapType::const_iterator itPcounterHitInThisEvent = pCounterMap.find(idetP_ID);
      if (itPcounterHitInThisEvent==pCounterMap.end()) {
	std::cout<<" ERROR!  musrAnalysis:  itPcounterHitInThisEvent not found ==> This should never happen!!!"<<std::endl;
	std::cout<<"      idetP="<<idetP<<" idetP_ID="<<idetP_ID<< std::endl;
      }
      else{
	(itPcounterHitInThisEvent->second)->FillTDChistogram(float(timeBin1-timeBin0+overallBinDelay)+0.5,wght);
      }
    }

    //  Check whether there is a good trigger in the next event
    //    std::cout<<" FillHistograms: timeBinOfThePreviouslyProcessedHit = "<<timeBinOfThePreviouslyProcessedHit<<std::endl;
    //    mCounterHitExistsForThisEventID = mCounter->GetNextGoodHitInThisEvent(eventID,timeBinOfThePreviouslyProcessedHit,0,'M',timeBin0,kEntry,idetM,idetM_ID,idetM_edep,doubleHitM);
    mCounterHitExistsForThisEventID = mCounter->GetNextGoodMuon(eventID,timeBinOfThePreviouslyProcessedHit,timeBin0,kEntry,idetM,idetM_ID,idetM_edep,doubleHitM);
    timeBinOfThePreviouslyProcessedHit = timeBin0;
    //    if (mCounterHitExistsForThisEventID) std::cout<<"   YYYYYYYYYYYYYYYYYYY check this : LOOOPING AGAIN"<<std::endl;
    oncePerEvent=false;				    
  }  while(mCounterHitExistsForThisEventID);

}

//================================================================
void musrAnalysis::InitialiseEvent() {
   runID_double        = runID;
   eventID_double      = eventID;
   muDecayDetID_double = muDecayDetID;
   det_n_double        = det_n;
   muDecayPosR         = sqrt(muDecayPosX*muDecayPosX+muDecayPosY*muDecayPosY);

   wght                = ( fChain->GetBranch("weight") ) ? weight : 1.;
}

//================================================================
Double_t musrAnalysis::PreprocessEvent(Long64_t iEn) {
  //  std::cout<<"musrAnalysis::PreprocessEvent()"<<std::endl;
  fChain->GetEntry(iEn); InitialiseEvent();
  //  std::cout<<"\n musrAnalysis::PreprocessEvent() Filling event "<<eventID<<std::endl;
  //  MyPrintTree();

  Double_t globTime = nextUnfilledEventTime;
  for (Int_t i=0; i<det_n; i++) {
    // //    Int_t detID=det_ID[i];
    std::map<int,musrCounter*>::iterator it;
    it = allCounterMap.find(det_ID[i]);
    if (it==allCounterMap.end()) {
      // uncomment      std::cout<<"Active detector with det_ID="<<det_ID[i]<<" not found !!!"<<std::endl;
    }
    else {
      //      Double_t omega = 851.372*fieldNomVal[0];
      Double_t dPhaseShift = (omega==0) ? 0:phaseShiftMap[det_ID[i]]/omega;
      //cDEL      if (det_ID[i]<20) std::cout<<"phaseShift = "<<phaseShiftMap[det_ID[i]]<<"    dPhaseShift="<<dPhaseShift<<"    tdcresolution="<<tdcresolution<<std::endl;
      Long64_t timeBin  = Long64_t( (globTime+det_time_start[i]            )/tdcresolution );
      Long64_t timeBin2 = Long64_t( (globTime+det_time_start[i]+dPhaseShift)/tdcresolution );
      //cDEL      if (det_ID[i]<20)       std::cout<<" timeBin = "<<timeBin<<"   ("<<globTime+det_time_start[i]<<")"<<std::endl;
      //cDEL      if (det_ID[i]<20)       std::cout<<" timeBin2 = "<<timeBin2<<"   ("<<globTime+det_time_start[i]+dPhaseShift<<")"<<std::endl;
      (*it).second->FillHitInCounter(det_edep[i],timeBin,timeBin2,iEn,eventID,i,det_ID[i]);
    }
  }

  //  std::cout<<"lastPreprocessedEntry+1="<<lastPreprocessedEntry+1<<"   iEn="<<iEn<<"   eventID="<<eventID<<std::endl;
  if ((lastPreprocessedEntry+1)!=iEn) {std::cout<<"ERROR PreprocessEvent - should never happen!!!"<<std::endl; }
  lastPreprocessedEntry = iEn;
  return timeToNextEvent*muonRateFactor;
}

//================================================================
Bool_t musrAnalysis::PositronCounterHit(Int_t evID, Long64_t dataBinMin, Long64_t dataBinMax, Long64_t& tBin1, Long64_t& tBin2, Int_t& kEntry, Int_t& idetP, Int_t& idetP_ID, Double_t& idetP_edep, Bool_t& doubleHit) {
  if (pCounterMap.empty()) return false;

  Bool_t positronHitFound = false;
  if (musrMode=='D') {
    // Loop over all positron counters
    for (counterMapType::const_iterator it = pCounterMap.begin(); it!=pCounterMap.end(); ++it) {
      //      std::cout<<" ===POSITRON==="<< pCounterMap.size()<<std::endl;
      //      Bool_t thereWasHit = (it->second)->GetNextGoodHitInThisEvent(evID,dataBinMin,dataBinMax,'P',tBin1,kEntry,idetP,idetP_ID,idetP_edep,doubleHit);
      Bool_t thereWasHit = (it->second)->GetNextGoodPositron(evID,dataBinMin,dataBinMax,tBin1,tBin2,kEntry,idetP,idetP_ID,idetP_edep,doubleHit);
      //      std::cout<<"000000000000  tBin1="<<tBin1<<"  tBin2="<<tBin2<<std::endl;
      if (doubleHit) {return false;}  // There were two hits in the same positron counter
      if (thereWasHit) {
	if (positronHitFound) {doubleHit = true; return false;}   // There were two hits in two different positron counters
	positronHitFound = true;
      }
    }
  }
  return positronHitFound;
}

//================================================================
void musrAnalysis::CopySubstring(char* inputChar, int iStart, int iEnd, char* outputChar) {
  int iiiEnd = std::min(iEnd,(int)strlen(inputChar));
  int lastChar=0;
  for (Int_t i=0; i<=(iiiEnd-iStart); i++) {
    outputChar[i]=inputChar[i+iStart];
    lastChar=i;
  }
  outputChar[lastChar+1]='\0';
  //  std::cout<<"iStart="<<iStart<<"  iiiEnd="<<iiiEnd<<std::endl;
}

//================================================================
void musrAnalysis::MyPrintTree() {
  std::cout<<"==========\n det_n="<<det_n<<std::endl;
  for (Int_t jj=0; jj<det_n; jj++) {
    std::cout<<"\t det_ID="<<det_ID[jj]<<"\t det_edep="<<det_edep[jj]<<"\t det_time_start="<<det_time_start[jj]<<std::endl;
  }
}

//================================================================
void musrAnalysis::MyPrintConditions() {
  std::cout<<"CONDITION NR = ";
  for (conditionMapType::const_iterator it = conditionMap.begin(); it!=conditionMap.end(); ++it) {
    std::cout<<"  "<<it->first;
  }
  std::cout<<std::endl<<"               ";
  for (conditionMapType::const_iterator it = conditionMap.begin(); it!=conditionMap.end(); ++it) {
    std::cout<<"  "<<*(it->second);
  }
  std::cout<<std::endl;
}
//================================================================
Double_t musrAnalysis::myAtan2(Double_t y, Double_t x) {
  if ( (y<-999.99999) && (y>-1000.00001) ) return -1000.;
  if ( (x<-999.99999) && (x>-1000.00001) ) return -1000.;
  if ( (x==0) && (y==0) ) return -1000.;
  return atan2(y,x);
}