EXP1803

//---------------------------------------------------------

Double_t fitf(Double_t *x,Double_t *par) {  // creating fit function

  Int_t A,Z;

  A = 4; Z = 2;

  Double_t fitval = par[4] + par[0]*( x[0]*(1 - (par[1]*A*Z*Z/x[0])*TMath::Log(1 + (1/(par[1]*A*Z*Z/x[0])) ) ) + par[3]*par[1]*A*Z*Z*TMath::Log( (x[0]+par[1]*A*Z*Z)/(A*par[2]+par[1]*A*Z*Z) ) );

  return fitval;

}

//---------------------------------------------------------

void CsI_fill() { 

    // Reading cal parameters 

  Float_t parXL1[32], parXL2[32],parYL1[16],parYL2[16],parXR1[32], parXR2[32],parYR1[16],parYR2[16];

  //------------------------------------------------------------------------------ 

  // for 1 mm Si detector

  TString line1;

  ifstream myfile1;

  Int_t count=-2;

  myfile1.open("/home/muzalevsky/work/exp1803/data/siCal/SQX_R.cal");

  while (! myfile1.eof() ){

    line1.ReadLine(myfile1);

    if(count < 0){

      count++;

      continue;

    }

    if(line1.IsNull()) break;

    sscanf(line1.Data(),"%g %g", parXR1+count,parXR2+count);

    count++;

  }

  cout << endl << " pars for XR strips" << endl;

  for(Int_t i=0;i<32;i++) cout << parXR1[i] << " " << parXR2[i] << endl;  

  ifstream myfile2;

  TString line2;

  count=-2;

  myfile2.open("/home/muzalevsky/work/exp1803/data/siCal/SQY_R.cal");

  while (! myfile2.eof() ){

    line2.ReadLine(myfile2);

    if(count < 0){

      count++;

      continue;

    }

    if(line2.IsNull()) break;

    sscanf(line2.Data(),"%g %g", parYR1+count,parYR2+count);

    count++;

  }

  cout << endl << " pars for YR strips" << endl;

  for(Int_t i=0;i<16;i++) cout << parYR1[i] << " " << parYR2[i] << endl;   

  ifstream myfile3;

  TString line3;

  count=-2;

  myfile3.open("/home/muzalevsky/work/exp1803/data/siCal/SQX_L.cal");

  while (! myfile3.eof() ){

    line3.ReadLine(myfile3);

    if(count < 0){

      count++;

      continue;

    }

    if(line3.IsNull()) break;

    sscanf(line3.Data(),"%g %g", parXL1+count,parXL2+count);

    count++;

  }

  cout << endl << " pars for XL strips" << endl;

  for(Int_t i=0;i<32;i++) cout << parXL1[i] << " " << parXL2[i] << endl;   

  ifstream myfile4;

  TString line4;

  count=-2;

  myfile4.open("/home/muzalevsky/work/exp1803/data/siCal/SQY_L.cal");

  while (! myfile4.eof() ){

    line4.ReadLine(myfile4);

    if(count < 0){

      count++;

      continue;

    }

    if(line4.IsNull()) break;

    sscanf(line4.Data(),"%g %g", parYL1+count,parYL2+count);

    count++;

  }

  cout << endl << " pars for YL strips" << endl;

  for(Int_t i=0;i<16;i++) cout << parYL1[i] << " " << parYL2[i] << endl;   

  // get pars for positions of times

  cout << endl << " time parameters " << endl;

  ifstream myfile5;

  TString line5;

  Float_t tPx_l1[32],tPx_l2[32];

  myfile5.open("/home/muzalevsky/work/exp1803/macro/dataMuzalevsky/exp201803/calib/Csi/tParX_L.clb");

  count = -2;

  while (! myfile5.eof() ){

    line5.ReadLine(myfile5);

    if(count < 0){

      count++;

      continue;

    }

    if(line5.IsNull()) break;

    sscanf(line5.Data(),"%g %g", tPx_l1+count,tPx_l2+count);

    count++;

  }

  cout << endl << " pars for X time positions" << endl;

  for(Int_t i=0;i<32;i++) cout << tPx_l1[i] << " " << tPx_l2[i] << endl;      

  ifstream myfile6;

  TString line6;

  Float_t tPy_l1[16],tPy_l2[16];

  myfile6.open("/home/muzalevsky/work/exp1803/macro/dataMuzalevsky/exp201803/calib/Csi/tParY_L.clb");

  count = -2;

  while (! myfile6.eof() ){

    line6.ReadLine(myfile6);

    if(count < 0){

      count++;

      continue;

    }

    if(line6.IsNull()) break;

    sscanf(line6.Data(),"%g %g", tPy_l1+count,tPy_l2+count);

    count++;

  }

  cout << endl << " pars for Y time positions" << endl;

  for(Int_t i=0;i<16;i++) cout << tPy_l1[i] << " " << tPy_l2[i] << endl;      

  //------------------------------------------------------------------------------ 

        TFile *f[100]; 

  TString input_file;

  Float_t CsI_R[16],F3[4],tF3[4],F5[4],tF5[4],SQX_R[32],SQY_R[16],SQX_L[32],SQY_L[16],tSQX_R[32],tSQY_R[16],tSQX_L[32],tSQY_L[16],zeroPos,shitf;

  TTree*                 t;

  UShort_t    NeEvent_CsI_R[16],NeEvent_F3[4],NeEvent_tF3[4],NeEvent_F5[4],NeEvent_tF5[4],NeEvent_SQX_R[32],NeEvent_SQY_R[16],NeEvent_SQX_L[32],NeEvent_SQY_L[16],NeEvent_tSQX_R[32],NeEvent_tSQY_R[16],NeEvent_tSQX_L[32],NeEvent_tSQY_L[16];

  TBranch    *b_NeEvent_CsI_R,*b_NeEvent_F3,*b_NeEvent_F5,*b_NeEvent_tF3,*b_NeEvent_tF5,*b_NeEvent_SQX_R,*b_NeEvent_SQY_R,*b_NeEvent_SQX_L,*b_NeEvent_SQY_L,*b_NeEvent_tSQX_R,*b_NeEvent_tSQY_R,*b_NeEvent_tSQX_L,*b_NeEvent_tSQY_L;

  Long64_t nentries1;

  Int_t maxE,multY_L,multX_L,multY_R,multX_R;

  // Creating outfile,outtree

  TFile *fw = new TFile("/home/muzalevsky/work/exp1803/data/exp1804/h5_14/out40cal.root", "RECREATE");

  TTree *tw = new TTree("tree", "data");

  tw->Branch("CsI_R",&CsI_R,"CsI_R[16]/F");

  tw->Branch("F3",&F3,"F3[4]/F");

  tw->Branch("F5",&F5,"F5[4]/F");

  tw->Branch("tF3",&tF3,"tF3[4]/F");

  tw->Branch("tF5",&tF5,"tF5[4]/F");

  tw->Branch("SQX_R",&SQX_R,"SQX_R[32]/F");

  tw->Branch("SQY_R",&SQY_R,"SQY_R[16]/F");

  tw->Branch("SQX_L",&SQX_L,"SQX_L[32]/F");

  tw->Branch("SQY_L",&SQY_L,"SQY_L[16]/F");

  tw->Branch("tSQX_R",&tSQX_R,"tSQX_R[32]/F");

  tw->Branch("tSQY_R",&tSQY_R,"tSQY_R[16]/F");

  tw->Branch("tSQX_L",&tSQX_L,"tSQX_L[32]/F");

  tw->Branch("tSQY_L",&tSQY_L,"tSQY_L[16]/F");

  tw->Branch("multY_L",&multY_L,"multY_L/I");

  tw->Branch("multX_L",&multX_L,"multX_L/I");

  for(Int_t n=10;n<50;n++) {

    input_file.Form("/media/analysis_nas/exp201804/rootfiles/h5_14_00%d.root",n);                

    f[n] = new TFile(input_file.Data());

    if (f[n]->IsZombie()) {

            cerr << "Input file was not opened " << input_file.Data() << endl;

            return 1;

    }

          f[n]->GetObject("AnalysisxTree",t);

          t->SetMakeClass(1);

          t->SetBranchAddress("NeEvent.CsI_R[16]", NeEvent_CsI_R, &b_NeEvent_CsI_R);

          t->SetBranchAddress("NeEvent.F3[4]", NeEvent_F3, &b_NeEvent_F3);

          t->SetBranchAddress("NeEvent.tF3[4]", NeEvent_tF3, &b_NeEvent_tF3);

          t->SetBranchAddress("NeEvent.F5[4]", NeEvent_F5, &b_NeEvent_F5);

          t->SetBranchAddress("NeEvent.tF5[4]", NeEvent_tF5, &b_NeEvent_tF5);

          t->SetBranchAddress("NeEvent.SQX_R[32]", NeEvent_SQX_R, &b_NeEvent_SQX_R);

          t->SetBranchAddress("NeEvent.SQY_R[16]", NeEvent_SQY_R, &b_NeEvent_SQY_R);

          t->SetBranchAddress("NeEvent.SQX_L[32]", NeEvent_SQX_L, &b_NeEvent_SQX_L);

          t->SetBranchAddress("NeEvent.SQY_L[16]", NeEvent_SQY_L, &b_NeEvent_SQY_L);

          t->SetBranchAddress("NeEvent.tSQX_R[32]", NeEvent_tSQX_R, &b_NeEvent_tSQX_R);

          t->SetBranchAddress("NeEvent.tSQY_R[16]", NeEvent_tSQY_R, &b_NeEvent_tSQY_R);

          t->SetBranchAddress("NeEvent.tSQX_L[32]", NeEvent_tSQX_L, &b_NeEvent_tSQX_L);

          t->SetBranchAddress("NeEvent.tSQY_L[16]", NeEvent_tSQY_L, &b_NeEvent_tSQY_L);

          nentries1 = t->GetEntries();

    cout << nentries1 << endl;

    maxE = nentries1;

    cout<<">>> filling TREE up to "<<maxE<< " event"<<endl;

          for (Long64_t jentry=0; jentry<maxE;jentry++) {

                  t->GetEntry(jentry);

      // обнуление

      for(Int_t i = 0; i<32;i++) {

        multY_L = 0;

        multX_L = 0;

        multY_R = 0;

        multX_R = 0;

        SQX_R[i]=0.;

        SQX_L[i]=0.;

        tSQX_R[i]=0.;

        tSQX_L[i]=0.;

        if(i<16) {

          CsI_R[i]=0.;

          SQY_R[i]=0.;

          SQY_L[i]=0.;

          tSQY_R[i]=0.;

          tSQY_L[i]=0.;

        }

        if(i<4){

          F3[i] = 0.;

          F5[i] = 0.;

          tF5[i] = 0.;

          tF3[i] = 0.;

        }

      }

       // selecting nice events

      //if( (NeEvent_F5[0]+NeEvent_F3[0]>400) && (NeEvent_F5[0]+NeEvent_F3[0]<750) && (NeEvent_tF5[0]-NeEvent_tF3[0]>640) && (NeEvent_tF5[0]-NeEvent_tF3[0]<750) && 

      //    (NeEvent_F5[1]+NeEvent_F3[1]>400) && (NeEvent_F5[1]+NeEvent_F3[1]<750) && (NeEvent_tF5[1]-NeEvent_tF3[1]>640) && (NeEvent_tF5[1]-NeEvent_tF3[1]<750) && 

      //    (NeEvent_F5[2]+NeEvent_F3[2]>400) && (NeEvent_F5[2]+NeEvent_F3[2]<750) && (NeEvent_tF5[2]-NeEvent_tF3[2]>640) && (NeEvent_tF5[2]-NeEvent_tF3[2]<750) && 

      //    (NeEvent_F5[3]+NeEvent_F3[3]>400) && (NeEvent_F5[3]+NeEvent_F3[3]<750) && (NeEvent_tF5[3]-NeEvent_tF3[3]>640) && (NeEvent_tF5[3]-NeEvent_tF3[3]<750) ) {

      //}

      for(Int_t i=0; i<32; i++) {

        SQX_R[i] = NeEvent_SQX_R[i]*parXR2[i] + parXR1[i];

        SQX_L[i] = NeEvent_SQX_L[i]*parXL2[i] + parXL1[i];

        if(SQX_L[i]>1.1) multX_L++; 

        if(SQX_R[i]>1.15) multX_R++; 

        tSQX_R[i] = NeEvent_tSQX_R[i];

        tSQX_L[i] = NeEvent_tSQX_L[i]*tPx_l2[i] + tPx_l1[i];

        if(i<16){

          CsI_R[i] = NeEvent_CsI_R[i];

          SQY_R[i] = NeEvent_SQY_R[i]*parYR2[i] + parYR1[i];

          SQY_L[i] = NeEvent_SQY_L[i]*parYL2[i] + parYL1[i];

          if(SQY_L[i]>0.85) multY_L++; 

          if(SQY_R[i]>1.25) multY_R++; 

          tSQY_R[i] = NeEvent_tSQY_R[i];

          tSQY_L[i] = NeEvent_tSQY_L[i]*tPy_l2[i] + tPy_l1[i];

        }

        if(i<4) {

          F3[i] = NeEvent_F3[i];

          F5[i] = NeEvent_F5[i];

          tF3[i] = NeEvent_tF3[i];

          tF5[i] = NeEvent_tF5[i];

        }

      }

                  tw->Fill();                        

          }//entries

  }//nfiles

  fw->cd();

        tw->Write();

        fw->Close();

}