Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 09:15:47

0001 #include "iostream"
0002 #include "fstream"
0003 #include "vector"
0004 #include "map"
0005 #include "utility"
0006 #include "TString.h"
0007 #include "TFile.h"
0008 #include "TCanvas.h"
0009 #include "TTree.h"
0010 #include "TF1.h"
0011 #include "TH1D.h"
0012 #include "TObjArray.h"
0013 #include "TObjString.h"
0014 
0015 using namespace std;
0016 
0017 struct MyBoard {
0018   int boardID;
0019   int LG[64];
0020   int HG[64];
0021   float time;
0022   float x[64];
0023   float y[64];
0024   float z[64];
0025 };
0026 struct LRC {
0027   int layer;
0028   int row;
0029   int col;
0030 };
0031 
0032 map<pair<int, int>, LRC> LoadMapping(TString);
0033 
0034 int main(int argc, char* argv[]){
0035 
0036   TF1* ped=new TF1("ped","gaus",0,600);
0037   TString InputName=Form("%s",argv[1]);
0038   TString OutputName=Form("%s.root",TString(InputName(0,InputName.Length()-4)).Data());
0039   cout<<OutputName.Data()<<endl;
0040   TString MappingName=Form("%s",argv[2]);
0041   cout<<MappingName<<endl;
0042   TString CalibName;
0043   int ReadWrite=0;
0044   fstream CalibFile;
0045   if(argc>3){
0046     if(TString(Form("%s",argv[4])).Atoi()==1) ReadWrite=1;
0047     cout<<"Taking calibration from "<<argv[3]<<endl;
0048     CalibName=Form("%s",argv[3]);
0049     if(ReadWrite) CalibFile.open(CalibName.Data(),ios::out);
0050     else CalibFile.open(CalibName.Data(),ios::in);
0051   }
0052   //  return 0;
0053   TFile* output=new TFile(OutputName.Data(),"RECREATE");
0054   int Trg;
0055   float time[8];
0056   int HGcaen[8][64];
0057   int LGcaen[8][64];
0058   float xcaen[8][64];
0059   float ycaen[8][64];
0060   float zcaen[8][64];
0061   //TString asscaen[8][64];
0062   int HGlrc[64][2][4];
0063   int LGlrc[64][2][4];
0064   float xlrc[64][2][4];
0065   float ylrc[64][2][4];
0066   float zlrc[64][2][4];
0067   //TString asslrc[64][2][4];
0068   TTree* tree = new TTree("tree","tree");
0069   tree->Branch("TriggerID",&Trg,"TriggerID/I");
0070   tree->Branch("time",time,"time[8]/F");
0071   tree->Branch("HGcaen",HGcaen,"HGcaen[8][64]/I");
0072   tree->Branch("LGcaen",LGcaen,"LGcaen[8][64]/I");
0073   tree->Branch("xcaen",xcaen,"xcaen[8][64]/F");
0074   tree->Branch("ycaen",ycaen,"ycaen[8][64]/F");
0075   tree->Branch("zcaen",zcaen,"zcaen[8][64]/F");
0076   //tree->Branch("asscaen",asscaen)
0077   tree->Branch("HGlrc",HGlrc,"HGlrc[64][2][4]/I");
0078   tree->Branch("LGlrc",LGlrc,"LGlrc[64][2][4]/I");
0079   tree->Branch("xlrc",xlrc,"xlrc[64][2][4]/F");
0080   tree->Branch("ylrc",ylrc,"ylrc[64][2][4]/F");
0081   tree->Branch("zlrc",zlrc,"zlrc[64][2][4]/F");
0082   fstream ASCIIinput;
0083   //ASCIIinput.open("/home/vandrieu/Documents/ePIC/LFHCAL_BeamTest2024/CAEN_RO/Run410_list.txt",ios::in);
0084   ASCIIinput.open(Form("%s",InputName.Data()),ios::in);
0085   TString aline;
0086   TObjArray* tokens;
0087   map<int, vector<MyBoard> > mEvents;
0088   map<int, vector<MyBoard> >::iterator itevent;
0089   map<pair<int,int>, LRC> mappingConversion=LoadMapping(Form("%s",MappingName.Data()));
0090   TH1D* hcalibHG[8][64];
0091   TH1D* hcalibLG[8][64];
0092   for(int iboard=0; iboard<8; iboard++){
0093     for(int ichannel=0; ichannel<64; ichannel++){
0094       hcalibHG[iboard][ichannel]=new TH1D(Form("hcalibHG%d_%d",iboard,ichannel),Form("HG Board %d Channel %d spectrum; ADC [counts];",iboard,ichannel),4000,0,4000);
0095       hcalibLG[iboard][ichannel]=new TH1D(Form("hcalibLG%d_%d",iboard,ichannel),Form("LG Board %d Channel %d spectrum; ADC [counts];",iboard,ichannel),4000,0,4000);
0096     }
0097   }
0098   while(!ASCIIinput.eof()){
0099     //getline(ASCIIinput,aline);
0100     aline.ReadLine(ASCIIinput);
0101     if(!ASCIIinput.good()) break;
0102     if(aline[0]=='/') continue;
0103     tokens=aline.Tokenize(" \t");
0104     tokens->SetOwner(true);
0105     int Nfields=tokens->GetEntries();
0106     if(((TObjString*)tokens->At(0))->String()=="Brd") {
0107       tokens->Clear();
0108       delete tokens;
0109       continue;
0110     }
0111     if(Nfields!=7){
0112       cout<<"Unexpected number of fields"<<endl;
0113       cout<<"Expected 7, got: "<<Nfields<<", exit"<<endl;
0114       for(int j=0; j<Nfields;j++) cout<<j<<"  "<<((TObjString*)tokens->At(j))->String()<<endl;
0115       tokens->Clear();
0116       delete tokens;
0117       return -1;
0118     }
0119     int TriggerID=((TObjString*)tokens->At(5))->String().Atoi();
0120     MyBoard aBoard;
0121     int NHits=((TObjString*)tokens->At(6))->String().Atoi();
0122     aBoard.time=((TObjString*)tokens->At(4))->String().Atof();
0123     aBoard.boardID=((TObjString*)tokens->At(0))->String().Atoi();
0124     int channel=((TObjString*)tokens->At(1))->String().Atoi();
0125     aBoard.HG[channel]=((TObjString*)tokens->At(3))->String().Atoi();
0126     aBoard.LG[channel]=((TObjString*)tokens->At(2))->String().Atoi();
0127     for(int ich=1; ich<NHits;ich++){
0128           aline.ReadLine(ASCIIinput);
0129       tokens=aline.Tokenize(" ");
0130       tokens->SetOwner(true);
0131       Nfields=tokens->GetEntries();
0132       if(Nfields!=4){
0133         cout<<"Expecting 4 fields but read "<<Nfields<<endl;
0134         return -1;
0135       }
0136       channel=((TObjString*)tokens->At(1))->String().Atoi();
0137       aBoard.HG[channel]=((TObjString*)tokens->At(3))->String().Atoi();
0138       aBoard.LG[channel]=((TObjString*)tokens->At(2))->String().Atoi();
0139     }
0140     itevent=mEvents.find(TriggerID);
0141     if(itevent!=mEvents.end()){
0142       itevent->second.push_back(aBoard);
0143       if(itevent->second.size()==8){
0144     //Fill the tree the event is complete and erase from the map
0145     Trg=itevent->first;
0146     for(int iboard=0; iboard<8; iboard++){
0147       int brd=itevent->second[iboard].boardID;
0148       time[brd]=itevent->second[iboard].time;
0149       for(int ich=0; ich<64; ich++){
0150         HGcaen[brd][ich]=itevent->second[iboard].HG[ich];
0151         LGcaen[brd][ich]=itevent->second[iboard].LG[ich];
0152         hcalibHG[brd][ich]->Fill(HGcaen[brd][ich]);
0153         hcalibLG[brd][ich]->Fill(LGcaen[brd][ich]);
0154         HGlrc[mappingConversion[make_pair(brd,ich)].layer][mappingConversion[make_pair(brd,ich)].row][mappingConversion[make_pair(brd,ich)].col]=HGcaen[brd][ich];
0155         LGlrc[mappingConversion[make_pair(brd,ich)].layer][mappingConversion[make_pair(brd,ich)].row][mappingConversion[make_pair(brd,ich)].col]=LGcaen[brd][ich];
0156         if(ich%8>3){
0157           xcaen[brd][ich]=-7.5+5*(3-(ich%8)%4);
0158           ycaen[brd][ich]=-2.5;
0159           //HGlrc[7-ich/8+8*brd][1][(ich%8)%4]=itevent->second[brd].HG[ich];
0160           //LGlrc[7-ich/8+8*brd][1][(ich%8)%4]=itevent->second[brd].LG[ich];
0161         }
0162         else{
0163           xcaen[brd][ich]=7.5-5*(3-(ich%8)%4);
0164           ycaen[brd][ich]=2.5;
0165           //HGlrc[7-ich/8+8*brd][0][3-(ich%8)%4]=itevent->second[brd].HG[ich];
0166           //LGlrc[7-ich/8+8*brd][0][3-(ich%8)%4]=itevent->second[brd].LG[ich];
0167         }
0168         zcaen[brd][ich]=((63-ich)/8)*(2.)+16*brd;
0169         xlrc[mappingConversion[make_pair(brd,ich)].layer][mappingConversion[make_pair(brd,ich)].row][mappingConversion[make_pair(brd,ich)].col]=xcaen[brd][ich];
0170         ylrc[mappingConversion[make_pair(brd,ich)].layer][mappingConversion[make_pair(brd,ich)].row][mappingConversion[make_pair(brd,ich)].col]=ycaen[brd][ich];
0171         zlrc[mappingConversion[make_pair(brd,ich)].layer][mappingConversion[make_pair(brd,ich)].row][mappingConversion[make_pair(brd,ich)].col]=zcaen[brd][ich];
0172       }
0173     }
0174     tree->Fill();
0175     mEvents.erase(itevent);
0176       }
0177     }
0178     else{//This is a new triggerID, ie new event
0179       vector<MyBoard> vboard;
0180       vboard.push_back(aBoard);
0181       mEvents[TriggerID]=vboard;  
0182     }
0183   }
0184   TCanvas* c = new TCanvas("c","c",800,400);
0185   c->Divide(2,1);
0186   c->cd(1)->SetLogy();
0187   c->cd(2)->SetLogy();
0188   if(CalibFile.is_open()&&(ReadWrite==1)){
0189     c->Print("TestDistributions.pdf[");
0190     for(int iboard=0; iboard<8; iboard++){
0191       for(int ichannel=0; ichannel<64; ichannel++){
0192     ped->SetParameters(1000.,70.,15.);
0193     ped->SetParLimits(1,1,500);
0194     ped->SetParLimits(2,1,60);
0195     c->cd(1);
0196     hcalibHG[iboard][ichannel]->GetXaxis()->SetRangeUser(0,600);
0197     hcalibHG[iboard][ichannel]->Draw();
0198     hcalibHG[iboard][ichannel]->Fit(ped);
0199     CalibFile<<iboard<<"\t"<<ichannel<<"\t"<<ped->GetParameter(1)<<"\t"<<ped->GetParameter(2);
0200     c->cd(2);
0201     hcalibLG[iboard][ichannel]->GetXaxis()->SetRangeUser(0,600);
0202     hcalibLG[iboard][ichannel]->Draw();
0203     hcalibLG[iboard][ichannel]->Fit(ped);
0204     CalibFile<<ped->GetParameter(1)<<"\t"<<ped->GetParameter(2)<<endl;
0205     c->Print("TestDistributions.pdf");
0206       }
0207     }
0208     c->Print("TestDistributions.pdf]");
0209     CalibFile.close();
0210   }
0211 
0212   cout<<"Number of partial events: "<<mEvents.size()<<endl;
0213   tree->Write();
0214   output->Close();
0215   return 0;
0216 }
0217 
0218 
0219 map<pair<int, int>, LRC> LoadMapping(TString name){
0220   cout<<"Loading the map from "<<name.Data()<<endl;
0221   map<pair<int, int>, LRC> m;
0222   fstream config;
0223   config.open(name.Data(),ios::in);
0224   if(!config.is_open()){
0225     cout<<"could not open the mapping file input"<<endl;
0226   }
0227   int board, channel, layer, chInLayer, row, column;
0228   string assemblyLayer;
0229   pair<int,int> mykey;
0230   LRC myvalue;
0231   while(!config.eof()){
0232     config>>board>>channel>>layer>>assemblyLayer>>chInLayer>>row>>column;
0233     if(!config.good()){
0234       break;
0235     }
0236     mykey.first=board;
0237     mykey.second=channel;
0238     myvalue.layer=layer;
0239     myvalue.row=row;
0240     myvalue.col=column;
0241     m[mykey]=myvalue;
0242   }
0243   return m;
0244 }