void tagghosts(Long64_t Nevts=10000)
{
//TString datadir="/data/glast/data/flight/ghostNumber/";
//TString baseFileName="periodic_nocaltrig_v19r4p1gr13_";

//TString datadir="/data/glast/data/flight/GRBs/";
//TString baseFileName="GRB090926_000015_000_";

//TString datadir="/data/glast/data/flight/SFRs/";
//TString baseFileName="SFR110924_2007_000025_012_";

TString datadir="/data47b/glast/mc/AG-GR-v19r4p1gr13-OVL/runs/";
TString baseFileName="AG-GR-v19r4p1gr13-OVL";

// Load libraries
gSystem->Load("libcommonRootData.so");
gSystem->Load("libdigiRootData.so");
gSystem->Load("libreconRootData.so");

class DigiEvent;
class ReconEvent;
class RelTable;

//TTree      *dgTree,*recTree, *relTree, *merit;
TChain      *dgTree,*recTree, *relTree, *merit;
DigiEvent  *dgEvent;
ReconEvent *recEvent;
RelTable   *relTable;

//TFile *dgFile=new TFile(datadir+baseFileName+"digi.root","READ");
//dgTree = (TTree*)dgFile->Get("Digi");
TChain *dgTree=new TChain("Digi");
dgTree->Add(datadir+"/digi/"+baseFileName+"*digi.root");
dgEvent = 0;
dgTree->SetBranchAddress("DigiEvent",&dgEvent);
cout<<"Digi Nentries "<<dgTree->GetEntries()<<endl;

//TFile *recFile=new TFile(datadir+baseFileName+"recon.root","READ");
//recTree = (TTree*)recFile->Get("Recon");
TChain *recTree=new TChain("Recon");
recTree->Add(datadir+"/recon/"+baseFileName+"*recon.root");
recEvent = 0;
recTree->SetBranchAddress("ReconEvent",&recEvent);
cout<<"Recon Nentries "<<recTree->GetEntries()<<endl;

//TFile *relFile=new TFile(datadir+baseFileName+"rel.root","READ");
//relTree = (TTree*)relFile->Get("Relations");
TChain *relTree=new TChain("Relations");
relTree->Add(datadir+"/relation/"+baseFileName+"*relation.root");
relTable = 0;
relTree->SetBranchAddress("RelTable",&relTable);
cout<<"Relations Nentries "<<relTree->GetEntries()<<endl;

//TFile *meritFile=new TFile(datadir+baseFileName+"merit.root","READ");
//merit = (TTree*)meritFile->Get("MeritTuple");
TChain *merit=new TChain("MeritTuple");
merit->Add(datadir+"/merit/"+baseFileName+"*merit.root");
cout<<"Merit Nentries "<<merit->GetEntries()<<endl;

if(merit->GetEntries()==0) break;

// Variables
Float_t GltGemSummary,CalUberEnergy;
UInt_t EvtRun,EvtEventId,CalXtalMaxEne;
merit->SetBranchAddress("GltGemSummary",&GltGemSummary);
merit->SetBranchAddress("CalUberEnergy",&CalUberEnergy);
merit->SetBranchAddress("EvtRun",&EvtRun);
merit->SetBranchAddress("EvtEventId",&EvtEventId);
merit->SetBranchAddress("CalXtalMaxEne",&CalXtalMaxEne);

Int_t calLo[16][8][2];
Int_t calHi[16][8][2];

Int_t calLoGhost[16][8][2];
Int_t calHiGhost[16][8][2];

Int_t ghostNumber=0;

// Initialize arrays
for(Int_t tower=0;tower<16;tower++){
  for(Int_t layer=0;layer<8;layer++){
    for(Int_t fc=0;fc<2;fc++){
      calLo[tower][layer][fc]=0;
      calHi[tower][layer][fc]=0;
      calLoGhost[tower][layer][fc]=0;
      calHiGhost[tower][layer][fc]=0;
      }
    }
  }

//-------------------------------
// Open output file and Define histograms
TFile *outF = new TFile(baseFileName+"_ghost.root","recreate");
TH2F *hGhostProb = new TH2F("hGhostProb","Ghost Prob vs Ghost Number",
                             40, -1., 1., 100., 0., 100.);
hGhostProb->SetXTitle("Ghost Prob");
hGhostProb->SetYTitle("Ghost Number");

TH2F *hHadProb = new TH2F("hHadProb","Had Prob vs Ghost Number",
                             40, -1., 1., 100., 0., 100.);
hHadProb->SetXTitle("Had Prob");
hHadProb->SetYTitle("Ghost Number");

TH2F *hHadGhostProb = new TH2F("hHadGhostProb","HadGhost Prob vs Ghost Number",
                             20, 0., 1., 100., 0., 100.);
hHadGhostProb->SetXTitle("Had+Ghos Prob");
hHadGhostProb->SetYTitle("Ghost Number");


//-------------------------------
// Event loop
for (Long64_t ievt = 0; ievt < Nevts; ievt++)
  {
   if(ievt%100==0)cout<<"Event "<<ievt<<endl;
   if(dgEvent)  dgEvent->Clear();
   if(recEvent) recEvent->Clear();
   if(relTable) relTable->Clear();

   dgTree->GetEvent(ievt);
   recTree->GetEvent(ievt);
   relTree->GetEvent(ievt);
   merit->GetEvent(ievt);
   
   if(CalXtalMaxEne>50){//(GltGemSummary&4)>0){
     //cout<<"Event "<<ievt<<" EvtRun "<<EvtRun<<" EvtEventId "<<EvtEventId<<" Energy "<<CalUberEnergy<<" GltGemSummary "<<GltGemSummary<<endl;
     //--------------------------------
     //-- read fle diagnostic info
     TClonesArray const*const calDiagCol = dgEvent->getCalDiagnosticCol();
     if (!calDiagCol) {
       cout << "no calDiagCol found for event# " << dgEvent->getEventId() << endl;
       return;
       }
     TIter calDiagIter(calDiagCol);
     const CalDiagnosticData  *pdiag = 0;
     // CalDiagnosticData loop
     while ((pdiag = (CalDiagnosticData*)calDiagIter.Next())) {
//         cout<<"Tower "<<pdiag->tower()<<"\tLayer "<<pdiag->layer()//<< endl;
//             <<"\tFLE "<<pdiag->low(0)<<" "<<pdiag->low(1)
//             <<"\tFHE "<<pdiag->high(0)<<" "<<pdiag->high(1)
//             <<" "<< GltGemSummary <<endl;
        calLo[pdiag->tower()][pdiag->layer()][CalXtalId::POS]=pdiag->low(CalXtalId::POS);
        calLo[pdiag->tower()][pdiag->layer()][CalXtalId::NEG]=pdiag->low(CalXtalId::NEG);
        calHi[pdiag->tower()][pdiag->layer()][CalXtalId::POS]=pdiag->high(CalXtalId::POS);
        calHi[pdiag->tower()][pdiag->layer()][CalXtalId::NEG]=pdiag->high(CalXtalId::NEG);
//    pdiag->Print();
       }

     // Get Relation Tables
     TObjArray const*const tables = relTable->getRelationTable();
     if (!tables) {
        cout << "no relation tables found for event# " << relTable->getEventId() << endl;
        return;
        }
//     TIter relIter(tables);
//     const Relation *rel = 0;
//     while ((rel = (Relation*)relIter.Next())) {
//          cout<<rel->getFirst()<<endl;
//      }
     
     // Get CalRecon cluster information
     TObjArray  const*const calCluCol=recEvent->getCalRecon()->getCalClusterCol();
     if (!calCluCol) {
       cout << "no Cal Cluster Collection found for event# " << recEvent->getEventId() << endl;
       return;
       }

     Int_t printFlag=1, clusterNum=0;
     
     TIter calCluIter(calCluCol);
     const CalCluster *cluster=0;
     while ((cluster = (CalCluster*)calCluIter.Next())) {
         //cout<<"New cluster"<<endl;
         //Initialize cluster ghost bits array
         ghostNumber=0;
         for(Int_t tower=0;tower<16;tower++){
           for(Int_t layer=0;layer<8;layer++){
             for(Int_t fc=0;fc<2;fc++){
               calLoGhost[tower][layer][fc]=0;
               calHiGhost[tower][layer][fc]=0;
               }
             }
           }
           
//         cout<<"Max Ene "<<cluster->getXtalsParams().getXtalEneMax()<<endl;        
         TIter relIter(tables);
         const Relation *rel = 0;
         while ((rel = (Relation*)relIter.Next())) {
           if(rel->getFirst()==cluster || rel->getSecond()==cluster){
             //cout<<rel->getFirst()<<" "<<rel->getSecond()<<" "<<cluster<<endl;
             CalXtalRecData *xtal=(CalXtalRecData*)rel->getFirst();
             Int_t tower=xtal->getPackedId().getTower();
             Int_t layer=xtal->getPackedId().getLayer();
             Float_t faceEP=xtal->getEnergy(0, CalXtalId::POS);
             Float_t faceEN=xtal->getEnergy(0, CalXtalId::NEG);
             if(faceEP>120 && !calLo[tower][layer][CalXtalId::POS]){
                 calLoGhost[tower][layer][CalXtalId::POS]+=1;
                 ghostNumber+=1;
                 //cout<<tower<<" "<<layer<<" "<<faceEP<<" "
                 //    <<calLo[tower][layer][CalXtalId::POS]<<endl;
                 }
             if(faceEN>120 && !calLo[tower][layer][CalXtalId::NEG]){
                 calLoGhost[tower][layer][CalXtalId::NEG]+=1;
                 ghostNumber+=1;
                 //cout<<tower<<" "<<layer<<" "<<faceEN<<" "
                 //    <<calLo[tower][layer][CalXtalId::NEG]<<endl;                 
                 }
             }
        } // loop on CalXtalRecData--cluster relations
        Float_t ghostProb = cluster->getClassParams().getProb("ghost");
        Float_t hadProb = cluster->getClassParams().getProb("had");
        if(ghostNumber){
            if(printFlag==1){
              printFlag=0;
              cout<<"Event "<<ievt<<" EvtRun "<<EvtRun
                  <<" EvtEventId "<<EvtEventId
                  <<" Energy "<<CalUberEnergy<<" GltGemSummary"<<GltGemSummary<<endl;
                  }
          cout<<"Cal Ghost Number "<<ghostNumber
              <<" GhostProb "      <<ghostProb
              <<" HadronProb "     <<hadProb
              <<" StatusBit "      <<hex<<cluster->getStatusBits()<<dec<<endl;
          }
        hGhostProb->Fill(ghostProb, ghostNumber);
        hHadProb->Fill(hadProb, ghostNumber);
        hHadGhostProb->Fill(hadProb+ghostProb, ghostNumber);
        clusterNum++;
//         for(Int_t tower=0;tower<16;tower++){
//           for(Int_t layer=0;layer<8;layer++){
//             for(Int_t fc=0;fc<2;fc++){
//               if(calLoGhost[tower][layer][fc]>0){
//                  cout<<tower<<" "<<layer<<" "                         
//                      <<fc<<" "<<calLoGhost[tower][layer][fc]<<endl;
//                 }
//               }
//             }
//           }
           
      } // loop on cluster
       
    } // end of if calLo trigger
  } // end of event loop

  outF->Write();
  outF->Close();
}