fillChain.cxx
| 1 |
#include "TSystem.h" |
|---|---|
| 2 |
#include "TFile.h" |
| 3 |
#include "TTree.h" |
| 4 |
#include "../src/TNeEvent.h" |
| 5 |
//#include "../../AculUtils/TELoss/TELoss.h"
|
| 6 |
|
| 7 |
using std::cout;
|
| 8 |
using std::endl;
|
| 9 |
|
| 10 |
void fillTree(Int_t nofile = 0, const Int_t noevents = 0) { |
| 11 |
|
| 12 |
TString inFile; |
| 13 |
// inFile.Form("~/data/exp1804/h5_14_%02d.root", nofile);
|
| 14 |
inFile.Form("~/data/exp1804/be_01_%02d.root", nofile);
|
| 15 |
// inFile.Form("~/data/exp1804/calib/si_20_03.root");
|
| 16 |
|
| 17 |
TString outFile; |
| 18 |
// outFile.Form("~/data/exp1804/h5_14_%02d_calib.root", nofile);
|
| 19 |
outFile.Form("~/data/exp1804/be_01_%02d_calib.root", nofile);
|
| 20 |
// outFile.Form("~/data/exp1804/calib/si_20_03_calib.root");
|
| 21 |
|
| 22 |
cout << "Input file: " << inFile << endl;
|
| 23 |
cout << "Output file: " << outFile << endl;
|
| 24 |
|
| 25 |
TFile *fr = new TFile(inFile);
|
| 26 |
TTree *tr = (TTree*)fr->Get("AnalysisxTree");
|
| 27 |
|
| 28 |
TNeEvent *revent = new TNeEvent();
|
| 29 |
|
| 30 |
tr->SetBranchAddress("NeEvent.", &revent);
|
| 31 |
|
| 32 |
TFile *fw = new TFile(outFile, "RECREATE"); |
| 33 |
TTree *tw = new TTree("cal", "Calibrated information"); |
| 34 |
|
| 35 |
Float_t SQ20E[16];
|
| 36 |
Float_t SQ20Ecorr[16];
|
| 37 |
Float_t SQ20Esum; |
| 38 |
Float_t SQLXE[32];
|
| 39 |
Float_t SQLXEsum; |
| 40 |
Float_t SQLYE[16];
|
| 41 |
Float_t SQLYEsum; |
| 42 |
|
| 43 |
Int_t SQLXmult; |
| 44 |
Int_t SQLYmult; |
| 45 |
Int_t SQ20EcorrMult; |
| 46 |
|
| 47 |
Float_t SQRXE[32];
|
| 48 |
Float_t SQRXEsum; |
| 49 |
Int_t SQRXmult; |
| 50 |
Float_t TOF, dEbeam; |
| 51 |
Float_t x1p, y1p, x2p, y2p, xt, yt; |
| 52 |
|
| 53 |
tw->Branch("SQ20E",SQ20E,"SQ20E[16]/F"); |
| 54 |
tw->Branch("SQ20Ecorr",SQ20Ecorr,"SQ20Ecorr[16]/F"); |
| 55 |
tw->Branch("SQ20Esum",&SQ20Esum,"SQ20Esum/F"); |
| 56 |
tw->Branch("SQLXE",SQLXE,"SQLXE[32]/F"); |
| 57 |
tw->Branch("SQLXEsum",&SQLXEsum,"SQLXE/F"); |
| 58 |
tw->Branch("SQLYE",SQLYE,"SQLYE[16]/F"); |
| 59 |
tw->Branch("SQLYEsum",&SQLYEsum,"SQLYEsum/F"); |
| 60 |
|
| 61 |
tw->Branch("SQLXmult",&SQLXmult,"SQLXmult/I"); |
| 62 |
tw->Branch("SQLYmult",&SQLYmult,"SQLYmult/I"); |
| 63 |
tw->Branch("SQ20EcorrMult",&SQ20EcorrMult,"SQ20EcorrMult/I"); |
| 64 |
|
| 65 |
tw->Branch("SQRXE",SQRXE,"SQRXE[32]/F"); |
| 66 |
tw->Branch("SQRXEsum",&SQRXEsum,"SQRXE/F"); |
| 67 |
tw->Branch("SQRXmult",&SQRXmult,"SQRXmult/I"); |
| 68 |
|
| 69 |
tw->Branch("TOF", &TOF, "TOF/F"); |
| 70 |
tw->Branch("dEbeam", &dEbeam, "dEbeam/F"); |
| 71 |
|
| 72 |
tw->Branch("x1p", &x1p, "x1p/F"); |
| 73 |
tw->Branch("y1p", &y1p, "y1p/F"); |
| 74 |
tw->Branch("x2p", &x2p, "x2p/F"); |
| 75 |
tw->Branch("y2p", &y2p, "y2p/F"); |
| 76 |
|
| 77 |
tw->Branch("xt", &xt, "xt/F"); |
| 78 |
tw->Branch("yt", &yt, "yt/F"); |
| 79 |
|
| 80 |
|
| 81 |
Long64_t nevents; |
| 82 |
if (noevents == 0) nevents = tr->GetEntries(); |
| 83 |
else nevents = noevents;
|
| 84 |
if (nevents > tr->GetEntries()) nevents = tr->GetEntries();
|
| 85 |
|
| 86 |
// TNeDet16 *pSQX_L_EC = new TNeDet16("SQX_L_EC");
|
| 87 |
// TNeDet16 pSQX_L_EC("../SQX_L_EC");
|
| 88 |
TNeDet16 pSQX_L_EC("./SQX_L");
|
| 89 |
pSQX_L_EC.ReadData(); |
| 90 |
|
| 91 |
TNeDet16 pSQX_L_TC("./SQX_L_time");
|
| 92 |
pSQX_L_TC.ReadData(); |
| 93 |
|
| 94 |
TNeDet16 pSQY_L_EC("./SQY_L");
|
| 95 |
pSQY_L_EC.ReadData(); |
| 96 |
|
| 97 |
TNeDet16 pSQX_R_EC("../SQX_R_EC");
|
| 98 |
pSQX_R_EC.ReadData(); |
| 99 |
|
| 100 |
TNeDet16 pSQ20_EC("./sq20_58");
|
| 101 |
pSQ20_EC.ReadData(); |
| 102 |
|
| 103 |
// for (Int_t i = 0; i < 32; i++) {
|
| 104 |
// cout << pSQX_L_EC.Energy(1, i) << endl;
|
| 105 |
// }
|
| 106 |
|
| 107 |
Float_t energy = 0;
|
| 108 |
|
| 109 |
cout << nevents << " entries will be treated." << endl;
|
| 110 |
|
| 111 |
// TFile fThickness("thicknessPreliminary.root", "OPEN");
|
| 112 |
TFile fThickness("thicknessFinal.root", "OPEN"); |
| 113 |
fr->cd(); |
| 114 |
TH2F *hThickness = new TH2F(*(TH2F*)fThickness.Get("hTh")); |
| 115 |
hThickness->Draw("col");
|
| 116 |
|
| 117 |
// std::cout << std::setprecision(1) << std::fixed;
|
| 118 |
|
| 119 |
|
| 120 |
const Int_t kSQL_X_strips = 32; |
| 121 |
const Int_t kSQL_Y_strips = 16; |
| 122 |
const Int_t kSQL_20_strips = 16; |
| 123 |
|
| 124 |
const Double_t kSQLY_energy_thr = 1.; |
| 125 |
const Double_t kSQLX_energy_thr = 1.; |
| 126 |
const Double_t kSQL20_energy_thr = 1.2; |
| 127 |
|
| 128 |
const Double_t kSQ20_norm_thickness = 20.; |
| 129 |
|
| 130 |
// for (Int_t yi = 0; yi < kSQL_Y_strips; yi++) {
|
| 131 |
// cout << "y bin: " << yi+1 << "\t\t";
|
| 132 |
// for (Int_t xi = 0; xi < kSQL_X_strips; xi++) {
|
| 133 |
// cout << hThickness->GetBinContent(xi+1, yi+1) << "\t";
|
| 134 |
// if (xi == kSQL_X_strips-1) cout << endl;
|
| 135 |
// }
|
| 136 |
// }
|
| 137 |
|
| 138 |
|
| 139 |
fw->cd(); |
| 140 |
|
| 141 |
|
| 142 |
// TELoss fAlphaSi;
|
| 143 |
// fAlphaSi.SetEL(1, 2.321); // density in g/cm3
|
| 144 |
// fAlphaSi.AddEL(14., 28.086, 1); //Z, mass
|
| 145 |
// // mSi.SetZP(1., 1.); //protons
|
| 146 |
// fAlphaSi.SetZP(2., 3.); //3He, Z, A
|
| 147 |
// fAlphaSi.SetEtab(100000, 100.); // ?, MeV calculate ranges
|
| 148 |
// fAlphaSi.SetDeltaEtab(300);
|
| 149 |
|
| 150 |
//////////////////////////////////
|
| 151 |
//event processing
|
| 152 |
//////////////////////////////////
|
| 153 |
|
| 154 |
for (Int_t eventNo = 0; eventNo < nevents; eventNo++) { |
| 155 |
tr->GetEvent(eventNo); |
| 156 |
// cout << revent->SQX_L[0] << endl;
|
| 157 |
// SQLXE[0] = revent->SQX_L[0];
|
| 158 |
|
| 159 |
SQ20Esum = 0.;
|
| 160 |
SQLXEsum = 0.;
|
| 161 |
SQLYEsum = 0.;
|
| 162 |
SQRXEsum = 0.;
|
| 163 |
SQRXmult = 0;
|
| 164 |
|
| 165 |
SQLXmult = 0;
|
| 166 |
SQLYmult = 0;
|
| 167 |
SQ20EcorrMult = 0;
|
| 168 |
|
| 169 |
|
| 170 |
for (Int_t i = 0; i < 32; i++) { |
| 171 |
// cout << pSQX_L_EC.Energy(1, i) << endl;
|
| 172 |
SQLXE[i] = 0;
|
| 173 |
SQRXE[i] = 0;
|
| 174 |
|
| 175 |
energy = pSQX_L_EC.Energy(revent->SQX_L[i], i); |
| 176 |
|
| 177 |
// cout << energy << "\t" << pSQX_L_EC.a[i][0] << "\t" << pSQX_L_EC.a[i][1] << "\t" << pSQX_L_EC.a[i][2] << endl;
|
| 178 |
// cout << energy << endl;
|
| 179 |
|
| 180 |
// if (i == 5 && i == 6) continue;
|
| 181 |
|
| 182 |
if (energy > kSQLX_energy_thr) {
|
| 183 |
SQLXE[i] = energy; |
| 184 |
SQLXEsum += SQLXE[i]; |
| 185 |
SQLXmult++; |
| 186 |
} |
| 187 |
|
| 188 |
energy = pSQX_R_EC.Energy(revent->SQX_R[i], i); |
| 189 |
|
| 190 |
if (energy > 1.0) { |
| 191 |
SQRXE[i] = energy; |
| 192 |
SQRXEsum += SQRXE[i]; |
| 193 |
SQRXmult++; |
| 194 |
} |
| 195 |
} |
| 196 |
|
| 197 |
for (Int_t i = 0; i < 16; i++) { |
| 198 |
// cout << pSQX_L_EC.Energy(1, i) << endl;
|
| 199 |
|
| 200 |
SQLYE[i] = 0.;
|
| 201 |
SQ20E[i] = 0.;
|
| 202 |
SQ20Ecorr[i] = 0.;
|
| 203 |
// SQRXE[i] = 0;
|
| 204 |
|
| 205 |
energy = pSQ20_EC.Energy(revent->SQ20[i], i); |
| 206 |
|
| 207 |
// if (energy > 1.) {
|
| 208 |
SQ20E[i] = energy; |
| 209 |
SQ20Esum += SQ20E[i]; |
| 210 |
// cout << energy << endl;
|
| 211 |
// cout << i << "\t" << energy << "\t" << pSQ20_EC.a[i][0] << "\t" << pSQ20_EC.a[i][1] << "\t" << pSQ20_EC.a[i][2] << endl;
|
| 212 |
// }
|
| 213 |
|
| 214 |
|
| 215 |
if (i==0) continue; |
| 216 |
|
| 217 |
energy = pSQY_L_EC.Energy(revent->SQY_L[i], i); |
| 218 |
|
| 219 |
// cout << energy << "\t" << pSQY_L_EC.a[i][0] << "\t" << pSQY_L_EC.a[i][1] << "\t" << pSQY_L_EC.a[i][2] << endl;
|
| 220 |
// cout << energy << endl;
|
| 221 |
|
| 222 |
if (energy > kSQLY_energy_thr) {
|
| 223 |
SQLYE[i] = energy; |
| 224 |
SQLYEsum += SQLYE[i]; |
| 225 |
SQLYmult++; |
| 226 |
} |
| 227 |
|
| 228 |
// energy = pSQX_R_EC.Energy(revent->SQX_R[i], i);
|
| 229 |
//
|
| 230 |
// if (energy > 1.0) {
|
| 231 |
// SQRXE[i] = energy;
|
| 232 |
// SQRXEsum += SQRXE[i];
|
| 233 |
// SQRXmult++;
|
| 234 |
// }
|
| 235 |
}//for 16
|
| 236 |
|
| 237 |
// cout << endl;
|
| 238 |
|
| 239 |
///////////////////////////////////////////////
|
| 240 |
//correction for thickness inhomogeneity
|
| 241 |
///////////////////////////////////////////////
|
| 242 |
|
| 243 |
Double_t currThickness; |
| 244 |
|
| 245 |
if (SQLXmult == 1 && SQLYmult == 1) { |
| 246 |
|
| 247 |
for (Int_t yi = 0; yi < kSQL_Y_strips; yi++) { |
| 248 |
// cout << "y bin: " << yi+1 << "\t\t";
|
| 249 |
if (SQLYE[yi]>kSQLY_energy_thr) {
|
| 250 |
|
| 251 |
for (Int_t xi = 0; xi < kSQL_X_strips; xi++) { |
| 252 |
|
| 253 |
if (SQLXE[xi] > kSQLX_energy_thr) {
|
| 254 |
currThickness = hThickness->GetBinContent(xi+1, yi+1); |
| 255 |
// if (currThickness<30.) {
|
| 256 |
// Int_t probableThinStrip = xi-2/2;
|
| 257 |
// SQ20Ecorr[0] =
|
| 258 |
// }
|
| 259 |
|
| 260 |
for (Int_t x20 = 0; x20 < kSQL_20_strips; x20++) { |
| 261 |
if (SQ20E[x20] > kSQL20_energy_thr && currThickness<30.) { |
| 262 |
// SQ20Ecorr[x20] = SQ20E[x20] + 1.;
|
| 263 |
|
| 264 |
// currThickness = kSQ20_norm_thickness + (currThickness - kSQ20_norm_thickness)*0.5;
|
| 265 |
|
| 266 |
SQ20Ecorr[x20] = (kSQ20_norm_thickness/currThickness)*SQ20E[x20]; |
| 267 |
SQ20EcorrMult++; |
| 268 |
|
| 269 |
// cout << xi << "\t" << yi << "\t" << x20
|
| 270 |
// << "\t" << currThickness
|
| 271 |
// << "\t" << SQ20E[x20] << "\t" << SQ20Ecorr[x20] << "\t" << endl;
|
| 272 |
} |
| 273 |
}//for*/
|
| 274 |
} |
| 275 |
|
| 276 |
|
| 277 |
|
| 278 |
// cout << hThickness->GetBinContent(xi+1, yi+1) << "\t";
|
| 279 |
// if (xi == kSQL_X_strips-1) cout << endl;
|
| 280 |
} |
| 281 |
} |
| 282 |
} |
| 283 |
|
| 284 |
} |
| 285 |
|
| 286 |
|
| 287 |
///////////////////////////////////////////////
|
| 288 |
//beam diagnostics
|
| 289 |
///////////////////////////////////////////////
|
| 290 |
|
| 291 |
dEbeam = 0.;
|
| 292 |
TOF = 0.;
|
| 293 |
if (revent->tF5[0]>0 && revent->tF5[1]>0 && revent->tF5[2]>0 && revent->tF5[3]>0 |
| 294 |
&& revent->tF3[0]>0 && revent->tF3[1]>0 && revent->tF3[2]>0 && revent->tF3[3]>0 |
| 295 |
&& revent->F5[0]>0 && revent->F5[1]>0 && revent->F5[2]>0 && revent->F5[3]>0 |
| 296 |
// && ( (revent->tF5[0]+revent->tF5[1]+revent->tF5[2]+revent->tF5[3]) - (revent->tF3[0]+revent->tF3[1]+revent->tF3[2]+revent->tF3[3]) )/4.*0.125+89.165<200
|
| 297 |
// && ( (tF5[0]+tF5[1]+tF5[2]+tF5[3]) - (tF3[0]+tF3[1]+tF3[2]+tF3[3]) )/4.*0.125+89.165>100
|
| 298 |
// && (F5[0]+F5[1]+F5[2]+F5[3])/4. < 2500
|
| 299 |
) {
|
| 300 |
dEbeam = (revent->F5[0]+revent->F5[1]+revent->F5[2]+revent->F5[3])/4.; |
| 301 |
TOF = ( (revent->tF5[0]+revent->tF5[1]+revent->tF5[2]+revent->tF5[3]) - (revent->tF3[0]+revent->tF3[1]+revent->tF3[2]+revent->tF3[3]) )/4.*0.125+89.165; |
| 302 |
} |
| 303 |
|
| 304 |
const Float_t l12 = 546.; |
| 305 |
const Float_t lt = 270.; |
| 306 |
|
| 307 |
x1p = -100.; |
| 308 |
y1p = -100.; |
| 309 |
x2p = -100.; |
| 310 |
y2p = -100.; |
| 311 |
xt = -100.; |
| 312 |
yt = -100.; |
| 313 |
|
| 314 |
if (revent->x1[0]<1000 && revent->y1[0]<1000 && revent->nx1==1 && revent->ny1==1) { |
| 315 |
x1p = revent->x1[0]*1.25-20.; |
| 316 |
y1p = revent->y1[0]*1.25-20.; |
| 317 |
} |
| 318 |
|
| 319 |
if (revent->x2[0]<1000 && revent->y2[0]<1000 && revent->nx2==1 && revent->ny2==1) { |
| 320 |
x2p = revent->x2[0]*1.25-20.; |
| 321 |
y2p = revent->y2[0]*1.25-20.; |
| 322 |
} |
| 323 |
|
| 324 |
if (x1p > -50. && y1p > -50. && x2p > -50. && y2p > -50.) { |
| 325 |
xt = x1p - (x1p - x2p)*(l12/lt); |
| 326 |
yt = y1p - (y1p - y2p)*(l12/lt); |
| 327 |
} |
| 328 |
|
| 329 |
|
| 330 |
// cout << endl;
|
| 331 |
|
| 332 |
tw->Fill(); |
| 333 |
|
| 334 |
// cout << eventNo << "\t" << eventNo%100 << endl;
|
| 335 |
if (eventNo%100000 == 0 && eventNo !=0) { |
| 336 |
cout << eventNo << " events processed." << endl;
|
| 337 |
} |
| 338 |
} |
| 339 |
|
| 340 |
cout << nevents << " entries processed." << endl;
|
| 341 |
|
| 342 |
fw->cd(); |
| 343 |
tw->Write(); |
| 344 |
fw->Close(); |
| 345 |
|
| 346 |
} |
| 347 |
|
| 348 |
|
| 349 |
void fillChain(const Int_t from = 0, const Int_t to = 9, const Int_t noevents = 0) { |
| 350 |
for (Int_t i = from; i <= to; i++) {
|
| 351 |
fillTree(i, noevents); |
| 352 |
} |
| 353 |
} |