X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=RICH%2FAliRICHDetect.cxx;h=ba7b1dcffd9b1620808499985d03ea414bdd1283;hb=88cb7938ca21d4a80991d4e7aa564008c29340f7;hp=b6e8f433cb9e1a6d3baf22b96af74f5cd5fe8d81;hpb=6e585aa20558e35335cd77fe14351976b68f4cfd;p=u%2Fmrichter%2FAliRoot.git diff --git a/RICH/AliRICHDetect.cxx b/RICH/AliRICHDetect.cxx index b6e8f433cb9..ba7b1dcffd9 100644 --- a/RICH/AliRICHDetect.cxx +++ b/RICH/AliRICHDetect.cxx @@ -13,54 +13,7 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -/* - $Log$ - Revision 1.15 2001/10/21 18:31:23 hristov - Several pointers were set to zero in the default constructors to avoid memory management problems - - Revision 1.14 2001/05/14 13:25:54 hristov - stdlib.h included (for Alpha) - - Revision 1.13 2001/05/10 12:26:31 jbarbosa - Totally reworked version of reconstruction algorithm. - - Revision 1.12 2001/02/27 22:15:03 jbarbosa - Removed compiler warning. - - Revision 1.11 2001/02/27 15:21:46 jbarbosa - Transition to SDigits. - - Revision 1.10 2001/02/13 20:39:06 jbarbosa - Changes to make it work with new IO. - - Revision 1.9 2001/01/22 21:39:11 jbarbosa - Several tune-ups - - Revision 1.8 2000/11/15 15:52:53 jbarbosa - Turned on spot algorithm. - - Revision 1.7 2000/11/01 15:37:05 jbarbosa - Updated to use its own rec. point object. - - Revision 1.6 2000/10/02 21:28:12 fca - Removal of useless dependecies via forward declarations - - Revision 1.5 2000/06/30 16:30:28 dibari - Disabled writing to rechits. - - Revision 1.4 2000/06/15 15:46:59 jbarbosa - Corrected compilation errors on HP-UX (replaced pow with TMath::Power) - - Revision 1.3 2000/06/13 13:15:41 jbarbosa - Still some code cleanup done (variable names) - - Revision 1.2 2000/06/12 15:19:30 jbarbosa - Cleaned up version. - - Revision 1.1 2000/04/19 13:05:14 morsch - J. Barbosa's spot reconstruction algorithm. - -*/ +/* $Id$ */ #include @@ -70,6 +23,8 @@ #include "AliRICHDetect.h" #include "AliRICHHit.h" #include "AliRICHDigit.h" +#include "AliRICHRawCluster.h" +#include "AliRICHCerenkov.h" #include "AliRICHSegmentationV0.h" #include "AliRun.h" #include "TParticle.h" @@ -79,8 +34,8 @@ #include "TH3.h" #include "TH2.h" #include "TCanvas.h" +#include -#include "malloc.h" ClassImp(AliRICHDetect) @@ -101,13 +56,19 @@ AliRICHDetect::AliRICHDetect(const char *name, const char *title) : TObject() { + TStyle *mystyle=new TStyle("Plain","mystyle"); + mystyle->SetPalette(1,0); + mystyle->cd(); + fc1= new TCanvas("c1","Reconstructed points",50,50,300,350); fc1->Divide(2,2); - fc2= new TCanvas("c2","Reconstructed points after SPOT",50,50,300,350); + fc2= new TCanvas("c2","Reconstructed points after SPOT",370,50,300,350); fc2->Divide(2,2); - fc3= new TCanvas("c3","Used Digits",50,50,300,350); - //fc3->Divide(2,1); + fc3= new TCanvas("c3","Used Digits",690,50,300,350); + fc4= new TCanvas("c4","Mesh activation data",50,430,600,350); + fc4->Divide(2,1); + } @@ -120,7 +81,7 @@ AliRICHDetect::~AliRICHDetect() } -void AliRICHDetect::Detect(Int_t nev) +void AliRICHDetect::Detect(Int_t nev, Int_t type) { // @@ -128,9 +89,9 @@ void AliRICHDetect::Detect(Int_t nev) //printf("Detection started!\n"); - Float_t omega,omega1,theta1,steptheta,stepphi,x,y,z,cx,cy,l,aux1,aux2,aux3,max,radius=0,meanradius=0; + Float_t omega,omega1,theta1,phi_relative,steptheta,stepphi,x,y,q=0,z,cx,cy,l,aux1,aux2,aux3,max,radius=0,meanradius=0; Int_t maxi,maxj,maxk; - //Float_t theta,phi,realomega,realtheta; + Float_t originalOmega, originalPhi, originalTheta; Float_t binomega, bintheta, binphi; Int_t intomega, inttheta, intphi; Int_t i,j,k; @@ -155,41 +116,23 @@ void AliRICHDetect::Detect(Int_t nev) const Int_t kSpot=0; //number of passes with spot algorithm - const Int_t kDimensionTheta=30; //Matrix dimension for angle Detection - const Int_t kDimensionPhi=45; + const Int_t kDimensionTheta=100; //Matrix dimension for angle Detection + const Int_t kDimensionPhi=100; const Int_t kDimensionOmega=100; - const Float_t SPOTp=1; //Percentage of spot action - const Float_t kMinOmega=20*kPi/180; - const Float_t kMaxOmega=70*kPi/180; //Maximum Cherenkov angle to identify + const Float_t SPOTp=.25; //Percentage of spot action + const Float_t kMinOmega=.4; + const Float_t kMaxOmega=.7; //Maximum Cherenkov angle to identify const Float_t kMinTheta=0; - const Float_t kMaxTheta=15*kPi/180; - //const Float_t kMaxTheta=0.1; + const Float_t kMaxTheta=10*kPi/180; const Float_t kMinPhi=0; const Float_t kMaxPhi=360*kPi/180; - - Float_t kCorr=0.61; //Correction factor, accounting for aberration, refractive index, etc. - //const Float_t kCorr=.9369; //from 0 incidence - //const Float_t kCorr=1; - - //TRandom* random=0; - Float_t rechit[6]; //Reconstructed point data - - - //printf("Creating matrices\n"); - //Float_t point[kDimensionTheta][kDimensionPhi][kDimensionOmega]; - //Float_t point1[kDimensionTheta][kDimensionPhi][kDimensionOmega]; - //printf("Created matrices\n"); - Int_t ***point = i3tensor(0,kDimensionTheta,0,kDimensionPhi,0,kDimensionOmega); Int_t ***point1 = i3tensor(0,kDimensionTheta,0,kDimensionPhi,0,kDimensionOmega); - //Int_t **point = new Int_t[kDimensionTheta][kDimensionPhi][kDimensionOmega]; - //Int_t **point1 = new Int_t[kDimensionTheta][kDimensionPhi][kDimensionOmega]; - steptheta=(kMaxTheta-kMinTheta)/kDimensionTheta; stepphi=(kMaxPhi-kMinPhi)/kDimensionPhi; @@ -202,6 +145,8 @@ void AliRICHDetect::Detect(Int_t nev) static TH2F *SpotOmegaTheta = new TH2F("OmegaTheta","Omega-Theta projection, spot",kDimensionTheta,0,kDimensionTheta,kDimensionOmega,0,kDimensionOmega); static TH2F *SpotOmegaPhi = new TH2F("OmegaPhi","Omega-Phi projection, spot",kDimensionPhi,0,kDimensionPhi,kDimensionOmega,0,kDimensionOmega); static TH2F *DigitsXY = new TH2F("DigitsXY","Pads used for reconstruction",150,-25,25,150,-25,25); + static TH1F *AngleAct = new TH1F("AngleAct","Activation per angle",100,.45,1); + static TH1F *Activation = new TH1F("Activation","Activation per ring",100,0,25); Points->SetXTitle("theta"); Points->SetYTitle("phi"); Points->SetZTitle("omega"); @@ -220,28 +165,30 @@ void AliRICHDetect::Detect(Int_t nev) SpotOmegaTheta->SetYTitle("omega"); SpotOmegaPhi->SetXTitle("phi"); SpotOmegaPhi->SetYTitle("omega"); + AngleAct->SetFillColor(5); + AngleAct->SetXTitle("rad"); + AngleAct->SetYTitle("activation"); + Activation->SetFillColor(5); + Activation->SetXTitle("activation"); + + Int_t ntracks = (Int_t)pRICH->TreeH()->GetEntries(); - Int_t ntracks = (Int_t)gAlice->TreeH()->GetEntries(); - //Int_t ntrks = gAlice->GetNtrack(); - Float_t trackglob[3]; Float_t trackloc[3]; //printf("Area de uma elipse com teta 0 e Omega 45:%f",Area(0,45)); - + Int_t track; for (track=0; trackResetHits(); - gAlice->TreeH()->GetEvent(track); + pRICH->TreeH()->GetEvent(track); TClonesArray *pHits = pRICH->Hits(); if (pHits == 0) return; Int_t nhits = pHits->GetEntriesFast(); if (nhits == 0) continue; //Int_t nent=(Int_t)gAlice->TreeD()->GetEntries(); - gAlice->TreeD()->GetEvent(0); AliRICHHit *mHit = 0; - AliRICHDigit *points = 0; //Int_t npoints=0; Int_t counter=0, counter1=0; @@ -261,219 +208,272 @@ void AliRICHDetect::Detect(Int_t nev) } } } - mHit = (AliRICHHit*) pHits->UncheckedAt(0); - //printf("Aqui vou eu\n"); - Int_t nch = mHit->Chamber(); - //printf("Aqui fui eu\n"); - trackglob[0] = mHit->X(); - trackglob[1] = mHit->Y(); - trackglob[2] = mHit->Z(); - printf("Chamber processed:%d\n",nch); + Int_t ncerenkovs = pRICH->Cerenkovs()->GetEntriesFast(); + + + originalOmega = 0; + counter = 0; + + for (Int_t hit=0;hitCerenkovs()->UncheckedAt(hit); + Float_t loss = cHit->fLoss; //did it hit the CsI? + Float_t production = cHit->fProduction; //was it produced in freon? + Float_t cherenkov = cHit->fCerenkovAngle; //production cerenkov angle + if (loss == 4 && production == 1) + { + counter +=1; + originalOmega += cherenkov; + //printf("%f\n",cherenkov); + } + } + + printf("Cerenkovs : %d\n",counter); + + if(counter != 0) //if there are cerenkovs + { + originalOmega = originalOmega/counter; + printf("Original omega: %f\n",originalOmega); + - printf("Reconstructing particle at (global coordinates): %3.1f %3.1f %3.1f,\n",trackglob[0],trackglob[1],trackglob[2]); - iChamber = &(pRICH->Chamber(nch-1)); - - //printf("Nch:%d\n",nch); + mHit = (AliRICHHit*) pHits->UncheckedAt(0); + Int_t nch = mHit->Chamber(); + originalTheta = mHit->Theta(); + originalPhi = mHit->Phi(); + trackglob[0] = mHit->X(); + trackglob[1] = mHit->Y(); + trackglob[2] = mHit->Z(); + + printf("\n--------------------------------------\n"); + printf("Chamber %d, track %d\n", nch, track); - iChamber->GlobaltoLocal(trackglob,trackloc); + + iChamber = &(pRICH->Chamber(nch-1)); - printf("Reconstructing particle at (local coordinates) : %3.1f %3.1f %3.1f\n",trackloc[0],trackloc[1],trackloc[2]); + //printf("Nch:%d\n",nch); - - iChamber->LocaltoGlobal(trackloc,trackglob); - - //printf("Transformation 2: %3.1f %3.1f %3.1f\n",trackglob[0],trackglob[1],trackglob[2]); + iChamber->GlobaltoLocal(trackglob,trackloc); - cx=trackloc[0]; - cy=trackloc[2]; + iChamber->LocaltoGlobal(trackloc,trackglob); + + + cx=trackloc[0]; + cy=trackloc[2]; + AliRICHDigit *points = 0; + TClonesArray *pDigits = pRICH->DigitsAddress(nch-1); + + AliRICHRawCluster *cluster =0; + TClonesArray *pClusters = pRICH->RawClustAddress(nch-1); - TClonesArray *pDigits = pRICH->DigitsAddress(nch-1); - Int_t ndigits = pDigits->GetEntriesFast(); - - //printf("Got %d digits\n",ndigits); - - counter=0; - printf("Starting calculations\n"); - for(Float_t theta=0;thetaTreeD()->GetEvent(0); + Int_t ndigits = pDigits->GetEntriesFast(); + maxcycle=ndigits; + printf("Got %d digits\n",ndigits); + } + + //digitize from clusters + + if(type==1) + { + Int_t nent=(Int_t)gAlice->TreeR()->GetEntries(); + gAlice->TreeR()->GetEvent(nent-1); + Int_t nclusters = pClusters->GetEntriesFast(); + maxcycle=nclusters; + printf("Got %d clusters\n",nclusters); + } + + + + + counter=0; + printf("Starting calculations\n"); + printf(" Start Finish\n"); + printf("Progress: "); + for(Float_t theta=0;thetaUncheckedAt(dig); - segmentation->GetPadC(points->PadX(), points->PadY(),x, y, z); - x=x-cx; - y=y-cy; - radius=TMath::Sqrt(TMath::Power(x,2)+TMath::Power(y,2)); - - if(radius>4) - { - //if(theta==0 && phi==0) - //{ - //printf("Radius: %f, Max Radius: %f\n",radius,kCorr*kHeight*tan(theta+kMaxOmega)*3/4); - meanradius+=radius; - counter++; - //} + printf("."); + for(Float_t phi=0;phi<=kMaxPhi;phi+=stepphi) + { + //printf("Phi:%3.1f\n", phi*180/kPi); + counter1=0; + for (Int_t cycle=0;cycleUncheckedAt(cycle); + segmentation->GetPadC(points->PadX(), points->PadY(),x, y, z); + } + + if(type==1) + { + cluster=(AliRICHRawCluster*) pClusters->UncheckedAt(cycle); + x=cluster->fX; + y=cluster->fY; + q=cluster->fQ; + } + + if(type ==0 || q > 100) + { - if(phi==0) - { - //printf("Radius: %f, Max Radius: %f\n",radius,2*kHeight*tan(theta+kMaxOmega)*3/4); - //printf("Loaded digit %d with coordinates x:%f, y%f\n",dig,x,y); - //printf("Using digit %d, for theta %f\n",dig,theta); - } - - counter1++; - - l=kHeight/cos(theta); + x=x-cx; + y=y-cy; + radius=TMath::Sqrt(TMath::Power(x,2)+TMath::Power(y,2)); - //x=x*kCorr; - //y=y*kCorr; - /*if(SnellAngle(theta+omega)<999) - { - //printf("(Angle)/(Snell angle):%f\n",(theta+omega)/SnellAngle(theta+omega)); - x=x*(theta+omega)/SnellAngle(theta+omega); - y=y*(theta+omega)/SnellAngle(theta+omega); - } - else - { - x=0; - y=0; - }*/ - - //main calculation - - DigitsXY->Fill(x,y,(float) 1); - - theta1=SnellAngle(theta)*1.5; - - aux1=-y*sin(phi)+x*cos(phi); - aux2=y*cos(phi)+x*sin(phi); - aux3=( TMath::Power(aux1,2)+TMath::Power(cos(theta1)*aux2 ,2))/TMath::Power(sin(theta1)*aux2+l,2); - omega=atan(sqrt(aux3)); + //calculation of relative phi angle of digit - //omega is distorted, theta1 is distorted - - if(InvSnellAngle(theta+omega)<999) - { - omega1=InvSnellAngle(omega+theta1) - theta; - //theta1=InvSnellAngle(omega+theta) - omega1; - //omega1=kCorr*omega; - - kCorr=InvSnellAngle(omega+theta)/(omega+theta); - theta1=kCorr*theta/1.4; - //if(phi==0) - //printf("Omega:%f Theta:%f Omega1:%f Theta1:%f ISA(o+t):%f ISA(t):%f\n",omega*180/kPi,theta*180/kPi,omega1*180/kPi,theta1*180/kPi,InvSnellAngle(omega+theta)*180/kPi,InvSnellAngle(theta)*180/kPi); - } - else - { - omega1=0; - theta1=0; - } + phi_relative = acos(y/radius); + phi_relative = TMath::Abs(phi_relative - phi); - //printf("Omega:%f\n",omega); - - - //if(SnellAngle(theta+omega)<999) - //printf("(Angle)/(Snell angle):%f\n",(theta+omega)/SnellAngle(theta+omega)); - if(theta==0 && phi==0) - { - //printf("Omega: %f Corrected Omega: %f\n",omega, omega/kCorr); - //omega=omega/kCorr; - } - //cout<<"\ni="<Fill(inttheta,intphi,intomega,(float) 1); - ThetaPhi->Fill(inttheta,intphi,(float) 1); - OmegaTheta->Fill(inttheta,intomega,(float) 1); - OmegaPhi->Fill(intphi,intomega,(float) 1); - //printf("Filling at %d %d %d\n",Int_t(theta*kDimensionTheta/kMaxTheta),Int_t(phi*kDimensionPhi/kMaxPhi),Int_t(omega*kDimensionOmega/kMaxOmega)); + //if (radius < (2*kHeight*tan(theta+kMaxOmega))) + if (radius < (2*kHeight*tan(kMaxOmega))) + //if(Fiducial(x,y,theta,phi,kHeight,kMaxOmega,kMinOmega)) + { + + if(phi==0) + { + //printf("Radius: %f, Max Radius: %f\n",radius,2*kHeight*tan(theta+kMaxOmega)*3/4); + //printf("Loaded digit %d with coordinates x:%f, y%f\n",dig,x,y); + //printf("Using digit %d, for theta %f\n",dig,theta); + } + + counter1++; + + l=kHeight/cos(theta); + + //main calculation + + DigitsXY->Fill(x,y,(float) 1); + + theta1=SnellAngle(theta); + + aux1=-y*sin(phi)+x*cos(phi); + aux2=y*cos(phi)+x*sin(phi); + aux3=( TMath::Power(aux1,2)+TMath::Power(cos(theta1)*aux2 ,2))/TMath::Power(sin(theta1)*aux2+l,2); + omega=atan(sqrt(aux3)); + + //omega is distorted, theta1 is distorted + + if(InvSnellAngle(omega+TMath::Abs(cos(phi_relative))*theta1)<999) + { + omega1=InvSnellAngle(omega+TMath::Abs(cos(phi_relative))*theta); + theta1=InvSnellAngle(theta1); + + } + else + { + omega1=0; + theta1=0; + } + + if(omega1kMinOmega) + { + //printf("Omega found:%f\n",omega); + omega1=omega1-kMinOmega; + + //printf("Omega: %f Theta: %3.1f Phi:%3.1f\n",omega, theta*180/kPi, phi*180/kPi); + + bintheta=theta1*kDimensionTheta/kMaxTheta; + binphi=phi*kDimensionPhi/kMaxPhi; + binomega=omega1*kDimensionOmega/(kMaxOmega-kMinOmega); + + if(Int_t(bintheta+0.5)==Int_t(bintheta)) + inttheta=Int_t(bintheta); + else + inttheta=Int_t(bintheta+0.5); + + if(Int_t(binomega+0.5)==Int_t(binomega)) + intomega=Int_t(binomega); + else + intomega=Int_t(binomega+0.5); + + if(Int_t(binphi+0.5)==Int_t(binphi)) + intphi=Int_t(binphi); + else + intphi=Int_t(binphi+0.5); + + //printf("Point added at %d %d %d\n",inttheta,intphi,intomega); + + if(type==0) + point[inttheta][intphi][intomega]+=1; + if(type==1) + point[inttheta][intphi][intomega]+=(Int_t)(q); + + //printf("Omega stored:%d\n",intomega); + Points->Fill(inttheta,intphi,intomega,(float) 1); + ThetaPhi->Fill(inttheta,intphi,(float) 1); + OmegaTheta->Fill(inttheta,intomega,(float) 1); + OmegaPhi->Fill(intphi,intomega,(float) 1); + //printf("Filling at %d %d %d\n",Int_t(theta*kDimensionTheta/kMaxTheta),Int_t(phi*kDimensionPhi/kMaxPhi),Int_t(omega*kDimensionOmega/kMaxOmega)); + } + //if(omegacd(1); - Points->Draw(); - fc1->cd(2); - ThetaPhi->Draw(); - fc1->cd(3); - OmegaTheta->Draw(); - fc1->cd(4); - OmegaPhi->Draw(); - fc3->cd(); - DigitsXY->Draw(); - } - else + meanradius=meanradius/counter; + //printf("Mean radius:%f, counter:%d\n",meanradius,counter); + rechit[5]=meanradius; + //printf("Used %d digits\n",counter1); + //printf("\n"); + + if(nev<2) { - //fc1->cd(1); - //Points->Draw("same"); - //fc1->cd(2); - //ThetaPhi->Draw("same"); - //fc1->cd(3); - //OmegaTheta->Draw("same"); - //fc1->cd(4); - //OmegaPhi->Draw("same"); + if(nev==0) + { + fc1->cd(1); + Points->Draw("colz"); + fc1->cd(2); + ThetaPhi->Draw("colz"); + fc1->cd(3); + OmegaTheta->Draw("colz"); + fc1->cd(4); + OmegaPhi->Draw("colz"); + fc3->cd(); + DigitsXY->Draw("colz"); + } + else + { + //fc1->cd(1); + //Points->Draw("same"); + //fc1->cd(2); + //ThetaPhi->Draw("same"); + //fc1->cd(3); + //OmegaTheta->Draw("same"); + //fc1->cd(4); + //OmegaPhi->Draw("same"); + } } - } - //SPOT execute twice - for(Int_t s=0;scd(1); - SpotPoints->Draw(); - fc2->cd(2); - SpotThetaPhi->Draw(); - fc2->cd(3); - SpotOmegaTheta->Draw(); - fc2->cd(4); - SpotOmegaPhi->Draw(); - } - else - { - //fc2->cd(1); - //SpotPoints->Draw("same"); - //fc2->cd(2); - //SpotThetaPhi->Draw("same"); - //fc2->cd(3); - //SpotOmegaTheta->Draw("same"); - //fc2->cd(4); - //SpotOmegaPhi->Draw("same"); + if(nev==0) + { + fc2->cd(1); + SpotPoints->Draw("colz"); + fc2->cd(2); + SpotThetaPhi->Draw("colz"); + fc2->cd(3); + SpotOmegaTheta->Draw("colz"); + fc2->cd(4); + SpotOmegaPhi->Draw("colz"); + } + else + { + //fc2->cd(1); + //SpotPoints->Draw("same"); + //fc2->cd(2); + //SpotThetaPhi->Draw("same"); + //fc2->cd(3); + //SpotOmegaTheta->Draw("same"); + //fc2->cd(4); + //SpotOmegaPhi->Draw("same"); + } } - } - - //Identification is equivalent to maximum determination - max=0;maxi=0;maxj=0;maxk=0; - printf(" Proceeding to identification"); - - for(i=0;imax) - { - //cout<<"maxi="<max) + { + //cout<<"maxi="<Fill(FinalOmega, (float) max); + Activation->Fill(max, (float) 1); + + if(nev==0) + { + fc4->cd(1); + AngleAct->Draw(); + fc4->cd(2); + Activation->Draw(); + } + else + { + fc4->cd(1); + AngleAct->Draw("same"); + fc4->cd(2); + Activation->Draw("same"); + } + + + //fscanf(omegas,"%f",&realomega); + //fscanf(thetas,"%f",&realtheta); + //printf("Real Omega: %f",realomega); + //cout<<"Detected:theta="<AddRecHit3D(nch-1,rechit); + pRICH->AddRecHit3D(nch-1,rechit,originalOmega, originalTheta, originalPhi); //printf("rechit %f %f %f %f %f\n",rechit[0],rechit[1],rechit[2],rechit[3],rechit[4]); //printf("Chamber:%d",nch); - } - //printf("\n\n\n\n"); - gAlice->TreeR()->Fill(); + } + + else //if no cerenkovs + + { + + rechit[0] = 0; + rechit[1] = 0; + rechit[2] = 0; + rechit[3] = 0; + rechit[4] = 0; + + } + + } + + if(type==1) //reco from clusters + { + pRICH->ResetRawClusters(); + //Int_t nent=(Int_t)gAlice->TreeR()->GetEntries(); + //gAlice->TreeR()->GetEvent(track); + //printf("Going to branch %d\n",track); + //gAlice->GetEvent(nev); + } + + + //printf("\n\n\n\n"); + gAlice->TreeR()->Fill(); TClonesArray *fRec; for (i=0;iRecHitsAddress3D(i); @@ -832,3 +876,47 @@ void AliRICHDetect::CreatePoints(Float_t theta, Float_t phi, Float_t omega, Floa fc3->cd(2); REllipse->Draw(); } + +Int_t AliRICHDetect::Fiducial(Float_t rotx, Float_t roty, Float_t theta, Float_t phi, Float_t height, Float_t maxOmega, Float_t minOmega) +{ + + Float_t x,y,y1,h,omega1,omega2; + + Float_t a,b, offset; + Float_t a1,b1, offset1; + + h = height; + + //refraction calculation + + theta = SnellAngle(theta); + phi = phi - TMath::Pi(); + omega1 = SnellAngle(maxOmega); + omega2 = SnellAngle(minOmega); + + //maximum zone + a = h*(tan(omega1+theta)+tan(omega1-theta))/2; + b = h*tan(omega1); + + offset = h*(tan(omega1+theta)-tan(omega1-theta))/2; + + + //minimum zone + a1 = h*(tan(omega2+theta)+tan(omega2-theta))/2; + b1 = h*tan(omega2); + + offset1 = h*(tan(omega2+theta)-tan(omega2-theta))/2; + + + //rotation to phi=0 + x = rotx*cos(phi)+roty*sin(phi); + y = -rotx*sin(phi)+roty*cos(phi) - offset; + y1 = -rotx*sin(phi)+roty*cos(phi) - offset1; + + + if(x*x/a+y*y/b<1 && x*x/a1+y1*y1/b1>1) + return 1; + else + return 0; + +}