Back to home page

Project CMSSW displayed by LXR
 
 

    

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

0001 // Lateral shower development studies - without integration along phi

0002 // To be compiled using the makefile

0003 // 2 loops to get the maximum containment point, 

0004 // 3rd loop to do the analysis (in +- 2mm arount it)

0005 
0006 // to run on calibrated rechits reduced trees

0007 
0008 /*

0009   dwjang recipe to run this file

0010 

0011 ###### data section

0012 #

0013 #foreach imom ($mom)

0014 #    echo "lateralAllCalib input_tb_${imom}.list tb_${imom} $imom 0 0"

0015 #    ./lateralAllCalib input_tb_${imom}.list tb_${imom} $imom 0 0 >&! log.tb_${imom}

0016 #    echo "Done."

0017 #end

0018 #

0019 ###### geant4

0020 #

0021 #foreach imom ($mom)

0022 #    echo "lateralAllCalib input_g4_${imom}.list g4_${imom} $imom 0 1"

0023 #    ./lateralAllCalib input_g4_${imom}.list g4_${imom} $imom 0 1 >&! log.g4_${imom}

0024 #    echo "Done."

0025 #end

0026 #

0027 ##### gflash

0028 

0029 #foreach imom ($mom)

0030 #    echo "lateralAllCalib input_gf_${imom}.list gf_${imom} $imom 0 1"

0031 #    ./lateralAllCalib input_gf_${imom}.list gf_${imom} $imom 0 1 >&! log.gf_${imom}

0032 #    echo "Done."

0033 #end

0034 

0035 

0036 */
0037 
0038 
0039 //! c++ includes              

0040 #include <string>
0041 #include <iostream>
0042 #include <unistd.h>
0043 #include <cmath>
0044 #include <sstream>
0045 #include <fstream>
0046 #include <map>
0047 
0048 //! ROOT includes        

0049 #include "TROOT.h"
0050 #include "TFile.h"
0051 #include "TSystem.h"
0052 #include "TChain.h"
0053 #include "TBranch.h"
0054 #include "TH1.h"
0055 #include "TH2.h"
0056 #include "TF1.h"
0057 #include "TProfile.h"
0058 #include "TMinuit.h"
0059 #include "TString.h"
0060 
0061 using namespace std;
0062 
0063 //! my includes

0064 #include "fitMaximumPoint.cc"
0065 #include "myFunctions.cc"
0066 
0067 // useful functions

0068 bool selections(int xnum, double eneCent, double qualX, double qualY, int table, int simul );
0069 bool atBorder(int xtalNum);
0070 
0071 // Main program

0072 int main ( int argc, char **argv)
0073 { 
0074   // --------------------------------------------

0075   // input parameters

0076   char inputFileName[500], outputFileName[150];
0077   int beamEne, appendFlag, simul;
0078   if (argc==6) {
0079     strcpy(inputFileName,argv[1]); 
0080     strcpy(outputFileName,argv[2]);
0081     beamEne    = atoi(argv[3]);
0082     appendFlag = atoi(argv[4]);
0083     simul      = atoi(argv[5]);
0084   }
0085   else { 
0086     cout << " " << endl;
0087     cout << " Input: 1) input file name  "                           << endl;
0088     cout << "        2) output file name "                           << endl;
0089     cout << "        3) beam energy"                                 << endl;
0090     cout << "        4) appendflag"                                  << endl;
0091     cout << "        5) simulation (1) or data (0)                 " << endl; 
0092     cout << endl;
0093     return 1; 
0094    }
0095 
0096   // to know what I'm doing:

0097   cout << endl << endl;
0098   cout << "transverse shape analysis - running on calibrated rechits!" << endl;
0099   cout << "beam energy = " << beamEne << endl;
0100   if (simul == 0){ cout << "run on testbeam data" << endl;  }
0101   if (simul == 1){ cout << "run on simulated data" << endl; }
0102   cout << endl << endl;
0103   
0104 
0105 
0106 
0107 
0108   // -----------------------------------------------------------------------------------------

0109   //

0110   //              INITIALIZATIONS    

0111   //

0112   // -----------------------------------------------------------------------------------------

0113   
0114   
0115   // --------------------------------------

0116   // output files

0117   char outputDir[150];  
0118   strcpy(outputDir,"./");
0119   char outputFileTxtRatio[150];  
0120   strcpy(outputFileTxtRatio,outputDir);
0121   strcat(outputFileTxtRatio,outputFileName);
0122   strcat(outputFileTxtRatio,"_Ratios.txt");
0123   char outputFileTxtMatrix[150];  
0124   strcpy(outputFileTxtMatrix,outputDir);
0125   strcat(outputFileTxtMatrix,outputFileName);
0126   strcat(outputFileTxtMatrix,"_Matrix.txt");
0127   char outputFileMcpTxt[150];
0128   strcpy(outputFileMcpTxt,outputDir);
0129   strcat(outputFileMcpTxt,outputFileName);
0130   strcat(outputFileMcpTxt,"_mcp.txt");
0131   
0132 
0133   // --------------------------------------

0134   // variables

0135   char outfile[200];
0136   char name[200];
0137   char title[200], titleXax[200], titleXay[200];
0138   TFile *file;
0139   
0140   // histo range

0141   double infY = 0.;
0142   double supY = 200.;
0143   int nbinY = 100;
0144   if(beamEne == 20){
0145     infY = 0.;  supY = 20.;  nbinY = (int)((supY - infY)/0.20);  
0146     cout << "ene = 20 --> histo range from " << infY << " to " << supY << " with " << nbinY << " bins" << endl;
0147   }
0148   if(beamEne == 30){
0149     infY = 10.;  supY = 30.;  nbinY = (int)((supY - infY)/0.20);  
0150     cout << "ene = 30 --> histo range from " << infY << " to " << supY << " with " << nbinY << " bins" << endl;
0151   }
0152   if(beamEne == 40){
0153     infY = 10.;  supY = 50.;  nbinY = (int)((supY - infY)/0.20);  
0154     cout << "ene = 40 --> histo range from " << infY << " to " << supY << " with " << nbinY << " bins" << endl;
0155   }
0156   if(beamEne == 50){
0157     infY = 10.;  supY = 50.;  nbinY = (int)((supY - infY)/0.20);  
0158     cout << "ene = 50 --> histo range from " << infY << " to " << supY << " with " << nbinY << " bins" << endl;
0159   }
0160   if(beamEne == 80){
0161     infY = 20.;  supY = 80.;  nbinY = (int)((supY - infY)/0.20);  
0162     cout << "ene = 80 --> histo range from " << infY << " to " << supY << " with " << nbinY << " bins" << endl;
0163   }
0164   if(beamEne == 100){
0165     infY = 20.;  supY = 100.;  nbinY = (int)((supY - infY)/0.20);  
0166     cout << "ene = 100 --> histo range from " << infY << " to " << supY << " with " << nbinY << " bins" << endl;
0167   }
0168   if(beamEne == 120){
0169     infY = 30.;  supY = 120.;  nbinY = (int)((supY - infY)/0.20);  
0170     cout << "ene = 120 --> histo range from " << infY << " to " << supY << " with " << nbinY << " bins" << endl;
0171   }
0172   if(beamEne == 150){
0173     infY = 40.;  supY = 150.;  nbinY = (int)((supY - infY)/0.20);  
0174     cout << "ene = 150 --> histo range from " << infY << " to " << supY << " with " << nbinY << " bins" << endl;
0175   }
0176   if(beamEne == 500){
0177     infY = 200.;  supY = 450.;  nbinY = (int)((supY - infY)/0.20);  
0178     cout << "ene = 150 --> histo range from " << infY << " to " << supY << " with " << nbinY << " bins" << endl;
0179   }
0180   
0181 
0182   // ------------------------------------------

0183   // reading the list of input files with trees and building the chain

0184   std::vector<int> treeEntries;    // progressive number of entries when adding files

0185   treeEntries.push_back(0);  
0186   
0187   TChain *T = new TChain("T1");
0188   char Buffer[500];
0189   char MyRootFile[2000];  // [max length filename]

0190   cout << "input: " << inputFileName << endl; 
0191   ifstream *inputFile = new ifstream(inputFileName);
0192   while( !(inputFile->eof()) ){
0193     inputFile->getline(Buffer,500);
0194     if (!strstr(Buffer,"#") && !(strspn(Buffer," ") == strlen(Buffer))){
0195       sscanf(Buffer,"%s",MyRootFile);
0196       T->Add(MyRootFile);     
0197       treeEntries.push_back( T->GetEntries() );  
0198       cout << "chaining " << MyRootFile << endl;
0199     }}
0200   inputFile->close();
0201   delete inputFile;
0202 
0203 
0204   const static int numTrees = treeEntries.size();
0205   cout << endl;  
0206   cout << "in total " << (numTrees-1) << " files have been chained" << endl;
0207   for (int ii=0; ii< (numTrees-1); ii++){cout << "file " << ii << ", from entry " << treeEntries[ii] << ", to entry " << treeEntries[ii+1] << endl; }
0208   cout << endl;
0209   
0210 
0211 
0212   // ------------------------------------------

0213   // Tree construction

0214   int run, event;
0215   int xtalSM;
0216   int xtalEta, xtalPhi;
0217   int tbMoving;
0218   int crystal[49];
0219   double amplit[49];
0220   double hodoX, hodoY;
0221   double hodoQualityX, hodoQualityY;
0222   double hodoSlopeX, hodoSlopeY;
0223   //

0224   T->SetMakeClass(1);  
0225   T->SetBranchStatus("*",0);
0226   T->SetBranchStatus("run",1);
0227   T->SetBranchStatus("event",1);
0228   T->SetBranchStatus("xtalSM",1);
0229   T->SetBranchStatus("xtalEta",1);
0230   T->SetBranchStatus("xtalPhi",1);
0231   T->SetBranchStatus("amplit",1);
0232   T->SetBranchStatus("tbMoving",1);
0233   T->SetBranchStatus("hodoX",1);
0234   T->SetBranchStatus("hodoY",1);
0235   T->SetBranchStatus("hodoQualityX",1);
0236   T->SetBranchStatus("hodoQualityY",1);
0237   T->SetBranchStatus("hodoSlopeX",1);
0238   T->SetBranchStatus("hodoSlopeY",1);
0239   T->SetBranchStatus("crystal",1);
0240   //

0241   T->SetBranchAddress("run",&run);
0242   T->SetBranchAddress("event",&event);
0243   T->SetBranchAddress("xtalSM",&xtalSM);
0244   T->SetBranchAddress("xtalEta",&xtalEta);
0245   T->SetBranchAddress("xtalPhi",&xtalPhi);
0246   T->SetBranchAddress("amplit",amplit);
0247   T->SetBranchAddress("tbMoving",&tbMoving);
0248   T->SetBranchAddress("hodoX",&hodoX);
0249   T->SetBranchAddress("hodoY",&hodoY);
0250   T->SetBranchAddress("hodoQualityX",&hodoQualityX);
0251   T->SetBranchAddress("hodoQualityY",&hodoQualityY);
0252   T->SetBranchAddress("hodoSlopeX",&hodoSlopeX);
0253   T->SetBranchAddress("hodoSlopeY",&hodoSlopeY);
0254   T->SetBranchAddress("crystal",crystal);
0255   
0256   // entries

0257   float nEnt=T->GetEntries();
0258   cout << "Total number of events in loop is " << nEnt << endl; 
0259   
0260 
0261 
0262 
0263 
0264   // -----------------------------------------------------------------------------------------

0265   //

0266   //  FIRST LOOP ON EVENTS: registering the number of crystals             

0267   //

0268   // -----------------------------------------------------------------------------------------

0269   
0270   // first loop on the events: registering the number of crystals 

0271   int thisXinBeam = 0;
0272   int prevXinBeam = 0;
0273   std::vector<int> xInBeam, eInBeam;
0274   for (int entry0=0; entry0<nEnt; entry0++) { 
0275     
0276     T -> GetEntry(entry0);
0277     
0278     // registering the first crystal

0279     if (entry0 == 0){ 
0280       xInBeam.push_back(xtalSM); 
0281       eInBeam.push_back(xtalEta); 
0282     }        
0283     
0284     // did the crystal change?

0285     prevXinBeam = thisXinBeam;
0286     if ( (entry0 == 0) || (prevXinBeam == 0) ){ prevXinBeam = xtalSM; }
0287     thisXinBeam = xtalSM;
0288     if ( thisXinBeam != prevXinBeam ){ 
0289       cout << "entry " << entry0 << ", event = " << event 
0290        << ": crystal changed--> before " << prevXinBeam << ", now " << thisXinBeam << endl; } 
0291     
0292     // checking we are not coming back

0293     if ( thisXinBeam != prevXinBeam ){ 
0294       bool itisnew = true; 
0295       for(std::vector<int>::const_iterator myXib = xInBeam.begin(); myXib != xInBeam.end(); ++myXib){
0296     if (thisXinBeam == *myXib){ itisnew = false; }
0297       }   
0298       if ( itisnew){ xInBeam.push_back(thisXinBeam); eInBeam.push_back(xtalEta); }  
0299     }
0300   }
0301   
0302   // number of crystals in beam in the whole scan and corresponding eta

0303   int xibC = 0;
0304   int eibC = 0;
0305   const static int numXinB = xInBeam.size();
0306   const static int numEinB = eInBeam.size();
0307   int myXinB[numXinB];
0308   int myEinB[numXinB];
0309   for(std::vector<int>::const_iterator myXib = xInBeam.begin(); myXib != xInBeam.end(); ++myXib){
0310     myXinB[xibC] = *myXib;
0311     xibC++;
0312   }
0313   for(std::vector<int>::const_iterator myEib = eInBeam.begin(); myEib != eInBeam.end(); ++myEib){
0314     myEinB[eibC] = *myEib;
0315     eibC++;
0316   }
0317   
0318   
0319   // to summarize

0320   cout << endl; 
0321   cout << "in total: " << xInBeam.size() << " crystals" << endl;
0322   cout << "Xtals in beam:" << endl;
0323   for (int myXib=0; myXib<numXinB; myXib++){ 
0324     cout << "number " << myXib << ": xtal " << myXinB[myXib] << ", eta = " << myEinB[myXib] << endl; }
0325   cout << endl;
0326   cout << endl;
0327   cout << "NOW start the analysis" << endl;
0328   cout << endl;
0329   cout << endl;
0330   
0331   
0332 
0333   // variables for different xtals  -----------------------------

0334   int totalInTheRun = 0;
0335   int wrongXtalInTheRun = 0;
0336   int wrongQuality[numXinB];
0337   int movingTable[numXinB];       
0338   int highAmpl[numXinB]; 
0339   int negaAmpl[numXinB]; 
0340   int good0[numXinB];     
0341   int good1[numXinB];     
0342   int good2[numXinB];     
0343   int good3[numXinB];
0344   int wrongMaxAmp[numXinB];
0345   for(int myXib=0; myXib<numXinB; myXib++){
0346     movingTable[myXib]  = 0;
0347     wrongQuality[myXib] = 0;
0348     highAmpl[myXib]     = 0;
0349     negaAmpl[myXib]     = 0;
0350     good0[myXib]        = 0;
0351     good1[myXib]        = 0;
0352     good2[myXib]        = 0;
0353     good3[myXib]        = 0;
0354     wrongMaxAmp[myXib]  = 0;
0355   }
0356   
0357   
0358   
0359   
0360  
0361   // -----------------------------------------------------------------------------------------

0362   //

0363   //  SECOND LOOP ON EVENTS: checking everything is fine

0364   //

0365   // -----------------------------------------------------------------------------------------

0366   TH2F *H_ampVsX_bef[numXinB];     
0367   TH2F *H_ampVsX_aft[numXinB];     
0368   TH2F *H_ampVsY_bef[numXinB];     
0369   TH2F *H_ampVsY_aft[numXinB];     
0370   TH2F *H_MovVsEve_bef[numXinB];
0371   TH2F *H_MovVsEve_aft[numXinB];
0372   TH1F *H_hodoX_bef[numXinB];     
0373   TH1F *H_hodoX_aft[numXinB];     
0374   TH1F *H_hodoY_bef[numXinB];     
0375   TH1F *H_hodoY_aft[numXinB];     
0376   TH1F *H_hodoQualityX_bef[numXinB];     
0377   TH1F *H_hodoQualityY_bef[numXinB];     
0378   
0379   for (int myXib=0; myXib<numXinB; myXib++){
0380     
0381     sprintf (name,  "H_ampVsX_bef[%d]", myXib);
0382     sprintf (title, "E1 vs X - before cuts");    
0383     H_ampVsX_bef[myXib] = new TH2F(name, title, 80, -20., 20., 150, 0., 150.); 
0384     
0385     sprintf (name,  "H_ampVsX_aft[%d]", myXib);
0386     sprintf (title, "E1 vs X - after cuts");    
0387     H_ampVsX_aft[myXib] = new TH2F(name, title, 80, -20., 20., 150, 0., 150.); 
0388     
0389     sprintf (name,  "H_ampVsY_bef[%d]", myXib);
0390     sprintf (title, "E1 vs Y - before cuts");    
0391     H_ampVsY_bef[myXib] = new TH2F(name, title, 80, -20., 20., 150, 0., 150.); 
0392     
0393     sprintf (name,  "H_ampVsY_aft[%d]", myXib);
0394     sprintf (title, "E1 vs Y - after cuts");    
0395     H_ampVsY_aft[myXib] = new TH2F(name, title, 80, -20., 20., 150, 0., 150.); 
0396     
0397     sprintf (name,  "H_MovVsEve_bef[%d]", myXib);
0398     sprintf (title, "Table checking, before cuts");    
0399     H_MovVsEve_bef[myXib] = new TH2F(name, title, 80000, 0., 800000., 2, 0., 2.);
0400     
0401     sprintf (name,  "H_MovVsEve_aft[%d]", myXib);
0402     sprintf (title, "Table checking, after cuts");    
0403     H_MovVsEve_aft[myXib] = new TH2F(name, title, 80000, 0., 800000., 2, 0., 2.);
0404     
0405     sprintf (name,  "H_hodoX_bef[%d]", myXib);
0406     sprintf (title, "HodoX, before cuts");    
0407     H_hodoX_bef[myXib] = new TH1F(name, title, 80, -20., 20.);
0408     
0409     sprintf (name,  "H_hodoX_aft[%d]", myXib);
0410     sprintf (title, "HodoX, after cuts");    
0411     H_hodoX_aft[myXib] = new TH1F(name, title, 80, -20., 20.);
0412     
0413     sprintf (name,  "H_hodoY_bef[%d]", myXib);
0414     sprintf (title, "HodoY, before cuts");    
0415     H_hodoY_bef[myXib] = new TH1F(name, title, 80, -20., 20.);
0416     
0417     sprintf (name,  "H_hodoY_aft[%d]", myXib);
0418     sprintf (title, "HodoY, after cuts");    
0419     H_hodoY_aft[myXib] = new TH1F(name, title, 80, -20., 20.);
0420     
0421     sprintf (name,  "H_hodoQualityX_bef[%d]", myXib);
0422     sprintf (title, "Hodo quality X, before cuts");    
0423     H_hodoQualityX_bef[myXib] = new TH1F(name, title, 100, 0., 6.);
0424     
0425     sprintf (name,  "H_hodoQualityY_bef[%d]", myXib);
0426     sprintf (title, "Hodo quality Y, before cuts");    
0427     H_hodoQualityY_bef[myXib] = new TH1F(name, title, 100, 0., 6.);  
0428   }  
0429   
0430   
0431   
0432 
0433   // ----------------------------------------

0434   for (int entry1=0; entry1<nEnt; entry1++) { 
0435     
0436     // charging the value in the branches

0437     T -> GetEntry(entry1);
0438     if (entry1%100000 == 0){ cout << "Check loop: entry = " << entry1 << endl; }
0439     
0440     // table movement

0441     int table = tbMoving;
0442     
0443     // xtal in beam

0444     int xtalInBeam = xtalSM;
0445     int myHere = 2000;
0446     for (int myXib=0; myXib<numXinB; myXib++){ if (xtalInBeam == myXinB[myXib]){myHere = myXib;} }
0447     
0448     // to know which file I'm analizing

0449     int thisII = 0;
0450     for (int ii=0; ii< (numTrees-1); ii++){if ( (entry1 >= treeEntries[ii]) && (entry1 < treeEntries[ii+1]) ){ thisII = ii; } }
0451     int globalEve = event + treeEntries[thisII];   // number of the event in the chain

0452     
0453     // histos before selections

0454     H_MovVsEve_bef[myHere]     -> Fill(globalEve, table); 
0455     H_ampVsX_bef[myHere]       -> Fill(hodoX, amplit[24]);
0456     H_ampVsY_bef[myHere]       -> Fill(hodoY, amplit[24]);
0457     H_hodoX_bef[myHere]        -> Fill(hodoX);     
0458     H_hodoY_bef[myHere]        -> Fill(hodoY);
0459     H_hodoQualityX_bef[myHere] -> Fill(hodoQualityX);
0460     H_hodoQualityY_bef[myHere] -> Fill(hodoQualityY);
0461     
0462     // selections (no cut on max amplit xtal here) 

0463     bool select = selections(myHere, amplit[24], hodoQualityX, hodoQualityY, table, simul);
0464     if (!select){ continue; }
0465     good0[myHere]++;       
0466     
0467     // histos after selections

0468     H_MovVsEve_aft[myHere]     -> Fill(globalEve, table); 
0469     H_ampVsX_aft[myHere]       -> Fill(hodoX, amplit[24]);
0470     H_ampVsY_aft[myHere]       -> Fill(hodoY, amplit[24]);
0471     H_hodoX_aft[myHere]        -> Fill(hodoX);     
0472     H_hodoY_aft[myHere]        -> Fill(hodoY);
0473   }
0474   
0475   
0476   // saving histos

0477   strcpy(outfile,outputDir);
0478   strcat(outfile,outputFileName);
0479   strcat(outfile,"_checks.root");    
0480   if(appendFlag==0){ file = new TFile (outfile,"RECREATE","istogrammi"); }
0481   else             { file = new TFile (outfile,"UPDATE","istogrammi");   }
0482   
0483   for (int myXib=0; myXib<numXinB; myXib++){ 
0484     
0485     sprintf (name, "H_ampVsX_bef_x%d", myXinB[myXib]); 
0486     H_ampVsX_bef[myXib]->SetTitle(name);
0487     H_ampVsX_bef[myXib]->Write(name);
0488     
0489     sprintf (name, "H_ampVsY_bef_x%d", myXinB[myXib]); 
0490     H_ampVsY_bef[myXib]->SetTitle(name);
0491     H_ampVsY_bef[myXib]->Write(name);
0492     
0493     sprintf (name, "H_ampVsX_aft_x%d", myXinB[myXib]); 
0494     H_ampVsX_aft[myXib]->SetTitle(name);
0495     H_ampVsX_aft[myXib]->Write(name);
0496     
0497     sprintf (name, "H_ampVsY_aft_x%d", myXinB[myXib]);
0498     H_ampVsY_aft[myXib]->SetTitle(name);
0499     H_ampVsY_aft[myXib]->Write(name);
0500     
0501     sprintf (name, "H_MovVsEve_bef_x%d", myXinB[myXib]); 
0502     H_MovVsEve_bef[myXib]->SetTitle(name);
0503     H_MovVsEve_bef[myXib]->Write(name);
0504     
0505     sprintf (name, "H_MovVsEve_aft_x%d", myXinB[myXib]); 
0506     H_MovVsEve_aft[myXib]->SetTitle(name);
0507     H_MovVsEve_aft[myXib]->Write(name);
0508     
0509     sprintf (name, "H_hodoX_bef_x%d", myXinB[myXib]); 
0510     H_hodoX_bef[myXib]->SetTitle(name);
0511     H_hodoX_bef[myXib]->Write(name);
0512     
0513     sprintf (name, "H_hodoX_aft_x%d", myXinB[myXib]); 
0514     H_hodoX_aft[myXib]->SetTitle(name);
0515     H_hodoX_aft[myXib]->Write(name);
0516     
0517     sprintf (name, "H_hodoY_bef_x%d", myXinB[myXib]); 
0518     H_hodoY_bef[myXib]->SetTitle(name);
0519     H_hodoY_bef[myXib]->Write(name);
0520     
0521     sprintf (name, "H_hodoY_aft_x%d", myXinB[myXib]); 
0522     H_hodoY_aft[myXib]->SetTitle(name);
0523     H_hodoY_aft[myXib]->Write(name);
0524     
0525     sprintf (name, "H_hodoQualityX_bef_x%d", myXinB[myXib]); 
0526     H_hodoQualityX_bef[myXib]->SetTitle(name);
0527     H_hodoQualityX_bef[myXib]->Write(name);
0528     
0529     sprintf (name, "H_hodoQualityY_bef_x%d", myXinB[myXib]); 
0530     H_hodoQualityY_bef[myXib]->SetTitle(name);
0531     H_hodoQualityY_bef[myXib]->Write(name);
0532   }
0533   
0534   // deleting

0535   for (int myXib=0; myXib<numXinB; myXib++){ 
0536     delete H_ampVsX_bef[myXib];     
0537     delete H_ampVsX_aft[myXib];     
0538     delete H_ampVsY_bef[myXib];    
0539     delete H_ampVsY_aft[myXib];     
0540     delete H_MovVsEve_bef[myXib];
0541     delete H_MovVsEve_aft[myXib];
0542     delete H_hodoX_bef[myXib];     
0543     delete H_hodoX_aft[myXib];     
0544     delete H_hodoY_bef[myXib];     
0545     delete H_hodoY_aft[myXib];     
0546     delete H_hodoQualityX_bef[myXib];     
0547     delete H_hodoQualityY_bef[myXib];     
0548   }
0549   
0550   file->Close();
0551   delete file;
0552   
0553   
0554 
0555 
0556 
0557   // -----------------------------------------------------------------------------------------

0558   //

0559   //  THIRD LOOP ON EVENTS: first search loop for maximum containment point

0560   //

0561   // -----------------------------------------------------------------------------------------

0562 
0563   // region definition: CUT on the impinging particles position       ------------------

0564   double UpxFitX   =  10.;       // limits for fit X       

0565   double DownxFitX = -10.;
0566   double UpyFitX   =   4.;
0567   double DownyFitX =  -4.;
0568   double UpxFitY   =   4.;       // limits for fit Y       

0569   double DownxFitY =  -4.;
0570   double UpyFitY   =  10.;
0571   double DownyFitY = -10.;
0572 
0573   // preparing histos for each crystal  

0574   TH2F *H_mapXhodoVsE_1[numXinB];     
0575   TH2F *H_mapYhodoVsE_1[numXinB];     
0576   for (int myXib=0; myXib<numXinB; myXib++){
0577     sprintf (name,  "H_mapXhodoVsE_1[%d]", myXib);
0578     sprintf (title, "Pulse maximum vs X");    
0579     H_mapXhodoVsE_1[myXib] = new TH2F(name, title, 80, -20., 20., nbinY, infY, supY);
0580     H_mapXhodoVsE_1[myXib] -> GetXaxis()->SetTitle("X (mm)");
0581     H_mapXhodoVsE_1[myXib] -> GetYaxis()->SetTitle("amplit (adc counts)");
0582     
0583     sprintf (name,  "H_mapYhodoVsE_1[%d]", myXib);
0584     sprintf (title, "Pulse maximum vs Y");    
0585     H_mapYhodoVsE_1[myXib] = new TH2F(name, title, 80, -20., 20., nbinY, infY, supY);
0586     H_mapYhodoVsE_1[myXib] -> GetXaxis()->SetTitle("Y (mm)");
0587     H_mapYhodoVsE_1[myXib] -> GetYaxis()->SetTitle("amplit (adc counts)");
0588    }
0589   
0590   
0591   // analysis --------------- 

0592   for (int entry2=0; entry2<nEnt; entry2++) { 
0593     
0594     // counting the events

0595     totalInTheRun++; 
0596      
0597     // charging the value in the branches

0598     T -> GetEntry(entry2);
0599     if (entry2%100000 == 0){ cout << "First analysis loop: entry = " << entry2 << endl; }
0600     
0601     // table movement

0602     int table = tbMoving;
0603     
0604     // xtal in beam

0605     int xtalInBeam = xtalSM;
0606     int myHere = 2000;
0607     for (int myXib=0; myXib<numXinB; myXib++){ if (xtalInBeam == myXinB[myXib]){myHere = myXib;} }
0608     
0609     // to know which file I'm analizing

0610     int thisII = 0;
0611     for (int ii=0; ii< (numTrees-1); ii++){ if ( (entry2 >= treeEntries[ii]) && (entry2 < treeEntries[ii+1]) ){ thisII = ii; } }
0612     
0613     // selections (no cut on max amplit xtal here)  

0614     bool select = true;
0615     if ( myHere >1700)       { wrongXtalInTheRun++; select = false; }
0616     if ( amplit[24] > 999.)  { highAmpl[myHere]++;  select = false; }
0617     if ( amplit[24] < 0.5)   { negaAmpl[myHere]++;  select = false; }
0618     if ( !simul)             { 
0619       if (table)                                    { movingTable[myHere]++;  select = false; }
0620       if ((hodoQualityX>2.0) || (hodoQualityX<0.0)) { wrongQuality[myHere]++; select = false; }
0621       if ((hodoQualityY>2.0) || (hodoQualityY<0.0)) { wrongQuality[myHere]++; select = false; }
0622     }
0623     if (!select){ continue; }
0624     good1[myHere]++;       
0625     
0626     // basic cuts on the position

0627     if ((hodoX < DownxFitX) || (hodoX > UpxFitX)){continue;}
0628     if ((hodoY < DownyFitY) || (hodoY > UpyFitY)){continue;}
0629     
0630     // maximum containment point: filling histos, after cuts on the impact point

0631     if ((DownyFitX<hodoY) && (hodoY<UpyFitX)){ H_mapXhodoVsE_1[myHere] -> Fill(hodoX, amplit[24]); }
0632     if ((DownxFitY<hodoX) && (hodoX<UpxFitY)){ H_mapYhodoVsE_1[myHere] -> Fill(hodoY, amplit[24]); }
0633     
0634   } // end loop over entries

0635   
0636 
0637   // First fit with the polynomial function + save histos

0638   double p0x_1[numXinB], p0y_1[numXinB];
0639   for (int myXib=0; myXib<numXinB; myXib++){ 
0640     p0x_1[myXib] = 1000.;
0641     p0y_1[myXib] = 1000.;
0642   }   
0643   
0644   strcpy(outfile,outputDir);
0645   strcat(outfile,outputFileName);
0646   strcat(outfile,"_histos.root");    
0647   if(appendFlag==0){ file = new TFile (outfile,"RECREATE","istogrammi"); }
0648   else             { file = new TFile (outfile,"UPDATE","istogrammi");   }
0649   
0650   for (int myXib=0; myXib<numXinB; myXib++){ 
0651     
0652     // fitting

0653     p0x_1[myXib] = Fit_MaximumPoint( H_mapXhodoVsE_1[myXib], 10);
0654     p0y_1[myXib] = Fit_MaximumPoint( H_mapYhodoVsE_1[myXib], 10);
0655 
0656     // saving

0657     sprintf (name, "H_mapXhodoVsE_1loop_x%d", myXinB[myXib]); 
0658     H_mapXhodoVsE_1[myXib]->SetTitle(name);
0659     H_mapXhodoVsE_1[myXib]->Write(name);
0660 
0661     sprintf (name, "H_mapYhodoVsE_1loop_x%d", myXinB[myXib]); 
0662     H_mapYhodoVsE_1[myXib]->SetTitle(name);
0663     H_mapYhodoVsE_1[myXib]->Write(name);
0664   }
0665   file->Close();
0666   delete file;
0667   
0668   
0669   // deleting

0670   for (int myXib=0; myXib<numXinB; myXib++){ 
0671     delete H_mapXhodoVsE_1[myXib];     
0672     delete H_mapYhodoVsE_1[myXib];     
0673   }
0674   
0675   
0676   
0677   // -----------------------------------------------------------------------------------------

0678   //

0679   //  FOURTH LOOP ON EVENTS: second search loop for maximum containment point

0680   //

0681   // -----------------------------------------------------------------------------------------

0682 
0683   // New limits ------------------------------------

0684   double DownyFitX_v[numXinB];  
0685   double DownxFitY_v[numXinB];  
0686   double UpyFitX_v[numXinB];  
0687   double UpxFitY_v[numXinB];  
0688   for (int myXib=0; myXib<numXinB; myXib++){ 
0689     DownyFitX_v[myXib] = p0y_1[myXib]-2;
0690     DownxFitY_v[myXib] = p0x_1[myXib]-2;
0691     UpyFitX_v[myXib]   = p0y_1[myXib]+2;
0692     UpxFitY_v[myXib]   = p0x_1[myXib]+2;
0693   }
0694   
0695   // preparing histos for each crystal  

0696   TH2F *H_mapXhodoVsE_2[numXinB];     
0697   TH2F *H_mapYhodoVsE_2[numXinB];     
0698   for (int myXib=0; myXib<numXinB; myXib++){
0699     
0700     sprintf (name,  "H_mapXhodoVsE_2[%d]", myXib);
0701     sprintf (title, "Pulse maximum vs X");    
0702     H_mapXhodoVsE_2[myXib] = new TH2F(name, title, 80, -20., 20., nbinY, infY, supY);
0703     H_mapXhodoVsE_2[myXib] -> GetXaxis()->SetTitle("X (mm)");
0704     H_mapXhodoVsE_2[myXib] -> GetYaxis()->SetTitle("amplit (adc counts)");
0705     
0706     sprintf (name,  "H_mapYhodoVsE_2[%d]", myXib);
0707     sprintf (title, "Pulse maximum vs Y");    
0708     H_mapYhodoVsE_2[myXib] = new TH2F(name, title, 80, -20., 20., nbinY, infY, supY);
0709     H_mapYhodoVsE_2[myXib] -> GetXaxis()->SetTitle("Y (mm)");
0710     H_mapYhodoVsE_2[myXib] -> GetYaxis()->SetTitle("amplit (adc counts)");
0711   }
0712   
0713   for (int entry3=0; entry3<nEnt; entry3++) { 
0714     
0715     // charging the value in the branches

0716     T -> GetEntry(entry3);
0717     if (entry3%100000 == 0){ cout << "Second analysis loop: entry = " << entry3 << endl; }
0718     
0719     // xtal in beam

0720     int xtalInBeam = xtalSM;
0721     int myHere = 2000;
0722     for (int myXib=0; myXib<numXinB; myXib++){ if (xtalInBeam == myXinB[myXib]){myHere = myXib;} }
0723     
0724     // table movement

0725     int table = tbMoving;
0726     
0727     // to know which file I'm analizing

0728     int thisII = 0;
0729     for (int ii=0; ii< (numTrees-1); ii++){if ( (entry3 >= treeEntries[ii]) && (entry3 < treeEntries[ii+1]) ){ thisII = ii; } }
0730     
0731     // selections (no cut on max amplit xtal here) 

0732     bool select = selections(myHere, amplit[24], hodoQualityX, hodoQualityY, table, simul); 
0733     if (!select){ continue; }
0734     good2[myHere]++;       
0735     
0736     // basic cuts on the position

0737     if ((hodoX < DownxFitX) || (hodoX > UpxFitX)){continue;}
0738     if ((hodoY < DownyFitY) || (hodoY > UpyFitY)){continue;}
0739 
0740     // maximum containment point: filling histos, after cuts on the impact point

0741     if ((DownyFitX_v[myHere]<hodoY) && (hodoY<UpyFitX_v[myHere])){ H_mapXhodoVsE_2[myHere] -> Fill(hodoX, amplit[24]); }
0742     if ((DownxFitY_v[myHere]<hodoX) && (hodoX<UpxFitY_v[myHere])){ H_mapYhodoVsE_2[myHere] -> Fill(hodoY, amplit[24]); }
0743      
0744   } // end loop over entries

0745   
0746 
0747   
0748   // ---------------------------------------------------------

0749   // Second fit with the polynomial function 

0750   strcpy(outfile,outputDir);
0751   strcat(outfile,outputFileName);
0752   strcat(outfile,"_histos.root");    
0753   file = new TFile (outfile,"UPDATE","istogrammi");   
0754   
0755   
0756   double p0x_2[numXinB], p0y_2[numXinB];
0757   double p0x_c22[numXinB],  p0y_c22[numXinB];
0758   Stat_t p0x_ent2[numXinB], p0y_ent2[numXinB];
0759   
0760   for (int myXib=0; myXib<numXinB; myXib++){ 
0761     
0762     // initializing

0763     p0x_2[myXib]    =  1000.;
0764     p0y_2[myXib]    =  1000.;
0765     p0x_c22[myXib]  =  1000.;
0766     p0y_c22[myXib]  =  1000.;
0767     p0x_ent2[myXib] = -15000;
0768     p0y_ent2[myXib] = -15000;
0769     
0770     // fitting: value

0771     p0x_2[myXib] = Fit_MaximumPoint( H_mapXhodoVsE_2[myXib], 10);
0772     p0y_2[myXib] = Fit_MaximumPoint( H_mapYhodoVsE_2[myXib], 10);
0773     
0774     // fitting: chi2s

0775     p0x_c22[myXib]  = Fit_MaximumPoint( H_mapXhodoVsE_2[myXib], 30);
0776     p0y_c22[myXib]  = Fit_MaximumPoint( H_mapYhodoVsE_2[myXib], 30);
0777     
0778     // histo entries

0779     p0x_ent2[myXib] = H_mapXhodoVsE_2[myXib]->GetEntries();
0780     p0y_ent2[myXib] = H_mapYhodoVsE_2[myXib]->GetEntries();
0781 
0782     // saving

0783     sprintf (name, "H_mapXhodoVsE_2loop_x%d", myXinB[myXib]); 
0784     H_mapXhodoVsE_2[myXib]->SetTitle(name);
0785     H_mapXhodoVsE_2[myXib]->Write(name);
0786 
0787     sprintf (name, "H_mapYhodoVsE_2loop_x%d", myXinB[myXib]); 
0788     H_mapYhodoVsE_2[myXib]->SetTitle(name);
0789     H_mapYhodoVsE_2[myXib]->Write(name);
0790   }
0791   file->Close();
0792   delete file;
0793   
0794   // deleting

0795   for (int myXib=0; myXib<numXinB; myXib++){ 
0796     delete H_mapXhodoVsE_2[myXib];     
0797     delete H_mapYhodoVsE_2[myXib];     
0798   }
0799   
0800   
0801   
0802   
0803   // --------------------------------------------------------------

0804   // now we have the maximum containment point for each xtal 

0805   ofstream *ResFileMc;
0806   if(appendFlag==0){
0807     ResFileMc = new ofstream(outputFileMcpTxt,ios::out);
0808     *ResFileMc << "#Xtal"    << "\t"
0809            << "entriesX" << "\t" << "maxX" << "\t" << "chi2X" << "\t"
0810            << "entriesY" << "\t" << "maxY" << "\t" << "chi2Y" << endl; }
0811   else { ResFileMc = new ofstream(outputFileMcpTxt,ios::app); }
0812   for (int myXib=0; myXib<numXinB; myXib++){
0813     *ResFileMc << myXinB[myXib]   << "\t"  
0814            << p0x_ent2[myXib] << "\t" << p0x_2[myXib] << "\t" << p0x_c22[myXib]  << "\t"
0815            << p0y_ent2[myXib] << "\t" << p0y_2[myXib] << "\t" << p0y_c22[myXib]  << endl;
0816   }
0817   cout << endl;
0818   
0819 
0820   // to skip xtals where the mcp is not known

0821   ofstream *badMcPFile = new ofstream("badMcpFile.txt",ios::out);
0822   map<int, bool> badMcp;
0823   for (int myXib=0; myXib<numXinB; myXib++){
0824     bool isItOk = true;
0825     if ( (p0x_c22[myXib] > 3.5) || (p0y_c22[myXib] > 3.5) ){ isItOk = false; }
0826     badMcp.insert(pair<int, bool>(myXinB[myXib], isItOk));
0827   }
0828   
0829   
0830 
0831 
0832   // -----------------------------------------------------------------------------------------

0833   //

0834   //  FIFTH LOOP ON EVENTS: final analysis

0835   //

0836   // -----------------------------------------------------------------------------------------

0837   
0838   // New limits ------------------------------------

0839   double downy[numXinB];  
0840   double downx[numXinB];  
0841   double upy[numXinB];  
0842   double upx[numXinB];  
0843   
0844   for (int myXib=0; myXib<numXinB; myXib++){ 
0845     downy[myXib] = p0y_2[myXib]-2;
0846     downx[myXib] = p0x_2[myXib]-2;
0847     upy[myXib]   = p0y_2[myXib]+2;
0848     upx[myXib]   = p0x_2[myXib]+2;
0849   }
0850   
0851   // histos --------------------------------

0852   TH1F *H_eMatrix[3][numXinB];         // [0] = 1x1         [1] = 3x3         [2] = 5x5    

0853   TH1F *H_eRatio[3][numXinB];          // [0] = e1/e9       [1] = e1/e25      [2] = e9/e25

0854 
0855   TH1F *H_energyNxN[3];
0856   TProfile *H_energyieta;
0857   TProfile *H_energyiphi;
0858 
0859   H_energyieta = new TProfile("H_energyieta","energy by ieta", 5, -2.5, 2.5,0.0,1.0);
0860   H_energyiphi = new TProfile("H_energyiphi","energy by iphi", 5, -2.5, 2.5,0.0,1.0);
0861 
0862   for (int ii=0; ii<3; ii++){     
0863     sprintf(name, "H_energy%dx%d",2*ii+1,2*ii+1);
0864     sprintf(title, "%dx%d crystal energy for beam energy = %d",2*ii+1,2*ii+1,beamEne);
0865     sprintf(titleXax, "E%dx%d (GeV)",2*ii+1,2*ii+1);
0866     
0867     H_energyNxN[ii] = new TH1F(name, title, 280, 0.5, 1.1);
0868     H_energyNxN[ii]-> GetXaxis()->SetTitle(titleXax);
0869   }
0870   TH1F* H_energyE1E9 = new TH1F("H_energyE1E9","Energy ratio of 1x1 to 3x3 crystals;E_{1}/E_{9}",200,0.7,0.9);
0871   TH1F* H_energyE1E25 = new TH1F("H_energyE1E25","Energy ratio of 1x1 to 5x5 crystals;E_{1}/E_{25}",200,0.7,0.9);
0872   TH1F* H_energyE9E25 = new TH1F("H_energyE9E25","Energy ratio of 3x3 to 5x5 crystals;E_{9}/E_{25}",200,0.9,1.0);
0873   
0874   for (int myXib=0; myXib<numXinB; myXib++){
0875     for (int ii=0; ii<3; ii++){     
0876       sprintf (name, "H_eMatrix[%d][%d]", ii, myXib);
0877       if(ii==0){ sprintf(title, "1x1 energy");  sprintf(titleXax, "E1 (GeV)"); }  
0878       if(ii==1){ sprintf(title, "3x3 energy");  sprintf(titleXax, "E9 (GeV)"); }   
0879       if(ii==2){ sprintf(title, "5x5 energy");  sprintf(titleXax, "E25 (GeV)"); }   
0880 
0881       if (beamEne==20){ 
0882     if (!simul){
0883       if (ii==0)         { H_eMatrix[ii][myXib] = new TH1F(name, title,  100,  0., 30.); }
0884       if (ii==1 || ii==2){ H_eMatrix[ii][myXib] = new TH1F(name, title,  100,  0., 30.); }
0885     }
0886     if (simul){
0887       if (ii==0)         { H_eMatrix[ii][myXib] = new TH1F(name, title,  100,  10., 20.); }
0888       if (ii==1)         { H_eMatrix[ii][myXib] = new TH1F(name, title,  100,  15., 22.); }
0889       if (ii==2)         { H_eMatrix[ii][myXib] = new TH1F(name, title,  100,  15., 22.); }
0890     }
0891       }
0892       else if (beamEne==30) { 
0893     if (!simul){
0894       if (ii==0)          { H_eMatrix[ii][myXib] = new TH1F(name, title,  100, 0., 40.); } 
0895       if (ii==1 || ii==2) { H_eMatrix[ii][myXib] = new TH1F(name, title,  100, 0., 40.); } 
0896     }
0897     if (simul){
0898       if (ii==0)          { H_eMatrix[ii][myXib] = new TH1F(name, title,  100, 15., 27.); }  
0899       if (ii==1)          { H_eMatrix[ii][myXib] = new TH1F(name, title,  100, 21., 32.); }  
0900       if (ii==2)          { H_eMatrix[ii][myXib] = new TH1F(name, title,  100, 24., 35.); }  
0901     }
0902       }
0903       else if (beamEne==40) { 
0904     if (!simul){
0905       if (ii==0)          { H_eMatrix[ii][myXib] = new TH1F(name, title,  100, 0., 60.); }  
0906       if (ii==1 || ii==2) { H_eMatrix[ii][myXib] = new TH1F(name, title,  100, 0., 60.); }  
0907     }
0908     if (simul){
0909       if (ii==0)          { H_eMatrix[ii][myXib] = new TH1F(name, title,  100, 25., 36.); }  
0910       if (ii==1)          { H_eMatrix[ii][myXib] = new TH1F(name, title,  100, 31., 42.); }  
0911       if (ii==2)          { H_eMatrix[ii][myXib] = new TH1F(name, title,  100, 33., 44.); }  
0912     }
0913       }
0914       else if (beamEne==50) { 
0915     if (!simul){
0916       if (ii==0)          { H_eMatrix[ii][myXib] = new TH1F(name, title,  100, 0., 60.); }  
0917       if (ii==1 || ii==2) { H_eMatrix[ii][myXib] = new TH1F(name, title,  100, 0., 60.); }  
0918     }
0919     if (simul){
0920       // electrons

0921       if (ii==0)          { H_eMatrix[ii][myXib] = new TH1F(name, title,  100, 26., 43.); }  
0922       if (ii==1)          { H_eMatrix[ii][myXib] = new TH1F(name, title,  100, 42., 50.); }  
0923       if (ii==2)          { H_eMatrix[ii][myXib] = new TH1F(name, title,  100, 45., 52.); }  
0924       // pions

0925       // if (ii==0)          { H_eMatrix[ii][myXib] = new TH1F(name, title,  100, 0., 43.); }  

0926       // if (ii==1)          { H_eMatrix[ii][myXib] = new TH1F(name, title,  100, 0., 50.); }  

0927       // if (ii==2)          { H_eMatrix[ii][myXib] = new TH1F(name, title,  100, 0., 52.); }  

0928     }
0929       }
0930       else if (beamEne==80) { 
0931     if (!simul){
0932       if (ii==0)          { H_eMatrix[ii][myXib] = new TH1F(name, title,  100, 0., 90.); }  
0933       if (ii==1)          { H_eMatrix[ii][myXib] = new TH1F(name, title,  100, 0., 90.); }  
0934       if (ii==2)          { H_eMatrix[ii][myXib] = new TH1F(name, title,  100, 0., 90.); }  
0935     }                                
0936     if (simul){                                                          
0937       if (ii==0)          { H_eMatrix[ii][myXib] = new TH1F(name, title,  100, 50., 70.); }  
0938       if (ii==1)          { H_eMatrix[ii][myXib] = new TH1F(name, title,  100, 65., 85.); }  
0939       if (ii==2)          { H_eMatrix[ii][myXib] = new TH1F(name, title,  100, 68., 88.); }  
0940     }
0941       }
0942       else if (beamEne==100){ 
0943     if (simul){                                              
0944       if (ii==0)          { H_eMatrix[ii][myXib] = new TH1F(name, title, 60,  60., 90.); } 
0945       if (ii==1)          { H_eMatrix[ii][myXib] = new TH1F(name, title, 100, 80., 105.); } 
0946       if (ii==2)          { H_eMatrix[ii][myXib] = new TH1F(name, title, 100, 90., 110.); } 
0947     }
0948       }
0949       else if (beamEne==120){ 
0950     if (!simul){
0951       // if (ii==0)          { H_eMatrix[ii][myXib] = new TH1F(name, title, 100,  95., 110.); }  

0952       // if (ii==1)          { H_eMatrix[ii][myXib] = new TH1F(name, title, 100, 118., 128.); }   

0953       // if (ii==2)          { H_eMatrix[ii][myXib] = new TH1F(name, title, 100, 123., 133.); }  

0954       if (ii==0)          { H_eMatrix[ii][myXib] = new TH1F(name, title, 100, 85., 100.); }  
0955       if (ii==1)          { H_eMatrix[ii][myXib] = new TH1F(name, title, 100, 90., 120.); }   
0956       // if (ii==2)          { H_eMatrix[ii][myXib] = new TH1F(name, title, 100, 0.95, 1.05); }  

0957       if (ii==2)          { H_eMatrix[ii][myXib] = new TH1F(name, title, 100, 100., 130.); }  
0958     }                                        
0959     if (simul){                                              
0960       // electrons

0961       if (ii==0)          { H_eMatrix[ii][myXib] = new TH1F(name, title, 100,  85., 110.); } 
0962       if (ii==1)          { H_eMatrix[ii][myXib] = new TH1F(name, title, 100, 106., 126.); } 
0963       if (ii==2)          { H_eMatrix[ii][myXib] = new TH1F(name, title, 100, 110., 130.); } 
0964       // pions

0965       // if (ii==0)          { H_eMatrix[ii][myXib] = new TH1F(name, title, 100, 0., 105.); } 

0966       // if (ii==1)          { H_eMatrix[ii][myXib] = new TH1F(name, title, 100, 0., 116.); } 

0967       // if (ii==2)          { H_eMatrix[ii][myXib] = new TH1F(name, title, 100, 0., 120.); } 

0968     }
0969       }
0970       else if (beamEne==150){ 
0971     if (!simul){
0972       if (ii==0)          { H_eMatrix[ii][myXib] = new TH1F(name, title, 100, 110., 130.); }  
0973       if (ii==1)          { H_eMatrix[ii][myXib] = new TH1F(name, title, 100, 130., 160.); }   
0974       if (ii==2)          { H_eMatrix[ii][myXib] = new TH1F(name, title, 100, 140., 160.); }  
0975     }                                        
0976     if (simul){                                              
0977       // electrons

0978       if (ii==0)          { H_eMatrix[ii][myXib] = new TH1F(name, title, 100, 110., 130.); } 
0979       if (ii==1)          { H_eMatrix[ii][myXib] = new TH1F(name, title, 100, 130., 160.); } 
0980       if (ii==2)          { H_eMatrix[ii][myXib] = new TH1F(name, title, 100, 140., 160.); } 
0981     }
0982       }
0983       else if (beamEne==500){ 
0984     if (simul){                                              
0985       //if (ii==0)          { H_eMatrix[ii][myXib] = new TH1F(name, title, 100, 380., 440.); } 

0986       //if (ii==1)          { H_eMatrix[ii][myXib] = new TH1F(name, title, 100, 400., 510.); } 

0987       //if (ii==2)          { H_eMatrix[ii][myXib] = new TH1F(name, title, 100, 400., 520.); } 

0988       //

0989       if (ii==0)          { H_eMatrix[ii][myXib] = new TH1F(name, title, 100, 380., 440.); } 
0990       if (ii==1)          { H_eMatrix[ii][myXib] = new TH1F(name, title, 101, 440., 500.); } 
0991       if (ii==2)          { H_eMatrix[ii][myXib] = new TH1F(name, title, 101, 440., 520.); } 
0992     }
0993       }
0994       H_eMatrix[ii][myXib] -> GetXaxis()->SetTitle(titleXax);
0995     }
0996 
0997     for (int ii=0; ii<3; ii++){      
0998       sprintf (name, "H_eRatio[%d][%d]", ii, myXib);
0999       if(ii==0){ sprintf(title, "E1/E9");   sprintf(titleXax, "E1/E9");  H_eRatio[ii][myXib] = new TH1F(name, title, 80, 0.80, 0.88); }
1000       if(ii==1){ sprintf(title, "E1/E25");  sprintf(titleXax, "E1/E25"); H_eRatio[ii][myXib] = new TH1F(name, title, 80, 0.75, 0.85);}
1001       if(ii==2){ sprintf(title, "E9/E25");  sprintf(titleXax, "E9/E25"); H_eRatio[ii][myXib] = new TH1F(name, title, 60, 0.94, 0.98);}
1002       H_eRatio[ii][myXib] -> GetXaxis()->SetTitle(titleXax);
1003     }
1004   }
1005   
1006 
1007 
1008   // -----------------------------------------------

1009 
1010   for (int entry3=0; entry3<nEnt; entry3++) { 
1011     
1012     // charging the value in the branches

1013     T -> GetEntry(entry3);
1014     if (entry3%100000 == 0){ cout << "Final analysis loop: entry = " << entry3 << endl; }
1015 
1016     // xtal in beam

1017     int xtalInBeam = xtalSM;
1018     int myHere = 2000;
1019     for (int myXib=0; myXib<numXinB; myXib++){ if (xtalInBeam == myXinB[myXib]){myHere = myXib;} }
1020     
1021     // table movement

1022     int table = tbMoving;
1023     
1024     // to know which file I'm analizing

1025     int thisII = 0;
1026     for (int ii=0; ii< (numTrees-1); ii++){
1027       if ( (entry3 >= treeEntries[ii]) && (entry3 < treeEntries[ii+1]) ){ thisII = ii; } 
1028     }
1029 
1030     // basic cuts (no cut on max amplit xtal here) 

1031     bool select = selections(myHere, amplit[24], hodoQualityX, hodoQualityY, table, simul); 
1032     if (!select){ continue; }
1033     good3[myHere]++;      
1034 
1035     // no xtals with bad determination of the max conteinment point

1036     if(!simul){
1037       bool mcpFound = true;
1038       map<int, bool>::const_iterator mcpIt;
1039       mcpIt = badMcp.find(xtalInBeam);
1040       if (mcpIt != badMcp.end()){ mcpFound = mcpIt->second; }
1041       if (!mcpFound){ *badMcPFile << mcpIt->first << endl; continue; }
1042     }
1043 
1044     // selections: cuts on the impact point 

1045     if ((hodoX < downx[myHere]) || (hodoX > upx[myHere])){continue;}
1046     if ((hodoY < downy[myHere]) || (hodoY > upy[myHere])){continue;}
1047 
1048     // no xtals at borders, basic cut to remove what enter at least 3x3

1049     bool amIatBorder = atBorder(xtalInBeam);
1050     if (amIatBorder){ continue; }      
1051 
1052     // further selection: only if xtal-in-beam = max-amplitude-xtal

1053     int maxAmpInMtx = maxAmplitInMatrix(amplit);
1054     if (maxAmpInMtx != 24){ wrongMaxAmp[myHere]++;  continue; }
1055 
1056     // computing variables

1057     double e1e9  = -50.;
1058     double e1e25 = -50.;
1059     double e9e25 = -50.;
1060     double e1    = ene1x1_xtal(amplit,24); 
1061     double e9    = ene3x3_xtal(amplit,24); 
1062     double e25   = ene5x5_xtal(amplit,24); 
1063 
1064     if( (e1>-20) && (e9>-20) ) { e1e9  = e1/e9;  }
1065     if( (e1>-20) && (e25>-20) ){ e1e25 = e1/e25; }
1066     if( (e9>-20) && (e25>-20) ){ e9e25 = e9/e25; }
1067     
1068     // filling histos 

1069     H_eMatrix[0][myHere] -> Fill(e1);
1070     H_eMatrix[1][myHere] -> Fill(e9);
1071     H_eMatrix[2][myHere] -> Fill(e25);
1072     H_eRatio[0][myHere]  -> Fill(e1e9);
1073     H_eRatio[1][myHere]  -> Fill(e1e25);
1074     H_eRatio[2][myHere]  -> Fill(e9e25);
1075 
1076     //@@@@syjun

1077     double Eieta[5] = {0};
1078     double Eiphi[5] = {0};
1079     if(beamEne!=0) {
1080       H_energyNxN[0]->Fill(e1/(1.0*beamEne));
1081       H_energyNxN[1]->Fill(e9/(1.0*beamEne));
1082       H_energyNxN[2]->Fill(e25/(1.0*beamEne));
1083       H_energyE1E9->Fill(e1/e9);
1084       H_energyE1E25->Fill(e1/e25);
1085       H_energyE9E25->Fill(e9/e25);
1086 
1087       energy_ieta(amplit,Eieta);
1088       energy_iphi(amplit,Eiphi);
1089 
1090       for(int ii = 0 ; ii < 5 ; ii++) {
1091     H_energyieta->Fill(-2.0+1.0*ii,Eieta[ii]/(1.0*beamEne));
1092     H_energyiphi->Fill(-2.0+1.0*ii,Eiphi[ii]/(1.0*beamEne));
1093       }
1094 
1095     }
1096 
1097   } // end loop on the events

1098 
1099 
1100   cout << endl;
1101   cout << "End Third Analysis Loop" << endl;
1102   cout << endl;
1103 
1104 
1105   // --------------------------------------------------------

1106   // Saving histos into files before the fit

1107   strcpy(outfile,outputDir);
1108   strcat(outfile,outputFileName);
1109   strcat(outfile,"_befFit.root");
1110   if(appendFlag==0){ file = new TFile (outfile,"RECREATE","istogrammi"); }
1111   else             { file = new TFile (outfile,"UPDATE","istogrammi"); }
1112 
1113   for (int ii=0; ii<3; ii++){     
1114     sprintf(name, "H_energy%dx%d",2*ii+1,2*ii+1);
1115     H_energyNxN[ii]->SetTitle(name);
1116     H_energyNxN[ii]->Write(name);
1117   }
1118   H_energyE1E9->Write();
1119   H_energyE1E25->Write();
1120   H_energyE9E25->Write();
1121   H_energyieta->Write();  
1122   H_energyiphi->Write();  
1123 
1124   for (int myXib=0; myXib<numXinB; myXib++){ 
1125 
1126     for (int ii=0; ii<3; ii++){      
1127       if(ii==0){ sprintf (name, "H_e1_x%d_ene%d",  myXinB[myXib], beamEne); }
1128       if(ii==1){ sprintf (name, "H_e9_x%d_ene%d",  myXinB[myXib], beamEne); }
1129       if(ii==2){ sprintf (name, "H_e25_x%d_ene%d", myXinB[myXib], beamEne); }
1130       H_eMatrix[ii][myXib] -> SetTitle(name);
1131       H_eMatrix[ii][myXib] -> Write(name);
1132     }
1133     
1134     for (int ii=0; ii<3; ii++){      
1135       if(ii==0){ sprintf (name, "H_e1e9_x%d_ene%d",  myXinB[myXib], beamEne); }
1136       if(ii==1){ sprintf (name, "H_e1e25_x%d_ene%d", myXinB[myXib], beamEne); }
1137       if(ii==2){ sprintf (name, "H_e9e25_x%d_ene%d", myXinB[myXib], beamEne); }
1138       H_eRatio[ii][myXib] -> SetTitle(name);
1139       H_eRatio[ii][myXib] -> Write(name);
1140     }
1141   }
1142 
1143   file->Close();
1144   delete file;
1145 
1146 
1147   // --------------------------------------------------------

1148   // -----------------------------------------------------------------------------------------

1149   // -----------------------------------------------------------------------------------------

1150   //

1151   //  FITTING 

1152   //

1153   // -----------------------------------------------------------------------------------------

1154   // -----------------------------------------------------------------------------------------

1155   // -----------------------------------------------------------------------------------------

1156 
1157   // NxN matrices

1158   ofstream *ResFileMatrix;
1159   if(appendFlag==0){ 
1160     ResFileMatrix = new ofstream(outputFileTxtMatrix,ios::out);    
1161     *ResFileMatrix << "#Energy"  << "\t" << "xtal"        << "\t" << "eta" << "\t"  
1162            << "sample"  << "\t" << "entries"     << "\t" 
1163            << "H_mean"  << "\t" << "H_rms"       << "\t" 
1164            << "G_mean"  << "\t" << "G_mean err"  << "\t" 
1165            << "G_sigma" << "\t" << "G_sigma err" << endl;} 
1166   else { ResFileMatrix = new ofstream(outputFileTxtMatrix,ios::app); }  
1167 
1168 
1169   // fitting 

1170   for (int myXib=0; myXib<numXinB; myXib++){ 
1171         
1172     for (int ii=0; ii<3; ii++){ 
1173 
1174       // histo parameters

1175       Stat_t h_entries = H_eMatrix[ii][myXib]->GetEntries();      
1176       int peakBin      = H_eMatrix[ii][myXib]->GetMaximumBin();
1177       double h_norm    = H_eMatrix[ii][myXib]->GetMaximum();
1178       double h_rms     = H_eMatrix[ii][myXib]->GetRMS();
1179       double h_mean    = H_eMatrix[ii][myXib]->GetMean();
1180       double h_peak    = H_eMatrix[ii][myXib]->GetBinCenter(peakBin);
1181 
1182 
1183       // gaussian fit to initialize

1184       TF1 *gausa = new TF1 ("gausa","[0]*exp(-1*(x-[1])*(x-[1])/2/[2]/[2])",h_peak-10*h_rms,h_peak+10*h_rms);
1185       gausa->SetParameters(h_norm,h_peak,h_rms);
1186       H_eMatrix[ii][myXib]->Fit(gausa,"","",h_peak-3*h_rms,h_peak+3*h_rms);
1187       double gausNorm  = gausa->GetParameter(0);
1188       double gausMean  = gausa->GetParameter(1);
1189       double gausSigma = fabs(gausa->GetParameter(2));
1190       double gausChi2  = gausa->GetChisquare()/gausa->GetNDF();
1191       if (gausChi2>100){ gausMean = h_peak; gausSigma = h_rms; }
1192 
1193       // crystalball limits

1194       double myXmin = gausMean - 3.*gausSigma;
1195       double myXmax = gausMean + 2.*gausSigma;
1196       
1197       // crystalball fit

1198       TF1 *cb_p = new TF1 ("cb_p",crystalball,myXmin,myXmax, 5) ;
1199       cb_p->SetParNames ("Mean","Sigma","alpha","n","Norm","Constant");
1200       cb_p->SetParameter(0, gausMean);
1201       cb_p->SetParameter(1, gausSigma);
1202       cb_p->FixParameter(3, 5.);    
1203       // cb_p->SetParameter(3, 5.);     // solo per caso 500 con BGO e 10% BGO birk        

1204       // cb_p->FixParameter(3, 3.);     // solo x caso 500 no Birk

1205       cb_p->SetParameter(4, gausNorm);
1206       cb_p->SetParLimits(2, 0.1, 5.);
1207       H_eMatrix[ii][myXib]->Fit(cb_p,"lR","",myXmin,myXmax);
1208 
1209       double matrix_gmean      = cb_p->GetParameter(0);
1210       double matrix_gsigma     = cb_p->GetParameter(1); 
1211       double matrix_gmean_err  = cb_p->GetParError(0); 
1212       double matrix_gsigma_err = cb_p->GetParError(1); 
1213 
1214       delete cb_p;
1215       delete gausa;
1216 
1217       // writing into the text file

1218       *ResFileMatrix << beamEne       << "\t" << myXinB[myXib]     << "\t" << myEinB[myXib] << "\t" 
1219              << ii            << "\t" << h_entries         << "\t" 
1220              << h_mean        << "\t" << h_rms             << "\t" 
1221              << matrix_gmean  << "\t" << matrix_gmean_err  << "\t" 
1222              << matrix_gsigma << "\t" << matrix_gsigma_err << endl; 
1223     }
1224   }
1225   
1226   
1227   // Ratios: 1/9, 1/25, 9/25

1228   ofstream *ResFileRatio;
1229   if(appendFlag==0){ 
1230     ResFileRatio = new ofstream(outputFileTxtRatio,ios::out);    
1231     *ResFileRatio << "#Ene"    << "\t" << "xtal"        << "\t" << "eta" << "\t" 
1232           << "sample"  << "\t" << "entries"     << "\t" 
1233           << "H_mean"  << "\t" << "H_rms"       << "\t" 
1234           << "G_mean"  << "\t" << "G_mean err"  << "\t" 
1235           << "G_sigma" << "\t" << "G_sigma err" << endl;} 
1236   else { ResFileRatio = new ofstream(outputFileTxtRatio,ios::app); }    
1237 
1238 
1239   for (int myXib=0; myXib<numXinB; myXib++){ 
1240     for (int ii=0; ii<3; ii++){ 
1241       
1242       // histo parameters

1243       Stat_t h_entries = H_eRatio[ii][myXib]->GetEntries();      
1244       int peakBin      = H_eRatio[ii][myXib]->GetMaximumBin();
1245       double h_norm    = H_eRatio[ii][myXib]->GetMaximum();
1246       double h_rms     = H_eRatio[ii][myXib]->GetRMS();
1247       double h_mean    = H_eRatio[ii][myXib]->GetMean();
1248       double h_peak    = H_eRatio[ii][myXib]->GetBinCenter(peakBin);
1249       
1250       // gaussian fit to initialize

1251       TF1 *gausa = new TF1 ("gausa","[0]*exp(-1*(x-[1])*(x-[1])/2/[2]/[2])",h_peak-10*h_rms,h_peak+10*h_rms);
1252       gausa->SetParameters(h_norm,h_peak,h_rms);
1253       H_eRatio[ii][myXib]->Fit(gausa,"","",h_peak-3*h_rms,h_peak+3*h_rms);
1254       double gausNorm  = gausa->GetParameter(0);
1255       double gausMean  = gausa->GetParameter(1);
1256       double gausSigma = fabs(gausa->GetParameter(2));
1257       double gausChi2  = gausa->GetChisquare()/gausa->GetNDF();
1258       if (gausChi2>100){ gausMean = h_peak; gausSigma = h_rms; }
1259 
1260       // crystalball limits

1261       double myXmin = gausMean - 3.*gausSigma;
1262       double myXmax = gausMean + 2.*gausSigma;
1263       
1264       // crystalball fit

1265       TF1 *cb_p = new TF1 ("cb_p",crystalball,myXmin,myXmax, 5) ;
1266       cb_p->SetParNames ("Mean","Sigma","alpha","n","Norm","Constant");
1267       cb_p->SetParameter(0, gausMean);
1268       cb_p->SetParameter(1, gausSigma);
1269       cb_p->FixParameter(3, 5.);
1270       cb_p->SetParameter(4, gausNorm);
1271       cb_p->SetParLimits(2, 0.1, 5.);
1272       H_eRatio[ii][myXib]->Fit(cb_p,"lR","",myXmin,myXmax);
1273 
1274       double ratio_gmean      = cb_p->GetParameter(0);
1275       double ratio_gsigma     = cb_p->GetParameter(1); 
1276       double ratio_gmean_err  = cb_p->GetParError(0); 
1277       double ratio_gsigma_err = cb_p->GetParError(1); 
1278 
1279       delete cb_p;
1280       delete gausa;
1281 
1282 
1283       // writing into the text file

1284       *ResFileRatio << beamEne      << "\t" << myXinB[myXib]    << "\t" << myEinB[myXib] << "\t" 
1285             << ii           << "\t" << h_entries        << "\t" 
1286             << h_mean       << "\t" << h_rms            << "\t" 
1287             << ratio_gmean  << "\t" << ratio_gmean_err  << "\t" 
1288             << ratio_gsigma << "\t" << ratio_gsigma_err << endl; 
1289     }
1290   }
1291   
1292 
1293   // ---------------------------------------------------------

1294   // Saving histos into files after the fit

1295   strcpy(outfile,outputDir);
1296   strcat(outfile,outputFileName);
1297   strcat(outfile,"_aftFit.root");
1298   if(appendFlag==0){ file = new TFile (outfile,"RECREATE","istogrammi"); }
1299   else             { file = new TFile (outfile,"UPDATE","istogrammi"); }  
1300   for (int myXib=0; myXib<numXinB; myXib++){ 
1301     
1302     for (int ii=0; ii<3; ii++){      
1303       if(ii==0){ sprintf (name, "H_e1_x[%d]",  myXinB[myXib]); }
1304       if(ii==1){ sprintf (name, "H_e9_x[%d]",  myXinB[myXib]); }
1305       if(ii==2){ sprintf (name, "H_e25_x[%d]", myXinB[myXib]); }
1306       H_eMatrix[ii][myXib] -> Write(name);
1307     }
1308     
1309     for (int ii=0; ii<3; ii++){      
1310       if(ii==0){ sprintf (name, "H_e1e9_x[%d]",  myXinB[myXib]); }
1311       if(ii==1){ sprintf (name, "H_e1e25_x[%d]", myXinB[myXib]); }
1312       if(ii==2){ sprintf (name, "H_e9e25_x[%d]", myXinB[myXib]); }
1313       H_eRatio[ii][myXib] -> Write(name);
1314     }    
1315   }
1316   file->Close();
1317   delete file;
1318 
1319 
1320 
1321   
1322   // -----------------------------------------------------------------------------------------

1323   // -----------------------------------------------------------------------------------------

1324   // -----------------------------------------------------------------------------------------

1325   //

1326   //  CONCLUSIONS

1327   //

1328   // -----------------------------------------------------------------------------------------

1329   // -----------------------------------------------------------------------------------------

1330   // -----------------------------------------------------------------------------------------

1331 
1332   // deleting --------------------------

1333   for (int myXib=0; myXib<numXinB; myXib++){
1334     for (int ii=0; ii<3; ii++){       
1335       delete H_eMatrix[ii][myXib];
1336       delete H_eRatio[ii][myXib];
1337     }
1338   }  
1339 
1340   // statistics ----------------------------

1341   cout << endl;
1342   cout << endl;
1343   cout << totalInTheRun        << " total events" << endl;
1344   cout << wrongXtalInTheRun    << " not existing xtal" << endl;
1345   cout << endl;
1346   for (int myXib=0; myXib<numXinB; myXib++){ 
1347     cout << "xtal " << myXinB[myXib] << endl; 
1348     cout << movingTable[myXib]  << " moving table" << endl;
1349     cout << highAmpl[myXib]     << " too high amplitude" << endl;
1350     cout << negaAmpl[myXib]     << " negative amplitude" << endl;
1351     cout << wrongQuality[myXib] << " bad hodoscope quality" << endl;
1352     cout << endl;
1353     cout << good0[myXib]        << " good events at check loop" << endl;
1354     cout << good1[myXib]        << " good events at 1st loop"   << endl;
1355     cout << good2[myXib]        << " good events at 2nd loop"   << endl;
1356     cout << good3[myXib]        << " good events at 3rd loop"   << endl;
1357     cout << "after other selections: " << wrongMaxAmp[myXib] << " xib != max amplitude xtal" << endl; 
1358     cout << endl;
1359   }
1360 
1361 } // end main

1362 
1363 
1364 
1365 
1366 // functions

1367 
1368 // selections

1369 bool selections(int xnum, double eneCent, double qualX, double qualY, int table, int simul){
1370 
1371   bool select = true;
1372   if ( xnum > 1700 )                                    { select = false; }        // problems with the xtal-in-beam

1373   if ( eneCent > 999. )                                 { select = false; }        // problems with amplitudes        

1374   if ( eneCent < 0.5)                                   { select = false; }        // problems with amplitudes

1375   if (!simul){ 
1376     if (table)                                          { select = false; }        // moving table

1377     if ((qualX>2.0) || (qualX<0.0))                     { select = false; }        // quality cut 

1378     if ((qualY>2.0) || (qualY<0.0))                     { select = false; }
1379   }
1380   return select;
1381 }
1382 
1383 
1384 // to skip xtals without 3x3 around

1385 bool atBorder(int xtalNum){
1386 
1387   bool amIatBd = false;
1388   if ( ((xtalNum-1)/20 == 0)  || ((xtalNum-1)/20 == 1)  ){amIatBd = true;}
1389   if ( ((xtalNum-1)/20 == 83) || ((xtalNum-1)/20 == 84) ){amIatBd = true;}
1390   if ( ((xtalNum-1)%20 == 0)  || ((xtalNum-1)%20 == 1)  ){amIatBd = true;}
1391   if ( ((xtalNum-1)%20 == 18) || ((xtalNum-1)%20 == 19) ){amIatBd = true;}
1392   return amIatBd;
1393 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!