analyze_PositionData.cc 6.85 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
//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);
    }

Emilio Depero's avatar
Emilio Depero committed
68
69
    const TString inputFile(argv[1]);
    const TString outputFile(argv[2]);
70
71
72
73
74
75
76
77
78

    // 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;

Emilio Depero's avatar
Emilio Depero committed
79
80
    std::cout << "Input file: " << inputFile << " output file: " << outputFile << "\n";

81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
    // 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
        for (int mmIndex=0;mmIndex<numberOfMM;mmIndex++)
        {
            // this loop is over all strips of the event (is there a faster way than looping 4 times?)
            for(int j=0; j<apv_raw_obj->mmStrip->size(); j++)
            {
                // for this strip, decide whether it is on the current MM. If so, perform clustering
                if (apv_raw_obj->mmChamber->at(j)==input_Names[mmIndex])
                {
                    if (apv_raw_obj->mmReadout->at(j)==49)
                        {theMMvector[mmIndex].xplane.accumulateHits_forStrip(apv_raw_obj->mmStrip->at(j),apv_raw_obj->raw_q->at(j));}
                    else if (apv_raw_obj->mmReadout->at(j)==48)
                        {theMMvector[mmIndex].yplane.accumulateHits_forStrip(apv_raw_obj->mmStrip->at(j),apv_raw_obj->raw_q->at(j));}
                    //std::cout << "pos 2.1.3 " << std::endl;
                    //DoMMReco(theMMvector[mmIndex]);
                }
            }
            theMMvector[mmIndex].computeHit();

            if (theMMvector[mmIndex].hasHit())
            {
	        bool passedCuts = true;
		//if ( theMMvector[mmIndex].xplane.GetNumberOfClusters()>maxNClusters_X || theMMvector[mmIndex].xplane.GetClusterSize(theMMvector[mmIndex].xplane.findBestCluster())>maxClusterSize_X )
                   // {passedCuts = false;}
		//if (theMMvector[mmIndex].yplane.GetNumberOfClusters()>maxNClusters_Y || theMMvector[mmIndex].yplane.GetClusterSize(theMMvector[mmIndex].yplane.findBestCluster())>maxClusterSize_Y )
	          //  {passedCuts = false;}
                if (passedCuts) theNtuples[mmIndex]->Fill(iEntry,theMMvector[mmIndex].GetHitPositionX(),theMMvector[mmIndex].GetHitPositionY(),theMMvector[mmIndex].GetHitTimeRise(),theMMvector[mmIndex].GetHitTimeMax(),theMMvector[mmIndex].GetHitChargeMax(),theMMvector[mmIndex].GetHitChargeMin(),theMMvector[mmIndex].GetHitChargeTot(),float(theMMvector[mmIndex].GetHitNumberOfStrips()));
            }
        }
positron cool guys's avatar
positron cool guys committed
186
        if (iEntry%10000 == 0) {std::cout << "Analyzed " << iEntry << " events of total " << entries << std::endl;}
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
    	
    }

    // write out Ntuples to file
    outFile->cd();
    for (int mmIndex=0;mmIndex<numberOfMM;mmIndex++)
    {    
      gDirectory->WriteObject(theNtuples[mmIndex],TString(output_Names[mmIndex]));
    }
    //outFile->Write();

    delete inFile;
    delete outFile;

}