// -------------------------------------------------------------------------
// -----                 PndMCTestTimeWalk source file             -----
// -----                  Created 18/07/08  by S.Esch              -----
// -------------------------------------------------------------------------
// libc includes
#include <iostream>

// Root includes
#include "TROOT.h"
#include "TClonesArray.h"
#include "TF1.h"

// framework includes
#include "FairRootManager.h"
#include "PndMCTestTimeWalk.h"
#include "FairRun.h"
#include "FairRuntimeDb.h"
#include "FairHit.h"
//#include "FairLinkedData.h"

#include "PndSdsMCPoint.h"
#include "PndSdsHit.h"
#include "PndMCEntry.h"
#include "PndSdsDigiPixelMCInfo.h"

#include "PndDetectorList.h"




// -----   Default constructor   -------------------------------------------
PndMCTestTimeWalk::PndMCTestTimeWalk()
	: FairTask("Creates PndMC test")
{
}
// -------------------------------------------------------------------------


// -----   Destructor   ----------------------------------------------------
PndMCTestTimeWalk::~PndMCTestTimeWalk()
{
}

// -----   Public method Init   --------------------------------------------
InitStatus PndMCTestTimeWalk::Init()
{
 //Histogramme anlegen

	hisMCEnergyLossVHitCharge = new TH2F("hisMCEnergyLossVHitCharge", "MC Energyloss vs. Hit Charge; MC Energyloss [e]; Hit Charge [e] ",1000,0,70000,1000,0,70000 );
	hisMCEnergyLossVNDigis = new TH2F("hisMCEnergyLossVNDigis","MC Energyloss vs. Number of Digis; MC Energyloss [e]; # Digis",1000,0,70000,15,0.5,15.5 );
	hisMCEnergyLossVHitChargeError = new TH2F("hisMCEnergyVHitChargeError","MC Energyloss vs.Hit Charge Error; MC Energyloss [e]; #Delta HitCharge = MCEnergyloss - HitCharge [e]",1000,0,70000,600,-3000,3000 );
	hisRecoTimeWalk = new TH1F("hisRecoTimeWalk","TimeWalkReco; TimeWalkReco [ns]", 1050,0,105);


	hisMCEnergyLossVHitChargeRelError = new TH2F("hisMCEnergyVHitChargeRelError","MC Energyloss vs. Hit Charge RelError; MC Energyloss [e]; #Delta HitCharge = (MCEnergyloss - HitCharge)/MCEnergyloss [e]",1000,0,70000,600,-0.3,0.3 );

	hisTotVMCCharge = new TH2F("hisTotVMCCharge", "Tot vs. MCCharge; MC Charge [e]; ToT [ns]", 1000,0,70000, 100,0,200);

	hisMCChargeVRecoChargeError = new TH2F("hisMCChargeVRecoChargeError", "MCCharge vs. #Delta ChargeReco; MC Charge [e]; #Delta RecoCharge = | RecoCharge - MC Charge | [e]", 700,0,70000,2000,0,100);

	hisMCChargeVCorrTimeStampError = new TH2F("hisMCChargeVCorrTimeStampError","MCCharge vs #Delta CorrTimeStamp;MC Charge [e] ;#Delta CorrTimeStamp = CorrTS - Tof - ET [ns]",700,0,70000,500,-10,15);

	hisMCChargeVCorrTimeWalkError = new TH2F("hisMCChargeVCorrTimeWalkError","MCCharge vs. #Delta CorrTimeWalk;MCCharge [e] ;#Delta CorrTimeWalk = | CorrTW - DigiTW |[ns]",700,0,70000,400,-1,7);

	hisMCChargeVCorrTimeWalk = new TH2F("hisMCChargeVCorrTimeWalk","MCCharge vs. CorrTimeWalk; MCCharge [e]; CorrTimeWalk [ns]", 1000,0,70000,550,0,110);

	hisMCChargeVDigiTimeWalk = new TH2F("hisMCChargeVDigiTimeWalk","MCCharge vs. DigiTimeWalk; MC Charge [e]; Digi TimeWalk [ns]",1000,0,70000,550,0,110);


	hisDigiChargeVRecoChargeError = new TH2F("hisDigiChargeVRecoChargeError", "DigiCharge vs. #Delta ChargeReco; DigiCharge [e]; #Delta RecoCharge = | RecoCharge - DigiCharge | [e]", 700,0,70000,1000,0,500);

	hisDigiChargeVCorrTimeStampError = new TH2F("hisDigiChargeVCorrTimeStampError","DigiCharge vs #Delta CorrTimeStamp;DigiCharge [e] ;#Delta CorrTimeStamp = CorrTS - Tof - ET [ns]",700,0,70000,500,-10,15);

	hisDigiChargeVCorrTimeWalkError = new TH2F("hisDigiChargeVCorrTimeWalkError","DigiCharge vs. #Delta CorrTimeWalk;DigiCharge [e] ;#Delta CorrTimeWalk = | CorrTW - DigiTW |[ns]",700,0,70000,400,-1,7);

	hisDigiChargeVCorrTimeWalk = new TH2F("hisDigiChargeVCorrTimeWalk","DigiCharge vs. CorrTimeWalk; DigiCharge [e]; CorrTimeWalk [ns]", 700,0,70000,440,0,110);

	hisDigiChargeVDigiTimeWalk = new TH2F("hisDigiChargeVDigiTimeWalk","DigiCharge vs. DigiTimeWalk; DigiCharge [e]; DigiTimeWalk [ns]",700,0,70000,440,0,110);

	hisDigiChargeVDigiTimeStampWOEventTime =new TH2F("hisDigiChargeVDigiTimeStampWOEventTime", "DigiCharge vs. DigiTimeStamp w/o EventTime; DigiCharge [e] ; DigiTimeStamp - EventTime", 700,0,70000,100,0,100);

	hisPointTofvMCDigiTof = new TH2F("hisPointTofvMCDigiTof","MCPointTof vs. DigiMCTof;DigiMC;MCPoint",200,0,2,200,0,2);

	hisDigiChargeVMCCharge = new TH2F("hisDigiChargeVMCCharge", "DigiCharge vs. MCCharge; DigiCharge [e]; MCCharge [e]", 1000,0,70000,1000,0,70000);

	hisDigiChargeVMCChargeError = new TH2F("hisDigiChargeVMCChargeError", "DigiCharge vs.  #Delta MCCharge; DigiCharge [e]; # Delta MCCharge = | MCCharge - DigiCharge | [e]", 1000,0,70000,1600,-800,800);

	DeltaTot = new TF1("DeltaTot","(1/[1]+[2]*[3]/(x*x))*[4]",0,70000);

	DeltaTot->SetParameter(1,60); // const current
	DeltaTot->SetParameter(2,1100);  // threshold
	DeltaTot->SetParameter(3,100); // rising time
	DeltaTot->SetParameter(4,200); // sigma noise


	DeltaChargeReco =new  TF1("DeltaChargeReco","DeltaTot*[1]/2.*(1-(([1]*[0]-[2]-x*[1])*0.5)/(sqrt(pow(([0]*[1]-[2]-x*[1])*0.5,2)+[0]*[2]*[1])))",0,80000);

	DeltaChargeReco->SetParameter(0,100); // rising time
	DeltaChargeReco->SetParameter(1,60); // const current
	DeltaChargeReco->SetParameter(2,1100); // threshold

	Final_add = new TF1("Final_add", "[3]+[2]*[0]/(x*x)*[1]*[3]/2.*(1+(x/2.-[4]/(2*x))/sqrt(pow([4]/(2*x)-x/2.,2)+[4]))",0,80000 );

	 Final_add->SetParameter(0,100); // t rising 100 ns
	 Final_add->SetParameter(1,60); // const current 60 e/ns
	 Final_add->SetParameter(2,1100); // Q threshold 1100 e
	 Final_add->SetParameter(3,6.67); // delta tot 6,67 ns
	 Final_add->SetParameter(4,6600000);
	 Final_add->SetLineWidth(1);
	 Final_add->SetLineColor(2);

	 Final_add2 = new TF1("Final_add2", "-([2]*[0]/(x*x)*[1]*[3]/2.*(1+(x/2.-[4]/(2*x))/sqrt(pow([4]/(2*x)-x/2.,2)+[4])))",0,80000 );

	 Final_add2->SetParameter(0,100); // t rising 100 ns
	 Final_add2->SetParameter(1,60); // const current 60 e/ns
	 Final_add2->SetParameter(2,1100); // Q threshold 1100 e
	 Final_add2->SetParameter(3,6.67); // delta tot 6,67 ns
	 Final_add2->SetParameter(4,6600000);
	 Final_add2->SetLineWidth(1);
	 Final_add2->SetLineColor(2);



//  fMCMatch->InitStage(kMCTrack, "", "MCTrack");

  FairRootManager* ioman = FairRootManager::Instance();
  	if (!ioman) {
  		std::cout << "-E- PndMCTestTimeWalk::Init: "
  				<< "RootManager not instantiated!" << std::endl;
  		return kFATAL;
  	}

  	fMCMatch = (PndMCMatch*)ioman->GetObject("MCMatch");
  	fPixelHit = (TClonesArray*)ioman->GetObject("MVDHitsPixel");
  	fMCPoint = (TClonesArray*)ioman->GetObject("MVDPoint");
  	fPixelDigiMCInfo = (TClonesArray*)ioman->GetObject("MVDPixelDigisMCInfo");
  	fPixelDigiCorr = (TClonesArray*)ioman->GetObject("MVDDigisCorr");
  	fPixelMCHeader = (FairMCEventHeader*)ioman->GetObject("MCEventHeader.");
  	fPixelDigiAdditionalInfo = (TClonesArray*)ioman->GetObject("MVDPixelDigisAdditionalInfo");

	std::cout << "-I- PndMCTestTimeWalk::Init: Initialization successfull" << std::endl;


  return kSUCCESS;
}


// -------------------------------------------------------------------------
void PndMCTestTimeWalk::SetParContainers()
{
  // Get Base Container
//  FairRun* ana = FairRun::Instance();
//  FairRuntimeDb* rtdb=ana->GetRuntimeDb();

}


// -----   Public method Exec   --------------------------------------------
void PndMCTestTimeWalk::Exec(Option_t* opt)
{
std::cout << "============= Begin PndMCTestTimeWalk::Exec" << std::endl;
std::cout << std::endl;

		PndMCResult myResult = fMCMatch->GetMCInfo("MVDHitsPixel", "MVDPoint");

		PndMCResult myResult2 = fMCMatch->GetMCInfo("MVDDigisCorr","MVDPoint");

		PndMCResult myResult4 = fMCMatch->GetMCInfo("MVDDigisCorr","MVDPixelDigisMCInfo");

		PndMCResult myResult5 = fMCMatch->GetMCInfo("MVDDigisCorr","MVDPixelDigisAdditionalInfo");

		std::cout << "Hits --> Points" << myResult << std::endl;
		std::cout << "Digis --> Points" << myResult2 << std::endl;
		std::cout << "Digis --> MCInfo" << myResult4 << std::endl;

		FairRootManager* ioman = FairRootManager::Instance();

		for (int i = 0; i < fPixelHit->GetEntries(); i++){ // ich loope ueber fPixelHit, weil ich zu jedem Hit meinen MCPoint finden moechte
			PndMCEntry myLinks = myResult.GetMCLink(i);

			PndSdsHit* myHit = (PndSdsHit*)fPixelHit->At(i);

			std::cout << "fPixelHit->GetEntries() " << fPixelHit->GetEntries()<<  " myResult->GetNEntries() " <<  myResult.GetNEntries() <<std::endl;

			for (int j = 0; j < myLinks.GetNLinks(); j++){
				std::cout << "GetNLinks " << myLinks.GetNLinks() << std::endl;

				if (myLinks.GetLink(j).GetType() == ioman->GetBranchId("MVDPoint")){

					PndSdsMCPoint* myMCPoint = (PndSdsMCPoint*)fMCPoint->At(myLinks.GetLink(j).GetIndex());

					if(myLinks.GetNLinks()==1)
					{
						// Namenskonvention his-XAchse-YAchse
						// MCCharge - DigiCharge in e
						// MCEnergyLosse - EnergyLoss MCPoint

						hisMCEnergyLossVHitCharge->Fill(myMCPoint->GetEnergyLoss() / (3.61e-9),myHit->GetCharge());
						hisMCEnergyLossVNDigis->Fill(myMCPoint->GetEnergyLoss() / (3.61e-9),myHit->GetNDigiHits() );
						hisMCEnergyLossVHitChargeError->Fill(myMCPoint->GetEnergyLoss() / (3.61e-9), myMCPoint->GetEnergyLoss() / (3.61e-9)-myHit->GetCharge());
						hisRecoTimeWalk->Fill(5);
						hisMCEnergyLossVHitChargeRelError->Fill(myMCPoint->GetEnergyLoss() / (3.61e-9), (3.61e-9)*(myMCPoint->GetEnergyLoss() / (3.61e-9) - myHit->GetCharge())/ myMCPoint->GetEnergyLoss() );
					}
				}
			}
			std::cout << std::endl;
		}

		for (int k=0; k < fPixelDigiCorr->GetEntries(); k++){  // loop ueber alle korrigierten digis
			PndMCEntry mylinks2 = myResult2.GetMCLink(k);
			PndMCEntry mylinks4 = myResult4.GetMCLink(k);
			PndMCEntry mylinks5 = myResult5.GetMCLink(k);

			PndSdsDigiPixel* myDigiCorr = (PndSdsDigiPixel*)fPixelDigiCorr->At(k);

			std::cout << " GetNLnks():  " << mylinks2.GetNLinks() << std::endl;

			for (int l = 0; l < mylinks2.GetNLinks(); l++){

				std::cout << "GetType: " << mylinks5.GetLink(l).GetType() << " ioman->GetBranchID: " <<  ioman->GetBranchId("MVDPixelDigisAdditionalInfo") << std::endl;

				if (mylinks2.GetLink(l).GetType() == ioman->GetBranchId("MVDPoint") && mylinks4.GetLink(l).GetType() == ioman->GetBranchId("MVDPixelDigisMCInfo") &&  mylinks5.GetLink(l).GetType() == ioman->GetBranchId("MVDPixelDigisAdditionalInfo")){

					PndSdsDigiPixelMCInfo* myDigiMC = (PndSdsDigiPixelMCInfo*)fPixelDigiMCInfo->At(mylinks4.GetLink(l).GetIndex());
					PndSdsDigiPixelMCInfo* myDigiAddInfo = (PndSdsDigiPixelMCInfo*)fPixelDigiAdditionalInfo->At(mylinks5.GetLink(l).GetIndex());
					PndSdsMCPoint* myMCPoint = (PndSdsMCPoint*)fMCPoint->At(mylinks2.GetLink(l).GetIndex());



					if(mylinks2.GetNLinks()==1 and mylinks4.GetNLinks()==1)
					{
						// Namenskonvention his-XAchse-YAchse
						// MCCharge - DigiCharge in e
						// MCEnergyLosse - EnergyLoss MCPoint

						hisMCChargeVCorrTimeStampError->Fill(myDigiMC->GetMCCharge(),myDigiCorr->GetTimeStamp()-myMCPoint->GetTime()-fPixelMCHeader->GetT());
						hisMCChargeVRecoChargeError->Fill(myDigiMC->GetMCCharge(),TMath::Abs(myDigiCorr->GetCharge()- myDigiMC->GetMCCharge()));
						hisMCChargeVCorrTimeWalkError->Fill(myDigiMC->GetMCCharge(),TMath::Abs(myDigiAddInfo->GetTimeWalkCorrection()-myDigiAddInfo->GetTimeWalk()));
						hisMCChargeVCorrTimeWalk->Fill(myDigiMC->GetMCCharge(),myDigiAddInfo->GetTimeWalkCorrection());
						hisMCChargeVDigiTimeWalk->Fill(myDigiMC->GetMCCharge(),myDigiAddInfo->GetTimeWalk());

						hisDigiChargeVCorrTimeStampError->Fill(myDigiMC->GetDigiCharge(),myDigiCorr->GetTimeStamp()-myMCPoint->GetTime()-fPixelMCHeader->GetT());
						hisDigiChargeVRecoChargeError->Fill(myDigiMC->GetDigiCharge(),TMath::Abs(myDigiCorr->GetCharge()- myDigiMC->GetDigiCharge()));
						hisDigiChargeVCorrTimeWalkError->Fill(myDigiMC->GetDigiCharge(),TMath::Abs(myDigiAddInfo->GetTimeWalkCorrection()-myDigiAddInfo->GetTimeWalk()));
						hisDigiChargeVCorrTimeWalk->Fill(myDigiMC->GetDigiCharge(),myDigiAddInfo->GetTimeWalkCorrection());
						hisDigiChargeVDigiTimeWalk->Fill(myDigiMC->GetDigiCharge(),myDigiAddInfo->GetTimeWalk());

						hisDigiChargeVMCCharge->Fill(myDigiMC->GetDigiCharge(),myDigiMC->GetMCCharge());
						hisDigiChargeVMCChargeError->Fill(myDigiMC->GetDigiCharge(),myDigiMC->GetMCCharge()-myDigiMC->GetDigiCharge());
						hisPointTofvMCDigiTof->Fill(myDigiMC->GetTof(),myMCPoint->GetTime());

						std::cout << *myDigiMC << std::endl;
					}
				}
			}
		}

//		for (int k=0; k < fPixelDigiMCInfo->GetEntries(); k++){  // loop ueber alle MC digis
//			PndMCEntry mylinks3 = myResult3.GetMCLink(k);
//
//			PndSdsDigiPixelMCInfo* myDigiMC = (PndSdsDigiPixelMCInfo*)fPixelDigiMCInfo->At(k);
//
//			for (int l = 0; l < mylinks3.GetNLinks(); l++){
//
//				if (mylinks3.GetLink(l).GetType() == ioman->GetBranchId("MVDPoint")){
//
//					myMCPoint = (PndSdsMCPoint*)fMCPoint->At(mylinks3.GetLink(l).GetIndex());
//
//					if(mylinks2.GetNLinks()==1)
//					{
//						// Namenskonvention his-XAchse-YAchse
//						// MCCharge - DigiCharge in e
//						// MCEnergyLosse - EnergyLoss MCPoint
//
//						//hisMCChargeChargeRecoError->Fill(myDigiMC->GetDigiMCCharge(),TMath::Abs(myDigiCorr->GetCharge()- myDigiMC->GetDigiMCCharge()));
//						//hisChargeMCTimeStampRecoError->Fill(myDigiMC->GetDigiMCCharge(),myDigiCorr->GetTime()-myMCPoint->GetTime());
//					}
//					//std::cout << "Cluster No.  " <<  myHit->GetClusterIndex()<< "Gehoert zu MCPoint "<< myMCPoint->GetIndex()  << std::endl;
//				}
//
//
//
//			}
//		}
//
		std::cout << "============= End PndMCTestTimeWalk::Exec" << std::endl;
		std::cout << std::endl;
}

void PndMCTestTimeWalk::Finish()
{
	// Hier die Histogramme speichern
	hisMCEnergyLossVHitCharge->Write();
	hisMCEnergyLossVNDigis->Write();
	hisMCEnergyLossVHitChargeError->Write();
	hisMCEnergyLossVHitChargeRelError->Write();
	hisTotVMCCharge->Write();
	hisMCChargeVRecoChargeError->Write();
	hisMCChargeVCorrTimeStampError->Write();
	hisMCChargeVCorrTimeWalkError->Write();
	hisMCChargeVCorrTimeWalk->Write();
	hisMCChargeVDigiTimeWalk->Write();

	hisDigiChargeVRecoChargeError->Write();
	hisDigiChargeVCorrTimeStampError->Write();
	hisDigiChargeVCorrTimeWalkError->Write();
	hisDigiChargeVCorrTimeWalk->Write();
	hisDigiChargeVDigiTimeWalk->Write();
	DeltaChargeReco->Write();
	DeltaTot->Write();
	Final_add->Write();
	Final_add2->Write();
	hisPointTofvMCDigiTof->Write();

	hisDigiChargeVMCCharge->Write();
	hisDigiChargeVMCChargeError->Write();
}


ClassImp(PndMCTestTimeWalk);