#define Selector_cxx // The class definition in Selector.h has been generated automatically // by the ROOT utility TTree::MakeSelector(). This class is derived // from the ROOT class TSelector. For more information on the TSelector // framework see $ROOTSYS/README/README.SELECTOR or the ROOT User Manual. // The following methods are defined in this file: // Begin(): called every time a loop on the tree starts, // a convenient place to create your histograms. // SlaveBegin(): called after Begin(), when on PROOF called only on the // slave servers. // Process(): called for each event, in this function you decide what // to read and fill your histograms. // SlaveTerminate: called at the end of the loop on the tree, when on PROOF // called only on the slave servers. // Terminate(): called at the end of the loop on the tree, // a convenient place to draw/fit your histograms. // // To use this file, try the following session on your Tree T: // // root> T->Process("Selector.C") // root> T->Process("Selector.C","some options") // root> T->Process("Selector.C+") // #include "Selector.h" #include #include //================================SETTINGS ======== const int chanmax=8000; const int bins2d=800; //const long eventsmax=2e+10; const long eventsmax=1e+6; //calibration-------------EXL 19 - 1st is zero double a[19]={0., 0.668983,0.519247,0.5649, 0.600488,0.579928,0.598377, 0.540216,0.488096,0.610848, 0.560869, 0.503247,0.464832, 0.848081,0.792648,0.799, 0.426236,0.485055,0.492694}; double b[19]={0., -71.1150, 132.36, 71.2, -58.9582, 72.8547, -56.2703, 91.2014, 95.1007, 9.26351, -2.42136, 86.8976, 105.457, 45.8667, 68.9514,2.84729, 18.4952, 18.7380, 116.032}; //thresholds EXL double tr[19]={0, 140.,140.,140., 140.,140.,140., 140.,140.,140., 140.,140.,140., 140.,140.,140., 140.,140.,140. }; // Experimental factor R1ex from MTX: includes gain G1 // the theoretical factor R1 should be used in correction // higher positive slope line - like exl1,3,5 signals double R1ex[9]={ 0.14,0.18,0.15, 0.16,0.18,0.16, 0.15,0.42,0.17 }; // lower positive slope line - like exl2,4,6 signal double R2ex[9]={ 0.16,0.15,0.14, 0.15,0.15,0.16, 0.16,0.15,0.16 }; //--- we compute gain ratio later double G1comp[9]; // arbitrarym good around 1.0 // we recalculate later double G2gain[9]={ 2616., 2368., 2639., 2717., 2728., 4308., 1822., 3168., 2796. }; double G2gain_factor=3000.; //=============================SETTTING END _------- void Selector::Begin(TTree * /*tree*/) { // The Begin() function is called at the start of the query. // When running with PROOF Begin() is only called on the client. // The tree argument is deprecated (on PROOF 0 is passed). char ch1[100]; char ch2[100]; TString option = GetOption(); printf("beginning the process\n%s",""); // HISTO DEFINITIONS------------------- for (int i=0;i<18;i++){ // SINGLES sprintf(ch1, "exl%02d", i+1 ); sprintf(ch2, "%s SINGLES", ch1 ); exl[i]=new TH1F(ch1,ch2,chanmax,0,chanmax); } for (int i=0;i<18;i++){ // SINGLES CALIB sprintf(ch1, "exlcal%02d", i+1 ); sprintf(ch2, "%s SINGLES CALIB", ch1 ); exlcal[i]=new TH1F(ch1,ch2,chanmax,0,chanmax); } exl_singlesum=new TH1F("exlcalSUM","TOTAL SUM of SINGLES / useful4 online",chanmax,0,chanmax); exl_pairsum=new TH1F("pairSUM","TOTAL SUM of summed pairs/ wrong concept",chanmax,0,chanmax); complete_pairs=new TH1F("complpairs","pair hits in the matrix",18,0,18); incomplete_pairs=new TH1F("cincomppairs","single hits in pairs in the matrix",18,0,18); // BIDIM-------------------------------------- for (int i=0;i<9;i++){ // PAIRS int j1=i*2+1; int j2=i*2+2; sprintf(ch1, "pair%02d%02d", j1, j2); sprintf(ch2, "%s PAIRS", ch1 ); expair[i]=new TH2F(ch1,ch2,bins2d,0,chanmax,bins2d,0,chanmax); } for (int i=0;i<9;i++){ // PAIRS CALIB int j1=i*2+1; int j2=i*2+2; sprintf(ch1, "paircal%02d%02d", j1, j2); sprintf(ch2, "%s PAIRS CALIB", ch1 ); expaircal[i]=new TH2F(ch1,ch2,bins2d,0,chanmax,bins2d,0,chanmax); } for (int i=0;i<9;i++){ // SUMPAIRS SINGLE int j1=i*2+1; int j2=i*2+2; sprintf(ch1, "paircalsum%02d%02d", j1, j2); sprintf(ch2, "%s SINGLES = PAIRS CALIBrated AND SUMMED, energy restored", ch1 ); exlsumpair[i]=new TH1F(ch1,ch2,chanmax,0,chanmax); } for (int i=0;i<6;i++){ // HNEIGH int j1=i*2+2; if (j1>=6){ j1=j1+2;} if (j1>=12){ j1=j1+2;} sprintf(ch1, "hneig%02d%02d", j1,j1+1 ); sprintf(ch2, "%s HNEIGHBOR", ch1 ); exhneig[i]=new TH2F(ch1,ch2,bins2d,0,chanmax,bins2d,0,chanmax); } for (int i=0;i<12;i++){ // VNEIGH int j1=i+1; sprintf(ch1, "vneigh%02d%02d", j1,j1+6 ); sprintf(ch2, "%s VNEIGHBOR", ch1 ); exvneig[i]=new TH2F(ch1,ch2,bins2d,0,chanmax,bins2d,0,chanmax); } //compensation for GAIN1 - odd crystals for (int i=0;i<9;i++){ // 2/1 bad sums G1comp[i]=TMath::Sqrt( R2ex[i]/R1ex[i] )*1.; // G2==1 // G1comp[i]=TMath::Sqrt( R1ex[i]/R2ex[i] )*1.; // G2==1 wrong printf("%d %f\n",i,G1comp[i] ); } } //==========BEGIN END------------------------------------------- void Selector::SlaveBegin(TTree * /*tree*/) { // The SlaveBegin() function is called after the Begin() function. // When running with PROOF SlaveBegin() is called on each slave server. // The tree argument is deprecated (on PROOF 0 is passed). TString option = GetOption(); } // invert matrix // for every event from PAIR matrix : restore the energy //------------------------- // // Bool_t Selector::inverse_mtx( int i9 ){ // i take raw channel int ii=2*i9+1; int jj=2*i9+2; double e1=EXL->GetEnergyCrystalNumber( ii ); double e2=EXL->GetEnergyCrystalNumber( jj ); if (e1<=tr[ii]) return kFALSE; if (e2<=tr[jj]) return kFALSE; // we try to substract the offset first // a*k+b c*(k-d) => b/a=d e1=e1-b[ii]/a[ii]; e2=e2-b[jj]/a[jj]; // determine Gain1 from alpha=beta double G2=G2gain[i9]/G2gain_factor; double G1=TMath::Sqrt( R2ex[i9]/R1ex[i9] )* G2; double R= R2ex[i9]*G2/G1; // get to the R1==R2 (alpha==beta) double alpha=R/(R+1.0);// from Zamoram /easy double beta=alpha; // both axes have the same leak double fact=G1*(1-alpha)*G2*(1-beta)-G1*G2*alpha*beta; fact=1.0/fact; //matrix" double N1=fact*(G2*(1-beta)*e1 - G1*alpha *e2); double N2=fact*(G1*(1-alpha)*e2 - G2*beta *e1); expaircal[i9]->Fill( N1,N2 ); return kTRUE; } Bool_t Selector::is_in_pair( int ii, int jj ){ double e1=EXL->GetEnergyCrystalNumber( ii ); double e2=EXL->GetEnergyCrystalNumber( jj ); if (e1<=tr[ii]) return kFALSE; if (e2<=tr[jj]) return kFALSE; // e1=a[ii]+b[ii]; // e2=a[jj]+b[jj]; return kTRUE; } Bool_t Selector::Process(Long64_t entry) { // The Process() function is called for each entry in the tree (or possibly // keyed object in the case of PROOF) to be processed. The entry argument // specifies which entry in the currently loaded tree is to be processed. // When processing keyed objects with PROOF, the object is already loaded // and is available via the fObject pointer. // // This function should contain the \"body\" of the analysis. It can contain // simple or elaborate selection criteria, run algorithms on the data // of the event and typically fill histograms. // // The processing can be stopped by calling Abort(). // // Use fStatus to set the return value of TTree::Process(). // // The return value is currently not used. double e,e1,e2, f,g; fReader.SetEntry(entry); count=count+1; if (count> eventsmax){ if (count==eventsmax+1)printf("Free ride now%s\n",""); //Abort("just abort"); }else{ //================================== START========== //================================== START========== //================================== START========== //================================== START========== //================================== START========== if (count%10000==0){ printf("%7ld k %d %f\n",count/1000, *PPT_Q, EXL->GetEnergyCrystalNumber(1) ); for ( int ii=0;ii<9;ii++){ if (is_in_pair( ii*2+1,ii*2+2 )){ printf(" %d\n",ii*2+1 ); } } } //---------------------------------------------------------- for (int i=0;i<18;i++){ //SINGLES e=EXL->GetEnergyCrystalNumber(i+1); if ((e>0)&& (e>tr[i+1]) ){ exl[i]->Fill( e ); } }//for 18 for (int i=0;i<18;i++){ //SINGLES CALIB WORKS but concept is wrong e=EXL->GetEnergyCrystalNumber(i+1); if ((e>0)&& (e>tr[i+1]) ){ f=gRandom->Rndm(); exlcal[i]->Fill( (e+f)*a[i+1]+b[i+1] ); // Fill the calibrated SINGLE exl_singlesum->Fill( (e+f)*a[i+1]+b[i+1] ); //SUM of singles (online) } }//for 18 //--------------------------------------------// PAIRS for (int i=0;i<9;i++){ // PAIRS e1=EXL->GetEnergyCrystalNumber(i*2+1); e2=EXL->GetEnergyCrystalNumber(i*2+2); f=gRandom->Rndm(); g=gRandom->Rndm(); if ( (e1>tr[i*2+1])&&(e2>tr[i*2+2]) ){ expair[i]->Fill( e1,e2 ); // single /* expaircal[i]->Fill( e1*a[i*2+1]+b[i*2+1],e2*a[i*2+2]+b[i*2+2] ); //CALIB2D //pairs summed ==== CASE OF PAIRS CALIB SUM double esumpair=(e1+f)*a[i*2+1]+b[i*2+1]+(e2+g)*a[i*2+2]+b[i*2+2]; exlsumpair[i]->Fill( esumpair ); // sums of pairs exl_pairsum->Fill( esumpair ); // sums of all sum of pairs */ inverse_mtx(i); //expaircal[i]->Fill( e1*G1comp[i],e2 ); //CALIB2D; now G1 evald //pairs summed ==== CASE OF PAIRS CALIB SUM double esumpair=(e1+f)*G1comp[i]+e2; exlsumpair[i]->Fill( esumpair ); // sums of pairs exl_pairsum->Fill( esumpair ); // sums of all sum of pairs //Statistics---------- complete_pairs->Fill( 2*i+1 ); // stat - complete pairs complete_pairs->Fill( 2*i+2 ); //stat - complete pairs }else if(e1>tr[i*2+1]){ incomplete_pairs->Fill( 2*i+1 ); //stat - incomplete pairs }else if(e2>tr[i*2+2]){ incomplete_pairs->Fill( 2*i+2 ); //stat - incomplete pairs } }//for 18 for (int i=0;i<9;i++){ // PAIRS CALIB e1=EXL->GetEnergyCrystalNumber(i*2+1); e2=EXL->GetEnergyCrystalNumber(i*2+2); if ( (e1>tr[i*2+1]) && (e2>tr[i*2+2]) ){ } }//for 18 for (int i=0;i<6;i++){ //hneigh int j1=i*2+2; if (j1>=6){ j1=j1+2;} if (j1>=12){ j1=j1+2;} e1=EXL->GetEnergyCrystalNumber(j1); e2=EXL->GetEnergyCrystalNumber(j1+1); if ((e1>tr[j1] )&&(e2>tr[j1+1])){ exhneig[i]->Fill( e1,e2 ); } }//for 6 for (int i=0;i<12;i++){ //vneig e1=EXL->GetEnergyCrystalNumber(i+1); e2=EXL->GetEnergyCrystalNumber(i+1+6); if ((e1>tr[i+1])&&(e2>tr[i+1+6])){ exvneig[i]->Fill( e1,e2 ); } }//for 12 } //==================================END=============== //=============================================== //D=============================================END=============== //=============================================== //D=============================================END=============== //=============================================== //D=============================================END=============== return kTRUE; } void Selector::SlaveTerminate() { // The SlaveTerminate() function is called after all entries or objects // have been processed. When running with PROOF SlaveTerminate() is called // on each slave server. } void Selector::Terminate() { // The Terminate() function is the last function to be called during // a query. It always runs on the client, it can be used to present // the results graphically or save the results to file. printf("terminating the processm saving spe\n%s",""); TFile *fo=new TFile("whatever_output.root","RECREATE"); //NEW for (int i=0;i<18;i++){ exl[i]->Write(); exlcal[i]->Write(); } for (int i=0;i<9;i++){ expair[i]->Write(); expaircal[i]->Write(); exlsumpair[i]->Write(); // and now i am exactly where before } for (int i=0;i<6;i++){ exhneig[i]->Write(); } for (int i=0;i<12;i++){ exvneig[i]->Write(); } exl_singlesum->Write(); complete_pairs->Write(); incomplete_pairs->Write(); exl_pairsum->Write(); fo->Close(); }