Evio2rootAnalysis limits.C

#include <math.h>

void evio2rootAnalysis_Limits()
{

  double theta,phi,PhiLimit;
  double thetaWire,phiWire,vector;
  double length;

  int      gpart;
  int     gpid[20];
  double  gpx[20];
  double  gpy[20];
  double  gpz[20];
  double  gx[20],gy[20],gz[20];

  Float_t   theoryLeft;
  Float_t   theoryRight;

  int  nhits = 0;
  int  dpid[105000];
  int  mpid[105000];
  int  mtid[105000];
  int  otid[105000];
  double X[105000];
  double Y[105000];
  double Z[105000];
  int  procID[105000];
  int  sector[105000];
  int  superlayer[500];
  int  layer[105000];
  int  wire[105000];
  int  dnhits[105000];

  TFile *f = new TFile("LH2_0Sol_0Tor_11GeV_0Phi_ShieldOut.root");
  TTree *T = (TTree*) f->Get("clas12");

  T->GetBranch("GenPart")->GetLeaf("gpart")->SetAddress(&gpart);
  T->GetBranch("GenPart")->GetLeaf("pid")->SetAddress(gpid);
  T->GetBranch("GenPart")->GetLeaf("px")->SetAddress(gpx);
  T->GetBranch("GenPart")->GetLeaf("py")->SetAddress(gpy);
  T->GetBranch("GenPart")->GetLeaf("pz")->SetAddress(gpz);
  T->GetBranch("GenPart")->GetLeaf("x")->SetAddress(gx);
  T->GetBranch("GenPart")->GetLeaf("y")->SetAddress(gy);
  T->GetBranch("GenPart")->GetLeaf("z")->SetAddress(gz);
  T->Notify();

  T->GetLeaf("theoryLeft")->SetAddress(&theoryLeft);
  T->GetLeaf("theoryRight")->SetAddress(&theoryRight);
  T->Notify();

  T->GetBranch("Detector")->GetLeaf("nhits")->SetAddress(&nhits);
  T->GetBranch("Detector")->GetLeaf("dpid")->SetAddress(dpid);
  T->GetBranch("Detector")->GetLeaf("mpid")->SetAddress(mpid);
  T->GetBranch("Detector")->GetLeaf("mtid")->SetAddress(mtid);
  T->GetBranch("Detector")->GetLeaf("otid")->SetAddress(otid);
  T->GetBranch("Detector")->GetLeaf("X")->SetAddress(X);
  T->GetBranch("Detector")->GetLeaf("Y")->SetAddress(Y);
  T->GetBranch("Detector")->GetLeaf("Z")->SetAddress(Z);
  T->GetBranch("Detector")->GetLeaf("procID")->SetAddress(procID);
  T->GetBranch("Detector")->GetLeaf("sector")->SetAddress(sector);
  T->GetBranch("Detector")->GetLeaf("superlayer")->SetAddress(superlayer);
  T->GetBranch("Detector")->GetLeaf("layer")->SetAddress(layer);
  T->GetBranch("Detector")->GetLeaf("wire")->SetAddress(wire);
  T->Notify();

  wire_test=1;



  for(int i = 0; i < 33400; i++)
  {

        T->GetEntry(i);

                        for(int k = 0; k < nhits; k++)
                        {
                                vector=X[k]*X[k]+Y[k]*Y[k]+Z[k]*Z[k];
                                vector=sqrt(vector);
                                theta=acos(Z[k]/vector)*180/3.14159265359;

                                vector=X[k]*X[k]+Y[k]*Y[k];
                                vector=sqrt(vector);
                                phi=acos(X[k]/vector)*180/3.14159265359;
                                PhiLimit=7.39537e-12*pow(theta,9)-1.81135e-9*pow(theta,8)
                                        +1.92315e-7*pow(theta,7)-0.0000116066*pow(theta,6)
                                        +0.000438844*pow(theta,5)-0.0108036*pow(theta,4)
                                        +0.174341*pow(theta,3)-1.8094*pow(theta,2)
                                        +11.4873*pow(theta,1)-8.3604;

                                        if(phi<0)
                                                PhiLimit=-PhiLimit;

                                        if(( phi < 0 && phi < PhiLimit ) || ( phi >= 0 && phi > PhiLimit ) ) //Not Within sector limits of theta and phi
                                        {

                                                if(superlayer[k]==1 && layer[k]==2 && dpid[k]==11 && mpid[k]==0 && mtid[k]==0 && otid[k]==2 && procID[k]==90 && sector[k]==1)
                                                        printf("{%.3f , %.3f , %.3f , %i , Event=%i\t hit=%i},",X[k],Y[k],Z[k],wire[k],i,k);
                                        }
                        }

    }
}