Commit fa5c7f39 authored by Emilio Depero's avatar Emilio Depero
Browse files

first commit with important script

parents
#define APV_RAW_PED_cxx
#include "APV_RAW_PED.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
void APV_RAW_PED::Loop()
{
// In a ROOT session, you can do:
// Root > .L APV_RAW_PED.C
// Root > APV_RAW_PED 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;
Long64_t nentries = fChain->GetEntriesFast();
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); nbytes += nb;
// if (Cut(ientry) < 0) continue;
}
}
//////////////////////////////////////////////////////////
// This class has been automatically generated on
// Thu Aug 14 19:35:56 2014 by ROOT version 5.34/11
// from TTree apv_raw_ped/APVRawPedestals
// found on file: /Users/stefanos/Desktop/run11141.root
//////////////////////////////////////////////////////////
#ifndef APV_RAW_PED_h
#define APV_RAW_PED_h
#include <TROOT.h>
#include <TChain.h>
#include <TFile.h>
// Header file for the classes stored in the TTree if any.
#include <vector>
#include <vector>
#include <vector>
#include <vector>
#include <vector>
// Fixed size dimensions of array or collections stored in the TTree if any.
class APV_RAW_PED {
public :
TTree *fChain; //!pointer to the analyzed TTree or TChain
Int_t fCurrent; //!current Tree number in a TChain
// Declaration of leaf types
ULong64_t evt;
UInt_t error;
Int_t daqTimeSec;
Int_t daqTimeMicroSec;
Int_t srsTimeStamp;
UInt_t srsTrigger;
std::vector<unsigned int> *srsFec;
std::vector<unsigned int> *srsChip;
std::vector<unsigned int> *srsChan;
std::vector<std::string> *mmChamber;
std::vector<int> *mmLayer;
std::vector<char> *mmReadout;
std::vector<int> *mmStrip;
std::vector<double> *ped_mean;
std::vector<double> *ped_stdev;
std::vector<double> *ped_sigma;
// List of branches
TBranch *b_evt; //!
TBranch *b_error; //!
TBranch *b_daqTimeSec; //!
TBranch *b_daqTimeMicroSec; //!
TBranch *b_srsTimeStamp; //!
TBranch *b_srsTrigger; //!
TBranch *b_srsFec; //!
TBranch *b_srsChip; //!
TBranch *b_srsChan; //!
TBranch *b_mmChamber; //!
TBranch *b_mmLayer; //!
TBranch *b_mmReadout; //!
TBranch *b_mmStrip; //!
TBranch *b_ped_mean; //!
TBranch *b_ped_stdev; //!
TBranch *b_ped_sigma; //!
APV_RAW_PED(TTree *tree=0);
virtual ~APV_RAW_PED();
virtual Int_t Cut(Long64_t entry);
virtual Int_t GetEntry(Long64_t entry);
virtual Long64_t LoadTree(Long64_t entry);
virtual void Init(TTree *tree);
virtual void Loop();
virtual Bool_t Notify();
virtual void Show(Long64_t entry = -1);
};
#endif
#ifdef APV_RAW_PED_cxx
APV_RAW_PED::APV_RAW_PED(TTree *tree) : fChain(0)
{
// if parameter tree is not specified (or zero), connect the file
// used to generate this class and read the Tree.
if (tree == 0) {
TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("run691.root");
//TFile *f = root_file;
if (!f || !f->IsOpen()) {
f = new TFile("run691.root");
// f = root_file;
}
f->GetObject("apv_raw_ped",tree);
}
Init(tree);
}
APV_RAW_PED::~APV_RAW_PED()
{
if (!fChain) return;
delete fChain->GetCurrentFile();
}
Int_t APV_RAW_PED::GetEntry(Long64_t entry)
{
// Read contents of entry.
if (!fChain) return 0;
return fChain->GetEntry(entry);
}
Long64_t APV_RAW_PED::LoadTree(Long64_t entry)
{
// Set the environment to read one entry
if (!fChain) return -5;
Long64_t centry = fChain->LoadTree(entry);
if (centry < 0) return centry;
if (fChain->GetTreeNumber() != fCurrent) {
fCurrent = fChain->GetTreeNumber();
Notify();
}
return centry;
}
void APV_RAW_PED::Init(TTree *tree)
{
// The Init() function is called when the selector needs to initialize
// a new tree or chain. Typically here the branch addresses and branch
// pointers of the tree will be set.
// It is normally not necessary to make changes to the generated
// code, but the routine can be extended by the user if needed.
// Init() will be called many times when running on PROOF
// (once per file to be processed).
// Set object pointer
srsFec = 0;
srsChip = 0;
srsChan = 0;
mmChamber = 0;
mmLayer = 0;
mmReadout = 0;
mmStrip = 0;
ped_mean = 0;
ped_stdev = 0;
ped_sigma = 0;
// Set branch addresses and branch pointers
if (!tree) return;
fChain = tree;
fCurrent = -1;
fChain->SetMakeClass(1);
fChain->SetBranchAddress("evt", &evt, &b_evt);
fChain->SetBranchAddress("error", &error, &b_error);
fChain->SetBranchAddress("daqTimeSec", &daqTimeSec, &b_daqTimeSec);
fChain->SetBranchAddress("daqTimeMicroSec", &daqTimeMicroSec, &b_daqTimeMicroSec);
fChain->SetBranchAddress("srsTimeStamp", &srsTimeStamp, &b_srsTimeStamp);
fChain->SetBranchAddress("srsTrigger", &srsTrigger, &b_srsTrigger);
fChain->SetBranchAddress("srsFec", &srsFec, &b_srsFec);
fChain->SetBranchAddress("srsChip", &srsChip, &b_srsChip);
fChain->SetBranchAddress("srsChan", &srsChan, &b_srsChan);
fChain->SetBranchAddress("mmChamber", &mmChamber, &b_mmChamber);
fChain->SetBranchAddress("mmLayer", &mmLayer, &b_mmLayer);
fChain->SetBranchAddress("mmReadout", &mmReadout, &b_mmReadout);
fChain->SetBranchAddress("mmStrip", &mmStrip, &b_mmStrip);
fChain->SetBranchAddress("ped_mean", &ped_mean, &b_ped_mean);
fChain->SetBranchAddress("ped_stdev", &ped_stdev, &b_ped_stdev);
fChain->SetBranchAddress("ped_sigma", &ped_sigma, &b_ped_sigma);
Notify();
}
Bool_t APV_RAW_PED::Notify()
{
// The Notify() function is called when a new file is opened. This
// can be either for a new TTree in a TChain or when when a new TTree
// is started when using PROOF. It is normally not necessary to make changes
// to the generated code, but the routine can be extended by the
// user if needed. The return value is currently not used.
return kTRUE;
}
void APV_RAW_PED::Show(Long64_t entry)
{
// Print contents of entry.
// If entry is not specified, print current entry
if (!fChain) return;
fChain->Show(entry);
}
Int_t APV_RAW_PED::Cut(Long64_t entry)
{
(void) entry;
// This function may be called from Loop.
// returns 1 if entry is accepted.
// returns -1 otherwise.
return 1;
}
#endif // #ifdef APV_RAW_PED_cxx
// ROOT
#include "TNtuple.h"
#include "TVector3.h"
#include "TProfile.h"
#include "TFile.h"
#include "TCanvas.h"
#include "TGraph.h"
#include "TMultiGraph.h"
#include "TStyle.h"
#include "TH2.h"
#include "TTree.h"
#include <vector>
#include <vector>
#include <TH1F.h>
#include <TH2F.h>
#include <TCanvas.h>
#include "TBranch.h"
#include "TKey.h"
#include "TClass.h"
#include "TDirectory.h"
#include "Riostream.h"
#include "TH1.h"
#include "TF1.h"
#include "TH1D.h"
#include "TH2.h"
#include "TROOT.h"
//main parameters GLOBALS
const char* treename = "t1";
const int NMM = 2;
const TVector3 targetpos(0.,0.,154.65); //target position of target received by Irene
struct MMinfo
{
TString name;
int range[2];
bool hasHit;
double x,y,z;
MMinfo()
:
name("dummy")
{
range[0] = 0;
range[1] = 1e5;
hasHit = false;
x = -9999;
y = -9999;
z = -9999;
}
MMinfo(TString name_,
int rl,
int rh)
:
name(name_)
{
range[0] = rl;
range[1] = rh;
hasHit = false;
x = -9999;
y = -9999;
z = -9999;
}
TVector3 pos()
{
return TVector3(x,y,z);
}
bool IsInRange(int detid)
{
return range[0] < detid && detid < range[1];
}
void RecordHit(double x_, double y_, double z_)
{
x = x_;
y = y_;
z = z_;
hasHit = true;
}
void Reset()
{
hasHit = false;
x = -9999;
y = -9999;
z = -9999;
}
};
//line extrapolation in plane defined by Y coordinate
void ExtrapolateLine(MMinfo mm[NMM],
double surfY,
double& recox, double& recoz, double& recoangle)
{
//for now assuming just two MM
const TVector3 p1 = mm[0].pos();
const TVector3 p2 = mm[1].pos();
const TVector3 diff = p2 - p1;
const double t = (surfY - p2.Y()) / diff.Y();
const TVector3 crossp = p2 + t * diff;
//assign it
recox = crossp.X();
recoz = crossp.Z();
recoangle = (p2-p1).Angle(targetpos);
return;
}
int main(int argc, char* argv[])
{
//guard
if(argc != 2)
{
std::cerr << "Invalid number of argument\n";
std::cout << "Usage: ./ConvertOutputMM <filename> \n";
return 1;
}
TFile* file = new TFile(argv[1]);
TTree* tree = (TTree*)file->Get(treename);
//trees
const Int_t nMax=1000;
double X[nMax],Y[nMax],Z[nMax]; //coordinates of hit
int EventID,DetID[nMax];
int nsave;
double mumDX,mumDY,mumDZ;
double recoX,recoZ,recoAngle;
//set branch address
tree->SetBranchAddress("save_n", &nsave);
tree->SetBranchAddress("eventID", &EventID);
tree->SetBranchAddress("save_x", X);
tree->SetBranchAddress("save_y", Y);
tree->SetBranchAddress("save_z", Z);
tree->SetBranchAddress("save_detID", DetID);
//muonium
tree->SetBranchAddress("muDecayPosX", &mumDX);
tree->SetBranchAddress("muDecayPosY", &mumDY);
tree->SetBranchAddress("muDecayPosZ", &mumDZ);
//outname file
TString name(argv[1]);
name.ReplaceAll(".root","_mmout.root");
TFile* outfile = new TFile(name,"RECREATE");
TTree* mumtree = new TTree("mumtruth","mumtruth");
mumtree->Branch("MumEventID", &EventID, "EventID/I");
mumtree->Branch("MumDecayX", &mumDX, "MumDecayX/D");
mumtree->Branch("MumDecayY", &mumDY, "MumDecayY/D");
mumtree->Branch("MumDecayZ", &mumDZ, "MumDecayZ/D");
mumtree->Branch("RecoX", &recoX, "RecoX/D");
mumtree->Branch("RecoZ", &recoZ, "RecoZ/D");
mumtree->Branch("RecoAngle", &recoAngle, "RecoAngle/D");
//create MMinfo
MMinfo mminfo[NMM];
unsigned int StartID = 8000;
for(unsigned int n(0);n < NMM; ++n)
{
const int thisid = StartID + n*10;
name.Form("MM%i", n+1);
mminfo[n] = MMinfo(name,
thisid, thisid+10);
}
//create Ntuples
std::vector<TNtuple*> theNtuples;
for(unsigned int n(0);n < NMM; ++n)
{
name.Form("MM%i",n+1);
theNtuples.push_back(new TNtuple(name, //name
name, //title
"EventID:x:y:t_rise:t_max:q_max:q_min:q_tot:nStrips")
);
}
//loop over tree
for(unsigned int entry(0); entry < tree->GetEntries(); ++entry)
{
//fill tree
tree->GetEntry(entry);
bool MMisHit(false);
//check if an MM is in range
for(unsigned int n(0);n < NMM; ++n)
{
//loop over indices
for(unsigned int index(0);index<nsave;++index)
{
//check if they are in range
if(mminfo[n].IsInRange(DetID[index]))
{
mminfo[n].RecordHit(X[index], Y[index], Z[index]);
//fill the Ntuples
theNtuples[n]->Fill(EventID, // event ID
X[index]-targetpos.X(),Z[index]-targetpos.Z(), // position in the micromegas, Y and Z exchanged
0,0,0,0,0,0 // all other at zeros in the simulation
);
}
}
}
ExtrapolateLine(mminfo,
targetpos.Y(),
recoX, recoZ, recoAngle);
if(mminfo[0].hasHit)mumtree->Fill();
if(entry%10000==0)std::cout << "processed entry: " << entry << "\n";
for(unsigned int n(0);n < NMM; ++n)mminfo[n].Reset();
}
//return 0;
//write objects in directory
outfile->cd();
for(unsigned int n(0);n < NMM; ++n)
{
outfile->WriteObject(theNtuples[n], mminfo[n].name);
}
mumtree->Write();
outfile->Close();
file->Close();
}
//compile with g++ -std=c++11 -o analyze_PositionData analyze_PositionData.cc `root-config --cflags --glibs`
//
// Date: November 2020
//
// Author: Jonas Nuber
//
// Purpose: write out ntuples of 4 micromegas from root file of micromegas daq (prepare data for plotting)
#include <vector>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <algorithm>
// ROOT
#include "TNtuple.h"
#include "TVector3.h"
#include "TProfile.h"
#include "TFile.h"
#include "TCanvas.h"
#include "TGraph.h"
#include "TMultiGraph.h"
#include "TStyle.h"
#include "TH2.h"
#include "TTree.h"
#include <vector>
#include <vector>
#include <TH1F.h>
#include <TH2F.h>
#include <TCanvas.h>
#include "TBranch.h"
#include "TKey.h"
#include "TClass.h"
#include "TDirectory.h"
#include "Riostream.h"
#include "TH1.h"
#include "TF1.h"
#include "TH1D.h"
#include "TH2.h"
#include "TROOT.h"
#include <string>
// custom imports
#include "apv_raw.C"
#include "APV_RAW_PED.C"
#include "globals.h"
#include "Micromega.hh"
#include "MicromegasPlane.hh"
int main(int argc, char* argv[])
{
// takes 2 arguments: infile name and outfile name
if(argc != 3)
{
std::cout << "Usage: ./analyze_PositionData [inputRootFile] [outputRootFile]" << std::endl;
std::exit(-1);
}
std::string inputFile = argv[1];
std::string outputFile = argv[2];
// cut on maximal number of strips in cluster
int maxClusterSize_X = 15;
int maxClusterSize_Y = 40;
// cut on maximal number of clusters
int maxNClusters_X = 27;
int maxNClusters_Y = 27;
// the input file
TFile* inFile = new TFile(TString(inputFile),"READ");
// the output file
TFile* outFile = new TFile(TString(outputFile),"RECREATE");
// number of Micromegas
int numberOfMM = 2;
// names of the MM in the input file
std::vector<std::string> input_Names;
input_Names.push_back("MM1");
input_Names.push_back("MM2");
// names in the outputfile: two sides of detectors (up/down), close and far each
std::vector<std::string> output_Names;
output_Names.push_back("MM1");
output_Names.push_back("MM2");
// define all MM structs and a vector for them
Micromega MM1 = Micromega();
Micromega MM2 = Micromega();
std::vector<Micromega> theMMvector = {MM1,MM2};
// initialize default thresholds of planes in MM objects
for (int i=0;i<theMMvector.size();i++){theMMvector[i].initPlaneThresholds();}
/// define the Ntuples for output
TNtuple* MM1_tuple = new TNtuple(TString(output_Names[0]),TString(output_Names[0]),"EventID:x:y:t_rise:t_max:q_max:q_min:q_tot:nStrips");
TNtuple* MM2_tuple = new TNtuple(TString(output_Names[1]),TString(output_Names[1]),"EventID:x:y:t_rise:t_max:q_max:q_min:q_tot:nStrips");
std::vector<TNtuple*> theNtuples = {MM1_tuple,MM2_tuple};
/// now set individual thresholds for all strips depending on their stdev in the pedestal
APV_RAW_PED* apv_raw_ped_obj = new APV_RAW_PED((TTree*)inFile->Get("apv_raw_ped"));
apv_raw_ped_obj->fChain->GetEntry(0);
// loop over all MM and strips to set thresholds
for (int mmIndex=0;mmIndex<numberOfMM;mmIndex++)
{
for(int j=0; j<apv_raw_ped_obj->mmStrip->size(); j++)
{
// for this stripe, decide whether it is on the current MM. If so, define threshold
if (apv_raw_ped_obj->mmChamber->at(j)==input_Names[mmIndex])
{
if (apv_raw_ped_obj->mmReadout->at(j)==49)
{theMMvector[mmIndex].xplane.SetThreshold_forStrip(apv_raw_ped_obj->mmStrip->at(j),4.*apv_raw_ped_obj->ped_stdev->at(j));}
else if (apv_raw_ped_obj->mmReadout->at(j)==48)
{theMMvector[mmIndex].yplane.SetThreshold_forStrip(apv_raw_ped_obj->mmStrip->at(j),4.*apv_raw_ped_obj->ped_stdev->at(j));}
}
}
}
// get the tree from file
apv_raw* apv_raw_obj = new apv_raw((TTree*)inFile->Get("apv_raw")); // Get the tree
const int entries = apv_raw_obj->fChain->GetEntries();
std::cout << "Entering event loop in file " << inFile << std::endl;
// loop over all event entries to do do clustering and get the hits
for(int iEntry=0; iEntry<entries; iEntry++)
{
// go to next event
apv_raw_obj->fChain->GetEntry(iEntry);
// (re-) initialize all the MM structs
for (int i=0;i<theMMvector.size();i++)
{
theMMvector[i].initPlanes();
}
// loop over all MM objects, do the clustering and identify hits