tags:

views:

9

answers:

0

Hi, This is a Root question, but I thought I'd try my luck in-case anyone has experience!

I have created a piece of code that will filter out final particles from a sample into a new tree. I have succeeded in doing so by creating a tree with repeated fields ie. event number, but am trying to make it better for analysis by using vectors for the particle data (energy, momentum etc) so that I can use one event number easily for all of those particles (I guess its best to describe like a 2D array).

The tree is created and the none vector fields are filled correctly but the vector fields all contain nothing! Am I missing anything out in my declaration or do I have to loop over all particles within that vector?? I'm at a bit of a loss so any help would be greatly appriciated!

Thanks in advance!

Ben

Sorry if there are unused variables I have been playing around with this for a couple of weeks without a breakthrough!


#define PT_test_cxx


#include "PT_test.h"

using namespace std;


const double GeV =0.001;
int EventCounter=0;
int Sample;

vector<float> mytruthPt;
vector<float> mytruthEta;

bool electronmother;
int motherparticle=0;

int k;
int m;
int P2;

int partref=0;

TFile *newfile1;  TTree *FinalPs;

  UInt_t T_N;
  Int_t Event, EventId, ParticleId;
  vector<double> T_Eta, T_Phi, T_Pt, T_E, T_pdgId, T_status;

TString option;

void PT_test::Begin(TTree *tree)
{

  cout<<"Begining"<<endl;

  Init(tree);
  option = GetOption();
  EventCounter=0;

  //create new trees for filtered particles
  newfile1 = new TFile("/export/atlas/bjdj20/"+option+"_AllFinals.root","RECREATE");
  FinalPs = new TTree("FinalPs","All Final Particles");

  FinalPs->Branch("Event",&Event,"Event/I");
  FinalPs->Branch("EventId",&EventId,"EventId/I");
  FinalPs->Branch("ParticleId",&ParticleId,"ParticleId/I");
  FinalPs->Branch("T_N",&T_N,"T_N/I");
  FinalPs->Branch("T_Eta",&T_Eta,"T_Eta[EventId]/D");
  FinalPs->Branch("T_Phi",&T_Phi,"T_Phi[EventId]/D");
  FinalPs->Branch("T_Pt",&T_Pt,"T_Pt[EventId]/D");
  FinalPs->Branch("T_E",&T_E,"T_E[EventId]/D");
  FinalPs->Branch("T_pdgId",&T_pdgId,"T_pdgId[EventId]/D");
  FinalPs->Branch("T_status",&T_status,"T_status[EventId]/D");

}


void PT_test::SlaveBegin(TTree *tree)
{
   // The SlaveBegin() function is called after the Begin() function.
   // When running with PROOF SlaveBegin() is called on each slave server.
   // The tree argument is deprecated (on PROOF 0 is passed).
   TString option = GetOption();
   cout<<"Slave Begins-Option: "<<option<<endl;
}


//Main Process
Bool_t PT_test::Process(Long64_t entry)
{
  TString OpName = GetOption();
  fChain->GetTree()->GetEntry(entry);
  EventCounter++;

if(EventCounter<20){  

  if(EventCounter>1){
    T_N=Truth_N;
    cout<<"Event: "<<EventCounter<<endl;
    cout<<"Eta: "<<T_Eta[3]<<endl;
    Event=EventCounter;
    cout<<"Number of Events: "<<P2<<endl;
    EventId=P2;



      //Fill a filtered tree
      FinalPs->Fill();


      T_Eta.clear();
      T_Phi.clear();
      T_Pt.clear();
      T_E.clear();
      T_pdgId.clear();
      T_status.clear();
  }

  //if(EventCounter%10==0){cout<<"Event: "<<EventCounter<<endl;}
  P2=0;
  mytruthPt.clear();
  mytruthEta.clear();

  int motherparticle=0;

  // loop over all Particles
  for (unsigned int particle=0; particle<Truth_N; particle++){

      bool FinalFlag = true;
      bool MotherZ=false;

      //Pt cut
      if ((*Truth_Pt)[particle]>1000){

          //loop over all secondary particles
          for (unsigned int i=0; i<Truth_N; i++){

              //find final particles
              if ((*Truth_OrigVtx)[i]==(*Truth_EndVtx)[particle] || (*Truth_status)[particle]!=1){
                  FinalFlag = false;
              }
              //find final particles with 
              if (((*Truth_OrigVtx)[particle]==(*Truth_EndVtx)[i]) && ((*Truth_pdgId)[i] == 23)){
                  MotherZ = true;
              }
              //min bias mother particles
              if (((*Truth_OrigVtx)[particle]==(*Truth_EndVtx)[i])&&((*Truth_pdgId)[i] != 23)){
                  motherparticle=i;                  
              }

          }


 //processing for final particles
          if (FinalFlag==true){
              bool unique=true;
              m=0;

              //checks the vector array for duplicates
              for (unsigned int x=1; x<mytruthPt.size(); x++){
                 if((*Truth_Pt)[particle]==mytruthPt[x] && (*Truth_Eta)[particle]==mytruthEta[x]){
                     unique=false;
                 }
             }

             //processing for unique final particles
             if (unique==true) {

                //add unique entries to the vector
                mytruthPt.push_back((*Truth_Pt)[particle]);
                mytruthEta.push_back((*Truth_Eta)[particle]);


                EventId=P2;
                T_Eta.push_back((*Truth_Eta)[particle]);

                T_Phi.push_back((*Truth_Phi)[particle]);

                T_Pt.push_back((*Truth_Pt)[particle]);

                T_E.push_back((*Truth_E)[particle]);

                T_pdgId.push_back((*Truth_pdgId)[particle]);

                T_status.push_back((*Truth_status)[particle]);
                P2++;


             }    //unique
         }    //final flag
     }    //PT>1000MeV    
  }    //all particles

}//events<20
  return kTRUE;
}//end


void PT_test::Terminate()
{
   // The Terminate() function is the last function to be called during
   // a query. It always runs on the client, it can be used to present
   // the results graphically or save the results to file.

  cout<<"This is Terminate"<<endl;

   //Write new Filtered Trees
   newfile1->Write();
   newfile1->Close();
   delete newfile1;


}