HcalCalibAlgos/macros/MinBias1_Full_Mean_8.9mln.C

0001 #define MinBias1_cxx
0002 #include "MinBias1.h_8.9mln"
0003 #include <TH2.h>
0004 #include <TStyle.h>
0005 #include <TCanvas.h>
0006 #include <iostream.h>
0007 #include <fstream>
0008 #include <sstream>
0009
0010 void MinBias1::Loop()
0011 {
0012 //   In a ROOT session, you can do:
0013 //      Root > .L MinBias1.C
0014 //      Root > MinBias1 t
0015 //      Root > t.GetEntry(12); // Fill t data members with entry number 12
0016 //      Root > t.Show();       // Show values of entry 12
0017 //      Root > t.Show(16);     // Read and show values of entry 16
0018 //      Root > t.Loop();       // Loop on all entries
0019 //
0020
0021 //     This is the loop skeleton where:
0022 //    jentry is the global entry number in the chain
0023 //    ientry is the entry number in the current Tree
0024 //  Note that the argument to GetEntry must be:
0025 //    jentry for TChain::GetEntry
0026 //    ientry for TTree::GetEntry and TBranch::GetEntry
0027 //
0028 //       To read only selected branches, Insert statements like:
0029 // METHOD1:
0030 //    fChain->SetBranchStatus("*",0);  // disable all branches
0031 //    fChain->SetBranchStatus("branchname",1);  // activate branchname
0032 // METHOD2: replace line
0033 //    fChain->GetEntry(jentry);       //read all branches
0034 //by  b_branchname->GetEntry(ientry); //read only this branch
0035    if (fChain == 0) return;
0036
0037    Long64_t nentries = fChain->GetEntriesFast();
0038
0039    Int_t nbytes = 0, nb = 0;
0040
0041    FILE *Out1 = fopen("mean_minbias_8.9mln.txt", "w+");
0042    FILE *Out2 = fopen("mean_noise_8.9mln.txt", "w+");
0043
0044      Float_t minmean_MB[43][73][5][5];
0045      Float_t minmean_NS[43][73][5][5];
0046      Float_t minmean_Diff[43][73][5][5];
0047      Int_t minnevetaphi[43][73][5][5];
0048
0049          Float_t minerr_MB[43][73][5][5];
0050          Float_t minerr_NS[43][73][5][5];
0051          Float_t minerr_Diff[43][73][5][5];
0052
0053      Float_t plmean_MB[43][73][5][5];
0054      Float_t plmean_NS[43][73][5][5];
0055      Float_t plmean_Diff[43][73][5][5];
0056      Int_t plnevetaphi[43][73][5][5];
0057
0058          Float_t plerr_MB[43][73][5][5];
0059          Float_t plerr_NS[43][73][5][5];
0060          Float_t plerr_Diff[43][73][5][5];
0061          Int_t mysubdetpl0[43][73][5][5];
0062          Int_t mysubdetmin0[43][73][5][5];
0063
0064   for(Int_t ietak = 0; ietak < 43; ietak++ )
0065   {
0066    for(Int_t idep = 0; idep < 5; idep++ )
0067    {
0068     for(Int_t isub = 0; isub < 5; isub++ )
0069    {
0070      for(Int_t iphik = 0; iphik < 73; iphik++ )
0071     {
0072
0073      minmean_MB[ietak][iphik][idep][isub] = 0.;
0074      minmean_NS[ietak][iphik][idep][isub] = 0.;
0075      minmean_Diff[ietak][iphik][idep][isub] = 0.;
0076      minnevetaphi[ietak][iphik][idep][isub] = 0;
0077
0078          minerr_MB[ietak][iphik][idep][isub] = 0.;
0079          minerr_NS[ietak][iphik][idep][isub] = 0.;
0080          minerr_Diff[ietak][iphik][idep][isub] = 0.;
0081
0082      plmean_MB[ietak][iphik][idep][isub] = 0.;
0083      plmean_NS[ietak][iphik][idep][isub] = 0.;
0084      plmean_Diff[ietak][iphik][idep][isub] = 0.;
0085      plnevetaphi[ietak][iphik][idep][isub] = 0;
0086
0087          plerr_MB[ietak][iphik][idep][isub] = 0.;
0088          plerr_NS[ietak][iphik][idep][isub] = 0.;
0089          plerr_Diff[ietak][iphik][idep][isub] = 0.;
0090
0091     }
0092     }
0093    }
0094   }
0095
0096    TH1F* h1noise = new TH1F ("h1noise","h1noise",100, -0.1,0.4);
0097    TH1F* h1noise_mean_1 = new TH1F ("h1noise_mean_1","h1noise_mean_1",100, -0.1,0.4);
0098    TH1F* h1noise_mean_2 = new TH1F ("h1noise_mean_2","h1noise_mean_2",100, -0.1,0.4);
0099    TH1F* h1noise_mean_3 = new TH1F ("h1noise_mean_3","h1noise_mean_3",100, -0.1,0.4);
0100    TH1F* h1noise_mean_4 = new TH1F ("h1noise_mean_4","h1noise_mean_4",100, -0.1,0.4);
0101    TH1F* h1noise_mean_6 = new TH1F ("h1noise_mean_6","h1noise_mean_6",100, -0.1,0.4);
0102
0103    TH2F* hlnoise = new TH2F ("hlnoise","hlnoise",28, -14, 14, 72, 0.5, 72.5);
0104    TH2F* hlnoise1 = new TH2F ("hlnoise1","hlnoise1",28, -14, 14, 72, 0.5, 72.5);
0105    TH2F* hlnoise2 = new TH2F ("hlnoise2","hlnoise2",28, -14, 14, 72, 0.5, 72.5);
0106    TH2F* hlnoise3 = new TH2F ("hlnoise3","hlnoise3",28, -14, 14, 72, 0.5, 72.5);
0107    TH2F* hlnoise4 = new TH2F ("hlnoise4","hlnoise4",28, -14, 14, 72, 0.5, 72.5);
0108    TH2F* hlnoise_stab4 = new TH2F ("hlnoise_stab4","hlnoise_stab4",1501, -0.5, 1500.5, 100, -0.1, 0.4);
0109    TH2F* hlnoise_stab3 = new TH2F ("hlnoise_stab3","hlnoise_stab3",1501, -0.5, 1500.5, 100, -0.1, 0.4);
0110    TH2F* hlnoise_stab2 = new TH2F ("hlnoise_stab2","hlnoise_stab2",1501, -0.5, 1500.5, 100, -0.1, 0.4);
0111    TH2F* hlnoise_stab1 = new TH2F ("hlnoise_stab1","hlnoise_stab1",1501, -0.5, 1500.5, 100, -0.1, 0.4);
0112    TH2F* hlnoise_stab6 = new TH2F ("hlnoise_stab6","hlnoise_stab6",1501, -0.5, 1500.5, 100, -0.1, 0.4);
0113
0114    TH2F* hlnoise_stab_HF_d1 = new TH2F ("hlnoise_stab_HF_d1","hlnoise_stab_HF_d1",1501, -0.5, 1500.5, 100, -0.1, 0.4);
0115    TH2F* hlnoise_stab_HF_d2 = new TH2F ("hlnoise_stab_HF_d2","hlnoise_stab_HF_d2",1501, -0.5, 1500.5, 100, -0.1, 0.4);
0116
0117
0118
0119    TH1F* h1signal = new TH1F ("h1signal","h1signal",100, -0.1,0.4);
0120    TH2F* hlsignal = new TH2F ("hlsignal","hlsignal",28, -14, 14, 72, 0.5, 72.5);
0121    TH2F* hlsignal1 = new TH2F ("hlsignal1","hlsignal1",28, -14, 14, 72, 0.5, 72.5);
0122    TH2F* hlsignal2 = new TH2F ("hlsignal2","hlsignal2",28, -14, 14, 72, 0.5, 72.5);
0123    TH2F* hlsignal3 = new TH2F ("hlsignal3","hlsignal3",28, -14, 14, 72, 0.5, 72.5);
0124    TH2F* hlsignal4 = new TH2F ("hlsignal4","hlsignal4",28, -14, 14, 72, 0.5, 72.5);
0125
0126
0127    Float_t myfile = 0.;
0128
0129
0130    for (Long64_t jentry=0; jentry<nentries;jentry++) {
0131       Long64_t ientry = LoadTree(jentry);
0132       if (ientry < 0) break;
0133       nb = fChain->GetEntry(jentry);   nbytes += nb;
0134       // if (Cut(ientry) < 0) continue;
0135
0136       if(ientry == 0) {cout<<" New file "<<myfile<<endl; myfile++;}
0137
0138       if(ieta == 0) continue;
0139 //      if( ieta<0) {
0140 //      cout<<" Mysubd "<<mysubd<<" "<<ieta<<" "<<" "<<iphi<<" "<<depth<<" "<<" "<<mom0_MB<<endl;
0141 //      cout<<" Mom "<<mom2_Diff<<" "<<mom2_MB<<" "<<mom2_Noise<<endl;
0142 //      cout<<" Mom "<<mom1_Diff<<" "<<mom1_MB<<" "<<mom1_Noise<<endl;
0143 //      }
0144       if(mom0_MB == 0) continue;
0145       Float_t mean_MB = mom1_MB/mom0_MB;
0146       Float_t mean_NS = mom1_Noise/mom0_MB;
0147
0148       Float_t disp_MB = mom2_MB/mom0_MB - (mom1_MB/mom0_MB)*(mom1_MB/mom0_MB);
0149       Float_t disp_NS = mom2_Noise/mom0_MB - (mom1_Noise/mom0_MB)*(mom1_Noise/mom0_MB);
0150       Float_t disp_Diff = mom2_Diff/mom0_MB - (mom1_Diff/mom0_MB)*(mom1_Diff/mom0_MB);
0151
0152       if(mysubd == 1)
0153       {
0154       if(abs(ieta)<15) {
0155           h1noise->Fill(disp_NS);
0156           if(disp_NS>0.22 && disp_NS<0.26) hlnoise4->Fill((Float_t)ieta,(Float_t)iphi,1.);
0157           if(disp_NS>0.16 && disp_NS<0.18) hlnoise3->Fill((Float_t)ieta,(Float_t)iphi,1.);
0158           if(disp_NS>0.12 && disp_NS<0.14) hlnoise2->Fill((Float_t)ieta,(Float_t)iphi,1.);
0159           if(disp_NS>0.11 && disp_NS<0.12) hlnoise1->Fill((Float_t)ieta,(Float_t)iphi,1.);
0160       if(disp_NS>0.11) hlnoise->Fill((Float_t)ieta,(Float_t)iphi,disp_NS);
0161
0162           h1signal->Fill(disp_MB);
0163           if(disp_NS>0.22 && disp_NS<0.26) {hlsignal4->Fill((Float_t)ieta,(Float_t)iphi,1.);}
0164           if(disp_NS>0.16 && disp_NS<0.18) {hlsignal3->Fill((Float_t)ieta,(Float_t)iphi,1.);}
0165           if(disp_NS>0.12 && disp_NS<0.14) {hlsignal2->Fill((Float_t)ieta,(Float_t)iphi,1.);}
0166           if(disp_NS>0.11 && disp_NS<0.12) {hlsignal1->Fill((Float_t)ieta,(Float_t)iphi,1.);}
0167       if(disp_NS>0.11) hlsignal->Fill((Float_t)ieta,(Float_t)iphi,disp_MB);
0168
0169       }
0170       if( ieta == -1 && iphi == 70 ) {h1noise_mean_1->Fill(mean_NS);hlnoise_stab1->Fill(myfile,disp_NS,1.);}
0171       if( ieta == -2 && iphi == 70 ) {h1noise_mean_2->Fill(mean_NS);hlnoise_stab2->Fill(myfile,disp_NS,1.);}
0172       if( ieta == -3 && iphi == 70 ) {h1noise_mean_3->Fill(mean_NS);hlnoise_stab3->Fill(myfile,disp_NS,1.);}
0173       if( ieta == -4 && iphi == 70 ) {h1noise_mean_4->Fill(mean_NS);hlnoise_stab4->Fill(myfile,disp_NS,1.);}
0174       if( ieta == -6 && iphi == 70 ) {h1noise_mean_6->Fill(mean_NS);hlnoise_stab6->Fill(myfile,disp_NS,1.);}
0175       }
0176
0177       if( mysubd == 4 ){
0178        if(depth == 1) hlnoise_stab_HF_d1->Fill(myfile,disp_NS,1.);
0179        if(depth == 2) hlnoise_stab_HF_d2->Fill(myfile,disp_NS,1.);
0180       }
0181
0182        if(ieta<0) {
0183 // Calculation of dispersion ===============================
0184
0185          mysubdetmin0[abs(ieta)][iphi][depth][mysubd] = mysubd;
0186
0187      minmean_MB[abs(ieta)][iphi][depth][mysubd]=minmean_MB[abs(ieta)][iphi][depth][mysubd] + (float)mean_MB;
0188      minmean_NS[abs(ieta)][iphi][depth][mysubd]=minmean_NS[abs(ieta)][iphi][depth][mysubd] + (float)mean_NS;
0189      minmean_Diff[abs(ieta)][iphi][depth][mysubd]=minmean_Diff[abs(ieta)][iphi][depth][mysubd] + (float)disp_Diff;
0190      minnevetaphi[abs(ieta)][iphi][depth][mysubd]++;
0191
0192          minerr_MB[abs(ieta)][iphi][depth][mysubd] = minerr_MB[abs(ieta)][iphi][depth][mysubd] + mean_MB*mean_MB;
0193          minerr_NS[abs(ieta)][iphi][depth][mysubd] = minerr_NS[abs(ieta)][iphi][depth][mysubd] + mean_NS*mean_NS;
0194          minerr_Diff[abs(ieta)][iphi][depth][mysubd] = minerr_Diff[abs(ieta)][iphi][depth][mysubd] + disp_Diff*disp_Diff;
0195
0196 // ==========================================================
0197
0198       }
0199       if(ieta>=0) {
0200
0201 // Calculation of dispersion ===============================
0202          mysubdetpl0[abs(ieta)][iphi][depth][mysubd] = mysubd;
0203      plmean_MB[abs(ieta)][iphi][depth][mysubd]=plmean_MB[abs(ieta)][iphi][depth][mysubd] + (float)mean_MB;
0204      plmean_NS[abs(ieta)][iphi][depth][mysubd]=plmean_NS[abs(ieta)][iphi][depth][mysubd] + (float)mean_NS;
0205      plmean_Diff[abs(ieta)][iphi][depth][mysubd]=plmean_Diff[abs(ieta)][iphi][depth][mysubd] + (float)disp_Diff;
0206      plnevetaphi[abs(ieta)][iphi][depth][mysubd]++;
0207
0208          plerr_MB[abs(ieta)][iphi][depth][mysubd] = plerr_MB[abs(ieta)][iphi][depth][mysubd] + mean_MB*mean_MB;
0209          plerr_NS[abs(ieta)][iphi][depth][mysubd] = plerr_NS[abs(ieta)][iphi][depth][mysubd] + mean_NS*mean_NS;
0210          plerr_Diff[abs(ieta)][iphi][depth][mysubd] = plerr_Diff[abs(ieta)][iphi][depth][mysubd] + disp_Diff*disp_Diff;
0211
0212
0213 // ==========================================================
0214
0215
0216       }
0217
0218
0219    } // jentry
0220
0221    cout<<" Finish cycle "<<endl;
0222
0223   Double_t plmeang_mean_MB = 0.;
0224   Double_t plmeang_mean_NS = 0.;
0225   Double_t plmeang_mean_Diff = 0.;
0226   Double_t plmeang_mean_Diff_av = 0.;
0227
0228   Double_t minmeang_mean_MB = 0.;
0229   Double_t minmeang_mean_NS = 0.;
0230   Double_t minmeang_mean_Diff = 0.;
0231   Double_t minmeang_mean_Diff_av = 0.;
0232
0233 // ieta>0
0234   for (int idep = 1; idep <5; idep++ )
0235   {
0236   for(int ietak = 1; ietak < 43; ietak++ )
0237   {
0238   for(int isub = 1; isub < 5; isub++ )
0239   {
0240     for(Int_t iphik = 1; iphik < 73; iphik++ )
0241     {
0242          if(mysubdetpl0[ietak][iphik][idep][isub] == 0) continue;
0243           if(plnevetaphi[ietak][iphik][idep][isub] == 0) {
0244             Float_t rr = 1.;
0245         Float_t ss = 0.;
0246         fprintf(Out1,"%d %d %d %d %.8f\n",mysubdetpl0[ietak][iphik][idep][isub],ietak,iphik,idep,ss);
0247         fprintf(Out2,"%d %d %d %d %.8f\n",mysubdetpl0[ietak][iphik][idep][isub],ietak,iphik,idep,ss);
0248 //      fprintf(Out3,"%d %d %d %d %.8f\n",mysubdetpl0[ietak][iphik][idep][isub],ietak,iphik,idep,ss);
0249 //      fprintf(Out4,"%d %d %d %d %.8f\n",mysubdetpl0[ietak][iphik][idep][isub],ietak,iphik,idep,ss);
0250             continue;
0251            }
0252
0253 // Mean dispersion
0254        plmeang_mean_MB = plmean_MB[ietak][iphik][idep][isub]/plnevetaphi[ietak][iphik][idep][isub];
0255        plmeang_mean_NS = plmean_NS[ietak][iphik][idep][isub]/plnevetaphi[ietak][iphik][idep][isub];
0256        plmeang_mean_Diff = plmean_Diff[ietak][iphik][idep][isub]/plnevetaphi[ietak][iphik][idep][isub];
0257        plmeang_mean_Diff_av = plmeang_mean_MB-plmeang_mean_NS;
0258
0259        fprintf(Out1,"%d %d %d %d %.8f\n",mysubdetpl0[ietak][iphik][idep][isub],ietak,iphik,idep,plmeang_mean_MB);
0260        fprintf(Out2,"%d %d %d %d %.8f\n",mysubdetpl0[ietak][iphik][idep][isub],ietak,iphik,idep,plmeang_mean_NS);
0261 //       fprintf(Out3,"%d %d %d %d %.8f\n",mysubdetpl0[ietak][iphik][idep][isub],ietak,iphik,idep,plmeang_mean_Diff);
0262 //       fprintf(Out4,"%d %d %d %d %.8f\n",mysubdetpl0[ietak][iphik][idep][isub],ietak,iphik,idep,plmeang_mean_Diff_av);
0263
0264     }
0265   }
0266   }
0267   }
0268 // ieta<0
0269   for (int idep = 1; idep <5; idep++ )
0270   {
0271   for(int ietak = 1; ietak < 43; ietak++ )
0272   {
0273     Int_t iietak = -1.*ietak;
0274   for(int isub = 1; isub < 5; isub++ )
0275   {
0276     for(Int_t iphik = 1; iphik < 73; iphik++ )
0277     {
0278         if(mysubdetmin0[ietak][iphik][idep][isub] == 0) continue;
0279           if(minnevetaphi[ietak][iphik][idep][isub] == 0) {
0280             Float_t rr = 1.;
0281         Float_t ss = 0.;
0282         fprintf(Out1,"%d %d %d %d %.8f\n",mysubdetmin0[ietak][iphik][idep][isub],iietak,iphik,idep,ss);
0283         fprintf(Out2,"%d %d %d %d %.8f\n",mysubdetmin0[ietak][iphik][idep][isub],iietak,iphik,idep,ss);
0284 //      fprintf(Out3,"%d %d %d %d %.8f\n",mysubdetmin0[ietak][iphik][idep][isub],iietak,iphik,idep,ss);
0285 //      fprintf(Out4,"%d %d %d %d %.8f\n",mysubdetmin0[ietak][iphik][idep][isub],iietak,iphik,idep,ss);
0286             continue;
0287            }
0288
0289 // Mean dispersion
0290        minmeang_mean_MB = minmean_MB[ietak][iphik][idep][isub]/minnevetaphi[ietak][iphik][idep][isub];
0291        minmeang_mean_NS = minmean_NS[ietak][iphik][idep][isub]/minnevetaphi[ietak][iphik][idep][isub];
0292        minmeang_mean_Diff = minmean_Diff[ietak][iphik][idep][isub]/minnevetaphi[ietak][iphik][idep][isub];
0293
0294        fprintf(Out1,"%d %d %d %d %.8f\n",mysubdetmin0[ietak][iphik][idep][isub],iietak,iphik,idep,minmeang_mean_MB);
0295        fprintf(Out2,"%d %d %d %d %.8f\n",mysubdetmin0[ietak][iphik][idep][isub],iietak,iphik,idep,minmeang_mean_NS);
0296 //       fprintf(Out3,"%d %d %d %d %.8f\n",mysubdetmin0[ietak][iphik][idep][isub],iietak,iphik,idep,minmeang_mean_Diff);
0297 //       fprintf(Out4,"%d %d %d %d %.8f\n",mysubdetmin0[ietak][iphik][idep][isub],iietak,iphik,idep,minmeang_mean_Diff_av);
0298
0299     }
0300   }
0301   }
0302   }
0303
0304   TFile f("signal_noise_mean_8.9mln.root","recreate");
0305   h1noise->Write();
0306   hlnoise->Write();
0307   hlnoise1->Write();
0308   hlnoise2->Write();
0309   hlnoise3->Write();
0310   hlnoise4->Write();
0311
0312   h1signal->Write();
0313   hlsignal->Write();
0314   hlsignal1->Write();
0315   hlsignal2->Write();
0316   hlsignal3->Write();
0317   hlsignal4->Write();
0318
0319   h1noise_mean_1->Write();
0320   h1noise_mean_2->Write();
0321   h1noise_mean_3->Write();
0322   h1noise_mean_4->Write();
0323   h1noise_mean_6->Write();
0324   hlnoise_stab4->Write();
0325   hlnoise_stab3->Write();
0326   hlnoise_stab2->Write();
0327   hlnoise_stab1->Write();
0328   hlnoise_stab6->Write();
0329   hlnoise_stab_HF_d1->Write();
0330   hlnoise_stab_HF_d2->Write();
0331
0332   fclose(Out1);
0333   fclose(Out2);
0334 //  fclose(Out3);
0335 //  fclose(Out4);
0336 }