System of coordinates in MWPC's inverted

analysis/fillChain.cxx

@@ -24,7 +24,7 @@ Int_t GetClustersMWPC(unsigned short n, unsigned short *x)
 	return clusters;
 } //--------------------------------------------------------------------
 
-Double_t GetClusterPositionMWPC(unsigned short n, unsigned short *x, Float_t planeOffset=0.)
+Double_t GetClusterPositionMWPC(unsigned short n, unsigned short *x, Float_t wireStep, Float_t planeOffset=0.)
 {
 	Double_t position = -100.;
 
@@ -34,7 +34,9 @@ Double_t GetClusterPositionMWPC(unsigned short n, unsigned short *x, Float_t pla
 //	if(n == 0) clusters = 0.;
 	if(n > 0 && n<=32)
 	{
-		position = (x[0]+n/2.+0.5)*1.25-20. + planeOffset;
+
+//		x2p = (revent->x2[0] + 0.5 - 16)*kMWPCwireStepX2 + MWPC2_X_offset;
+		position = (x[0]+n/2. + 0.5 -16)*wireStep + planeOffset;
 //		cout << n << endl;
 	}
 	return position;
@@ -271,16 +273,24 @@ void fillTree(const TString beam = "he", Int_t nofile = 0, const Int_t noevents
 	// Parameters related to geometry
 	///////////////////////////////////////////////////
 
+	//all distances in mm
+
+	//MWPC
+	const Float_t kMWPCz1 = -815.;	//z coordinate of the center of MWPC1
+	const Float_t kMWPCz2 = -270.;	//z coordinate of the center of MWPC2
 
-	//MWPC //todo: convert to minus
-	const Float_t l12 = 546.;	//z coordinate of the center of MWPC1
-	const Float_t lt = 270.;	//z coordinate of the center of MWPC2
+	//step between two wires;
+	//sign mean onward (+) or rearward (-) to physical axis
+	const Float_t kMWPCwireStepX1 = -1.25;
+	const Float_t kMWPCwireStepY1 = 1.25;		//step between two wires
+	const Float_t kMWPCwireStepX2 = -1.25;		//step between two wires
+	const Float_t kMWPCwireStepY2 = 1.25;		//step between two wires
 
 	//offsets taken from S. Krupko
-	const Float_t MWPC1_X_offset = -1.0;
-	const Float_t MWPC1_Y_offset = -2.1375;
+	const Float_t MWPC1_X_offset = -1.19;
+	const Float_t MWPC1_Y_offset = -2.12;
 	const Float_t MWPC2_X_offset = 0.2;
-	const Float_t MWPC2_Y_offset = -1.125;
+	const Float_t MWPC2_Y_offset = -1.02;
 
 	//left telescope
 	const Float_t z1mm = 230.;
@@ -289,8 +299,8 @@ void fillTree(const TString beam = "he", Int_t nofile = 0, const Int_t noevents
 	const Float_t yThinOffset = -1.8;
 
 
-	//todo: add sign taking into account direction of numbering
-	const Int_t kSQL_X_strips = 32;
+	//TODO: add sign taking into account direction of numbering
+	const Int_t kSQL_X_strips = -32;
 	const Int_t kSQL_Y_strips = 16;
 	const Int_t kSQL_20_strips = 16;
 
@@ -393,6 +403,8 @@ void fillTree(const TString beam = "he", Int_t nofile = 0, const Int_t noevents
 
 		if (trigger != 3) continue;
 
+		//reset values to be written
+
 		SQ20Esum = 0.;
 		SQLXEsum = 0.;
 		SQLYEsum = 0.;
@@ -415,6 +427,25 @@ void fillTree(const TString beam = "he", Int_t nofile = 0, const Int_t noevents
 		CsI_L_veto = kFALSE;
 		expectedThinStrip = -1;
 
+		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;
+
 
 		///////////////////////////////////////////////
 		//beam diagnostics
@@ -444,83 +475,79 @@ void fillTree(const TString beam = "he", Int_t nofile = 0, const Int_t noevents
 			///////////////////////////////////////////////
 
 
-		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];
 
+		//check for MWPC timing
+		//TODO: parametrization of time thresholds
 		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;
 
+		if (flagMWPC==0) continue;
+
 		//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;
+//				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.;
+				x1p = (revent->x1[0] + 0.5 - 16)*kMWPCwireStepX1 + MWPC1_X_offset;
+				y1p = (revent->y1[0] + 0.5 - 16)*kMWPCwireStepY1 + MWPC1_Y_offset;
 			}
 
 			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;
+//				x2p = (revent->x2[0]+0.5)*1.25-20. + MWPC2_X_offset;
+//				y2p = (revent->y2[0]+0.5)*1.25-20. + MWPC2_Y_offset;
+
+				x2p = (revent->x2[0] + 0.5 - 16)*kMWPCwireStepX2 + MWPC2_X_offset;
+				y2p = (revent->y2[0] + 0.5 - 16)*kMWPCwireStepY2 + 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 = 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;
+
+				xt = x1p - (x2p -x1p)*kMWPCz1/(kMWPCz2-kMWPCz1);
+				yt = y1p - (y2p -y1p)*kMWPCz1/(kMWPCz2-kMWPCz1);
+
 			}
 		}//if wire multiplicity == 1
 
 //		cout << flagMWPC << endl;
 
-		//MWPC: work with cluster multiplicity equal to 2
+		//MWPC: work with cluster multiplicity equal to 1
 		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);
 
+			//set wire step as input parameter for the function
 			if (x1MultC==1 && y1MultC==1) {
-				x1c = GetClusterPositionMWPC(revent->nx1, revent->x1, MWPC1_X_offset);
-				y1c = GetClusterPositionMWPC(revent->ny1, revent->y1, MWPC1_Y_offset);
+				x1c = GetClusterPositionMWPC(revent->nx1, revent->x1, kMWPCwireStepX1, MWPC1_X_offset);
+				y1c = GetClusterPositionMWPC(revent->ny1, revent->y1, kMWPCwireStepY1, 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);
+				x2c = GetClusterPositionMWPC(revent->nx2, revent->x2, kMWPCwireStepX2, MWPC2_X_offset);
+				y2c = GetClusterPositionMWPC(revent->ny2, revent->y2, kMWPCwireStepY2, 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 = 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;
+
+				xt = x1c - (x2c -x1c)*kMWPCz1/(kMWPCz2-kMWPCz1);
+				yt = y1c - (y2c -y1c)*kMWPCz1/(kMWPCz2-kMWPCz1);
 			}
 		}
 
@@ -642,7 +669,9 @@ void fillTree(const TString beam = "he", Int_t nofile = 0, const Int_t noevents
 //		cout << endl;
 
 		///////////////////////////////////////////////
-		//correction for thickness inhomogeneity
+		// correction for thickness inhomogeneity
+		// when not taking into account hitting
+		// particle tracking
 		///////////////////////////////////////////////
 
 		Double_t currThickness;
@@ -653,35 +682,22 @@ void fillTree(const TString beam = "he", Int_t nofile = 0, const Int_t noevents
 				//			cout << "y bin: " << yi+1 << "\t\t";
 				if (SQLYE[yi]>0) {
 
-					for (Int_t xi = 0; xi < kSQL_X_strips; xi++) {
+					for (Int_t xi = 0; xi < TMath::Abs(kSQL_X_strips); xi++) {
 
 						if (SQLXE[xi] > 0) {
+							//TODO: check the mapping here at first
 							currThickness = hThickness->GetBinContent(xi+1, yi+1);
-//							if (currThickness<30.) {
-//								Int_t probableThinStrip = xi-2/2;
-//								SQ20Ecorr[0] =
-//							}
 
 							for (Int_t x20 = 0; x20 < kSQL_20_strips; x20++) {
 								if (SQ20E[x20] > 0/*SQ20timeMult==1*/ && currThickness<30.) {
-//									SQ20Ecorr[x20] = SQ20E[x20] + 1.;
-
-//									currThickness = kSQ20_norm_thickness + (currThickness - kSQ20_norm_thickness)*0.5;
 
 									SQ20Ecorr[x20] = (kSQ20_norm_thickness/currThickness)*SQ20E[x20];
 									SQ20EcorrMult++;
 
-//									cout << xi << "\t" << yi << "\t" << x20
-//											<< "\t" << currThickness
-//											<< "\t" << SQ20E[x20]  << "\t" << SQ20Ecorr[x20] << "\t" << endl;
 								}
 							}//for*/
 						}
 
-
-
-						//					cout << hThickness->GetBinContent(xi+1, yi+1) << "\t";
-						//					if (xi == kSQL_X_strips-1) cout << endl;
 					}
 				}
 			}
@@ -692,7 +708,7 @@ void fillTree(const TString beam = "he", Int_t nofile = 0, const Int_t noevents
 
 
 		///////////////////////////////////////////////
-		//position of hit in left telescope
+		// position of hit in left telescope
 		///////////////////////////////////////////////
 
 		//reset
@@ -707,14 +723,20 @@ void fillTree(const TString beam = "he", Int_t nofile = 0, const Int_t noevents
 
 		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) {
+					//todo: check the system of coordinates
+					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) {
+					//todo: check the system of coordinates
+					y1mm = -30. + (i+1./2.)*60./16.;
+				}
 			}
 		}
 
@@ -739,6 +761,7 @@ void fillTree(const TString beam = "he", Int_t nofile = 0, const Int_t noevents
 			yThin = vHit1mm.Y()+yt;
 
 			//calculation of corrected energy
+			//TODO: check the mapping here at second
 			currThickness = hThickness->GetBinContent(mapXbin+1, mapYbin+1);
 
 			for (Int_t x20 = 0; x20 < kSQL_20_strips; x20++) {

analysis/showBananasHe.cxx

@@ -546,6 +546,65 @@ void showBananasHe(Long64_t drawEntries = 100000, const TString beam = "he", TSt
 //	return;
 
 
+/////////////////////////////////////////////////////////////////////
+//	c17
+
+	TCanvas *c17 = new TCanvas("c17", "dE-E corrected (second approximation, clusters, cleaned, check for thin strip); target center", 1600, 800);
+	c17->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();
+
+	c17->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++) {
+		c17->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"
+				"&& abs(xt)<10. && abs(yt)<10."
+				"&& xThin<25. && yThin<24.",
+				i+firstThinStrip,
+				i+firstThinStrip,
+				max20deposit,
+				max1mmDeposit,
+				i+firstThinStrip);
+
+
+		drawnEventsCorr2nd[i] = trCal->Draw(varName, condition, "", drawEntries);
+		trCalBe->SetMarkerStyle(20);
+		trCalBe->SetMarkerColor(kRed);
+		trCalBe->SetMarkerSize(.4);
+		trCalBe->Draw(varName, condition, "same", drawEntries);
+
+		TH2F *h = (TH2F*)gPad->GetPrimitive("htemp");
+		h->SetStats();
+
+		c17->Update();
+
+	}//*/
+
+//	return;
+
 
 
 

analysis/showBeam.cxx

@@ -42,6 +42,27 @@ void showBeam(const Long64_t drawEntries = 100000, const TString beam = "he")
 
 	c1->cd(1);
 
+	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 "
+			"&& F5[0]>0 && F5[1]>0 && F5[2]>0 && F5[3]>0"
+			"&& ( (tF5[0]+tF5[1]+tF5[2]+tF5[3]) - (tF3[0]+tF3[1]+tF3[2]+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"
+			"&& trigger==1",
+			"col", drawEntries);
+	TH2F *hTOF1 = (TH2F*)gPad->GetPrimitive("htemp");
+//	hTOF->SetStats(0);
+	hTOF1->SetTitle("trigger 1");
+	hTOF1->SetXTitle("TOF (ns)");
+	hTOF1->SetYTitle("dE (a.u.)");
+	hTOF1->GetXaxis()->CenterTitle();
+	hTOF1->GetXaxis()->SetTitleSize(0.06);
+	hTOF1->GetXaxis()->SetTitleOffset(0.75);
+	hTOF1->GetYaxis()->CenterTitle();
+	hTOF1->GetYaxis()->SetTitleSize(0.06);
+	hTOF1->GetYaxis()->SetTitleOffset(0.75);
+
 	c1->cd(3);
 
 	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;",
@@ -51,19 +72,19 @@ void showBeam(const Long64_t drawEntries = 100000, const TString beam = "he")
 			"&& ( (tF5[0]+tF5[1]+tF5[2]+tF5[3]) - (tF3[0]+tF3[1]+tF3[2]+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"
-			"&& trigger==1",
+			"&& trigger==3",
 			"col", drawEntries);
-	TH2F *hTOF = (TH2F*)gPad->GetPrimitive("htemp");
-	hTOF->SetStats(0);
-	hTOF->SetTitle("");
-	hTOF->SetXTitle("TOF (ns)");
-	hTOF->SetYTitle("dE (a.u.)");
-	hTOF->GetXaxis()->CenterTitle();
-	hTOF->GetXaxis()->SetTitleSize(0.06);
-	hTOF->GetXaxis()->SetTitleOffset(0.75);
-	hTOF->GetYaxis()->CenterTitle();
-	hTOF->GetYaxis()->SetTitleSize(0.06);
-	hTOF->GetYaxis()->SetTitleOffset(0.75);
+	TH2F *hTOF3 = (TH2F*)gPad->GetPrimitive("htemp");
+//	hTOF->SetStats(0);
+	hTOF3->SetTitle("trigger 3");
+	hTOF3->SetXTitle("TOF (ns)");
+	hTOF3->SetYTitle("dE (a.u.)");
+	hTOF3->GetXaxis()->CenterTitle();
+	hTOF3->GetXaxis()->SetTitleSize(0.06);
+	hTOF3->GetXaxis()->SetTitleOffset(0.75);
+	hTOF3->GetYaxis()->CenterTitle();
+	hTOF3->GetYaxis()->SetTitleSize(0.06);
+	hTOF3->GetYaxis()->SetTitleOffset(0.75);
 
 	c1->cd(2);
 	tr->Draw("(y1[0]-16)*1.25:(x1[0]-16.)*1.25>>hXY1(32, -20., 20., 32, -20., 20.)", "x1[1]<1000 && y1[1]<1000 && trigger==1 && nx1==1 && ny1==1", "col", drawEntries);