1 /*----------------------
2 GATE version name: gate_v...
4 Copyright (C): OpenGATE Collaboration
6 This software is distributed under the terms
7 of the GNU Lesser General Public Licence (LGPL)
8 See GATE/LICENSE.txt for further details
9 ----------------------*/
17 #include <TIterator.h>
34 #include "GateMergeManager.hh"
38 void GateMergeManager::StartMerging(string splitfileName){
40 // get the files to merge
41 ReadSplitFile(splitfileName);
43 if (m_fastMerge==true) FastMergeRoot();
46 //if we are here the merging has been successful
47 //we mark the directory as ready for cleanup
48 //string dir = splitfileName.substr(0,splitfileName.rfind("/"));
49 //string ready = "touch "+dir+"/ready_for_delete";
50 //const int res=system(ready.c_str());
51 //if(res) cout<<"Strange?? Can't mark "<<dir<<" as done!"<<endl;
54 void GateMergeManager::StartMergingFromFilenames(Strings filenames, string outputfile)
56 for (Strings::const_iterator iter=filenames.begin(); iter!=filenames.end(); iter++)
58 m_vRootFileNames.push_back(*iter);
59 if(m_verboseLevel>2) cout<<"Root input file name: "<<m_vRootFileNames.back()<<endl;
62 m_Nfiles = m_vRootFileNames.size();
63 m_RootTargetName=m_outDir+outputfile;
64 if(m_verboseLevel>2) cout<<"Root output file name: "<<m_RootTargetName<<endl;
66 if (m_fastMerge==true) FastMergeRoot();
70 // to process the splitfile
71 void GateMergeManager::ReadSplitFile(string splitfileName){
74 splitfile.open(splitfileName.c_str());
76 cout<<"Can't open split file: "<<splitfileName<<endl;
80 //avoid multiple calls
81 if(m_vRootFileNames.size()>0) return;
83 stringstream ssnfiles;
85 while(splitfile){ //we allow for rubbish before "Number of files:"
86 splitfile.getline(cline,512);
87 if(!strncmp(cline,"Number of files:",16)){
88 ssnfiles.str(cline+16);
89 ssnfiles>>m_Nfiles; // the number of files to merge
94 // now we look for the file names
97 splitfile.getline(cline,512);
100 if(!strncmp(cline,"Root filename:",14)){
101 m_vRootFileNames.push_back(strtok(cline+15," "));
102 m_vRootFileNames[iRoot]+=".root";
103 if(m_verboseLevel>2) cout<<"Root input file name: "<<m_vRootFileNames[iRoot]<<endl;
107 else if(!strncmp(cline,"Original Root filename:",23)){
108 m_RootTargetName=strtok(cline+23," ");
109 m_RootTargetName+=".root";
110 if(m_outDir!=""){// output directory option used
111 //replace path with outDir
112 size_t pos=m_RootTargetName.rfind('/',m_RootTargetName.length());
113 if(pos==string::npos) pos=-1;
114 m_RootTargetName=m_outDir+m_RootTargetName.substr(pos+1);
116 if(m_verboseLevel>2) cout<<"Root output file name: "<<m_RootTargetName<<endl;
120 // check if number of root files correct
121 if(iRoot && iRoot!=m_Nfiles) {
122 cout<<"Inconsistent number of root file entries in split file!"<<endl;
127 /************************************************************************************/
128 void GateMergeManager::FastMergeRoot()
131 if(m_RootTargetName=="") {
132 cout<<"No root output file name given in split file"<<endl;
136 // number of files to merge
137 int nfiles=m_vRootFileNames.size();
139 cout<<"No files to merge!"<<endl;
142 //we try to recover the last_event_ID in all root files
143 for(int i=0;i<nfiles;i++) {
144 filearr.push_back(TFile::Open(m_vRootFileNames[i].c_str(),"OLD"));
145 if(filearr[i]==NULL){
146 cout<<"Not a readable file "<<m_vRootFileNames[i]<<" - exit!"<<endl;
150 m_lastEvents.resize(nfiles+1);
152 bool flgLstEvntID(false);
154 // latest_event_ID histogram
157 for(int i=0;i<nfiles;i++){
158 h = (TH1*)filearr[i]->Get("latest_event_ID");
159 m_lastEvents[i+1]=int(h->GetMean());
163 cout<<"FastMerge: No latest_event_ID histogram - exit!"<<endl;
167 /*for (int i=0;i<nfiles;i++)
169 cout<<"Last ID is: "<<m_lastEvents[i+1]<<" File " <<i+1<<endl;
173 vector<string> treeNames;
174 TFile* node = TFile::Open(m_vRootFileNames[0].c_str(),"OLD");
175 TIter nextkey(node->GetListOfKeys());
177 while ((key = (TKey*)nextkey()))
180 TObject* obj = (TObject*)key->ReadObj();
181 if( obj->IsA()->InheritsFrom("TTree"))
184 bool singleName=true;
185 for(unsigned int i=0;i<treeNames.size();i++)
187 if(treeNames[i].compare(obj->GetName())==0)
193 if(singleName) treeNames.push_back(obj->GetName());
196 //take care of all trees
197 for(unsigned int i=0;i<treeNames.size();i++)
199 if(!MergeTree(treeNames[i])) if(m_verboseLevel>1) cout<<"Problem with merging "<<treeNames[i]<<endl;
202 /************************************************************************************/
203 void GateMergeManager::MergeRoot(){
206 if(m_RootTargetName=="") {
207 cout<<"No root output file name given in split file"<<endl;
211 // number of files to merge
212 int nfiles=m_vRootFileNames.size();
214 // the root output file
216 m_RootTarget = TFile::Open(m_RootTargetName.c_str(),"RECREATE");
218 m_RootTarget = TFile::Open(m_RootTargetName.c_str(),"NEW");
220 cout<<"The root ouput file already exists! Try -f to overwrite and erase the file."<<endl;
226 cout<<"No files to merge!"<<endl;
230 TFile* filearr[nfiles];
231 for(int i=1;i<nfiles;i++) {
232 filearr[i] = TFile::Open(m_vRootFileNames[i].c_str(),"OLD");
233 if(filearr[i]==NULL){
234 cout<<"Not a readable file "<<m_vRootFileNames[i]<<" - exit!"<<endl;
239 // first we copy all histos (only top directory)
240 // and look for the latestEventID
241 // and find out which trees/ntuples exist
242 m_lastEvents.resize(nfiles+1);
245 bool flgLstEvntID(false);
247 if(m_verboseLevel>0) {
249 for(int i=0;i<nfiles;i++) cout<<m_vRootFileNames[i]<<" ";
250 cout<<"-> "<<m_RootTarget->GetName()<<endl;
253 vector<string> treeNames;
254 TFile* node = TFile::Open(m_vRootFileNames[0].c_str(),"OLD");
256 cout<<"Not a readable file "<<m_vRootFileNames[0]<<" - exit!"<<endl;
259 TIter nextkey(node->GetListOfKeys());
261 while ((key = (TKey*)nextkey())) {
263 TObject* obj = (TObject*)key->ReadObj();
264 if( obj->IsA()->InheritsFrom("TH1")){
266 TH1 *h1 = (TH1 *)obj->Clone();
267 if(strcmp(h1->GetName(),"latest_event_ID")){
268 //any other histogram
269 for(int i=1;i<nfiles;i++){
270 TH1 *h2 = (TH1*)filearr[i]->Get( h1->GetName() );
275 // latest_event_ID histogram
277 m_lastEvents[1]=int(h1->GetMean());
279 for(int i=1;i<nfiles;i++){
280 h = (TH1*)filearr[i]->Get("latest_event_ID");
281 m_lastEvents[i+1]=int(h->GetMean());
283 if(h!=NULL) h->Write(); // we keep the highest latest_event_ID to allow for successive merging
287 if( obj->IsA()->InheritsFrom("TH2")){
289 TH2 *h1 = (TH2 *)obj->Clone();
290 for(int i=1;i<nfiles;i++){
291 TH2 *h2 = (TH2*)filearr[i]->Get( h1->GetName() );
296 if( obj->IsA()->InheritsFrom("TTree")){
298 bool singleName=true;
299 for(unsigned int i=0;i<treeNames.size();i++) {
300 if(treeNames[i].compare(obj->GetName())==0){
305 if(singleName) treeNames.push_back(obj->GetName());
309 cout<<"Merge: No latest_event_ID histogram - exit!"<<endl;
313 //now we take care of the trees
314 for(unsigned int i=0;i<treeNames.size();i++)
315 if(!MergeTree(treeNames[i])) if(m_verboseLevel>1) cout<<"Problem with merging "<<treeNames[i]<<endl;
317 // everything is done
320 /*******************************************************************************************/
321 // cleanup after merging
322 void GateMergeManager::StartCleaning(string splitfileName,bool test){
324 // get the files to erase
325 ReadSplitFile(splitfileName);
327 string dir=splitfileName.substr(0,splitfileName.rfind("/"));
329 // test for the mark: ready_for_delete
330 string touched=dir+"/ready_for_delete";
332 ready.open(touched.c_str());
334 cout<<"Cannot do the cleanup - directory "<<dir<<" not marked as ready!"<<endl;
339 if(test) cout<<"I would like to erase:"<<endl;
340 if(!test&&m_verboseLevel>1) cout<<"Going to erase the following files:"<<endl;
341 if(m_verboseLevel>1||test){
342 cout<<" --> "<<dir<<endl;
344 for(unsigned int i=0;i<m_vRootFileNames.size();i++){
345 cout<<" --> "<<m_vRootFileNames[i]<<endl;
351 for(unsigned int i=0;i<m_vRootFileNames.size();i++){
352 const string rmfiles("rm -f "+m_vRootFileNames[i]);
353 const int res = system(rmfiles.c_str());
355 cout<<"Could not remove "<<m_vRootFileNames[i]<<endl;
356 cout<<"Please remove it manually!"<<endl;
359 const string rmfiles("rm -f "+dir+"/*");
360 system(rmfiles.c_str());
361 const string rmdir="rm -f -r "+dir;
362 const int res2 = system(rmdir.c_str());
364 cout<<"Could not remove "<<dir<<endl;
365 cout<<"Please remove it manually!"<<endl;
370 /*******************************************************************************************/
371 //find out what kind of tree we have and call the right merger
372 bool GateMergeManager::MergeTree(string chainName){
373 if (m_fastMerge==false)
375 TChain* chain = new TChain(chainName.c_str());
377 // number of files to merge
378 int nfiles=m_vRootFileNames.size();
380 for(int i=0;i<nfiles;i++) chain->Add(m_vRootFileNames[i].c_str());
381 int nentries=chain->GetEntries();
383 if(m_verboseLevel>1) cout<<chain->GetName()<<" is empty!"<<endl;
387 if(chainName=="Gate") MergeGate(chain);
389 if(chain->FindBranch("eventID1")!=NULL) {
390 if( (chain->FindBranch("runID")==NULL)
391 ||(chain->FindBranch("time1")==NULL)
392 ||(chain->FindBranch("time2")==NULL) ) {
393 cout<<"Cannot find one of: runID, time1, time2 in "<<chain->GetName()<<endl;
397 return MergeCoin(chain);
399 } else if( (chain->FindBranch("runID")==NULL)||(chain->FindBranch("eventID")==NULL)
400 ||(chain->FindBranch("time")==NULL) )
402 cout<<"Cannot find one of: runID, eventID, time in "<<chain->GetName()<<endl;
406 return MergeSing(chain);
410 if (chainName.find("Hits",0)!=string::npos) return FastMergeSing(chainName);
411 if (chainName.find("Singles",0)!=string::npos) return FastMergeSing(chainName);
412 if (chainName.find("Gate",0)!=string::npos) return FastMergeGate(chainName);
413 if (chainName.find("Coincidences",0)!=string::npos) return FastMergeCoin(chainName);
418 /*******************************************************************************************/
419 bool GateMergeManager::FastMergeGate(string name)
421 //create the output file
426 //float maxtime =-999999999;
427 string clusterName=name+"_cluster";
429 for(int j=0;j<m_Nfiles;j++)
431 filearr[j]->ReOpen("UPDATE");
432 TTree *oldTree = (TTree*)gDirectory->Get(name.c_str());
433 TBranch *branch = oldTree->GetBranch("event");
434 branch->SetAddress(&event);
435 TBranch* newBranch=oldTree->Branch("eventcluster",&event,"event");
437 Int_t nentries=(Int_t)oldTree->GetEntries();
440 else offset+=m_lastEvents[currentfile];
443 for(int i=0;i<nentries;i++)
454 if((i==nentries-1) && (j==m_Nfiles-1))
459 else if (i!=nentries-1)
474 /*******************************************************************************************/
475 bool GateMergeManager::FastMergeSing(string name)
482 string clusterName=name+"_cluster";
484 for(int j=0;j<m_Nfiles;j++)
486 filearr[j]->ReOpen("UPDATE");
487 TTree *oldTree = (TTree*)gDirectory->Get(name.c_str());
488 TBranch *branch = oldTree->GetBranch("eventID");
489 branch->SetAddress(&eventID);
490 TBranch *branch2 = oldTree->GetBranch("runID");
491 branch2->SetAddress(&runID);
492 TBranch* newBranch=oldTree->Branch("eventIDcluster",&eventID,"eventID/I");
494 Int_t nentries=(Int_t)oldTree->GetEntries();
497 else offset+=m_lastEvents[currentfile];
500 for(int i=0;i<nentries;i++)
516 /*******************************************************************************************/
517 bool GateMergeManager::FastMergeCoin(string name)
519 //create the output file
526 string clusterName=name+"_cluster";
527 cout<<"starting coincidence merging..."<<endl;
528 for(int j=0;j<m_Nfiles;j++)
530 cout<<"working on file..."<<j<<endl;
531 filearr[j]->ReOpen("UPDATE");
532 TTree *oldTree = (TTree*)gDirectory->Get(name.c_str());
533 TBranch *branch1 = oldTree->GetBranch("eventID1");
534 branch1->SetAddress(&eventID1);
535 TBranch *branch2 = oldTree->GetBranch("eventID2");
536 branch2->SetAddress(&eventID2);
537 TBranch *branch3 = oldTree->GetBranch("runID");
538 branch3->SetAddress(&runID);
540 TBranch* newBranch1=oldTree->Branch("eventID1cluster",&eventID1,"eventID1/I");
541 TBranch* newBranch2=oldTree->Branch("eventID2cluster",&eventID2,"eventID2/I");
543 Int_t nentries=(Int_t)oldTree->GetEntries();
546 else offset+=m_lastEvents[currentfile];
547 cout<<"the offset in this file is: "<<offset<<endl;
550 for(int i=0;i<nentries;i++){
551 branch1->GetEvent(i);
552 branch2->GetEvent(i);
553 branch3->GetEvent(i);
569 /*******************************************************************************************/
571 bool GateMergeManager::MergeGate(TChain* chainG) {
573 int nentries=chainG->GetEntries();
576 chainG->SetBranchAddress("event",&event);
578 chainG->SetBranchAddress("iontime",&iontime);
580 //create the new tree
581 //Allow for large files and do not write every bit separately
583 TTree * newTree = chainG->CloneTree(0);
584 newTree->SetAutoSave(2000000000);
585 if(m_maxRoot!=0) newTree->SetMaxTreeSize(m_maxRoot);
586 else newTree->SetMaxTreeSize(17179869184LL);
588 // changing the compression level everywhere
590 TIter next(newTree->GetListOfBranches());
591 while ((br=(TBranch*)next())) br->SetCompressionLevel(m_CompLevel);
593 // the main part: modify the eventID
597 float maxtime =-999999999;
598 for(int i=0;i<nentries;i++){
599 if (chainG->GetEntry(i) <= 0) break; //end of chain
600 if(chainG->GetTreeNumber()>currentTree) {
602 offset+=m_lastEvents[currentTree];
603 // check for overlaping time intervalls between different files
604 // (not within the same file i.e. no time order assumed)
607 cout<<"Warning - overlapping Gate iontime ("
608 <<iontime<<") in file: "<<m_vRootFileNames[currentTree].c_str()<<endl;
610 // run end is marked by repeating event
611 if(lastEvent!=event) {
613 // if (fpclassify(iontime)!=FP_NAN && fpclassify(iontime)!=FP_INFINITE) {
615 if(maxtime<iontime) maxtime=iontime;
616 // the offset to get a unique event numbering
620 cout<<"Warning - inf or NaN in Gate tree for iontime! "<<m_vRootFileNames[currentTree].c_str()<<endl;
626 if(chainG->GetEntry(i+1) <= 0) { //chain end fill the double event
629 } else if(chainG->GetTreeNumber()==currentTree //no new file or
630 ||(chainG->GetTreeNumber()>currentTree
631 &&m_lastEvents[chainG->GetTreeNumber()]==0) ) { //new file with no offset fill double event
635 lastEvent=0;// we may have triple 1 (1 event run)
645 /*******************************************************************************************/
646 // Singles and Hits tree merger
647 bool GateMergeManager::MergeSing(TChain* chainS){
649 int nentriesS=chainS->GetEntries();
655 chainS->SetBranchAddress("eventID",&eventID);
656 chainS->SetBranchAddress("runID",&runID);
657 chainS->SetBranchAddress("time",&time);
660 TTree * newSing = chainS->CloneTree(0);
661 newSing->SetAutoSave(2000000000);
662 if(m_maxRoot!=0) newSing->SetMaxTreeSize(m_maxRoot);
663 else newSing->SetMaxTreeSize(17179869184LL);
665 // changing CompLevel everywhere
667 TIter next(newSing->GetListOfBranches());
668 while ((br=(TBranch*)next())) br->SetCompressionLevel(m_CompLevel);
670 //the main part: changing eventID
674 float maxtime =-999999999;
675 for(int i=0;i<nentriesS;i++){
676 if(chainS->GetEntry(i) <= 0) break; //end of chain
677 if(chainS->GetTreeNumber()>currentTree) {
679 offset+=m_lastEvents[currentTree];
680 // check for overlaping time intervalls between different files
681 // (not within the same file i.e. no time order assumed)
684 cout<<"Warning - overlapping Singles time ("
685 <<time<<") in file: "<<m_vRootFileNames[currentTree].c_str()<<endl;
687 // the offset to get a unique event numbering
691 offset=0; //new run in file we must not change the eventID anymore
695 if(maxtime<time)maxtime=time;
702 /*******************************************************************************************/
703 // Coincidences tree merger
704 bool GateMergeManager::MergeCoin(TChain* chainC){
706 int nentriesC=chainC->GetEntries();
709 chainC->SetBranchAddress("eventID1",&eventID1);
711 chainC->SetBranchAddress("time1",&time1);
713 chainC->SetBranchAddress("eventID2",&eventID2);
715 chainC->SetBranchAddress("time2",&time2);
717 chainC->SetBranchAddress("runID",&runID);
720 TTree * newCoin = chainC->CloneTree(0);
721 newCoin->SetAutoSave(2000000000);
722 if(m_maxRoot!=0) newCoin->SetMaxTreeSize(m_maxRoot);
723 else newCoin->SetMaxTreeSize(17179869184LL);
725 // changing CompLevel everywhere
727 TIter next(newCoin->GetListOfBranches());
728 while ((br=(TBranch*)next())) br->SetCompressionLevel(m_CompLevel);
733 float maxtime =-999999999;
734 for(int i=0;i<nentriesC;i++){
735 if(chainC->GetEntry(i) <= 0) break; //end of chain
736 if(chainC->GetTreeNumber()>currentTree) {
738 offset+=m_lastEvents[currentTree];
739 // check for overlaping time intervalls between different files
740 // (not within the same file i.e. no time order assumed)
743 cout<<"Warning - overlapping Coincidences time1 ("
744 <<time1<<") in file: "<<m_vRootFileNames[currentTree].c_str()<<endl;
747 cout<<"Warning - overlapping Coincidences time2 ("
748 <<time2<<") in file: "<<m_vRootFileNames[currentTree].c_str()<<endl;
750 // the offset to get a unique event numbering
754 offset=0; //new run in file we must not change the eventID anymore
759 if(maxtime<time1)maxtime=time1;
760 if(maxtime<time2)maxtime=time2;
766 /*******************************************************************************************/