Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:30:02

0001 static const int lutIndxMax1 = 4294967295;  //32 bit integer
0002 static const int lutIndxMax2 = 4294967295;  //32 bit integer
0003 static const int binningIndxMax = 9;
0004 static const int maxIndxMax  = 9;
0005 static const int partLutIndxMax = 64;
0006 
0007 TH1I* lutPartIDLut = new TH1I("lutPartIDLut","lutPartIDLut",partLutIndxMax,0,partLutIndxMax);
0008 TH1I* binInfo = new TH1I("binInfo","binInfo",binningIndxMax,0,binningIndxMax);
0009 TH1I* binFactor = new TH1I("binFactor","binFactor",maxIndxMax,0,maxIndxMax);
0010 TH1I* binStep = new TH1I("binStep","binStep",maxIndxMax,0,maxIndxMax);
0011 TH1I* maxBitsInfo = new TH1I("maxBitsInfo","maxBitsInfo",maxIndxMax,0,maxIndxMax);
0012 TH2I* lutMatrixEAverage =new TH2I("lutMatrixEAverage","lutMatrixEAverage",lutIndxMax1,0,lutIndxMax1,lutIndxMax2,0,lutIndxMax2);
0013 TH2I* lutMatrixESigma   =new TH2I("lutMatrixESigma","lutMatrixESigma",lutIndxMax1,0,lutIndxMax1,lutIndxMax2,0,lutIndxMax2);
0014 TH2I* lutMatrixEDist    =new TH2I("lutMatrixEDist","lutMatrixEDist",lutIndxMax1,0,lutIndxMax1,lutIndxMax2,0,lutIndxMax2);
0015 
0016 
0017 int zdcLutTableGen(){
0018   
0019   TFile f("zdcLutTable.root","RECREATE");
0020   fillBinInfo();
0021   fillMaxBitsInfo();
0022   fillBinFactor();
0023   fillBinStep();
0024   fillParticleTableInfo();
0025   //fillLut();
0026   binInfo->Write();
0027   maxBitsInfo->Write();
0028   lutPartIDLut->Write();
0029   lutMatrixEAverage->Write();
0030   lutMatrixESigma->Write();
0031   lutMatrixEDist->Write();
0032   f.Close();
0033   return 0;
0034 }
0035 
0036 
0037 int fillLut(){
0038   double dE =0;
0039   double dEsigma =0;
0040   double dEdist =0;
0041   for(int ienergy= 0; ienergy <= int(maxBitsInfo->GetBinContent(1)*binFactor->GetBinContent(1)); ienergy+=binStep->GetBinContent(1))
0042     for(int itheta=  -int(maxBitsInfo->GetBinContent(2)*binFactor->GetBinContent(2)); 
0043     itheta <= int(maxBitsInfo->GetBinContent(2)*binFactor->GetBinContent(2)); itheta+=binStep->GetBinContent(2))
0044       for(int iphi =  - int(maxBitsInfo->GetBinContent(3)*binFactor->GetBinContent(3)); 
0045       iphi <= int(maxBitsInfo->GetBinContent(3)*binFactor->GetBinContent(3)); iphi+=binStep->GetBinContent(3))
0046     for(int iside= 0; iside <= int(maxBitsInfo->GetBinContent(4)*binFactor->GetBinContent(4)); iside+=binStep->GetBinContent(4))
0047       for(int isection= 0; isection <= int(maxBitsInfo->GetBinContent(5)*binFactor->GetBinContent(5)); isection+=binStep->GetBinContent(5))
0048         for(int channel = 0; channel <= int(maxBitsInfo->GetBinContent(6)*binFactor->GetBinContent(6)); channel+=binStep->GetBinContent(6))
0049           for(int ixin=  - int(maxBitsInfo->GetBinContent(7)*binFactor->GetBinContent(7)); 
0050           ixin <= int(maxBitsInfo->GetBinContent(7)*binFactor->GetBinContent(7)); ixin+=binStep->GetBinContent(7))
0051         for(int iyin=  - int(maxBitsInfo->GetBinContent(8)*binFactor->GetBinContent(8)); 
0052             iyin <= int(maxBitsInfo->GetBinContent(8)*binFactor->GetBinContent(8)); iyin+=binStep->GetBinContent(8))
0053           for(int izin= - int(maxBitsInfo->GetBinContent(9)*binFactor->GetBinContent(9)); 
0054               izin <= int(maxBitsInfo->GetBinContent(9)*binFactor->GetBinContent(9)); izin+=binStep->GetBinContent(9))
0055             for(int iparCode= 0; iparCode <= int(maxBitsInfo->GetBinContent(10)*binFactor->GetBinContent(10));iparCode+=binStep->GetBinContent(10)){
0056               int partID = lutPartIDLut->GetBinContent(iparCode);
0057               std::cout<<ienergy<<" "
0058                    <<itheta<<" "
0059                    <<iphi<<" "
0060                    <<iside<<" "
0061                    <<isection<<" "
0062                    <<channel<<" "
0063                    <<ixin<<" "
0064                    <<iyin<<" "
0065                    <<izin<<" "
0066                    <<iparCode<<" "
0067                    <<std::endl;
0068               dE = 11111.000;
0069               dEsigma = 11.00;
0070               dEdist = 1.0;
0071               long int iLutIndex1 = encode1(iphi,itheta,ixin,iyin,izin);
0072               long int iLutIndex2 = encode2(ienergy,isection,iside,channel,iparCode);
0073               lutMatrixEAverage->SetBinContent(iLutIndex1,iLutIndex1,dE);
0074               lutMatrixESigma->SetBinContent(iLutIndex1,iLutIndex1,dEsigma);
0075               lutMatrixEDist->SetBinContent(iLutIndex1,iLutIndex1,dEdist);      
0076             }
0077   return 0;
0078 }
0079 
0080 int fillMaxBitsInfo(){
0081   // maxBitsInfo*binFactor = max value of variable
0082   maxBitsInfo->SetBinContent(1,512); // energy
0083   maxBitsInfo->SetBinContent(2,64);  // theta
0084   maxBitsInfo->SetBinContent(3,64);  // phi
0085   maxBitsInfo->SetBinContent(4,2);   // detector side 
0086   maxBitsInfo->SetBinContent(5,4);   // detector section 
0087   maxBitsInfo->SetBinContent(6,8);   // detector channel
0088   maxBitsInfo->SetBinContent(7,16);  // X
0089   maxBitsInfo->SetBinContent(8,16);  // Y
0090   maxBitsInfo->SetBinContent(9,32);  // Z
0091   maxBitsInfo->SetBinContent(10,64);  // pid
0092   return 0;
0093   }
0094 
0095 int fillBinInfo(){
0096   binInfo->SetBinContent(1,50000); // 20 GeV bins 50,000 -- 100 GeV 10,000 -- 500 GeV 2,000, -- 5000 GeV 200
0097   binInfo->SetBinContent(2,64);  // 3 degree theta binning 
0098   binInfo->SetBinContent(3,64);  // 3 degree phi binning 
0099   binInfo->SetBinContent(4,1);  // side
0100   binInfo->SetBinContent(5,1);  // section 
0101   binInfo->SetBinContent(6,1);  // channel
0102   binInfo->SetBinContent(7,10); // x 1 cm bin
0103   binInfo->SetBinContent(8,10); // y 1 cm bin
0104   binInfo->SetBinContent(9,40); // z 4 cm bin
0105   binInfo->SetBinContent(10,1); // 1 bin part ID
0106   return 0;
0107 }
0108 
0109 int fillBinFactor(){
0110   binFactor->SetBinContent(1,20); // 20 GeV bins 50000, 100 GeV 10,000 , etc.
0111   binFactor->SetBinContent(2,2);  // 3 degree theta binning 
0112   binFactor->SetBinContent(3,3);  // 3 degree phi binning 
0113   binFactor->SetBinContent(4,1);  // side
0114   binFactor->SetBinContent(5,1);  // section 
0115   binFactor->SetBinContent(6,1);  // channel
0116   binFactor->SetBinContent(7,1); // x 1 cm bin
0117   binFactor->SetBinContent(8,1); // y 1 cm bin
0118   binFactor->SetBinContent(9,4); // z 4 cm bin
0119   binFactor->SetBinContent(10,1); // 1 bin part ID
0120   return 0;
0121 }
0122 
0123 int fillBinStep(){
0124   binStep->SetBinContent(1,195); // 20 GeV bins, 100 GeV 10,000 , 500, etc.
0125   binStep->SetBinContent(2,64);  // 3 degree theta binning 
0126   binStep->SetBinContent(3,64);  // 3 degree phi binning 
0127   binStep->SetBinContent(4,1);  // side
0128   binStep->SetBinContent(5,1);  // section 
0129   binStep->SetBinContent(6,1);  // channel
0130   binStep->SetBinContent(7,1); // x 1 cm bin
0131   binStep->SetBinContent(8,1); // y 1 cm bin
0132   binStep->SetBinContent(9,4); // z 4 cm bin
0133   binStep->SetBinContent(10,1); // 1 bin part ID
0134   return 0;
0135 }
0136 
0137 int fillParticleTableInfo(){
0138   //initialize the array with a non existing pid
0139   // fill by hand according to PGD convension in cmssw:
0140   // http://cmslxr.fnal.gov/lxr/source/SimGeneral/HepPDTESSource/data/particle.tbl?v=CMSSW_2_0_6
0141   lutPartIDLut->SetBinContent(1,11);  // e^-
0142   lutPartIDLut->SetBinContent(2,-11); // e^+
0143   lutPartIDLut->SetBinContent(3,13);   // mu^-
0144   lutPartIDLut->SetBinContent(4,-13); //mu^+
0145   lutPartIDLut->SetBinContent(5,21);  //g
0146   lutPartIDLut->SetBinContent(6,22); //gamma
0147   lutPartIDLut->SetBinContent(7,111); //pi^0
0148   lutPartIDLut->SetBinContent(8,211); //pi^+
0149   lutPartIDLut->SetBinContent(9,-211); // pi^-
0150   lutPartIDLut->SetBinContent(10,311); //K^0             
0151   lutPartIDLut->SetBinContent(11,321); //K^+           
0152   lutPartIDLut->SetBinContent(12,-321); //K^-          
0153   lutPartIDLut->SetBinContent(13,2112); //n^0 
0154   lutPartIDLut->SetBinContent(14,-2112); //n~^0
0155   lutPartIDLut->SetBinContent(15,2212); //p^+  
0156   lutPartIDLut->SetBinContent(16,-2212); //p~^-
0157   lutPartIDLut->SetBinContent(17,1000010020); //Deuterium    
0158   lutPartIDLut->SetBinContent(18,1000010030); //Tritium      
0159   lutPartIDLut->SetBinContent(19,1000020030); //He3          
0160   lutPartIDLut->SetBinContent(20,100002004); //Alpha-(He4)  
0161   lutPartIDLut->SetBinContent(partLutIndxMax,666);             
0162   return 0;
0163 }
0164 
0165 void decode1(const unsigned long & lutidx,int& iphi, int& itheta, int& ix,int& iy, int& iz){
0166   int iphisgn = (lutidx>>29)&1;
0167   int ithsgn  = (lutidx>>28)&1;
0168   int izsgn   = (lutidx>>27)&1;
0169   int iysgn   = (lutidx>>26)&1;
0170   int ixsgn   = (lutidx>>25)&1;
0171   itheta = (lutidx>>19)&63;
0172   iphi = (lutidx>>13)&63;
0173   iz = (lutidx>>8)&31;
0174   iy = (lutidx>>4)&15;
0175   ix = (lutidx)&15;
0176 
0177   if(ithsgn == 0)itheta*= -1;
0178   if(iphisgn == 0)iphi*= -1;
0179   if(izsgn == 0)iz*= -1;
0180   if(iysgn == 0)iy*= -1;
0181   if(ixsgn == 0)ix*= -1;
0182   return;
0183 }
0184 
0185 void decode2(const unsigned long & lutidx,int& ien, int& isec, int& isid, int& icha, int& iparID){
0186   ien = (lutidx>>12)&511;
0187   iparID = (lutidx>>6)&63;
0188   icha = (lutidx>>3)&7;
0189   isec = (lutidx>>1)&3;
0190   isid = 1 +(lutidx&1);
0191   return;
0192 }
0193 
0194 unsigned long encode1(int iphi, int itheta, int ix, int iy, int iz){
0195   int ixsgn = 1;
0196   if(ix<0){
0197     ix = -ix;
0198     ixsgn = 0;
0199   }
0200   int iysgn = 1;
0201   if(iy<0){
0202     iy = -iy;
0203     iysgn = 0;
0204   }
0205   int izsgn = 1;  
0206   if(iz<0){
0207     iz = -iz;
0208     izsgn = 0;
0209   }
0210   int ithsgn = 1;
0211   if(itheta<0){
0212     itheta = -itheta;
0213     ithsgn = 0;
0214   }
0215   int iphsgn = 1;
0216   if(iphi<0){
0217     iphi = -iphi;
0218     iphsgn = 0;
0219   }
0220 
0221   unsigned long lutindex = (iphsgn&1)<<29;
0222   lutindex += (ithsgn&1) <<28;
0223   lutindex += (izsgn&1)  <<27;
0224   lutindex += (iysgn&1)  <<26;
0225   lutindex += (ixsgn&1)  <<25;    //bits 25
0226   lutindex += (itheta&63)<<19;    //bits 19-24
0227   lutindex += (iphi&63)  <<13;    //bits 13-18
0228   lutindex += (iz&31)    <<8;     //bits  8-12
0229   lutindex += (iy&15)    <<4;     //bits  4- 7
0230   lutindex += (ix&15);            //bits  0- 3
0231 
0232   //  int newiphi, newitheta, newix, newiy, newiz; 
0233   //  decode1(lutindex, newiphi, newitheta, newix, newiy, newiz);    
0234   return lutindex;
0235 
0236 }
0237 
0238 unsigned long encode2(int ien, int isec, int isid, int icha, int iparID){
0239   unsigned long  lutindex = (ien&511)<<12;   //bits  12-20
0240   lutindex += (iparID&63)<<6;                //bits  6-11
0241   lutindex += (icha&7)   <<3;                //bits  3- 5
0242   lutindex += (isec&3)   <<1;                //bits  1- 2
0243   lutindex += ((isid-1)&1);                  //bits  0
0244   //int newien, newisec, newisid, newicha, newipar; 
0245   //decode2(lutindex, newien, newisec, newisid, newicha, newipar);    
0246   return lutindex;
0247 }
0248