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
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
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
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
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
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
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
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
0160
0161 }
0162 else{
0163 xcaen[brd][ich]=7.5-5*(3-(ich%8)%4);
0164 ycaen[brd][ich]=2.5;
0165
0166
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{
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 }