ReadMiniTree3.C

Go to the documentation of this file.

/**************************************************************************
 * Copyright(c) 1998, Agata Experiment, All rights reserved. *
 *                                                                        *
 * Author: The Agata Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/


//
//
#if !defined(__CINT__) || defined(__MAKECINT__)

#include <Riostream.h>
#include "TFile.h"
#include "TROOT.h"
#include "TCanvas.h"
#include "TFolder.h"
#include "TTree.h"
#include "TString.h"
#include "TH1F.h"
#include "TStopwatch.h"
#include "CorrelatedSpaceTree2I.h"
#include "CubicFilter.h"

#endif
Int_t dim = 256;
Int_t nDiv = 8;
Int_t nTree = 512;

//Int_t dim = 1024;
//Int_t nDiv = 2;
//Int_t nTree = 8;

std::vector<int> idxTree;;

void AddIdx(Int_t idx)
{
   Bool_t flag = false;
   for (UInt_t i = 0; i < idxTree.size(); ++i) {
          if (idx == idxTree[i])
                 flag = true;
   }  
   if (!flag)
          idxTree.push_back(idx);
}

void SetMiniCubeId(UShort_t* xCube, UShort_t* xEnergies, Int_t idx)
{
   Int_t i = idx/(nDiv*nDiv);
   Int_t j = (idx - i*(nDiv*nDiv))/nDiv;
   Int_t k = (idx - i*(nDiv*nDiv) - j*nDiv);
   
   xEnergies[0] = (xCube[0] + i*dim);
   xEnergies[1] = (xCube[1] + j*dim);
   xEnergies[2] = (xCube[2] + k*dim);   
}


//__________________________________________________________
Int_t IsInside(UShort_t* energies, Int_t mult, Short_t* inGate, Int_t* gateInf, Int_t* gateSup, Int_t multg)
{   
   Int_t nGate = 0;

   for(Int_t i = 0; i < mult; ++i) {
          
          inGate[i] = 0;
          Bool_t IsInGate = false;
          
          for(Int_t k = 0; k < multg; ++k) {
                 if( (energies[i] >= gateInf[k]) && (energies[i] <= gateSup[k]) ) {
                        IsInGate = true;
                 }
          }  
          
          if (IsInGate) {
                 inGate[i] = 1;
                 nGate++;
          }
   }
   
   return nGate;
}

//__________________________________________________________
Int_t IsInside(Int_t* energies, Int_t mult, Short_t* inGate, UShort_t* gateInf, UShort_t* gateSup, Int_t multg)
{   
   Int_t nGate = 0;
   Int_t usedGate[multg];
   
   for(Int_t k = 0; k < multg; ++k) 
          usedGate[k] = 0;
   
   for(Int_t i = 0; i < mult; ++i) {
          
          inGate[i] = 0;
          Bool_t IsInGate = false;
          
          for(Int_t k = 0; k < multg; ++k) {
                 if( (energies[i] >= gateInf[k]) && (energies[i] <= gateSup[k]) ) {
                        if (usedGate[k] == 0) {
                           IsInGate = true;
                           usedGate[k] = 1;
                           break;
                        }
                 }
          }  
          
          if (IsInGate) {
                 inGate[i] = 1;
                 nGate++;
          }
   }
   
   return nGate;
}

//__________________________________________________________
void ReadMiniTree3(Int_t gateInf1 = 195, Int_t gateSup1 = 205, Int_t gateInf2 = 292, Int_t gateSup2 = 300,
//void ReadMiniTree3(Int_t gateInf1 = 195, Int_t gateSup1 = 205, Int_t gateInf2 = 1322, Int_t gateSup2 = 1330,
                                   TString fileName = "test8.root")
{
  
   Int_t hDim = 2000;
   TStopwatch watch;
   UShort_t e[nTree][3];
   UShort_t er[3];
   UShort_t limInf[3] = {0, 0, 0};
   UShort_t limSup[3] = {1023, 1023, 1023};
   UShort_t limInfr[3];
   UShort_t limSupr[3];
   
   Int_t gateInf[2] = {gateInf1, gateInf2};
   Int_t gateSup[2] = {gateSup1, gateSup2};
   
   Int_t multg      = 2;
   Short_t inGate[3];
   Int_t nofGates  = 0;
   
   
   Int_t idx = 0;

   for (Int_t iTree = 0; iTree < nTree; ++iTree) {
          SetMiniCubeId(limInf, limInfr, iTree);
          SetMiniCubeId(limSup, limSupr, iTree);
          printf("%d %d %d %d\n", iTree, limInfr[0], limInfr[1], limInfr[2]);
          printf("%d %d %d %d\n\n", iTree, limSupr[0], limSupr[1], limSupr[2]);
          nofGates = IsInside(gateInf, 2,  inGate, limInfr, limSupr, 3);
          if (inGate[0]*inGate[1] == 1) {
                 AddIdx(iTree);
          }     
          nofGates = IsInside(gateSup, 2,  inGate, limInfr, limSupr, 3);

          if (inGate[0]*inGate[1] == 1) {
                 AddIdx(iTree);
          }       
          
   }
   
   
//   idxTree.push_back(0);
//   idxTree.push_back(1);
//   idxTree.push_back(2);
//   idxTree.push_back(3);
//   idxTree.push_back(4);
//   idxTree.push_back(5);
//   idxTree.push_back(6);
//   idxTree.push_back(7);

   TTree* oak[nTree];
   TBranch* b[nTree];
   TFile* fileIn = new TFile(fileName.Data(), "READ");
   for (Int_t i = 0; i < nTree; ++i) {
          oak[i] = (TTree*)fileIn->Get(Form("oak%d", i));
          b[i] = oak[i]->GetBranch("e");
          b[i]->SetAddress(&e[i]);
   }
   
   TH1F* h = new TH1F("h", "Projection", hDim, 0, hDim);
   
   watch.Start();
   for (UInt_t iTree = 0; iTree < idxTree.size(); ++iTree) {
          idx = idxTree[iTree];
          
          Int_t nEntries = b[idx]->GetEntries();
          printf("%d %d\n", idx, nEntries);
          
          for (Int_t i = 0; i < nEntries; ++i) {
                 b[idx]->GetEntry(i);
                 SetMiniCubeId(e[idx], er, idx);
                 
                 nofGates = IsInside(er, 3, inGate, gateInf, gateSup, multg);
                 
                 for (Int_t k= 0; k < 3; ++k) {
                        if (nofGates == 2) {
                           if (inGate[k] == 0) {
                                  h->Fill(er[k]); 
                           }
                        }
                        if (nofGates == 3) {
                           h->Fill(er[k]); 
                        }
                 }
          }
   }
   watch.Stop();
   h->Draw();
   printf("Real time %6.3f\n", watch.RealTime());
}