void RootTutorial1_2022_answer(){


	//->load data file
	TFile *f = new TFile("zjet_unrec.root", "read");

	//->create new file to write new tree
	TFile *f_save = new TFile("zjet_rec.root", "recreate");
	
	//->read data tree
	TTree* t = (TTree*)f->Get("Tdata");

	//->tree where reconstructed z boson will be saved
	TTree* t_save = new TTree("Tdata", "Tdata");

	//->histogram to verify boson reconstruction
	TH1F* h = new TH1F("h", "Z boson mass", 100, 0, 150);

	//->declare vectors to where event variables should be assigned
	vector<int> *id = 0;
	vector<double> *m = 0;
	vector<double> *px = 0;
	vector<double> *py = 0;
	vector<double> *pz = 0;
	vector<double> *E = 0;

	//->declare scalars to where you want to save the new ttree information
	Double_t m_save, px_save, py_save, pz_save, E_save;



	//->assign the branches for the data tree
	t->SetBranchAddress("id",&id);
	t->SetBranchAddress("m",&m);
	t->SetBranchAddress("px",&px);
	t->SetBranchAddress("py",&py);
	t->SetBranchAddress("pz",&pz);
	t->SetBranchAddress("En",&E);


	//->create new branches for the variables that you want to save in the new tree
	t_save->Branch("px", &px_save);
	t_save->Branch("py", &py_save);
	t_save->Branch("pz", &pz_save);
	t_save->Branch("En", &m_save);
	t_save->Branch("m", &m_save);


	//->get number of entries in data tree
	int NEntries = t->GetEntries();

	//loop through the events
	for(int i = 0; i< NEntries;i++){

		//update for the event i
		t->GetEntry(i);

		//auxiliary variables for the invariant mass calculation
		float muon[4] = {0.}, muon_[4] = {0.};//will save the muons/antimuons px, py, pz, E: muon - muon, muon_ - antiparticle muon


		//loop through each particle in the event i
		for ( int iObj = 0; iObj < (*id).size(); ++iObj ) {

		 //muon
		  if ( (*id)[iObj] == 13){
		  	muon[1] = (*px)[iObj];
		  	muon[2] = (*py)[iObj];
		  	muon[3] = (*pz)[iObj];
		  	muon[0] = (*E)[iObj];
		  }
		  //muon's antiparticle
		  else if ( (*id)[iObj] == -13){
		  	muon_[1] = (*px)[iObj];
		  	muon_[2] = (*py)[iObj];
		  	muon_[3] = (*pz)[iObj];
		  	muon_[0] = (*E)[iObj];
		  }
		}

		//check if we found the two, and in that case calculate invariant mass
		if(muon[0]!=0. && muon_[0]!=0.){

			//first term: quadratic sum of energies
			float sum_E = muon[0]+muon_[0];
			float Px = (muon[1] + muon_[1]);
			float Py = (muon[2] + muon_[2]);
			float Pz = (muon[3] + muon_[3]);

			//second term: quadrtic sum of momentums
			float Zp = sqrt(Px*Px + Py*Py + Pz*Pz);
			float im = sqrt(sum_E*sum_E - Zp*Zp);//GEv
			

			//fill histogram with invariant mass calculated
			h->Fill(im,1);

			//verify first cases' invariant masses
			if (i<10)
				printf("Invariant mass calculated: %f\n", im);


			px_save = Px;
			py_save = Py;
			pz_save = Pz;
			E_save = sum_E;
			m_save = im;
			t_save->Fill();
		}

	}

	//Verify if spectrum makes sense
	h->Draw();

	//save new TTree in file
	t_save->Write();

	
	return;
}