diff --git a/analysis/fillChain.cxx b/analysis/fillChain.cxx index 787ec59350d301eea41b1b4e14652063c20af74d..6521f7b180effabacbf7f9da329a0889844beb7a 100644 --- a/analysis/fillChain.cxx +++ b/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++) { diff --git a/analysis/showBananasHe.cxx b/analysis/showBananasHe.cxx index ff2b55291daed8daf5a31f58cbe2e1f5c8f9ac3c..7a6a4232b89705f87dff6e6a485ffe270daeeac1 100644 --- a/analysis/showBananasHe.cxx +++ b/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; + diff --git a/analysis/showBeam.cxx b/analysis/showBeam.cxx index bca881213cc006620e36dacdb8dbfa89f2c11057..a6b39deefc988cf77d15eccdad8349963ec91e66 100644 --- a/analysis/showBeam.cxx +++ b/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);