EXP1803

{

        TFile* f = new TFile("sim_digi.root");

           TTree* t = (TTree*)f->Get("er");

           t->AddFriend("er","reco.root");

        TCanvas *c1 = new TCanvas("c1", "Telescope Coordinate reconstruction", 1400, 1000);

        TClonesArray* T1_DoubleSi_points = new TClonesArray("ERQTelescopeSiPoint",100);

        TClonesArray* T1_SingleSi_digi = new TClonesArray("ERQTelescopeSiDigi",100);

        TClonesArray* T1_DoubleSi_digiX = new TClonesArray("ERQTelescopeSiDigi",100);

        TClonesArray* T1_DoubleSi_digiY = new TClonesArray("ERQTelescopeSiDigi",100);

        TClonesArray* T1_CsI_digi = new TClonesArray("ERQTelescopeCsIDigi",100);

        TClonesArray* T2_DoubleSi_points = new TClonesArray("ERQTelescopeSiPoint",100);

        TClonesArray* T2_DoubleSi_digiX = new TClonesArray("ERQTelescopeSiDigi",100);

        TClonesArray* T2_DoubleSi_digiY = new TClonesArray("ERQTelescopeSiDigi",100);

        TClonesArray* T2_CsI_digi = new TClonesArray("ERQTelescopeCsIDigi",100);

        t->SetBranchAddress("ERQTelescopeSiPoint_T1_DoubleSi_SD2_XY_0_X",&T1_DoubleSi_points);

        t->SetBranchAddress("ERQTelescopeSiDigi_T1_SingleSi_1_X_0",&T1_SingleSi_digi);

        t->SetBranchAddress("ERQTelescopeSiDigi_T1_DoubleSi_SD2_XY_0_X",&T1_DoubleSi_digiX);

        t->SetBranchAddress("ERQTelescopeSiDigi_T1_DoubleSi_SD2_XY_0_Y",&T1_DoubleSi_digiY);

        t->SetBranchAddress("ERQTelescopeCsIDigi_T1_CsI_1_0",&T1_CsI_digi);

        t->SetBranchAddress("ERQTelescopeSiPoint_T2_DoubleSi_SD2_XY_1_X",&T2_DoubleSi_points);

        t->SetBranchAddress("ERQTelescopeSiDigi_T2_DoubleSi_SD2_XY_1_X",&T2_DoubleSi_digiX);

        t->SetBranchAddress("ERQTelescopeSiDigi_T2_DoubleSi_SD2_XY_1_Y",&T2_DoubleSi_digiY);

        t->SetBranchAddress("ERQTelescopeCsIDigi_T2_CsI_1_1",&T2_CsI_digi);

        TH2F* T1deEHe3 = new TH2F("T1deEHe3","T1 de vs E He3",1000,0,0.25,1000,0,0.0014);

        T1deEHe3->GetXaxis()->SetTitle("E");

        T1deEHe3->GetYaxis()->SetTitle("de");

        TH2F* T1deEH3 = new TH2F("T1deEH3","T1 deE H3",1000,0,0.25,1000,0,0.0014);

        TH2F* T1deEOther = new TH2F("T1deEOther","T1 deE Other",1000,0,0.25,1000,0,0.0014);

        TH2F* T2deEHe3 = new TH2F("T2deEHe3","T2 de vs E He3",1000,0,0.25,1000,0,0.04);

        T2deEHe3->GetXaxis()->SetTitle("E");

        T2deEHe3->GetYaxis()->SetTitle("de");

        TH2F* T2deEH3 = new TH2F("T2deEH3","T2 deE H3",1000,0,0.25,1000,0,0.04);

        TH2F* T2deEOther = new TH2F("T2deEOther","T2 deE Other",1000,0,0.25,1000,0,0.04);

        TCanvas *c2 = new TCanvas("c2", "Telescope de/E on digi info", 1400, 1000);

        c2->Divide(2,1);

        Long64_t nentries = t->GetEntriesFast();

        //Long64_t nentries = 100;

        for (Long64_t i=0; i<nentries; i++) {

          t->GetEntry(i);

          if (T1_DoubleSi_digiX->GetEntriesFast() != 1)

                  continue;

          if (T1_SingleSi_digi->GetEntriesFast() != 1)

                  continue;

          Double_t de = ((ERQTelescopeSiDigi*)T1_SingleSi_digi->At(0))->GetEdep();

          Double_t E = ((ERQTelescopeSiDigi*)T1_SingleSi_digi->At(0))->GetEdep() + 

                                  (((ERQTelescopeSiDigi*)T1_DoubleSi_digiX->At(0))->GetEdep()        +

                                   ((ERQTelescopeSiDigi*)T1_DoubleSi_digiY->At(0))->GetEdep())/2.;

          for (Int_t iDigi = 0; iDigi < T1_CsI_digi->GetEntriesFast(); iDigi++){

                  E+= ((ERQTelescopeCsIDigi*)T1_CsI_digi->At(iDigi))->Edep();

          }

          Bool_t he3 = kFALSE;

          Bool_t h3 = kFALSE;

          for (Int_t iPoint = 0; iPoint < T1_DoubleSi_points->GetEntriesFast(); iPoint++){

                  ERQTelescopeSiPoint* p = (ERQTelescopeSiPoint*)T1_DoubleSi_points->At(iPoint);

                  if (p->GetPDG() == 1000020030)

                          he3 = kTRUE;

                  if (p->GetPDG() == 1000010030)

                          h3 = kTRUE;

          }

          if (he3)

                  T1deEHe3->Fill(E,de);

          else if (h3)

                  T1deEH3->Fill(E,de);

          else

                  T1deEOther->Fill(E,de);

        }

        c2->cd(1);

        T1deEHe3->SetMarkerColor(kBlue);

        T1deEHe3->SetMarkerStyle(21);

        T1deEHe3->Draw();

        T1deEH3->SetMarkerColor(kRed);

        T1deEH3->SetMarkerStyle(21);

        T1deEH3->Draw("same");

        T1deEOther->SetMarkerColor(kGreen);

        T1deEOther->SetMarkerStyle(21);

        T1deEOther->Draw("same");

        auto legend = new TLegend(0.7,0.5,0.95,0.7);

        legend->SetHeader("Ions");

        legend->AddEntry(T1deEHe3,"He3","p");

        legend->AddEntry(T1deEH3,"H3","p");

        legend->AddEntry(T1deEOther,"Others","p");

        legend->Draw();

        for (Long64_t i=0; i<nentries; i++) {

          t->GetEntry(i);

          if (T2_DoubleSi_digiX->GetEntriesFast() != 1)

                  continue;

          Double_t de = (((ERQTelescopeSiDigi*)T2_DoubleSi_digiX->At(0))->GetEdep()        +

                                   ((ERQTelescopeSiDigi*)T2_DoubleSi_digiY->At(0))->GetEdep())/2.;

          Double_t E = de;

          for (Int_t iDigi = 0; iDigi < T2_CsI_digi->GetEntriesFast(); iDigi++){

                  E+= ((ERQTelescopeCsIDigi*)T2_CsI_digi->At(iDigi))->Edep();

          }

          Bool_t he3 = kFALSE;

          Bool_t h3 = kFALSE;

          for (Int_t iPoint = 0; iPoint < T2_DoubleSi_points->GetEntriesFast(); iPoint++){

                  ERQTelescopeSiPoint* p = (ERQTelescopeSiPoint*)T2_DoubleSi_points->At(iPoint);

                  if (p->GetPDG() == 1000020030)

                          he3 = kTRUE;

                  if (p->GetPDG() == 1000010030)

                          h3 = kTRUE;

          }

          if (he3)

                  T2deEHe3->Fill(E,de);

          else if (h3)

                  T2deEH3->Fill(E,de);

          else

                  T2deEOther->Fill(E,de);

        }

        c2->cd(2);

        T2deEHe3->SetMarkerColor(kBlue);

        T2deEHe3->SetMarkerStyle(21);

        T2deEHe3->Draw();

        T2deEH3->SetMarkerColor(kRed);

        T2deEH3->SetMarkerStyle(21);

        T2deEH3->Draw("same");

        T2deEOther->SetMarkerColor(kGreen);

        T2deEOther->SetMarkerStyle(21);

        T2deEOther->Draw("same");

        legend = new TLegend(0.7,0.5,0.95,0.7);

        legend->SetHeader("Ions");

        legend->AddEntry(T2deEHe3,"He3","p");

        legend->AddEntry(T2deEH3,"H3","p");

        legend->AddEntry(T2deEOther,"Others","p");

        legend->Draw();

}