Many modifications made in long period.

analysis/fillChain.cxx

@@ -8,19 +8,52 @@
 using std::cout;
 using std::endl;
 
-void fillTree(Int_t nofile = 0, const Int_t noevents = 0) {
+Int_t GetClustersMWPC(unsigned short n, unsigned short *x)
+{
+	Int_t clusters = 1;
+
+	if(n == 0) clusters = 0.;
+	else if(n > 1 && n<=32)
+	{
+		for(Int_t i = 1; i < n; ++i)
+		{
+//			if (x[i]>50 || x[i-1]>50) continue;
+			if( (abs(x[i] - x[i-1]) > 1) && (abs(x[i] - x[i-1]) < 32) ) clusters += 1.;
+		}
+	}
+	return clusters;
+} //--------------------------------------------------------------------
+
+Double_t GetClusterPositionMWPC(unsigned short n, unsigned short *x, Float_t planeOffset=0.)
+{
+	Double_t position = -100.;
+
+
+
+
+//	if(n == 0) clusters = 0.;
+	if(n > 0 && n<=32)
+	{
+		position = (x[0]+n/2.+0.5)*1.25-20. + planeOffset;
+//		cout << n << endl;
+	}
+	return position;
+} //--------------------------------------------------------------------
+
+
+void fillTree(const TString beam = "he", Int_t nofile = 0, const Int_t noevents = 0) {
 
 	TString inFile;
-	inFile.Form("~/data/exp1804/h5_14_%02d.root", nofile);
-//	inFile.Form("~/data/exp1804/be10_03_%d0.root", nofile);		//files 00,10,...,90
+	if (beam.Contains("he")) inFile.Form("~/data/exp1804/h5_14_%02d.root", nofile);
+	if (beam.Contains("be")) inFile.Form("~/data/exp1804/be10_03_%d0.root", nofile);		//files 00,10,...,90
 																//where 70 is run 01
 																//		80 is run 02
 																//		90 is run 05
 //	inFile.Form("~/data/exp1804/calib/si_20_03.root");
 
 	TString outFile;
-	outFile.Form("~/data/exp1804/h5_14_%02d_calib.root", nofile);
-//	outFile.Form("~/data/exp1804/be10_03_%d0_calib.root", nofile);	//files 00,10,...,60
+	if (beam.Contains("he")) outFile.Form("~/data/exp1804/h5_14_%02d_calib.root", nofile);
+	if (beam.Contains("be")) outFile.Form("~/data/exp1804/be10_03_%d0_calib.root", nofile);	//files 00,10,...,60
 																	//where 70 is run 01
 																	//		80 is run 02
 																	//		90 is run 05
@@ -39,9 +72,12 @@ void fillTree(Int_t nofile = 0, const Int_t noevents = 0) {
 	TFile *fw = new TFile(outFile, "RECREATE");
 	TTree *tw = new TTree("cal", "Calibrated information");
 
+	Int_t trigger;
+
 	Float_t SQ20E[16];
 	Float_t SQ20Ecorr[16];
 	Float_t SQ20EcorrHit[16];
+	Float_t SQ20EcorrHitC[16];
 	Float_t SQ20Esum;
 
 	Float_t SQ20time[16];
@@ -68,13 +104,26 @@ void fillTree(Int_t nofile = 0, const Int_t noevents = 0) {
 	Int_t SQLYmult;
 	Int_t SQ20EcorrMult;
 	Int_t SQ20EcorrHitMult;
+	Int_t SQ20EcorrHitCMult;
+
+	Bool_t CsI_L_veto;
+	Int_t expectedThinStrip;
+
 
 	Float_t SQRXE[32];
 	Float_t SQRXEsum;
 	Int_t SQRXmult;
 	Float_t TOF, dEbeam;
+
+	Float_t tF5;
+
+	//MWPC, wire multiplicity
 	Float_t x1p, y1p, x2p, y2p, xt, yt;
 
+	//MWPC, cluster multiplicity
+	Int_t x1MultC, x2MultC, y1MultC, y2MultC;
+	Float_t x1c, y1c, x2c, y2c, xtc, ytc;
+
 
 	//left 1 mm position
 	Float_t x1mm, y1mm;
@@ -91,6 +140,8 @@ void fillTree(Int_t nofile = 0, const Int_t noevents = 0) {
 	TVector3 vNorm(0.,0.,1.);
 	Double_t angleLeft;
 
+	tw->Branch("trigger",&trigger,"trigger/I");
+
 	tw->Branch("angleLeft",&angleLeft,"angleLeft/D");
 	tw->Branch("x1mm",&x1mm,"x1mm/F");
 	tw->Branch("y1mm",&y1mm,"y1mm/F");
@@ -103,6 +154,7 @@ void fillTree(Int_t nofile = 0, const Int_t noevents = 0) {
 	tw->Branch("SQ20E",SQ20E,"SQ20E[16]/F");
 	tw->Branch("SQ20Ecorr",SQ20Ecorr,"SQ20Ecorr[16]/F");
 	tw->Branch("SQ20EcorrHit",SQ20EcorrHit,"SQ20EcorrHit[16]/F");
+	tw->Branch("SQ20EcorrHitC",SQ20EcorrHitC,"SQ20EcorrHitC[16]/F");
 	tw->Branch("SQ20Esum",&SQ20Esum,"SQ20Esum/F");
 
 	tw->Branch("SQ20time",SQ20time,"SQ20time[16]/F");
@@ -130,6 +182,11 @@ void fillTree(Int_t nofile = 0, const Int_t noevents = 0) {
 	tw->Branch("SQLXmult",&SQLXmult,"SQLXmult/I");
 	tw->Branch("SQLYmult",&SQLYmult,"SQLYmult/I");
 	tw->Branch("SQ20EcorrMult",&SQ20EcorrMult,"SQ20EcorrMult/I");
+	tw->Branch("SQ20EcorrHitCMult",&SQ20EcorrHitCMult,"SQ20EcorrHitCMult/I");
+//	tw->Branch("SQ20EcorrMult",&SQ20EcorrMult,"SQ20EcorrMult/I");
+
+	tw->Branch("CsI_L_veto", &CsI_L_veto, "CsI_L_veto/O");
+	tw->Branch("expectedThinStrip",&expectedThinStrip,"expectedThinStrip/I");
 
 	tw->Branch("SQRXE",SQRXE,"SQRXE[32]/F");
 	tw->Branch("SQRXEsum",&SQRXEsum,"SQRXE/F");
@@ -137,6 +194,7 @@ void fillTree(Int_t nofile = 0, const Int_t noevents = 0) {
 
 	tw->Branch("TOF", &TOF, "TOF/F");
 	tw->Branch("dEbeam", &dEbeam, "dEbeam/F");
+	tw->Branch("tF5", &tF5, "tF5/F");
 
 	tw->Branch("x1p", &x1p, "x1p/F");
 	tw->Branch("y1p", &y1p, "y1p/F");
@@ -146,6 +204,19 @@ void fillTree(Int_t nofile = 0, const Int_t noevents = 0) {
 	tw->Branch("xt", &xt, "xt/F");
 	tw->Branch("yt", &yt, "yt/F");
 
+	tw->Branch("x1MultC",&x1MultC,"x1MultC/I");
+	tw->Branch("x2MultC",&x2MultC,"x2MultC/I");
+	tw->Branch("y1MultC",&y1MultC,"y1MultC/I");
+	tw->Branch("y2MultC",&y2MultC,"y2MultC/I");
+
+	tw->Branch("x1c", &x1c, "x1c/F");
+	tw->Branch("y1c", &y1c, "y1c/F");
+	tw->Branch("x2c", &x2c, "x2c/F");
+	tw->Branch("y2c", &y2c, "y2c/F");
+
+	tw->Branch("xtc", &xtc, "xtc/F");
+	tw->Branch("ytc", &ytc, "ytc/F");
+
 
 	Long64_t nevents;
 	if (noevents == 0) nevents = tr->GetEntries();
@@ -190,9 +261,9 @@ void fillTree(Int_t nofile = 0, const Int_t noevents = 0) {
 	const Int_t kSQL_Y_strips = 16;
 	const Int_t kSQL_20_strips = 16;
 
-	const Double_t kSQLY_energy_thr = 1.;
-	const Double_t kSQLX_energy_thr = 1.;
-	const Double_t kSQL20_energy_thr = 1.2;
+//	const Double_t kSQLY_energy_thr = 1.;
+//	const Double_t kSQLX_energy_thr = 1.;
+//	const Double_t kSQL20_energy_thr = 1.2;
 
 	const Double_t kSQ20_norm_thickness = 20.;
 
@@ -210,6 +281,82 @@ void fillTree(Int_t nofile = 0, const Int_t noevents = 0) {
 							6, 4, -4, 13
 	};
 
+	Int_t thinChThresholds[16] = {
+		113, 121, 125, 125,
+		125, 132, 130, 115,
+		118, 128, 132, 130,
+		130, 130, 130, 112
+	};
+
+	Int_t thickLYChThresholds[16] = {
+		4000, 130, 82, 110,
+		87, 112, 82, 108,
+		82, 112, 82, 108,
+		80, 110, 90, 105,
+	};
+
+	Int_t thickLXChThresholds[32] = {
+		60, 45, 45, 45,
+		45, 4000, 4000, 45,
+		45, 45, 45, 45,
+		45, 45, 45, 45,
+		/////
+		45, 45, 45, 45,
+		45, 45, 45, 45,
+		45, 45, 45, 45,
+		45, 45, 45, 60
+	};
+
+	Int_t minTimeLthin[16] {
+		410, 410, 420, 415,
+		420, 410, 415, 415,
+		420, 420, 420, 420,
+		410, 425, 0, 0
+	};
+
+	Int_t maxTimeLthin[16] {
+		520, 525, 520, 530,
+		520, 530, 525, 525,
+		530, 530, 530, 520,
+		510, 525, 0, 0
+	};
+
+	Double_t minTimeLX[32] {
+		0, 0, 325, 328,
+		330, 0, 0, 330,
+		326, 328, 326, 335,
+		326, 328, 328, 330,
+		//
+		420, 425, 420, 420,
+		420, 420, 420, 420,
+		430, 425, 420, 420,
+		415, 415, 425, 415
+	};
+
+	Double_t maxTimeLX[32] {
+		0, 0, 345, 345,
+		350, 0, 0, 350,
+		345, 345, 345, 350,
+		340, 342, 345, 345,
+		//
+		460, 460, 460, 460,
+		460, 460, 460, 460,
+		460, 460, 460, 460,
+		460, 460, 460, 453
+	};
+
+	const UShort_t CsIleftThr = 180;
+
+	const Double_t thinXoffset = 1.;
+	const Double_t thinYoffset = -1.8;
+
+	const Float_t MWPC1_X_offset = -1.0;
+	const Float_t MWPC1_Y_offset = -2.1375;
+	const Float_t MWPC2_X_offset = 0.2;
+	const Float_t MWPC2_Y_offset = -1.125;
+
+
+
 	fw->cd();
 
 	//////////////////////////////////
@@ -222,6 +369,10 @@ void fillTree(Int_t nofile = 0, const Int_t noevents = 0) {
 
 		tr->GetEvent(eventNo);
 
+		trigger = revent->trigger;
+
+		if (trigger != 3) continue;
+
 //		cout << eventNo << endl;
 //		cout << revent->SQX_L[0] << endl;
 //		SQLXE[0] = revent->SQX_L[0];
@@ -239,16 +390,143 @@ void fillTree(Int_t nofile = 0, const Int_t noevents = 0) {
 		SQLYmult = 0;
 		SQ20EcorrMult = 0;
 		SQ20EcorrHitMult = 0;
+		SQ20EcorrHitCMult = 0;
 
 		SQ20timeMult = 0;
 		SQLXtimeMult = 0;
 //		SQLXtimeSum = 0;
 		SQLYtimeMult = 0;
 
+		CsI_L_veto = kFALSE;
+		expectedThinStrip = -1;
+
+
+		///////////////////////////////////////////////
+		//beam diagnostics
+		///////////////////////////////////////////////
+
+			///////////////////////////////////////////////
+			//TOF
+			///////////////////////////////////////////////
+
+		dEbeam = 0.;
+		TOF = 0.;
+		if (revent->tF5[0]>0 && revent->tF5[1]>0 && revent->tF5[2]>0 && revent->tF5[3]>0
+				&& revent->tF3[0]>0 && revent->tF3[1]>0 && revent->tF3[2]>0 && revent->tF3[3]>0
+				&& revent->F5[0]>0 && revent->F5[1]>0 && revent->F5[2]>0 && revent->F5[3]>0
+//				&& ( (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
+//				&& ( (tF5[0]+tF5[1]+tF5[2]+tF5[3]) - (tF3[0]+tF3[1]+tF3[2]+tF3[3]) )/4.*0.125+89.165>100
+//				&& (F5[0]+F5[1]+F5[2]+F5[3])/4. < 2500
+		) {
+			dEbeam = (revent->F5[0]+revent->F5[1]+revent->F5[2]+revent->F5[3])/4.;
+			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;
+		}
+
+		const Float_t l12 = 546.;
+		const Float_t lt = 270.;
+
+			///////////////////////////////////////////////
+			//MWPC's
+			///////////////////////////////////////////////
+
+
+		x1p = -100.;
+		y1p = -100.;
+		x2p = -100.;
+		y2p = -100.;
+		xt = -100.;
+		yt = -100.;
+
+		x1c = -100.;
+		y1c = -100.;
+		x2c = -100.;
+		y2c = -100.;
+		xtc = -100.;
+		ytc = -100.;
+
+		x1MultC = 0;
+		y1MultC = 0;
+		x2MultC = 0;
+		y2MultC = 0;
+
+
+		Bool_t flagMWPC = 1;
+
+		tF5 = 0.125*(revent->tF5[0]+revent->tF5[1]+revent->tF5[2]+revent->tF5[3])/4.;
+//		tF5 = 0.125*revent->tF5[0];
+
+		if(revent->nx1>10 || revent->ny1>10 || revent->nx2>10 || revent->ny2>10) flagMWPC=0;
+		if((0.125*revent->tMWPC[0]-tF5)<60 || (0.125*revent->tMWPC[0]-tF5)>77) flagMWPC=0;
+		if((0.125*revent->tMWPC[1]-tF5)<60 || (0.125*revent->tMWPC[1]-tF5)>80) flagMWPC=0;
+		if((0.125*revent->tMWPC[2]-tF5)<70 || (0.125*revent->tMWPC[2]-tF5)>90) flagMWPC=0;
+		if((0.125*revent->tMWPC[3]-tF5)<60 || (0.125*revent->tMWPC[3]-tF5)>80) flagMWPC=0;
+
+		//MWPC: work with wire multiplicity equal to 1
+		if (flagMWPC!=0) {
+			if (revent->x1[0]<1000 && revent->y1[0]<1000 && revent->nx1==1 && revent->ny1==1) {
+				x1p = (revent->x1[0]+0.5)*1.25-20. + MWPC1_X_offset;
+				y1p = (revent->y1[0]+0.5)*1.25-20. + MWPC1_Y_offset;
+
+//						(x[0]+n/2.+0.5)*1.25-20.;
+			}
+
+			if (revent->x2[0]<1000 && revent->y2[0]<1000 && revent->nx2==1 && revent->ny2==1) {
+				x2p = (revent->x2[0]+0.5)*1.25-20. + MWPC2_X_offset;
+				y2p = (revent->y2[0]+0.5)*1.25-20. + MWPC2_Y_offset;
+			}
+
+			if (x1p > -50. && y1p > -50. && x2p > -50. && y2p > -50.) {
+				xt = x1p - (x1p - x2p)*((l12+lt)/l12);
+				yt = y1p - (y1p - y2p)*((l12+lt)/l12);
+//				xt = (l12*x1p + (l12+lt)*(x2p-x1p))/(l12 - (x2p-x1p)*TMath::Tan(TMath::DegToRad()*12.));
+//				yt = (y2p-y1p)*(xt-x1p)/(x2p-x1p) + y1p;
+			}
+		}//if wire multiplicity == 1
+
+//		cout << flagMWPC << endl;
+
+		//MWPC: work with cluster multiplicity equal to 2
+		if (flagMWPC!=0) {
+			x1MultC = GetClustersMWPC(revent->nx1, revent->x1);
+			x2MultC = GetClustersMWPC(revent->nx2, revent->x2);
+			y1MultC = GetClustersMWPC(revent->ny1, revent->y1);
+			y2MultC = GetClustersMWPC(revent->ny2, revent->y2);
+
+			if (x1MultC==1 && y1MultC==1) {
+				x1c = GetClusterPositionMWPC(revent->nx1, revent->x1, MWPC1_X_offset);
+				y1c = GetClusterPositionMWPC(revent->ny1, revent->y1, MWPC1_Y_offset);
+			}
+
+			if (x2MultC==1 && y2MultC==1) {
+				x2c = GetClusterPositionMWPC(revent->nx2, revent->x2, MWPC2_X_offset);
+				y2c = GetClusterPositionMWPC(revent->ny2, revent->y2, MWPC2_Y_offset);
+			}
+
+			if (x1c > -50. && y1c > -50. && x2c > -50. && y2c > -50.) {
+				xtc = x1c - (x1c - x2c)*((l12+lt)/l12);
+				ytc = y1c - (y1c - y2c)*((l12+lt)/l12);
+//				xtc = (l12*x1c + (l12+lt)*(x2c-x1c))/(l12 - (x2c-x1c)*TMath::Tan(TMath::DegToRad()*12.));
+//				ytc = (y2c-y1c)*(xtc-x1c)/(x2c-x1c) + y1c;
+			}
+		}
+
 
+		if (x1MultC!=1 || y1MultC!=1 || x2MultC!=1 || y2MultC!=1) continue;
+
+
+
+		///////////////////////////////////////////////
+		//work with Si detectors
+		///////////////////////////////////////////////
 
 //cout << endl;
 
+		//CdI left -veto:
+		for(Int_t i = 0; i<16; i++) {
+			if (revent->CsI_L[i]>CsIleftThr) CsI_L_veto = kTRUE;
+		}
+
+
 		for (Int_t i = 0; i < 32; i++) {
 //			cout << pSQX_L_EC.Energy(1, i) << endl;
 			SQLXE[i] = 0;
@@ -256,19 +534,16 @@ void fillTree(Int_t nofile = 0, const Int_t noevents = 0) {
 			SQLXtime[i] = 0;
 			SQLXEtimeFiltered[i] = 0.;
 
-			energy = pSQX_L_EC.Energy(revent->SQX_L[i], i);
-
-//			if (i == 5 && i == 6) continue;
-
-			if (energy > kSQLX_energy_thr) {
+			if (revent->SQX_L[i] > thickLXChThresholds[i]) {
+				energy = pSQX_L_EC.Energy(revent->SQX_L[i], i);
 				SQLXE[i] = energy;
 				SQLXEsum += SQLXE[i];
 				SQLXmult++;
 			}
 
-			if (i<16 && i!=5 && i!=6 && i!=0 && i!=1
-					&& SQLXE[i]>kSQLX_energy_thr
-					&& revent->tSQX_L[i]*0.3>320 && revent->tSQX_L[i]*0.3<380) {
+			if (i<16 && i!=0 && i!=1
+					&& SQLXE[i]>0
+					&& revent->tSQX_L[i]*0.3>minTimeLX[i] && revent->tSQX_L[i]*0.3<maxTimeLX[i]) {
 				SQLXtime[i] = revent->tSQX_L[i]*0.3;
 				SQLXtimeMult++;
 //				SQLXtimeSum = SQLXtimeSum + SQLXtime[i];
@@ -278,8 +553,8 @@ void fillTree(Int_t nofile = 0, const Int_t noevents = 0) {
 //				cout << i << endl;
 			}
 			if (i>15
-					&& SQLXE[i]>kSQLX_energy_thr
-					&& revent->tSQX_L[i]*0.3>410 && revent->tSQX_L[i]*0.3<470) {
+					&& SQLXE[i]>0
+					&& revent->tSQX_L[i]*0.3>minTimeLX[i] && revent->tSQX_L[i]*0.3<maxTimeLX[i]) {
 				SQLXtime[i] = revent->tSQX_L[i]*0.3;
 				SQLXtimeMult++;
 //				SQLXtimeSum = SQLXtimeSum + SQLXtime[i];
@@ -295,7 +570,7 @@ void fillTree(Int_t nofile = 0, const Int_t noevents = 0) {
 				SQRXEsum += SQRXE[i];
 				SQRXmult++;
 			}
-		}
+		}// for 32
 
 		for (Int_t i = 0; i < 16; i++) {
 			//			cout << pSQX_L_EC.Energy(1, i) << endl;
@@ -305,32 +580,33 @@ void fillTree(Int_t nofile = 0, const Int_t noevents = 0) {
 			SQ20time[i] = 0.;
 			SQ20Ecorr[i] = 0.;
 			SQ20EcorrHit[i] = 0.;
+			SQ20EcorrHitC[i] = 0.;
 			//				SQRXE[i] = 0;
 
 			SQLYtime[i] = 0.;
 			SQLYEtimeFiltered[i] = 0.;
 
-			energy = pSQ20_EC.Energy(revent->SQ20[i], i);
-
-//			if (energy > 1.) {
+			if (revent->SQ20[i] > thinChThresholds[i]) {
+				energy = pSQ20_EC.Energy(revent->SQ20[i], i);
 				SQ20E[i] = energy;
 				SQ20Esum += SQ20E[i];
 //				cout << energy << endl;
 //				cout << i << "\t" << energy << "\t" << pSQ20_EC.a[i][0] << "\t" << pSQ20_EC.a[i][1] << "\t" << pSQ20_EC.a[i][2] << endl;
-//			}
+			}
 
 			SQ20time[i] = revent->tSQ20[i]*0.3;
-			if (SQ20E[i]>kSQL20_energy_thr && SQ20time[i]>400 && SQ20time[i]<500) SQ20timeMult++;
+			if (SQ20E[i]>0 && SQ20time[i]>minTimeLthin[i] && SQ20time[i]<maxTimeLthin[i]) SQ20timeMult++;
+
 
+//			if (i==0) continue;
 
-			if (i==0) continue;
 
-			energy = pSQY_L_EC.Energy(revent->SQY_L[i], i);
 
 			//				cout << energy << "\t" << pSQY_L_EC.a[i][0] << "\t" << pSQY_L_EC.a[i][1] << "\t" << pSQY_L_EC.a[i][2] << endl;
 			//			cout << energy << endl;
 
-			if (energy > kSQLY_energy_thr) {
+			if (revent->SQY_L[i] > thickLYChThresholds[i]) {
+				energy = pSQY_L_EC.Energy(revent->SQY_L[i], i);
 				SQLYE[i] = energy;
 				SQLYEsum += SQLYE[i];
 				SQLYmult++;
@@ -338,23 +614,13 @@ void fillTree(Int_t nofile = 0, const Int_t noevents = 0) {
 
 			Double_t time = revent->tSQY_L[i]*0.3 + timeCorr[i];
 
-			if (SQLYE[i]>kSQLY_energy_thr && time>325 && time<333) {
-//				cout << "!!!!!!!!!!!!" << endl;
+			if (SQLYE[i]>0 && time>325 && time<333) {
 				SQLYtime[i] = time;
 				SQLYtimeMult++;
 				SQLYEtimeFiltered[i]=energy;
 				SQLYEtimeFilteredSum = SQLYEtimeFilteredSum + SQLYEtimeFiltered[i];
 			}
 
-
-
-			//				energy = pSQX_R_EC.Energy(revent->SQX_R[i], i);
-			//
-			//				if (energy > 1.0) {
-			//					SQRXE[i] = energy;
-			//					SQRXEsum += SQRXE[i];
-			//					SQRXmult++;
-			//				}
 		}//for 16
 
 //		cout << SQLXtimeMult << endl;
@@ -371,11 +637,11 @@ void fillTree(Int_t nofile = 0, const Int_t noevents = 0) {
 
 			for (Int_t yi = 0; yi < kSQL_Y_strips; yi++) {
 				//			cout << "y bin: " << yi+1 << "\t\t";
-				if (SQLYE[yi]>kSQLY_energy_thr) {
+				if (SQLYE[yi]>0) {
 
 					for (Int_t xi = 0; xi < kSQL_X_strips; xi++) {
 
-						if (SQLXE[xi] > kSQLX_energy_thr) {
+						if (SQLXE[xi] > 0) {
 							currThickness = hThickness->GetBinContent(xi+1, yi+1);
 //							if (currThickness<30.) {
 //								Int_t probableThinStrip = xi-2/2;
@@ -383,7 +649,7 @@ void fillTree(Int_t nofile = 0, const Int_t noevents = 0) {
 //							}
 
 							for (Int_t x20 = 0; x20 < kSQL_20_strips; x20++) {
-								if (SQ20E[x20] > kSQL20_energy_thr/*SQ20timeMult==1*/ && currThickness<30.) {
+								if (SQ20E[x20] > 0/*SQ20timeMult==1*/ && currThickness<30.) {
 //									SQ20Ecorr[x20] = SQ20E[x20] + 1.;
 
 //									currThickness = kSQ20_norm_thickness + (currThickness - kSQ20_norm_thickness)*0.5;
@@ -409,48 +675,6 @@ void fillTree(Int_t nofile = 0, const Int_t noevents = 0) {
 		}//if for correction according energy multiplicity
 
 
-		///////////////////////////////////////////////
-		//beam diagnostics
-		///////////////////////////////////////////////
-
-		dEbeam = 0.;
-		TOF = 0.;
-		if (revent->tF5[0]>0 && revent->tF5[1]>0 && revent->tF5[2]>0 && revent->tF5[3]>0
-				&& revent->tF3[0]>0 && revent->tF3[1]>0 && revent->tF3[2]>0 && revent->tF3[3]>0
-				&& revent->F5[0]>0 && revent->F5[1]>0 && revent->F5[2]>0 && revent->F5[3]>0
-//				&& ( (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
-//				&& ( (tF5[0]+tF5[1]+tF5[2]+tF5[3]) - (tF3[0]+tF3[1]+tF3[2]+tF3[3]) )/4.*0.125+89.165>100
-//				&& (F5[0]+F5[1]+F5[2]+F5[3])/4. < 2500
-		) {
-			dEbeam = (revent->F5[0]+revent->F5[1]+revent->F5[2]+revent->F5[3])/4.;
-			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;
-		}
-
-		const Float_t l12 = 546.;
-		const Float_t lt = 270.;
-
-		x1p = -100.;
-		y1p = -100.;
-		x2p = -100.;
-		y2p = -100.;
-		xt = -100.;
-		yt = -100.;
-
-		if (revent->x1[0]<1000 && revent->y1[0]<1000 && revent->nx1==1 && revent->ny1==1) {
-			x1p = revent->x1[0]*1.25-20.;
-			y1p = revent->y1[0]*1.25-20.;
-		}
-
-		if (revent->x2[0]<1000 && revent->y2[0]<1000 && revent->nx2==1 && revent->ny2==1) {
-			x2p = revent->x2[0]*1.25-20.;
-			y2p = revent->y2[0]*1.25-20.;
-		}
-
-		if (x1p > -50. && y1p > -50. && x2p > -50. && y2p > -50.) {
-			xt = x1p - (x1p - x2p)*(l12/lt);
-			yt = y1p - (y1p - y2p)*(l12/lt);
-		}
-
 
 
 		///////////////////////////////////////////////
@@ -469,56 +693,42 @@ void fillTree(Int_t nofile = 0, const Int_t noevents = 0) {
 
 		if (SQLXtimeMult==1) {
 			for (Int_t i = 0; i<32; i++) {
-				if (SQLXEtimeFiltered[i]>0) x1mm = -30. + (i+1/2)*60./32.;
+				if (SQLXEtimeFiltered[i]>0) x1mm = -30. + (i+1./2.)*60./32.;
 			}
 		}
 //		cout << SQLYtimeMult << endl;
 
 		if (SQLYtimeMult==1) {
 			for (Int_t i = 0; i<16; i++) {
-				if (SQLYEtimeFiltered[i]>0) y1mm = -30. + (i+1/2)*60./16.;
+				if (SQLYEtimeFiltered[i]>0) y1mm = -30. + (i+1./2.)*60./16.;
 			}
 		}
 
+		//for MWPC wire multiplicity == 1
 
 		if (SQLXtimeMult==1 && SQLYtimeMult==1 && xt>-50. && yt>-50.) {
 			vHit1mm.SetXYZ(x1mm-xt, y1mm-yt, z1mm);
-//			cout << vHit1mm.X() << "\t" << vHit1mm.Y() << endl;
-//			cout << "\t" << zThin/z1mm*vHit1mm.X() << "\t" << zThin/z1mm*vHit1mm.Y() << endl;
 			vHit1mm *= zThin/z1mm;
-//			cout << "\t" << vHit1mm.X() << "\t" << vHit1mm.Y() << "\t" << vHit1mm.Z() << endl;
-//			cout << "\t" << xt << "\t" << yt << "\t" << 0 << endl;
-//			cout << zThin/z1mm << endl;
 			angleLeft = vHit1mm.Angle(vNorm);
-//			xThin = vHit1mm.X()+xt;
-//			yThin = vHit1mm.Y()+yt;
-
-//			cout << angleLeft << endl;
 		}
 
-//		cout << angleLeft << endl;
-
-//		cout << x1mm << "\t" << y1mm << endl;
-
-//		cout << vHit1mm.X() << "\t" << vHit1mm.Y() << endl;
-
-
 
 		if (SQLXtimeMult==1 && SQLYtimeMult==1
 				&& vHit1mm.X()+xt>-30. && vHit1mm.X()+xt<30.
 				&& vHit1mm.Y()+yt>-30. && vHit1mm.Y()+yt<30.) {
-			mapYbin = (vHit1mm.Y()+yt+30.-0.5*60./16.)*16./60.;
-			mapXbin = (vHit1mm.X()+xt+30.-0.5*60./32.)*32./60.;
+
+
+			mapYbin = (vHit1mm.Y()+yt+thinYoffset+30.-0.5*60./16.)*16./60.;
+			mapXbin = (vHit1mm.X()+xt+thinXoffset+30.-0.5*60./32.)*32./60.;
+			//coordinate in the thin detector plane
 			xThin = vHit1mm.X()+xt;
 			yThin = vHit1mm.Y()+yt;
-//			cout << "\t\t" << xThin << "\t" << yThin << endl;
-
-//			cout << mapYbin << "\t" << "\t" << y1mm << "\t" << xt << endl;
 
 			//calculation of corrected energy
 			currThickness = hThickness->GetBinContent(mapXbin+1, mapYbin+1);
+
 			for (Int_t x20 = 0; x20 < kSQL_20_strips; x20++) {
-				if (SQ20E[x20] > kSQL20_energy_thr/*SQ20timeMult==1*/ && currThickness<30.) {
+				if (SQ20E[x20] > 0/*SQ20timeMult==1*/ && currThickness<30.) {
 					SQ20EcorrHit[x20] = (kSQ20_norm_thickness/currThickness)*SQ20E[x20];
 					SQ20EcorrHitMult++;
 //					cout << mapXbin << "\t" << mapYbin << "\t" << x20
@@ -530,13 +740,58 @@ void fillTree(Int_t nofile = 0, const Int_t noevents = 0) {
 
 		}//if inside the map
 
-//		cout << mapYbin << "\t" << "\t" << y1mm << endl;
-//		cout << mapXbin << "\t" << xThin << "\t" << x1mm << endl;
+		//end of MWPC wire multiplicity == 1
 
+		//for MWPC cluster multiplicity == 1
 
-//		cout << mapXbin << "\t" << mapYbin << endl;
+		if (SQLXtimeMult==1 && SQLYtimeMult==1 && xtc>-50. && ytc>-50.) {
+			vHit1mm.SetXYZ(x1mm-xtc, y1mm-ytc, z1mm);
+			vHit1mm *= zThin/z1mm;
+			angleLeft = vHit1mm.Angle(vNorm);
+		}
 
-		tw->Fill();
+
+		if (SQLXtimeMult==1 && SQLYtimeMult==1
+				&& vHit1mm.X()+xtc>-30. && vHit1mm.X()+xtc<30.
+				&& vHit1mm.Y()+ytc>-30. && vHit1mm.Y()+ytc<30.) {
+
+
+			mapYbin = (vHit1mm.Y()+ytc+thinYoffset+30.-0.5*60./16.)*16./60.;
+			mapXbin = (vHit1mm.X()+xtc+thinXoffset+30.-0.5*60./32.)*32./60.;
+			//coordinate in the thin detector plane
+			xThin = vHit1mm.X()+xtc;
+			yThin = vHit1mm.Y()+ytc;
+
+			//calculation of corrected energy
+			currThickness = hThickness->GetBinContent(mapXbin+1, mapYbin+1);
+
+			for (Int_t x20 = 0; x20 < kSQL_20_strips; x20++) {
+				if (SQ20E[x20] > 0/*SQ20timeMult==1*/ && currThickness<30.) {
+					SQ20EcorrHitC[x20] = (kSQ20_norm_thickness/currThickness)*SQ20E[x20];
+					SQ20EcorrHitCMult++;
+					//					cout << mapXbin << "\t" << mapYbin << "\t" << x20
+					//						<< "\t" << currThickness
+					//						<< "\t" << SQ20E[x20]  << "\t" << SQ20EcorrHit[x20] << "\t" << endl;
+
+				}
+			}//for*/
+
+			if (SQ20EcorrHitCMult==1) {
+				expectedThinStrip = (xThin + 25. - thinXoffset)*(16./50.)-0.5;
+			}
+
+		}//if inside the map
+
+		//end of MWPC wire multiplicity == 1
+
+		if (
+				SQLXtimeMult>0 && SQLXEtimeFilteredSum>0
+				&& trigger==3
+				&& x1MultC==1 && y1MultC==1 && x2MultC==1 && y2MultC==1
+		) {
+
+			tw->Fill();
+		}//if TTree::Fill
 
 		if (eventNo%100000 == 0 && eventNo !=0) {
 			cout << eventNo << " events processed." << endl;
@@ -552,8 +807,15 @@ void fillTree(Int_t nofile = 0, const Int_t noevents = 0) {
 }
 
 
-void fillChain(const Int_t from = 0, const Int_t to = 9, const Int_t noevents = 0) {
+void fillChain(const TString beam = "he", const Int_t from = 0, const Int_t to = 9, const Int_t noevents = 0) {
+
+	TStopwatch stopwatch;
+	stopwatch.Start();
+
 	for (Int_t i = from; i <= to; i++) {
-		fillTree(i, noevents);
+		fillTree(beam, i, noevents);
 	}
+
+	cout << "Finished in "<< stopwatch.RealTime() << " seconds" << endl;
+	cout << "Finished in "<< stopwatch.RealTime()/60. << " minutes" << endl;
 }

analysis/rootlogon.C

@@ -4,4 +4,6 @@
 //	gInterpreter->AddIncludePath("/home/vratik/workspace/AculUtils/TELoss");
 	gSystem->Load("../libUserAnalysis.so");
 //	gSystem->Load("../../AculUtils/libTELoss.so");
+
+	gStyle->SetOptStat(1);
 }

analysis/searchOfThresholds.cxx

+#include "TSystem.h"
+#include "TFile.h"
+#include "TTree.h"
+#include "TChain.h"
+#include "TCanvas.h"
+#include "TBox.h"
+#include "TCut.h"
+#include "TCutG.h"
+#include "TStopwatch.h"
+
+using std::cout;
+using std::endl;
+
+void searchOfThresholds() {
+
+	Long64_t drawEntries = 10000000;
+
+//	TFile *fr = new TFile("~/data/exp1804/h5_11_00.root");
+	TChain *tr = new TChain("AnalysisxTree");
+	tr->Add("~/data/exp1804/h5_14_0?.root");
+//	tr->Add("~/data/exp1804/be10_0?_?0.root");
+
+	TChain *trCal = new TChain("cal");
+//	trCal->Add("~/data/exp1804/h5_11_0?_calib.root");
+//	trCal->Add("~/data/exp1804/h5_14_0?_calib.root");
+//	trCal->Add("~/data/exp1804/be10_0?_?0_calib.root");
+
+
+
+	//calibration file for thin detector
+	TFile *fcal = new TFile("~/data/exp1804/calib/si_20_03.root");
+	TTree *tAlphaCal = (TTree*)fcal->Get("AnalysisxTree");
+
+	//calibration file for 1mm left detector
+	TFile *fcal1000 = new TFile("~/data/exp1804/calib/si_1000_LR_02.root");
+	TTree *tAlphaCal1000 = (TTree*)fcal1000->Get("AnalysisxTree");
+
+	tr->GetListOfFiles()->Print();
+	trCal->GetListOfFiles()->Print();
+
+	cout << tr->GetEntries() << " events." << endl;
+	cout << trCal->GetEntries() << " calibrated events." << endl;
+
+	TString varName;
+	TString condition;
+
+	Int_t firstThinStrip = 5;
+
+
+	TStopwatch stopwatch;
+	stopwatch.Start();
+
+
+//return;
+
+
+/////////////////////////////////////////////////////////////////////
+//	c0_1: thresholds search for SQ20 in channels:
+
+
+	/*TCanvas *c0_1 = new TCanvas("c0_1", "amplitudes in ADC channel unit", 1600, 1200);
+	c0_1->Divide(4,4);
+
+	firstThinStrip = 0;
+
+	//according to calibration data
+//	Int_t thinChThresholds[16] = {
+//		110, 112, 113, 115,
+//		112, 118, 113, 110,
+//		104, 115, 120, 113,
+//		115, 119, 113, 102
+//	};
+
+	//according to experimental data
+	Int_t thinChThresholds[16] = {
+		113, 121, 125, 125,
+		125, 132, 130, 115,
+		118, 128, 132, 130,
+		130, 130, 130, 112
+	};
+
+
+	cout << tr->GetEntries() << " events." << endl;
+	cout << trCal->GetEntries() << " calibrated events." << endl;
+	drawEntries = 10000000;
+	cout << drawEntries << " entries in input chain will be processed." << endl;
+
+	for (Int_t i = 0; i < 16; i++) {
+		c0_1->cd(i+1);
+		gPad->SetLogy();
+
+		varName.Form("SQ20[%d]", i+firstThinStrip);
+		condition.Form("SQ20[%d]>0 && SQ20[%d]<1200", i+firstThinStrip, i+firstThinStrip);
+
+		tAlphaCal->SetLineColor(kBlack);
+		tAlphaCal->Draw(varName, condition, "", drawEntries);
+		c0_1->Update();
+
+		condition.Form("SQ20[%d]>%d && SQ20[%d]<1200", i+firstThinStrip, thinChThresholds[i], i+firstThinStrip);
+		cout << condition << endl;
+
+		tAlphaCal->SetLineColor(kRed);
+		tAlphaCal->Draw(varName, condition, "same", drawEntries);
+		c0_1->Update();
+
+
+	}//*/
+
+//	return;
+
+
+/////////////////////////////////////////////////////////////////////
+//	c0_1_check: check of thresholds for SQ20 in channels
+//	on real experimental data:
+
+
+	/*TCanvas *c0_1_check = new TCanvas("c0_1_check", "amplitudes in ADC channel unit", 1600, 1200);
+	c0_1_check->Divide(4,4);
+
+	firstThinStrip = 0;
+
+	cout << tr->GetEntries() << " events." << endl;
+	cout << trCal->GetEntries() << " calibrated events." << endl;
+	drawEntries = 1000000;
+	cout << drawEntries << " entries in input chain will be processed." << endl;
+
+	for (Int_t i = 0; i < 16; i++) {
+		c0_1_check->cd(i+1);
+		gPad->SetLogy();
+
+		varName.Form("SQ20[%d]", i+firstThinStrip);
+		condition.Form("SQ20[%d]>0 && SQ20[%d]<200", i+firstThinStrip, i+firstThinStrip);
+
+		tr->SetLineColor(kBlack);
+		tr->Draw(varName, condition, "", drawEntries);
+		c0_1_check->Update();
+
+		condition.Form("SQ20[%d]>%d && SQ20[%d]<200", i+firstThinStrip, thinChThresholds[i], i+firstThinStrip);
+		cout << condition << endl;
+
+		tr->SetLineColor(kRed);
+		tr->Draw(varName, condition, "same", drawEntries);
+		c0_1_check->Update();
+
+
+	}//*/
+
+//	return;
+
+
+/////////////////////////////////////////////////////////////////////
+//	c0_2: thresholds search for SQX_L in channels:
+
+
+	/*TCanvas *c0_2 = new TCanvas("c0_2", "amplitudes in ADC channel unit (0-15)", 1600, 1200);
+	c0_2->Divide(4,4);
+	TCanvas *c0_2_2 = new TCanvas("c0_2_2", "amplitudes in ADC channel unit (15-32)", 1600, 1200);
+	c0_2_2->Divide(4,4);
+
+	firstThinStrip = 0;
+
+	Int_t thickChThresholds[32] = {
+		60, 45, 45, 45,
+		45, 4000, 4000, 45,
+		45, 45, 45, 45,
+		45, 45, 45, 45,
+		/////
+		45, 45, 45, 45,
+		45, 45, 45, 45,
+		45, 45, 45, 45,
+		45, 45, 45, 60
+	};
+
+	cout << tr->GetEntries() << " events." << endl;
+	cout << trCal->GetEntries() << " calibrated events." << endl;
+	drawEntries = 10000000;
+	cout << drawEntries << " entries in input chain will be processed." << endl;
+
+	for (Int_t i = 0; i < 32; i++) {
+		if (i<16) c0_2->cd(i+1);
+		if (i>15) c0_2_2->cd(i-16+1);
+		gPad->SetLogy();
+
+		varName.Form("SQX_L[%d]", i+firstThinStrip);
+		condition.Form("SQX_L[%d]>200 && SQX_L[%d]<2000", i+firstThinStrip, i+firstThinStrip);
+
+		tAlphaCal1000->SetLineColor(kBlack);
+		tAlphaCal1000->Draw(varName, condition, "", drawEntries);
+		c0_2->Update();
+		c0_2_2->Update();
+
+		condition.Form("SQX_L[%d]>%d && SQX_L[%d]<2000", i+firstThinStrip, thickChThresholds[i], i+firstThinStrip);
+		cout << condition << endl;
+
+		tAlphaCal1000->SetLineColor(kRed);
+		tAlphaCal1000->Draw(varName, condition, "same", drawEntries);
+		c0_2->Update();
+		c0_2_2->Update();
+
+
+	}//*/
+
+//	return;
+
+/////////////////////////////////////////////////////////////////////
+//	c0_2_check: check of thresholds for SQX_L in channels
+//	on real experimental data:
+
+
+	/*TCanvas *c0_2_check = new TCanvas("c0_2_check", "amplitudes in ADC channel unit (0-15); exp", 1600, 1200);
+	c0_2_check->Divide(4,4);
+	TCanvas *c0_2_2_check = new TCanvas("c0_2_2_check", "amplitudes in ADC channel unit (15-32); exp", 1600, 1200);
+	c0_2_2_check->Divide(4,4);
+
+	firstThinStrip = 0;
+
+	cout << tr->GetEntries() << " events." << endl;
+	cout << trCal->GetEntries() << " calibrated events." << endl;
+	drawEntries = 1000000;
+	cout << drawEntries << " entries in input chain will be processed." << endl;
+
+	for (Int_t i = 0; i < 32; i++) {
+		if (i<16) c0_2_check->cd(i+1);
+		if (i>15) c0_2_2_check->cd(i-16+1);
+		gPad->SetLogy();
+
+		varName.Form("SQX_L[%d]", i+firstThinStrip);
+		condition.Form("SQX_L[%d]>0 && SQX_L[%d]<400", i+firstThinStrip, i+firstThinStrip);
+
+		tr->SetLineColor(kBlack);
+		tr->Draw(varName, condition, "", drawEntries);
+		c0_2_check->Update();
+		c0_2_2_check->Update();
+
+		condition.Form("SQX_L[%d]>%d && SQX_L[%d]<400", i+firstThinStrip, thickChThresholds[i], i+firstThinStrip);
+		cout << condition << endl;
+
+		tr->SetLineColor(kRed);
+		tr->Draw(varName, condition, "same", drawEntries);
+		c0_2_check->Update();
+		c0_2_2_check->Update();
+
+
+	}//*/
+
+//	return;
+
+
+/////////////////////////////////////////////////////////////////////
+//	c0_3: thresholds search for SQY_L in channels:
+
+
+	/*TCanvas *c0_3 = new TCanvas("c0_3", "amplitudes in ADC channel unit (0-15); Y-strips", 1600, 1200);
+	c0_3->Divide(4,4);
+
+	firstThinStrip = 0;
+
+	//according to calibration data
+//	Int_t thickYChThresholds[16] = {
+//		4000, 130, 80, 110,
+//		90, 120, 80, 110,
+//		80, 110, 80, 110,
+//		80, 110, 90, 115,
+//	};
+
+	//according to experimental data
+	Int_t thickYChThresholds[16] = {
+		4000, 130, 82, 110,
+		87, 112, 82, 108,
+		82, 112, 82, 108,
+		80, 110, 90, 105,
+	};
+
+
+	cout << tr->GetEntries() << " events." << endl;
+	cout << trCal->GetEntries() << " calibrated events." << endl;
+	drawEntries = 10000000;
+	cout << drawEntries << " entries in input chain will be processed." << endl;
+
+	for (Int_t i = 0; i < 16; i++) {
+		c0_3->cd(i+1);
+
+		gPad->SetLogy();
+
+		varName.Form("SQY_L[%d]", i+firstThinStrip);
+		condition.Form("SQY_L[%d]>0 && SQY_L[%d]<2000", i+firstThinStrip, i+firstThinStrip);
+
+		tAlphaCal1000->SetLineColor(kBlack);
+		tAlphaCal1000->Draw(varName, condition, "", drawEntries);
+		c0_3->Update();
+
+		condition.Form("SQY_L[%d]>%d && SQY_L[%d]<2000", i+firstThinStrip, thickYChThresholds[i], i+firstThinStrip);
+		cout << condition << endl;
+
+		tAlphaCal1000->SetLineColor(kRed);
+		tAlphaCal1000->Draw(varName, condition, "same", drawEntries);
+		c0_3->Update();
+
+
+	}//*/
+
+//	return;
+
+/////////////////////////////////////////////////////////////////////
+//	c0_3_check: check of thresholds for SQY_L in channels
+//	on real experimental data:
+
+
+	/*TCanvas *c0_3_check = new TCanvas("c0_3_check", "amplitudes in ADC channel unit (0-15); Y-strips; exp", 1600, 1200);
+	c0_3_check->Divide(4,4);
+
+	firstThinStrip = 0;
+
+	cout << tr->GetEntries() << " events." << endl;
+	cout << trCal->GetEntries() << " calibrated events." << endl;
+	drawEntries = 1000000;
+	cout << drawEntries << " entries in input chain will be processed." << endl;
+
+	for (Int_t i = 0; i < 16; i++) {
+		c0_3_check->cd(i+1);
+
+		gPad->SetLogy();
+
+		varName.Form("SQY_L[%d]", i+firstThinStrip);
+		condition.Form("SQY_L[%d]>0 && SQY_L[%d]<200", i+firstThinStrip, i+firstThinStrip);
+
+		tr->SetLineColor(kBlack);
+		tr->Draw(varName, condition, "", drawEntries);
+		c0_3_check->Update();
+
+		condition.Form("SQY_L[%d]>%d && SQY_L[%d]<200", i+firstThinStrip, thickYChThresholds[i], i+firstThinStrip);
+		cout << condition << endl;
+
+		tr->SetLineColor(kRed);
+		tr->Draw(varName, condition, "same", drawEntries);
+		c0_3_check->Update();
+
+	}//*/
+
+
+
+/////////////////////////////////////////////////////////////////////
+//	c0_3_check: check of thresholds for SQY_L in channels
+//	on real experimental data:
+
+
+	TCanvas *c2_check = new TCanvas("c2_check", "amplitudes in ADC channel unit; CsI_left; exp", 1200, 1100);
+	c2_check->Divide(4,4);
+
+	firstThinStrip = 0;
+
+	cout << tr->GetEntries() << " events." << endl;
+	cout << trCal->GetEntries() << " calibrated events." << endl;
+	drawEntries = 1000000;
+	cout << drawEntries << " entries in input chain will be processed." << endl;
+
+	for (Int_t i = 0; i < 16; i++) {
+		c2_check->cd(i+1);
+
+//		gPad->SetLogy();
+
+		varName.Form("CsI_L[%d]", i+firstThinStrip);
+		condition.Form("CsI_L[%d]>0 && CsI_L[%d]<2000", i+firstThinStrip, i+firstThinStrip);
+
+		tr->SetLineColor(kBlack);
+		tr->Draw(varName, condition, "", drawEntries);
+		c2_check->Update();
+
+		condition.Form("CsI_L[%d]>%d && CsI_L[%d]<2000", i+firstThinStrip, 180, i+firstThinStrip);
+		cout << condition << endl;
+
+		tr->SetLineColor(kRed);
+		tr->Draw(varName, condition, "same", drawEntries);
+		c2_check->Update();
+
+	}//*/
+
+
+	return;
+
+
+}

analysis/searchTimeConditions.cxx

+#include "TSystem.h"
+#include "TFile.h"
+#include "TTree.h"
+#include "TChain.h"
+#include "TCanvas.h"
+#include "TBox.h"
+#include "TCut.h"
+#include "TCutG.h"
+#include "TStopwatch.h"
+
+using std::cout;
+using std::endl;
+
+void searchTimeConditions() {
+
+	Long64_t drawEntries = 10000000;
+
+//	TFile *fr = new TFile("~/data/exp1804/h5_11_00.root");
+	TChain *tr = new TChain("AnalysisxTree");
+	tr->Add("~/data/exp1804/h5_14_0?.root");
+//	tr->Add("~/data/exp1804/be10_0?_?0.root");
+
+	tr->SetMarkerStyle(20);
+	tr->SetMarkerSize(0.4);
+
+
+	TChain *trCal = new TChain("cal");
+//	trCal->Add("~/data/exp1804/h5_11_0?_calib.root");
+	trCal->Add("~/data/exp1804/h5_14_0?_calib.root");
+//	trCal->Add("~/data/exp1804/be10_0?_?0_calib.root");
+
+	tr->AddFriend(trCal);
+
+	TFile *fcal = new TFile("~/data/exp1804/calib/si_20_03_calib.root");
+	TTree *tAlphaCal = (TTree*)fcal->Get("cal");
+	tAlphaCal->SetMarkerColor(kRed);
+	tAlphaCal->SetMarkerStyle(20);
+	tAlphaCal->SetMarkerSize(.6);
+
+//	TFile *fc = new TFile("cutsIdentification.root", "READ");
+//	TCutG *cUS = (TCutG*)fc->Get("c7UpperShadow");
+//	TCutG *cLS = (TCutG*)fc->Get("c7LowerShadow");
+//	cLS->SetLineColor(kMagenta);
+//	TCutG *cMA = (TCutG*)fc->Get("c7MainAlpha");
+//	TCutG *cA = (TCutG*)fc->Get("c7All");
+//	TCutG *cBL = (TCutG*)fc->Get("cBeamLeft");
+//	TCutG *cBR = (TCutG*)fc->Get("cBeamRight");
+
+
+//	const Int_t drawEntries = 6000000;
+//	const Int_t drawEntries = tr->GetEntries();
+
+	tr->GetListOfFiles()->Print();
+	trCal->GetListOfFiles()->Print();
+
+	cout << tr->GetEntries() << " events." << endl;
+	cout << trCal->GetEntries() << " calibrated events." << endl;
+
+	TString varName;
+	TString condition;
+
+//	cout << drawEntries << " entries in input chain will be processed." << endl;
+
+	TString timeSiLeftCondition;
+	timeSiLeftCondition.Form("(SQLXtime[0]>0");
+	for (Int_t j = 1; j<32; j++) {
+		condition.Form("|| SQLXtime[%d]>0", j);
+		timeSiLeftCondition.Append(condition);
+	}
+	timeSiLeftCondition.Append(")");
+
+	Int_t firstThinStrip = 5;
+
+	Int_t conTimeMult = 1;
+
+	TStopwatch stopwatch;
+	stopwatch.Start();
+
+/////////////////////////////////////////////////////////////////////
+//	c1
+
+
+	TCanvas *c1 = new TCanvas("c1", "time in thin detector", 1600, 1200);
+	c1->Divide(4,4);
+
+	firstThinStrip = 0;
+
+	Int_t minTimeLthin[16] {
+		410, 410, 420, 415,
+		420, 410, 415, 415,
+		420, 420, 420, 420,
+		410, 425, 0, 0
+	};
+
+	Int_t maxTimeLthin[16] {
+		520, 525, 520, 530,
+		520, 530, 525, 525,
+		530, 530, 530, 520,
+		510, 525, 0, 0
+	};
+
+	cout << tr->GetEntries() << " events." << endl;
+	cout << trCal->GetEntries() << " calibrated events." << endl;
+	drawEntries = 10000000;
+	cout << drawEntries << " entries in input chain will be processed." << endl;
+
+	for (Int_t i = 0; i < 15; i++) {
+		c1->cd(i+1);
+
+		varName.Form("SQ20E[%d]:tSQ20[%d]*0.3", i+firstThinStrip, i+firstThinStrip);
+		condition.Form("tSQ20[%d]*0.3>300 && tSQ20[%d]*0.3<700 && SQ20E[%d]>0 && SQ20E[%d]<4",
+				i+firstThinStrip, i+firstThinStrip, i+firstThinStrip, i+firstThinStrip);
+
+		tr->SetMarkerColor(kBlack);
+		tr->Draw(varName, condition, "col", drawEntries);
+		c1->Update();
+
+		condition.Form("tSQ20[%d]>0 && SQ20E[%d]>0 && SQ20E[%d]<40 && SQ20timeMult==%d",
+				i+firstThinStrip, i+firstThinStrip, i+firstThinStrip, conTimeMult);
+		tr->SetMarkerColor(kRed);
+		tr->Draw(varName, condition, "same", drawEntries);
+		c1->Update();
+
+	}//*/
+
+//	return;
+
+/////////////////////////////////////////////////////////////////////
+//	c2
+
+	TCanvas *c2 = new TCanvas("c2", "time issues in 32 X-strips: 0-15", 1600, 1200);
+	c2->Divide(4,4);
+
+	firstThinStrip = 0;
+
+	Double_t minTimeLX[32] {
+		0, 0, 325, 328,
+		330, 0, 0, 330,
+		326, 328, 326, 335,
+		326, 328, 328, 330,
+		//
+		420, 425, 420, 420,
+		420, 420, 420, 420,
+		430, 425, 420, 420,
+		415, 415, 425, 415
+	};
+
+	Double_t maxTimeLX[32] {
+		0, 0, 345, 345,
+		350, 0, 0, 350,
+		345, 345, 345, 350,
+		340, 342, 345, 345,
+		//
+		460, 460, 460, 460,
+		460, 460, 460, 460,
+		460, 460, 460, 460,
+		460, 460, 460, 453
+	};
+
+
+	cout << tr->GetEntries() << " events." << endl;
+	cout << trCal->GetEntries() << " calibrated events." << endl;
+	drawEntries = 1000000;
+	cout << drawEntries << " entries in input chain will be processed." << endl;
+
+	for (Int_t i = 0; i < 16; i++) {
+		c2->cd(i+1);
+
+		varName.Form("SQLXE[%d]:0.3*tSQX_L[%d]", i+firstThinStrip, i+firstThinStrip);
+		condition.Form("0.3*tSQX_L[%d]>310 && 0.3*tSQX_L[%d]<380 && SQLXE[%d]>0 && trigger==3",
+				i+firstThinStrip, i+firstThinStrip, i+firstThinStrip);
+
+		tr->SetMarkerColor(kRed);
+		tr->SetMarkerStyle(20);
+		tr->SetMarkerSize(.4);
+		tr->Draw(varName, condition, "col", drawEntries);
+
+		varName.Form("SQLXE[%d]:SQLXtime[%d]", i+firstThinStrip, i+firstThinStrip);
+		condition.Form("SQLXtime[%d]>0 && SQLXE[%d]>0 && SQLXtimeMult==%d",
+				i+firstThinStrip, i+firstThinStrip, conTimeMult);
+
+		varName.Form("SQLXEtimeFiltered[%d]:0.3*tSQX_L[%d]", i+firstThinStrip, i+firstThinStrip);
+		condition.Form("0.3*tSQX_L[%d]>200 && 0.3*tSQX_L[%d]<500 && SQLXEtimeFilteredSum>0 && SQLXtimeMult==%d",
+				i+firstThinStrip, i+firstThinStrip, conTimeMult);
+
+		cout << condition << endl;
+
+		tr->Draw(varName, condition, "same", drawEntries);
+
+
+		c2->Update();
+	}//for */
+
+
+
+//	return;
+
+
+/////////////////////////////////////////////////////////////////////
+//	c2_1
+
+	TCanvas *c2_1 = new TCanvas("c2_1", "time issues in 32 X-strips: 16-31", 1600, 1200);
+	c2_1->Divide(4,4);
+
+	firstThinStrip = 16;
+
+	cout << tr->GetEntries() << " events." << endl;
+	cout << trCal->GetEntries() << " calibrated events." << endl;
+	drawEntries = 1000000;
+	cout << drawEntries << " entries in input chain will be processed." << endl;
+
+	for (Int_t i = 0; i < 16; i++) {
+		c2_1->cd(i+1);
+
+		varName.Form("SQLXE[%d]:0.3*tSQX_L[%d]", i+firstThinStrip, i+firstThinStrip);
+		condition.Form("0.3*tSQX_L[%d]>390 && 0.3*tSQX_L[%d]<480 && SQLXE[%d]>0", i+firstThinStrip, i+firstThinStrip, i+firstThinStrip);
+
+		tr->Draw(varName, condition, "col", drawEntries);
+
+		varName.Form("SQLXE[%d]:SQLXtime[%d]", i+firstThinStrip, i+firstThinStrip);
+		condition.Form("SQLXtime[%d]>0 && SQLXE[%d]>0 && SQLXtimeMult==%d",
+				i+firstThinStrip, i+firstThinStrip, conTimeMult);
+		tr->Draw(varName, condition, "same", drawEntries);
+
+		c2_1->Update();
+	}//for */
+
+//	return;
+
+/////////////////////////////////////////////////////////////////////
+//	c3
+
+	TCanvas *c3 = new TCanvas("c3", "time issues in 16 Y-strips", 1600, 1200);
+	c3->Divide(4,4);
+
+	firstThinStrip = 0;
+
+	cout << tr->GetEntries() << " events." << endl;
+	cout << trCal->GetEntries() << " calibrated events." << endl;
+	drawEntries = 1000000;
+	cout << drawEntries << " entries in input chain will be processed." << endl;
+
+	Float_t timeCorr[16] = {0, -2, 0, 2,
+							-2.5, 5, -6, 4,
+							7, 3, 3, 4,
+							6, 4, -4, 13
+	};
+
+
+	for (Int_t i = 0; i < 16; i++) {
+		c3->cd(i+1);
+
+
+		varName.Form("SQLYE[%d]:0.3*tSQY_L[%d]+%f", i+firstThinStrip, i+firstThinStrip, timeCorr[i]);
+		condition.Form("0.3*tSQY_L[%d]>300 && 0.3*tSQY_L[%d]<360 && SQLYE[%d]>0", i+firstThinStrip, i+firstThinStrip, i+firstThinStrip);
+
+		tr->SetMarkerColor(kRed);
+		tr->SetMarkerStyle(20);
+		tr->SetMarkerSize(.4);
+		tr->Draw(varName, condition, "col", drawEntries);
+//		tr->Draw(varName, condition, "", drawEntries);
+
+		varName.Form("SQLYEtimeFiltered[%d]:SQLYtime[%d]", i+firstThinStrip, i+firstThinStrip);
+		condition.Form("0.3*tSQY_L[%d]>200 && 0.3*tSQY_L[%d]<500 && SQLYEtimeFilteredSum>0 && SQLYtimeMult==%d",
+				i+firstThinStrip, i+firstThinStrip, conTimeMult);
+
+		cout << condition << endl;
+
+		tr->Draw(varName, condition, "same", drawEntries);
+		c3->Update();
+
+		condition.Form("0.3*tSQY_L[%d]+%f>325 && 0.3*tSQY_L[%d]+%f<333 && SQLYE[%d]>0",
+				i+firstThinStrip, timeCorr[i], i+firstThinStrip, timeCorr[i], i+firstThinStrip);
+
+		tr->SetMarkerColor(kRed);
+		tr->Draw(varName, condition, "same", drawEntries);
+
+		//		tAlphaCal->Draw(varName, "", "same");
+
+		c3->Update();
+	}//for */
+
+//	return;
+
+/////////////////////////////////////////////////////////////////////
+//	c4
+
+	TCanvas *c4 = new TCanvas("c4", "Multiplicity", 1800, 900);
+	c4->Divide(3,2);
+
+	cout << tr->GetEntries() << " events." << endl;
+	cout << trCal->GetEntries() << " calibrated events." << endl;
+	drawEntries = 1000000;
+	cout << drawEntries << " entries in input chain will be processed." << endl;
+
+
+	c4->cd(1);
+	varName.Form("SQLXmult");
+	condition.Form("trigger==3 && SQLXmult>0 && SQLXmult<8");
+	tr->Draw(varName, condition, "", drawEntries);
+	c4->Update();
+
+	c4->cd(4);
+	varName.Form("SQLXtimeMult");
+	condition.Form("trigger==3 && SQLXtimeMult>0 && SQLXtimeMult<8");
+	tr->Draw(varName, condition, "", drawEntries);
+	c4->Update();
+
+	c4->cd(2);
+	varName.Form("SQLYmult");
+	condition.Form("trigger==3 && SQLYmult>0 && SQLYmult<8");
+	tr->Draw(varName, condition, "", drawEntries);
+	c4->Update();
+
+	c4->cd(5);
+	varName.Form("SQLYtimeMult");
+	condition.Form("trigger==3 && SQLYtimeMult>0 && SQLYtimeMult<18");
+	tr->Draw(varName, condition, "", drawEntries);
+	c4->Update();
+
+	c4->cd(3);
+	varName.Form("SQ20EcorrMult");
+	condition.Form("trigger==3 && SQ20EcorrMult>0");
+	tr->Draw(varName, condition, "", drawEntries);
+	c4->Update();
+
+	c4->cd(6);
+	varName.Form("SQ20timeMult");
+	condition.Form("trigger==3 && SQ20timeMult>0");
+	tr->Draw(varName, condition, "", drawEntries);
+	c4->Update();
+
+	//*/
+
+	return;
+
+}

analysis/showBananas2.cxx

@@ -76,7 +76,7 @@ void showBananas2() {
 //	c1
 
 
-	/*TCanvas *c1 = new TCanvas("c1", "time in thin detector", 1600, 1200);
+	TCanvas *c1 = new TCanvas("c1", "time in thin detector", 1600, 1200);
 	c1->Divide(4,4);
 
 	firstThinStrip = 0;
@@ -104,7 +104,7 @@ void showBananas2() {
 /////////////////////////////////////////////////////////////////////
 //	c2
 
-	/*TCanvas *c2 = new TCanvas("c2", "time issues in 32 X-strips: 0-15", 1600, 1200);
+	TCanvas *c2 = new TCanvas("c2", "time issues in 32 X-strips: 0-15", 1600, 1200);
 	c2->Divide(4,4);
 
 	firstThinStrip = 0;
@@ -147,7 +147,7 @@ void showBananas2() {
 
 
 
-//	return;
+	return;
 
 
 /////////////////////////////////////////////////////////////////////

analysis/showBananasHe.cxx

+#include "TSystem.h"
+#include "TFile.h"
+#include "TTree.h"
+#include "TChain.h"
+#include "TCanvas.h"
+#include "TBox.h"
+#include "TCut.h"
+#include "TCutG.h"
+#include "TStopwatch.h"
+#include "TStyle.h"
+
+using std::cout;
+using std::endl;
+
+void showBananasHe(Long64_t drawEntries = 100000, const TString beam = "he", TString opt = "") {
+
+//	Long64_t drawEntries = 10000000;
+
+	TChain *trCal = new TChain("cal");
+	TChain *trCalBe = new TChain("cal");
+	if (beam.EqualTo("he")) trCal->Add("~/data/exp1804/h5_14_0?_calib.root");
+	if (beam.EqualTo("be")) trCal->Add("~/data/exp1804/be10_03_?0_calib.root");
+	if (beam.Contains("he") && beam.Contains("be")) {
+		trCalBe->Add("~/data/exp1804/be10_03_?0_calib.root");
+		trCal->Add("~/data/exp1804/h5_14_0?_calib.root");
+	}
+
+	trCal->GetListOfFiles()->Print();
+
+	cout << trCal->GetEntries() << " calibrated events." << endl;
+
+	TString varName;
+	TString condition;
+
+	Int_t firstThinStrip = 5;
+	Double_t max1mmDeposit = 31.;
+//	Double_t max1mmDeposit = 65.;
+
+//	Double_t max20deposit = 5.;
+	Double_t max20deposit = 4.;
+
+	Long64_t drawnEventsCorr2nd[16];
+
+
+	TStopwatch stopwatch;
+	stopwatch.Start();
+
+//return;
+
+
+/////////////////////////////////////////////////////////////////////
+//	amplitudes and their cuts for left telescope
+//	may be found in searchOfThresholds.cxx
+//
+//	time conditions for better choice of events
+//	may be found in searchTimeConditions.cxx
+//
+/////////////////////////////////////////////////////////////////////
+
+
+/////////////////////////////////////////////////////////////////////
+//	c1
+
+	/*TCanvas *c1 = new TCanvas("c1", "1D energy deposits in left Si detectors: different corrections", 1400, 700);
+	c1->Divide(3,2);
+
+
+
+	c1->cd(1);
+	trCal->Draw("SQLXEsum", "SQLXEsum>0 && trigger==3", "", 1000000);
+	c1->Update();
+
+	c1->cd(2);
+	trCal->Draw("SQLXEtimeFilteredSum", "SQLXEtimeFilteredSum>0 && trigger==3", "", 1000000);
+	c1->Update();
+
+	c1->cd(4);
+	trCal->Draw("SQ20E", "SQ20E>0 && trigger==3", "", 1000000);
+	c1->Update();
+
+	c1->cd(5);
+	trCal->Draw("SQ20Ecorr", "SQ20Ecorr>0 && trigger==3", "", 1000000);
+	c1->Update();
+
+	c1->cd(6);
+	trCal->Draw("SQ20EcorrHit", "SQ20EcorrHit>0 && trigger==3", "", 1000000);
+	c1->Update();
+
+	//*/
+
+
+//	return;
+
+/////////////////////////////////////////////////////////////////////
+//	c2
+
+	/*TCanvas *c2 = new TCanvas("c5", "Hits of reaction products in left telescope", 1400, 900);
+	c2->Divide(3,3);
+
+//	cout << tr->GetEntries() << " events." << endl;
+	cout << trCal->GetEntries() << " calibrated events." << endl;
+//	drawEntries = 10000000;
+	cout << drawEntries << " entries in input chain will be processed." << endl;
+	c2->Update();
+
+	c2->cd(1);
+	varName.Form("yt:xt>>hTarget(32, -40., 40., 32, -40., 40.)");
+	//	condition.Form("yThin>-80. && xThin>-80.");
+	condition.Form("trigger==3");
+	trCal->Draw(varName, condition, "col", drawEntries);
+	c2->Update();
+
+	c2->cd(2);
+	varName.Form("yThin:xThin>>hThin(32, -40., 40., 16, -40., 40.)");
+	condition.Form("trigger==3 && yThin>-80. && xThin>-80.");
+	//	condition.Form("yThin>-80. && xThin>-80.");
+	trCal->Draw(varName, condition, "col", drawEntries);
+	c2->Update();
+
+
+	c2->cd(3);
+	varName.Form("y1mm:x1mm>>h(32, -30., 30., 16, -30., 30.)");
+	condition.Form("trigger==3 && y1mm>-80. && x1mm>-80.");
+	trCal->Draw(varName, condition, "col", drawEntries);
+
+	TBox *boxThinDet = new TBox(25., -25., -25., 25.);
+	boxThinDet->SetLineColor(kRed);
+	boxThinDet->SetLineWidth(3);
+	boxThinDet->SetFillStyle(0);
+	boxThinDet->Draw("same");
+
+	Double_t thinXoffset = 1.;
+	Double_t thinYoffset = -1.8;
+
+	TBox *boxThinDetCorrPosition
+		= new TBox(25.+thinXoffset, -25.+thinYoffset, -25.+thinXoffset, 25.+thinYoffset);
+	boxThinDetCorrPosition->SetLineColor(kGreen);
+	boxThinDetCorrPosition->SetLineWidth(3);
+	boxThinDetCorrPosition->SetFillStyle(0);
+	boxThinDetCorrPosition->Draw("same");
+	c2->Update();
+
+	c2->cd(4);
+	varName.Form("xThin");
+	condition.Form("trigger==3 && yThin>-80. && xThin>10. && xThin<20.");
+	trCal->Draw(varName, condition, "col", drawEntries);
+
+	TLine *lx1 = new TLine(16.67, 0., 16.67, 100000.);
+	lx1->SetLineColor(kRed);
+	lx1->SetLineWidth(3);
+	lx1->Draw("same");
+	TLine *lx2 = new TLine(18.542, 0., 18.542, 100000.);
+	lx2->SetLineColor(kRed);
+	lx2->SetLineWidth(3);
+	lx2->Draw("same");
+	c2->Update();
+
+	c2->cd(5);
+	varName.Form("yThin");
+	condition.Form("trigger==3 && yThin>-10. && yThin>10. && xThin>-80.");
+	trCal->Draw(varName, condition, "col", drawEntries);
+
+	TLine *ly1 = new TLine(12.869, 0., 12.869, 15000.);
+	ly1->SetLineColor(kRed);
+	ly1->SetLineWidth(3);
+	ly1->Draw("same");
+	TLine *ly2 = new TLine(16.305, 0., 16.305, 15000.);
+	ly2->SetLineColor(kRed);
+	ly2->SetLineWidth(3);
+	ly2->Draw("same");
+	c2->Update();
+
+	c2->cd(6);
+	varName.Form("y1mm:x1mm");
+	condition.Form("trigger==3 && y1mm>-80. && x1mm>-80.");
+	condition.Form("trigger==3 && y1mm>-80. && x1mm>-10. && x1mm<10.");
+	trCal->Draw(varName, condition, "", drawEntries);
+
+
+
+	c2->cd(7);
+	varName.Form("mapXbin");
+	condition.Form("trigger==3 && mapXbin>-1");
+	trCal->Draw(varName, condition, "col", drawEntries);
+
+	TLine *lxMap = new TLine(25.5, 0., 25.5, 745252.);
+	lxMap->SetLineColor(kRed);
+	lxMap->SetLineWidth(3);
+	lxMap->Draw("same");
+	c2->Update();
+
+	c2->cd(8);
+	varName.Form("mapYbin");
+	condition.Form("trigger==3 && mapYbin>-1");
+	trCal->Draw(varName, condition, "col", drawEntries);
+	TLine *lyMap = new TLine(7.5, 0., 7.5, 268000.);
+	lyMap->SetLineColor(kRed);
+	lyMap->SetLineWidth(3);
+	lyMap->Draw("same");
+	c2->Update();
+
+	c2->cd(9);
+	//	varName.Form("SQ20timeMult");
+	//	condition.Form("trigger==3 && SQ20timeMult>0");
+	trCal->Draw("angleLeft", "angleLeft<3", "", drawEntries);
+	c2->Update();
+
+	//*/
+
+//	return;
+
+/////////////////////////////////////////////////////////////////////
+//	c3
+
+
+	/*TCanvas *c3 = new TCanvas("c3", "Multiplicity - clusters", 1400, 900);
+	c3->Divide(3,3);
+
+//	cout << tr->GetEntries() << " events." << endl;
+	cout << trCal->GetEntries() << " calibrated events." << endl;
+
+	c3->cd(1);
+	trCal->Draw("x1MultC", "x1MultC<4", "", drawEntries);
+	c3->Update();
+
+	c3->cd(2);
+	trCal->Draw("y1MultC", "y1MultC<4", "", drawEntries);
+	c3->Update();
+
+	c3->cd(4);
+	trCal->Draw("x2MultC", "x2MultC<4", "", drawEntries);
+	c3->Update();
+
+	c3->cd(5);
+	trCal->Draw("y2MultC", "y2MultC<4", "", drawEntries);
+	c3->Update();
+	//*/
+
+//	return;
+
+
+/////////////////////////////////////////////////////////////////////
+//	c11
+
+	/*TCanvas *c11 = new TCanvas("c11", "dE-E uncorrected", 1600, 800);
+	c11->Divide(4,2);
+
+	firstThinStrip = 5;
+
+//	Int_t drawnEvents1[8];
+
+//	cout << tr->GetEntries() << " events." << endl;
+	cout << trCal->GetEntries() << " calibrated events." << endl;
+//	drawEntries = 1000000000;
+	cout << drawEntries << " entries in input chain will be processed." << endl;
+
+	c11->cd(1);
+	trCal->Draw("SQ20E[0]:SQ20E[0]+SQLXEsum>>h11(150, 0, 35, 150, 0, 5)", "", "col", 1000);
+
+	for (Int_t i = 0; i < 8; i++) {
+		c11->cd(i+1);
+
+		if (opt.Contains("logz")) gPad->SetLogz();
+
+		varName.Form("SQ20E[%d]:SQ20E[%d]+SQLXEsum", i+firstThinStrip, i+firstThinStrip);
+		condition.Form("SQ20E[%d]>0 && SQ20E[%d]<%f "
+//				"&& SQLXEsum>1.1 "
+				"&& SQLXEsum<%f "
+				//				"&& SQLXtimeMult==1"
+				"&& trigger==3",
+				i+firstThinStrip, i+firstThinStrip, max20deposit, max1mmDeposit);
+
+		trCal->Draw(varName, condition, "col", drawEntries);
+
+		c11->Update();
+
+		trCalBe->SetMarkerStyle(20);
+		trCalBe->SetMarkerColor(kRed);
+		trCalBe->SetMarkerSize(.4);
+		trCalBe->Draw(varName, condition, "same", drawEntries);
+
+		TH2F *h = (TH2F*)gPad->GetPrimitive("htemp");
+		h->SetStats();
+
+		c11->Update();
+	}//*/
+
+//	return;
+
+
+
+
+/////////////////////////////////////////////////////////////////////
+//	c12
+
+	/*TCanvas *c12 = new TCanvas("c12", "dE-E corrected (first approximation)", 1600, 800);
+	c12->Divide(4,2);
+
+//	Int_t drawnEvents2[8];
+
+	firstThinStrip = 5;
+
+	cout << tr->GetEntries() << " events." << endl;
+	cout << trCal->GetEntries() << " calibrated events." << endl;
+	drawEntries = 1000000000;
+	cout << drawEntries << " entries in input chain will be processed." << endl;
+
+	c12->cd(1);
+	trCal->Draw("SQ20Ecorr[0]:SQ20Ecorr[0]+SQLXEsum>>h12(100, 0, 35, 100, 0, 5)", "", "col", 1000);
+
+	for (Int_t i = 0; i < 8; i++) {
+		c12->cd(i+1);
+
+		if (opt.Contains("logz")) gPad->SetLogz();
+
+		varName.Form("SQ20Ecorr[%d]:SQ20Ecorr[%d]+SQLXEsum",
+				i+firstThinStrip, i+firstThinStrip);
+		condition.Form(""
+				"SQ20Ecorr[%d]>0 "
+				"&& SQ20Ecorr[%d]<5 "
+//				"&& SQLXEsum>1.1 "
+				"&& SQLXEsum<%f "
+				"&& SQLXtimeMult==1"
+				"&& trigger==3",
+				i+firstThinStrip,
+				i+firstThinStrip,
+				max1mmDeposit);
+
+		trCal->Draw(varName, condition, "col", drawEntries);
+		c12->Update();
+	}//*/
+
+//	return;
+
+/////////////////////////////////////////////////////////////////////
+//	c13
+
+	/*TCanvas *c13 = new TCanvas("c13", "dE-E corrected (second approximation)", 1600, 800);
+	c13->Divide(4,2);
+
+	firstThinStrip = 5;
+
+//	cout << tr->GetEntries() << " events." << endl;
+	cout << trCal->GetEntries() << " calibrated events." << endl;
+//	drawEntries = 1000000000;
+	cout << drawEntries << " entries in input chain will be processed." << endl;
+
+	c13->cd(1);
+//	trCal->Draw("SQ20Ecorr[0]:SQ20Ecorr[0]+SQLXEsum>>h13(150, 0, 35, 150, 0, 5)", "", "col", 1000);
+	trCal->Draw("SQ20Ecorr[0]:SQ20Ecorr[0]+SQLXEsum>>h13(150, 0, 45, 150, 0, 5)", "", "col", 1000);
+
+	for (Int_t i = 0; i < 8; i++) {
+		c13->cd(i+1);
+
+		if (opt.Contains("logz")) gPad->SetLogz();
+
+		drawnEventsCorr2nd[i] = 0;
+
+		varName.Form("SQ20EcorrHit[%d]:SQ20EcorrHit[%d]+SQLXEtimeFilteredSum", i+firstThinStrip, i+firstThinStrip);
+		condition.Form(""
+				"SQ20EcorrHit[%d]>0"
+				" && SQ20EcorrHit[%d]<%f "
+//				"&& SQLXEsum>1.1"
+				" && SQLXEsum<%f "
+				"&& SQLXtimeMult==1"
+				"&& SQLXEtimeFilteredSum>0.5"
+				"&& trigger==3",
+				i+firstThinStrip,
+				i+firstThinStrip,
+				max20deposit,
+				max1mmDeposit);
+
+
+		drawnEventsCorr2nd[i] = trCal->Draw(varName, condition, "col", drawEntries);
+
+		TH2F *h = (TH2F*)gPad->GetPrimitive("htemp");
+		h->SetStats();
+
+		c13->Update();
+
+	}//*/
+
+//	return;
+
+
+/////////////////////////////////////////////////////////////////////
+//	c14
+
+	/*TCanvas *c14 = new TCanvas("c14", "dE-E corrected (second approximation, clusters)", 1600, 800);
+	c14->Divide(4,2);
+
+	firstThinStrip = 5;
+
+//	cout << tr->GetEntries() << " events." << endl;
+	cout << trCal->GetEntries() << " calibrated events." << endl;
+//	drawEntries = 1000000000;
+	cout << drawEntries << " entries in input chain will be processed." << endl;
+
+	c14->cd(1);
+	trCal->Draw("SQ20Ecorr[0]:SQ20Ecorr[0]+SQLXEsum>>h13(150, 0, 35, 150, 0, 5)", "", "col", 1000);
+
+	for (Int_t i = 0; i < 8; i++) {
+		c14->cd(i+1);
+		if (opt.Contains("logz")) gPad->SetLogz();
+		drawnEventsCorr2nd[i] = 0;
+
+		varName.Form("SQ20EcorrHitC[%d]:SQ20EcorrHitC[%d]+SQLXEtimeFilteredSum", i+firstThinStrip, i+firstThinStrip);
+		condition.Form(""
+				"SQ20EcorrHitC[%d]>0"
+				" && SQ20EcorrHitC[%d]<%f "
+//				"&& SQLXEsum>1.1"
+				" && SQLXEsum<%f "
+				"&& SQLXtimeMult==1"
+				"&& SQLXEtimeFilteredSum>0.5"
+				"&& trigger==3",
+				i+firstThinStrip,
+				i+firstThinStrip,
+				max20deposit,
+				max1mmDeposit);
+
+
+		drawnEventsCorr2nd[i] = trCal->Draw(varName, condition, "col", drawEntries);
+		trCalBe->SetMarkerStyle(20);
+		trCalBe->SetMarkerColor(kRed);
+		trCalBe->SetMarkerSize(.4);
+		trCalBe->Draw(varName, condition, "same", drawEntries);
+
+		TH2F *h = (TH2F*)gPad->GetPrimitive("htemp");
+		h->SetStats();
+
+		c14->Update();
+
+	}//*/
+
+//	return;
+
+/////////////////////////////////////////////////////////////////////
+//	c15
+
+	/*TCanvas *c15 = new TCanvas("c15", "dE-E corrected (second approximation, clusters, cleaned)", 1600, 800);
+	c15->Divide(4,2);
+
+	firstThinStrip = 5;
+
+//	cout << tr->GetEntries() << " events." << endl;
+	cout << trCal->GetEntries() << " calibrated events." << endl;
+//	drawEntries = 1000000000;
+	cout << drawEntries << " entries in input chain will be processed." << endl;
+
+	c15->cd(1);
+	trCal->Draw("SQ20Ecorr[0]:SQ20Ecorr[0]+SQLXEsum>>h13(150, 0, 35, 150, 0, 5)", "", "col", 1000);
+
+	for (Int_t i = 0; i < 8; i++) {
+		c15->cd(i+1);
+		if (opt.Contains("logz")) gPad->SetLogz();
+		drawnEventsCorr2nd[i] = 0;
+
+		varName.Form("SQ20EcorrHitC[%d]:SQ20EcorrHitC[%d]+SQLXEtimeFilteredSum", i+firstThinStrip, i+firstThinStrip);
+		condition.Form(""
+				"SQ20EcorrHitC[%d]>0"
+				" && SQ20EcorrHitC[%d]<%f "
+//				"&& SQLXEsum>1.1"
+				" && SQLXEsum<%f "
+				"&& SQLXtimeMult==1"
+				"&& SQLXEtimeFilteredSum>0.5"
+				"&& SQ20EcorrHitCMult"
+				"&& CsI_L_veto==0"
+				"&& trigger==3",
+				i+firstThinStrip,
+				i+firstThinStrip,
+				max20deposit,
+				max1mmDeposit);
+
+
+		drawnEventsCorr2nd[i] = trCal->Draw(varName, condition, "col", drawEntries);
+		trCalBe->SetMarkerStyle(20);
+		trCalBe->SetMarkerColor(kRed);
+		trCalBe->SetMarkerSize(.4);
+		trCalBe->Draw(varName, condition, "same", drawEntries);
+
+		TH2F *h = (TH2F*)gPad->GetPrimitive("htemp");
+		h->SetStats();
+
+		c15->Update();
+
+	}//*/
+
+//	return;
+
+
+/////////////////////////////////////////////////////////////////////
+//	c16
+
+	TCanvas *c16 = new TCanvas("c16", "dE-E corrected (second approximation, clusters, cleaned, check for thin strip)", 1600, 800);
+	c16->Divide(4,2);
+
+	firstThinStrip = 5;
+
+//	cout << tr->GetEntries() << " events." << endl;
+	cout << trCal->GetEntries() << " calibrated events." << endl;
+//	drawEntries = 1000000000;
+	cout << drawEntries << " entries in input chain will be processed." << endl;
+
+	gStyle->SetOptStat();
+
+	c16->cd(1);
+	trCal->Draw("SQ20Ecorr[0]:SQ20Ecorr[0]+SQLXEsum>>h13(150, 0, 35, 150, 0, 5)", "", "col", 1000);
+
+	for (Int_t i = 0; i < 8; i++) {
+		c16->cd(i+1);
+		if (opt.Contains("logz")) gPad->SetLogz();
+		drawnEventsCorr2nd[i] = 0;
+
+		varName.Form("SQ20EcorrHitC[%d]:SQ20EcorrHitC[%d]+SQLXEtimeFilteredSum", i+firstThinStrip, i+firstThinStrip);
+		condition.Form(""
+				"SQ20EcorrHitC[%d]>0"
+				" && SQ20EcorrHitC[%d]<%f "
+//				"&& SQLXEsum>1.1"
+				" && SQLXEsum<%f "
+				"&& SQLXtimeMult==1"
+				"&& SQLXEtimeFilteredSum>0.5"
+				"&& SQ20EcorrHitCMult"
+				"&& CsI_L_veto==0"
+				"&& abs(expectedThinStrip-%d)<=1"
+				"&& trigger==3",
+				i+firstThinStrip,
+				i+firstThinStrip,
+				max20deposit,
+				max1mmDeposit,
+				i+firstThinStrip);
+
+
+		drawnEventsCorr2nd[i] = trCal->Draw(varName, condition, "col", drawEntries);
+		trCalBe->SetMarkerStyle(20);
+		trCalBe->SetMarkerColor(kRed);
+		trCalBe->SetMarkerSize(.4);
+		trCalBe->Draw(varName, condition, "same", drawEntries);
+
+		TH2F *h = (TH2F*)gPad->GetPrimitive("htemp");
+		h->SetStats();
+
+		c16->Update();
+
+	}//*/
+
+//	return;
+
+
+
+
+
+//	TCanvas *c = new TCanvas();
+//	trCal->Draw("expectedThinStrip", "", "", drawEntries);
+
+
+
+	for (Int_t i = 0; i<8; i++) {
+		cout << "Events after 2nd correction in strip No." << i << ": " << drawnEventsCorr2nd[i] << endl;
+	}
+
+	cout << "Finished in "<< stopwatch.RealTime() << " seconds" << endl;
+	cout << "Finished in "<< stopwatch.RealTime()/60. << " minutes" << endl;
+
+	return;
+
+}

analysis/showBeam.cxx

@@ -7,21 +7,17 @@
 using std::cout;
 using std::endl;
 
-void showBeam(const Long64_t drawEntries = 100000)
+void showBeam(const Long64_t drawEntries = 100000, const TString beam = "he")
 {
 
 	TChain *tr = new TChain("AnalysisxTree");
-//	tr->Add("~/data/exp1804/h5_14_0?.root");
-//	tr->Add("~/data/exp1804/output.root");
-	tr->Add("~/data/exp1804/h5_12.root");
+	if (beam.Contains("he")) tr->Add("~/data/exp1804/h5_14_0?.root");
+	if (beam.Contains("be")) tr->Add("~/data/exp1804/be10_03_?0.root");
 
 
 	TChain *trCal = new TChain("cal");
-//	trCal->Add("~/data/exp1804/h5_11_0?_calib.root");
-//	trCal->Add("~/data/exp1804/h5_14_0?_calib.root");
-//	trCal->Add("~/data/exp1804/h5_14_0?_calib.root");
-
-//	tr->AddFriend(trCal);
+	if (beam.Contains("he")) trCal->Add("~/data/exp1804/h5_14_0?_calib.root");
+	if (beam.Contains("be")) trCal->Add("~/data/exp1804/be10_03_?0_calib.root");
 
 //	const Int_t drawEntries = 6000000;
 //	const Int_t drawEntries = tr->GetEntries();
@@ -37,27 +33,17 @@ void showBeam(const Long64_t drawEntries = 100000)
 
 //	cout << drawEntries << " entries in input chain will be processed." << endl;
 
+/////////////////////////////////////////////////////////////////////
+//	c1
+
 
-//	TCanvas *c1 = new TCanvas("c1", "R: 0-7", 1600, 800);
-//	c1->Divide(4,2);
 	TCanvas *c1 = new TCanvas("beam", "beam info", 1200, 1200);
 	c1->Divide(2,2);
 
 	c1->cd(1);
-//	tr->Draw("trigger", "", "", drawEntries);
 
 	c1->cd(3);
 
-//	af5 = (NeEvent->F5[0]+NeEvent->F5[1]+NeEvent->F5[2]+NeEvent->F5[3]+4.*gRandom->Uniform())/4.
-//	tf5 = (NeEvent->tF5[0]+NeEvent->tF5[1]+NeEvent->tF5[2]+NeEvent->tF5[3]+4.*gRandom->Uniform())/4.;
-//	tf3 = (NeEvent->tF3[0]+NeEvent->tF3[1]+NeEvent->tF3[2]+NeEvent->tF3[3]+4.*gRandom->Uniform())/4.;
-//	ToF = (tf5 - tf3)*0.125+89.165;
-//	af5:TOF
-//	tr->Draw("(NeEvent->F5[0]+NeEvent->F5[1]+NeEvent->F5[2]+NeEvent->F5[3])/4.", "", "", drawEntries);
-//	tr->Draw("(F5[0]+F5[1]+F5[2]+F5[3])/4.", "", "", drawEntries);
-//	tr->Draw("(tF5[0]+tF5[1]+tF5[2]+tF5[3])/4.;", "", "", drawEntries);
-//	tr->Draw("( (tF5[0]+tF5[1]+tF5[2]+tF5[3]) - (tF3[0]+tF3[1]+tF3[2]+tF3[3]) )/4.*0.125+89.165;",
-//			"tF5[0]>0 && tF5[1]>0 && tF5[2]>0 && tF5[3]>0 && tF3[0]>0 && tF3[1]>0 && tF3[2]>0 && tF3[3]>0", "", drawEntries);
 	tr->Draw("(F5[0]+F5[1]+F5[2]+F5[3])/4.:( (tF5[0]+tF5[1]+tF5[2]+tF5[3]) - (tF3[0]+tF3[1]+tF3[2]+tF3[3]) )/4.*0.125+89.165;",
 			"tF5[0]>0 && tF5[1]>0 && tF5[2]>0 && tF5[3]>0 "
 			"&& tF3[0]>0 && tF3[1]>0 && tF3[2]>0 && tF3[3]>0 "
@@ -108,6 +94,9 @@ void showBeam(const Long64_t drawEntries = 100000)
 	h2->GetYaxis()->SetTitleSize(0.06);
 	h2->GetYaxis()->SetTitleOffset(0.75);
 
+/////////////////////////////////////////////////////////////////////
+//	c2
+
 
 	TCanvas *c2 = new TCanvas("c2", "Bugs in MWPC", 1200, 1200);
 	c2->Divide(2,2);
@@ -126,6 +115,128 @@ void showBeam(const Long64_t drawEntries = 100000)
 	c2->cd(4);
 	tr->Draw("y2[0]", "ny2==0", "", drawEntries);
 
+	c2->Update();
+
+/////////////////////////////////////////////////////////////////////
+//	c3
+
+
+	TCanvas *c3 = new TCanvas("c3", "Wire multiplicity", 1200, 1200);
+	c3->Divide(2,2);
+
+	c3->cd(1);
+//	tr->Draw("x1[0]", "x2[0]<1000 && nx2==1", "", drawEntries);
+	tr->Draw("nx1", "nx1<30000", "", drawEntries);
+
+	c3->cd(2);
+	tr->Draw("ny1", "ny1<30000", "", drawEntries);
+
+	c3->cd(3);
+	tr->Draw("nx2", "nx2<30000", "", drawEntries);
+
+	c3->cd(4);
+	tr->Draw("ny2", "ny2<30000", "", drawEntries);
+
+
+
+/////////////////////////////////////////////////////////////////////
+//	c4
+
+
+	TCanvas *c4 = new TCanvas("c4", "MWPC time", 1200, 1200);
+	c4->Divide(2,2);
+
+	c4->cd(1);
+//	tr->Draw("x1[0]", "x2[0]<1000 && nx2==1", "", drawEntries);
+	tr->Draw("tMWPC[0]*0.125", "nx1<30000", "", drawEntries);
+	tr->SetLineColor(kRed);
+	tr->Draw("tMWPC[0]*0.125", "nx1<30000 && trigger==1", "same", drawEntries);
+	tr->SetLineColor(kBlue);
+	tr->Draw("tMWPC[0]*0.125", "nx1<30000 && trigger==3", "same", drawEntries);
+
+
+	tr->SetLineColor(kBlack);
+	c4->cd(2);
+	tr->Draw("tMWPC[1]*0.125", "ny1<30000", "", drawEntries);
+	tr->Draw("tMWPC[1]*0.125-tF5[0]*0.125", "ny1<30000", "", drawEntries);
+
+	c4->cd(3);
+	tr->Draw("tMWPC[2]*0.125", "nx2<30000", "", drawEntries);
+
+	c4->cd(4);
+	tr->Draw("tMWPC[3]*0.125", "ny2<30000", "", drawEntries);
+
+	c4->Update();
+
+/////////////////////////////////////////////////////////////////////
+//	c5
+
+
+	TCanvas *c5 = new TCanvas("c5", "Multiplicity filtered for time", 1200, 1200);
+	c5->Divide(2,2);
+
+	c5->cd(1);
+	tr->Draw("nx1", "nx1<30000 && (tMWPC[0]*0.125-tF5[0]*0.125)>60 && (tMWPC[0]*0.125-tF5[0]*0.125)<77", "", drawEntries);
+
+	c5->cd(2);
+	tr->Draw("ny1", "ny1<30000 && (tMWPC[1]*0.125-tF5[0]*0.125)>60 && (tMWPC[1]*0.125-tF5[0]*0.125)<80", "", drawEntries);
+
+	c5->cd(3);
+	tr->Draw("nx2", "nx2<30000 && (tMWPC[2]*0.125-tF5[0]*0.125)>70 && (tMWPC[2]*0.125-tF5[0]*0.125)<90", "", drawEntries);
+
+	c5->cd(4);
+	tr->Draw("ny2", "ny2<30000 && (tMWPC[3]*0.125-tF5[0]*0.125)>60 && (tMWPC[3]*0.125-tF5[0]*0.125)<80", "", drawEntries);
+
+	c5->Update();
+
+/////////////////////////////////////////////////////////////////////
+//	c6
+
+
+	TCanvas *c6 = new TCanvas("c6", "Cluster multiplicity", 1500, 1000);
+	c6->Divide(3,2);
+
+	c6->cd(1);
+	trCal->Draw("x1MultC", "x1MultC>0", "", drawEntries);
+
+	c6->cd(2);
+	trCal->Draw("y1MultC", "y1MultC>0", "", drawEntries);
+
+	c6->cd(4);
+	trCal->Draw("x2MultC", "x2MultC>0", "", drawEntries);
+
+	c6->cd(5);
+	trCal->Draw("y2MultC", "y2MultC>0", "", drawEntries);
+
+	c6->cd(3);
+	tr->Draw("nx1", "ny2==1 && nx2==1 && ny1==1 && nx1==1", "", drawEntries);
+
+	c6->cd(6);
+	trCal->Draw("x1MultC", "x1MultC==1 && y1MultC==1 && x2MultC==1 && y2MultC==1", "", drawEntries);
+//	tr->Draw("x1", "nx1<1", "", drawEntries);
+
+	c6->Update();
+
+
+/////////////////////////////////////////////////////////////////////
+//	c1
+
+
+	TCanvas *c7 = new TCanvas("c7", "beam info - clusters", 1200, 1200);
+	c7->Divide(2,2);
+
+	c7->cd(1);
+	trCal->Draw("y1p:x1p", "y1p>-50. && x1p>-50.", "col", drawEntries);
+
+	c7->cd(3);
+	trCal->Draw("y2p:x2p", "y2p>-50. && x2p>-50.", "col", drawEntries);
+
+	c7->cd(2);
+	trCal->Draw("y1c:x1c", "y1c>-50. && x1c>-50.", "col", drawEntries);
+
+	c7->cd(4);
+	trCal->Draw("y2c:x2c", "y2c>-50. && x2c>-50.", "col", drawEntries);
+
 
 	return;