Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-06 08:14:00 

0001 #define Analysis_cxx
0002 #include "Analysis.hh"
0003 
0004 //////////////////////////////////////////////////////////////////////
0005 // Constructor                                                      //
0006 //////////////////////////////////////////////////////////////////////
0007 Analysis::Analysis( int run ) :
0008   fChain(0),
0009   run_( run )
0010 {
0011 // if parameter tree is not specified (or zero), connect the file
0012 // used to generate this class and read the Tree.
0013   TTree* tree = 0;
0014   //if (tree == 0)
0015     {
0016       string file_name = "results/tracking_run" + to_string( run_ ) + ".root";
0017       cout << file_name << endl;
0018       // TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject( file_name.c_str() );
0019       TFile *f = new TFile( file_name.c_str() );
0020       f->GetObject("clus_tree",tree);
0021       
0022     }
0023     
0024   Init( tree );
0025 }
0026 
0027 //////////////////////////////////////////////////////////////////////
0028 // Destructor                                                       //
0029 //////////////////////////////////////////////////////////////////////
0030 Analysis::~Analysis()
0031 {
0032    if (!fChain)
0033      return;
0034 
0035    delete fChain->GetCurrentFile();
0036 }
0037 
0038 //////////////////////////////////////////////////////////////////////
0039 // Init                                                             //
0040 //////////////////////////////////////////////////////////////////////
0041 void Analysis::Init(TTree *tree)
0042 {
0043    // The Init() function is called when the selector needs to initialize
0044    // a new tree or chain. Typically here the branch addresses and branch
0045    // pointers of the tree will be set.
0046    // It is normally not necessary to make changes to the generated
0047    // code, but the routine can be extended by the user if needed.
0048    // Init() will be called many times when running on PROOF
0049    // (once per file to be processed).
0050 
0051    // Set object pointer
0052    x_in = 0;
0053    y_in = 0;
0054    z_in = 0;
0055    r_in = 0;
0056    size_in = 0;
0057    phi_in = 0;
0058    theta_in = 0;
0059    adc_in = 0;
0060    is_associated_in = 0;
0061    track_incoming_theta_in = 0;
0062    x_out = 0;
0063    y_out = 0;
0064    z_out = 0;
0065    r_out = 0;
0066    size_out = 0;
0067    phi_out = 0;
0068    theta_out = 0;
0069    adc_out = 0;
0070    is_associated_out = 0;
0071    track_incoming_theta_out = 0;
0072 
0073    // Set branch addresses and branch pointers
0074    if( !tree )
0075      return;
0076    
0077    fChain = tree;
0078    fCurrent = -1;
0079    fChain->SetMakeClass(1);
0080 
0081    fChain->SetBranchAddress("evt_clus"          , &evt_clus, &b_evt_clus                      );
0082    fChain->SetBranchAddress("x_in"          , &x_in, &b_x_in                          );
0083    fChain->SetBranchAddress("y_in"          , &y_in, &b_y_in                          );
0084    fChain->SetBranchAddress("z_in"          , &z_in, &b_z_in                          );
0085    fChain->SetBranchAddress("r_in"          , &r_in, &b_r_in                          );
0086    fChain->SetBranchAddress("size_in"           , &size_in, &b_size_in                        );
0087    fChain->SetBranchAddress("phi_in"            , &phi_in, &b_phi_in                          );
0088    fChain->SetBranchAddress("theta_in"          , &theta_in, &b_theta_in                      );
0089    fChain->SetBranchAddress("adc_in"            , &adc_in, &b_adc_in                          );
0090    fChain->SetBranchAddress("is_associated_in"      , &is_associated_in, &b_is_associated_in              );
0091    fChain->SetBranchAddress("track_incoming_theta_in"   , &track_incoming_theta_in, &b_track_incoming_theta_in        );
0092    fChain->SetBranchAddress("x_out"         , &x_out, &b_x_out                        );
0093    fChain->SetBranchAddress("y_out"         , &y_out, &b_y_out                        );
0094    fChain->SetBranchAddress("z_out"         , &z_out, &b_z_out                        );
0095    fChain->SetBranchAddress("r_out"         , &r_out, &b_r_out                        );
0096    fChain->SetBranchAddress("size_out"          , &size_out, &b_size_out                      );
0097    fChain->SetBranchAddress("phi_out"           , &phi_out, &b_phi_out                        );
0098    fChain->SetBranchAddress("theta_out"         , &theta_out, &b_theta_out                    );
0099    fChain->SetBranchAddress("adc_out"           , &adc_out, &b_adc_out                        );
0100    fChain->SetBranchAddress("is_associated_out"     , &is_associated_out, &b_is_associated_out            );
0101    fChain->SetBranchAddress("track_incoming_theta_out"  , &track_incoming_theta_out, &b_track_incoming_theta_out      );
0102    fChain->SetBranchAddress("z_vertex"          , &z_vertex, &b_z_vertex                      );
0103    Notify();
0104    
0105   string title = "ADC distribution;ADC [arb. units];Counts [arb. units]";
0106 
0107   hist_all_ = new MipHist( "Raw", title );  
0108   hist_all_->SetColorAlpha( kGreen+2, 0.2 );
0109   hist_all_->SetTag( "Raw ADC" );
0110 
0111   hist_aso_ = new MipHist( "Track_assoiation", title );
0112   hist_aso_->SetColorAlpha( kRed, 0.1 );
0113   hist_aso_->SetTag( "Tracklet asso." );
0114 
0115   hist_no_aso_ = new MipHist( "Track_not_assoiated", title );
0116   hist_no_aso_->SetColorAlpha( kBlack, 0.1 );
0117   hist_no_aso_->SetTag( "No tracklet asso." );
0118 
0119   ///////////////////////////////////////////////////////////////
0120   // Misaki's config
0121   ///////////////////////////////////////////////////////////////
0122   hist_ang85_ = new MipHist( "Ang85_90", title );
0123   hist_ang85_->SetColorAlpha( kBlue, 0.1 );
0124   hist_ang85_->SetTag( "85#circ < |#theta^{INTT}| #leq 90#circ" );
0125 
0126   hist_ang45_ = new MipHist( "Ang40_45", title );
0127   hist_ang45_->SetColorAlpha( kRed, 0.1 );
0128   hist_ang45_->SetTag( "40#circ < |#theta^{INTT}| #leq 45#circ" );
0129 
0130   hist_ang35_ = new MipHist( "Ang30_35", title );
0131   hist_ang35_->SetColorAlpha( kSpring-1, 0.1 );
0132   hist_ang35_->SetTag( "30#circ < |#theta^{INTT}| #leq 35#circ" );
0133 
0134   hist_ang25_ = new MipHist( "Ang20_25", title );
0135   hist_ang25_->SetColorAlpha( kOrange+1, 0.1 );
0136   hist_ang25_->SetTag( "20#circ < |#theta^{INTT}| #leq 25#circ" );
0137 
0138   ///////////////////////////////////////////////////////////////
0139   // my config                                                 //
0140   ///////////////////////////////////////////////////////////////
0141   hist_ang80_ = new MipHist( "Ang80_90", title );
0142   hist_ang80_->SetColorAlpha( kRed, 0.1 );
0143   hist_ang80_->SetTag( "80#circ < |#theta^{INTT}| #leq 90#circ" );
0144 
0145   hist_ang70_ = new MipHist( "Ang70_80", title );
0146   hist_ang70_->SetColorAlpha( kRed, 0.1 );
0147   hist_ang70_->SetTag( "70#circ < |#theta^{INTT}| #leq 80#circ" );
0148 
0149   hist_ang60_ = new MipHist( "Ang60_70", title );
0150   hist_ang60_->SetColorAlpha( kRed, 0.1 );
0151   hist_ang60_->SetTag( "60#circ < |#theta^{INTT}| #leq 70#circ" );
0152 
0153   hist_ang50_ = new MipHist( "Ang50_60", title );
0154   hist_ang50_->SetColorAlpha( kRed, 0.1 );
0155   hist_ang50_->SetTag( "50#circ < |#theta^{INTT}| #leq 60#circ" );
0156 
0157   hist_ang40_ = new MipHist( "Ang40_50", title );
0158   hist_ang40_->SetColorAlpha( kRed, 0.1 );
0159   hist_ang40_->SetTag( "40#circ < |#theta^{INTT}| #leq 50#circ" );
0160 
0161   hist_ang30_ = new MipHist( "Ang30_40", title );
0162   hist_ang30_->SetColorAlpha( kRed, 0.1 );
0163   hist_ang30_->SetTag( "30#circ < |#theta^{INTT}| #leq 40#circ" );
0164 
0165   hist_ang20_ = new MipHist( "Ang20_30", title );
0166   hist_ang20_->SetColorAlpha( kRed, 0.1 );
0167   hist_ang20_->SetTag( "20#circ < |#theta^{INTT}| #leq 30#circ" );
0168 
0169   hist_ang10_ = new MipHist( "Ang10_20", title );
0170   hist_ang10_->SetColorAlpha( kRed, 0.1 );
0171   hist_ang10_->SetTag( "10#circ < |#theta^{INTT}| #leq 20#circ" );
0172 
0173   hist_ang0_ = new MipHist( "Ang0_10", title );
0174   hist_ang0_->SetColorAlpha( kRed, 0.1 );
0175   hist_ang0_->SetTag( "0#circ < |#theta^{INTT}| #leq 10#circ" );
0176 
0177   ///////////////////////////////////////////////////////////////
0178   // fine selection                                            //
0179   ///////////////////////////////////////////////////////////////
0180   hist_ang10_11_ = new MipHist( "Ang10_11", title );
0181   hist_ang10_11_->SetColorAlpha( kRed, 0.1 );
0182   hist_ang10_11_->SetTag( "10#circ < |#theta^{INTT}| #leq 11#circ" );
0183 
0184   hist_ang20_21_ = new MipHist( "Ang20_21", title );
0185   hist_ang20_21_->SetColorAlpha( kRed, 0.1 );
0186   hist_ang20_21_->SetTag( "20#circ < |#theta^{INTT}| #leq 21#circ" );
0187 
0188   hist_ang30_31_ = new MipHist( "Ang30_31", title );
0189   hist_ang30_31_->SetColorAlpha( kRed, 0.1 );
0190   hist_ang30_31_->SetTag( "30#circ < |#theta^{INTT}| #leq 31#circ" );
0191 
0192   ///////////////////////////////////////////////////////////////
0193   // put all into a vector                                     //
0194   ///////////////////////////////////////////////////////////////
0195   hists_.push_back( hist_all_ );
0196   hists_.push_back( hist_aso_ );
0197   hists_.push_back( hist_no_aso_ );
0198   hists_.push_back( hist_ang85_ );
0199   hists_.push_back( hist_ang45_ );
0200   hists_.push_back( hist_ang35_ );
0201   hists_.push_back( hist_ang25_ );
0202 
0203   hists_.push_back( hist_ang80_ );
0204   hists_.push_back( hist_ang70_ );
0205   hists_.push_back( hist_ang60_ );
0206   hists_.push_back( hist_ang50_ );
0207   hists_.push_back( hist_ang40_ );
0208   hists_.push_back( hist_ang30_ );
0209   hists_.push_back( hist_ang20_ );
0210   hists_.push_back( hist_ang10_ );
0211   hists_.push_back( hist_ang0_ );
0212   
0213   hists_.push_back( hist_ang10_11_ );
0214   hists_.push_back( hist_ang20_21_ );
0215   hists_.push_back( hist_ang30_31_ );
0216 
0217   hist_association_ = new TH1D( "association_ratio", "Ratio of associated hits", 100, 0, 1 );
0218   HistSetting( hist_association_ );
0219 
0220   hist_association_in_ = new TH1D( "association_ratio_in", "Ratio of associated hits", 100, 0, 1 );
0221   HistSetting( hist_association_in_ );
0222 
0223   hist_association_out_ = new TH1D( "association_ratio_out", "Ratio of associated hits", 100, 0, 1 );
0224   HistSetting( hist_association_out_ );
0225 
0226   hist_association_in_out_ = new TH2D( "association_ratio_in_out", "Ratio of associated hits", 100, 0, 1, 100, 0, 1 );
0227   HistSetting( hist_association_in_out_ );
0228 
0229   hist_correlation_ = new TH2D( "nhit_correlation_barrel",
0230                    "Cluster correlation between barrels;Number of clusters at inner barrel;Number of clusters at outer barrel;",
0231                    100, 0, 100, 
0232                    100, 0, 100 ); 
0233   HistSetting( hist_correlation_ );
0234 }
0235 
0236 //////////////////////////////////////////////////////////////////////
0237 // private functions                                                //
0238 //////////////////////////////////////////////////////////////////////
0239 string Analysis::GetDate()
0240 {
0241   return "9/13/2024";
0242   
0243   TDatime dt;
0244   int year  = dt.GetYear();
0245   int month = dt.GetMonth();
0246   int day   = dt.GetDay();
0247 
0248   // format: mm/dd/yyyy
0249   std::stringstream ss;
0250   ss << month << "/" << day << "/" << year;
0251 
0252   return ss.str();
0253 }
0254 
0255 void Analysis::DrawWords()
0256 {
0257 
0258   ///////////////////////////////////////////////////////////////////////////////
0259   // Writting words in the canvas                                              //
0260   ///////////////////////////////////////////////////////////////////////////////
0261   TLatex* tex = new TLatex();
0262 
0263   double line_height = 0.05;
0264   double first_margin = 0.005;  
0265   double pos_y = 1.0 - top_margin_ + first_margin;// - line_height;
0266 
0267   // Date
0268   tex->DrawLatexNDC( 0.7, pos_y,
0269              string("#it{" + GetDate() + "}").c_str() );
0270 
0271   // sPHENIX Internal or sPHENIX Prelimnary
0272   pos_y -= line_height - first_margin + 0.025;
0273   double pos_x = 0.45 - right_margin_;
0274   if( is_preliminary_ == false )
0275     {
0276       tex->DrawLatexNDC( pos_x, pos_y, "#it{#bf{sPHENIX}} Internal" );
0277     }
0278   else
0279     {
0280       //pos_x = 0.6;
0281       tex->DrawLatexNDC( pos_x, pos_y, "#it{#bf{sPHENIX}} Preliminary" );
0282     }
0283 
0284   // p+p 200 GeV
0285   pos_y -= line_height;
0286   tex->DrawLatexNDC( pos_x, pos_y, "Run-24 #it{p+p} 200 GeV" );
0287 
0288   pos_y -= line_height;
0289   tex->DrawLatexNDC( pos_x, pos_y, "Run 50889 " );
0290 
0291   pos_y -= line_height;
0292   tex->DrawLatexNDC( pos_x, pos_y, "INTT Streaming readout " );
0293 }
0294 
0295 Int_t Analysis::GetEntry(Long64_t entry)
0296 {
0297 // Read contents of entry.
0298    if (!fChain)
0299      return 0;
0300    
0301    return fChain->GetEntry(entry);
0302 }
0303 
0304 Long64_t Analysis::LoadTree(Long64_t entry)
0305 {
0306 // Set the environment to read one entry
0307    if (!fChain)
0308      return -5;
0309    
0310    Long64_t centry = fChain->LoadTree(entry);
0311    if (centry < 0)
0312      return centry;
0313    
0314    if (fChain->GetTreeNumber() != fCurrent)
0315      {
0316        fCurrent = fChain->GetTreeNumber();
0317        Notify();
0318      }
0319    
0320    return centry;
0321 }
0322 
0323 
0324 Bool_t Analysis::Notify()
0325 {
0326    // The Notify() function is called when a new file is opened. This
0327    // can be either for a new TTree in a TChain or when when a new TTree
0328    // is started when using PROOF. It is normally not necessary to make changes
0329    // to the generated code, but the routine can be extended by the
0330    // user if needed. The return value is currently not used.
0331 
0332    return kTRUE;
0333 }
0334 
0335 Int_t Analysis::Cut(Long64_t entry)
0336 {
0337 // This function may be called from Loop.
0338 // returns  1 if entry is accepted.
0339 // returns -1 otherwise.
0340    return 1;
0341 }
0342 
0343 TBox* Analysis::DrawExpectedPeakRegion( double theta_min, double theta_max )
0344 {
0345   double pos_min = GetExpectedPeakPosition( theta_max );
0346   double pos_max = GetExpectedPeakPosition( theta_min );
0347 
0348   double alpha = 0.3;
0349   if( pos_max - pos_min < 3 )
0350     {
0351       pos_max = pos_min + 3;
0352       alpha = 0.8;
0353     }
0354   else if( pos_max - pos_min < 10 )
0355     {
0356       alpha = 0.8;
0357     }
0358 
0359   TBox* box = new TBox( pos_min, 0, pos_max, 0.3 );
0360   box->SetFillColorAlpha( kAzure-1, alpha );
0361   box->Draw();
0362   return box;
0363 }
0364 
0365 void Analysis::DrawSingle( MipHist* mip_hist )
0366 {
0367   int color = mip_hist->GetColor();
0368   double alpha = mip_hist->GetAlpha();
0369 
0370   mip_hist->SetColorAlpha( kBlack, 0.1, true );
0371   
0372   HistSetting1D( mip_hist->GetHist() );
0373   mip_hist->GetHist()->Draw();
0374   this->DrawWords();
0375 
0376   string page_title = string("Title: ") + mip_hist->GetHist()->GetName();
0377   c_->Print( c_->GetName(), page_title.c_str() );
0378 
0379   string output = c_->GetName();
0380   output = output.substr( 0, output.find_last_of("_") ) + "_" + mip_hist->GetHist()->GetName();
0381   if( is_preliminary_ == true )
0382     output += "_preliminary.pdf";
0383   else
0384     output += "_internal.pdf";
0385     
0386   c_->Print( output.c_str() );
0387 
0388   // restore color and alpha
0389   mip_hist->SetColorAlpha( color, alpha, true );
0390 }
0391 
0392 void Analysis::DrawMultiple( vector < MipHist* >& mip_hists, string output_tag, bool use_color_palette )
0393 {
0394 
0395   TLegend* leg = new TLegend( 0.5, 0.65 - mip_hists.size() * 0.05, 0.7, 0.7);
0396   int counter = 0;
0397   for( auto& mip_hist : mip_hists )
0398     {
0399       //double top_val = 0.05 * ++counter + 1;
0400       double top_val = ++counter ;
0401       
0402       auto hist = mip_hist->GetNormalizedHist( top_val );
0403       HistSetting1D( hist );
0404 
0405       hist->GetYaxis()->SetRangeUser( 0, 0.3 );
0406 
0407       if( use_color_palette )
0408     {
0409       hist->SetFillStyle( 0 );
0410       hist->SetFillColorAlpha( kWhite, 0 );
0411       hist->Draw( (mip_hist == mip_hists[0] ? "HISTE PLC" : "HISTE PLC same" ) );
0412     }
0413       else 
0414     hist->Draw( (mip_hist == mip_hists[0] ? "HISTE" : "HISTE same" ) );
0415 
0416       leg->AddEntry( hist, mip_hist->GetTag().c_str() );
0417       // auto f = mip_hist->GetNormalizedFunction( top_val );
0418       // f->Draw( "same" );
0419       // mip_hist->DrawLine( f );
0420       
0421     }
0422 
0423   leg->Draw();
0424   this->DrawWords();
0425   c_->Print( c_->GetName() );  
0426 
0427   string output = c_->GetName();
0428   output = output.substr( 0, output.find_last_of("_") ) + "_" + output_tag;
0429   if( is_preliminary_ == true )
0430     output += "_preliminary.pdf";
0431   else
0432     output += "_internal.pdf";
0433 
0434   c_->Print( output.c_str() );
0435 }
0436 
0437 void Analysis::DrawMultiPanel( vector < MipHist* >& mip_hists, string output_tag )
0438 {
0439   TCanvas* c = new TCanvas( "aaa", "title", 2500, 1000);
0440   c->Clear();
0441   c->Divide( 4, 2 );
0442   
0443   double top_val = 0.4;
0444 
0445   TBox* box;
0446   for( int i=0; i<mip_hists.size(); i++ )
0447     {
0448       c->cd(i+1);
0449       auto mip_hist = mip_hists[i];
0450 
0451       // values to restore configuration
0452       int color = mip_hist->GetColor();
0453       double alpha = mip_hist->GetAlpha();
0454 
0455       mip_hist->SetColorAlpha( kBlack, 0.1, true );
0456       auto hist = mip_hist->GetNormalizedHist( top_val );
0457       hist->GetYaxis()->SetRangeUser( 0, 0.3 );
0458       
0459       HistSetting1D( hist );
0460       hist->GetXaxis()->SetLabelSize( 0.07 );
0461       hist->GetXaxis()->SetTitleSize( 0.07 );
0462       hist->GetXaxis()->SetTitleOffset( 1 );
0463 
0464       hist->GetYaxis()->SetLabelSize( 0.07 );
0465       hist->GetYaxis()->SetTitleSize( 0.07 );
0466       hist->GetYaxis()->SetTitleOffset( 1.2 );
0467       
0468       hist->Draw( "HISTE" );
0469       // TLegend* leg = new TLegend( 0.5, 0.65, 0.92, 0.92);
0470       // leg->AddEntry( hist, mip_hist->GetTag().c_str() );
0471       // leg->Draw();
0472 
0473       TLatex* tex = new TLatex();
0474       tex->SetTextSize( 0.1 );
0475       tex->DrawLatexNDC( 0.35, 0.8, mip_hist->GetTag().c_str() );
0476       
0477       double theta_min = stoi( string(hist->GetName()).substr( 3, 2 ) );
0478       double theta_max = stoi( string(hist->GetName()).substr( 6, 2 ) );
0479       box = DrawExpectedPeakRegion( theta_min, theta_max );
0480       
0481       // restore color and alpha
0482       mip_hist->SetColorAlpha( color, alpha, true );
0483     }
0484 
0485   {
0486     c->cd(8);
0487     ///////////////////////////////////////////////////////////////////////////////
0488     // Writting words in the canvas                                              //
0489     ///////////////////////////////////////////////////////////////////////////////
0490     TLatex* tex = new TLatex();
0491     tex->SetTextSize( 0.1 );
0492     double line_height = 0.1;
0493     double first_margin = 0.005;  
0494     double pos_y = 0.85; //1.0 - top_margin_ + first_margin;// - line_height;
0495     double pos_x = 0.0; // 0.45 - right_margin_;
0496 
0497     // Date
0498     tex->DrawLatexNDC( pos_x, pos_y,
0499                string("#it{" + GetDate() + "}").c_str() );
0500 
0501     // sPHENIX Internal or sPHENIX Prelimnary
0502     pos_y -= line_height - first_margin + 0.025;
0503     if( is_preliminary_ == false )
0504       {
0505     tex->DrawLatexNDC( pos_x, pos_y, "#it{#bf{sPHENIX}} Internal" );
0506       }
0507     else
0508       {
0509     //pos_x = 0.6;
0510     tex->DrawLatexNDC( pos_x, pos_y, "#it{#bf{sPHENIX}} Preliminary" );
0511       }
0512 
0513     // p+p 200 GeV
0514     pos_y -= line_height;
0515     tex->DrawLatexNDC( pos_x, pos_y, "Run-24, #it{p+p} 200 GeV" );
0516   
0517     pos_y -= line_height;
0518     tex->DrawLatexNDC( pos_x, pos_y, "Run 50889 " );
0519   
0520     pos_y -= line_height;
0521     tex->DrawLatexNDC( pos_x, pos_y, "INTT Streaming readout " );
0522 
0523     pos_y -= line_height;
0524     //    TLegend* leg = new TLegend( pos_x, pos_y );
0525     TLegend* leg = new TLegend( pos_x, 0.1, pos_x + 0.6, pos_y + 0.01 );
0526     leg->SetTextSize( 0.1 );
0527     leg->AddEntry( box, "Expected peak region", "f" );
0528     leg->Draw();
0529   }
0530 
0531   
0532   string page_title = string("Title: Multi panel");// + mip_hist->GetHist()->GetName();
0533   c_->Print( c_->GetName(), page_title.c_str() );
0534   string output = c_->GetName();
0535   output = output.substr( 0, output.find_last_of("_") ) + "_" + output_tag;
0536   if( is_preliminary_ == true )
0537     output += "_preliminary.pdf";
0538   else
0539     output += "_internal.pdf";
0540 
0541   c->Print( output.c_str() );
0542 
0543 }
0544 
0545 void Analysis::FillClusterInfo( int mode )
0546 {
0547 
0548   vector<double>  *x;
0549   vector<double>  *y;
0550   vector<double>  *z;
0551   vector<double>  *r;
0552   vector<int>     *size;
0553   vector<double>  *phi;
0554   vector<double>  *theta;
0555   vector<double>  *adcs;
0556   vector<bool>    *is_associated;
0557   vector<double>  *track_incoming_theta;
0558 
0559   if( mode == 0 )
0560     {
0561       x = x_in;
0562       y = y_in;
0563       z = z_in;
0564       r = r_in;
0565       size = size_in;
0566       phi = phi_in;
0567       theta = theta_in;
0568       adcs = adc_in;
0569       is_associated = is_associated_in;
0570       track_incoming_theta = track_incoming_theta_in;
0571     }
0572   else if( mode == 1 )
0573     {
0574       x = x_out;
0575       y = y_out;
0576       z = z_out;
0577       r = r_out;
0578       size = size_out;
0579       phi = phi_out;
0580       theta = theta_out;
0581       adcs = adc_out;
0582       is_associated = is_associated_out;
0583       track_incoming_theta = track_incoming_theta_out;
0584     }
0585 
0586   for( int i=0; i<adcs->size(); i++ )
0587     {
0588       
0589       if( (*size)[i] > cluster_size_max_ )
0590     {
0591       counter_too_large_cluster_++;
0592       continue;
0593     }
0594 
0595       /* if( fabs( z_vertex - (*z)[i] ) > 5 ) */ // no need...
0596       /*   continue; */
0597              
0598       int adc = (*adcs)[i];
0599       bool is_single_hit_cluster_adc7  = ( (*size)[i] == 1 && adc == adc7_ );
0600 
0601       // Modification of bins containig 1 hit with ADC7, No need to use it because it's not fair modification.
0602       // bool is_double_hit_cluster_adc7  = ( (*size)[i] == 2 
0603       //                       && (adc == 245 || adc == 255 || adc == 390 )
0604       //                       );
0605       // is_single_hit_cluster_adc7 = is_single_hit_cluster_adc7 || is_double_hit_cluster_adc7;
0606       
0607       bool is_double_hit_cluster_adc14 = ( (*size)[i] == 2 && adc == adc7_ * 2 );
0608       hist_all_->FillAll( adc, is_single_hit_cluster_adc7, is_double_hit_cluster_adc14 );
0609     
0610       // if this cluster is not associated with a tracklet, store info here and skip the rest parts
0611       if( (*is_associated)[i] == false )
0612     {
0613 
0614       hist_no_aso_->FillAll( adc, is_single_hit_cluster_adc7, is_double_hit_cluster_adc14 );
0615       continue;
0616     }
0617        
0618       hist_aso_->FillAll( adc, is_single_hit_cluster_adc7, is_double_hit_cluster_adc14 );
0619       ///////////////////////////////////////////////////////////////
0620       // Misaki's config                                           //
0621       ///////////////////////////////////////////////////////////////
0622       double angle_abs = fabs( (*track_incoming_theta)[i] );
0623       
0624       if( 85 <= angle_abs )
0625     {
0626       hist_ang85_->FillAll( adc, is_single_hit_cluster_adc7, is_double_hit_cluster_adc14 );
0627     }
0628       else if( 40 < angle_abs && angle_abs <= 45 )
0629     {
0630       hist_ang45_->FillAll( adc, is_single_hit_cluster_adc7, is_double_hit_cluster_adc14 );
0631     }
0632       else if( 30 < angle_abs && angle_abs <= 35 )
0633     {
0634       hist_ang35_->FillAll( adc, is_single_hit_cluster_adc7, is_double_hit_cluster_adc14 );
0635     }
0636       else if( 20 < angle_abs && angle_abs <= 25 )
0637     {
0638       hist_ang25_->FillAll( adc, is_single_hit_cluster_adc7, is_double_hit_cluster_adc14 );
0639     }
0640   
0641       ///////////////////////////////////////////////////////////////
0642       // My config                                                 //
0643       ///////////////////////////////////////////////////////////////
0644       if( 80 < angle_abs && angle_abs <= 90 )
0645     {
0646       hist_ang80_->FillAll( adc, is_single_hit_cluster_adc7, is_double_hit_cluster_adc14 );
0647     }
0648       else if( 70 < angle_abs && angle_abs <= 80 )
0649     {
0650       hist_ang70_->FillAll( adc, is_single_hit_cluster_adc7, is_double_hit_cluster_adc14 );
0651     }
0652       else if( 60 < angle_abs && angle_abs <= 70 )
0653     {
0654       hist_ang60_->FillAll( adc, is_single_hit_cluster_adc7, is_double_hit_cluster_adc14 );
0655     }
0656       else if( 50 < angle_abs && angle_abs <= 60 )
0657     {
0658       hist_ang50_->FillAll( adc, is_single_hit_cluster_adc7, is_double_hit_cluster_adc14 );
0659     }
0660       else if( 40 < angle_abs && angle_abs <= 50 )
0661     {
0662       hist_ang40_->FillAll( adc, is_single_hit_cluster_adc7, is_double_hit_cluster_adc14 );
0663     }
0664       else if( 30 < angle_abs && angle_abs <= 40 )
0665     {
0666       hist_ang30_->FillAll( adc, is_single_hit_cluster_adc7, is_double_hit_cluster_adc14 );
0667     }
0668       else if( 20 < angle_abs && angle_abs <= 30 )
0669     {
0670       hist_ang20_->FillAll( adc, is_single_hit_cluster_adc7, is_double_hit_cluster_adc14 );
0671     }
0672       else if( 10 < angle_abs && angle_abs <= 20 )
0673     {
0674       hist_ang10_->FillAll( adc, is_single_hit_cluster_adc7, is_double_hit_cluster_adc14 );
0675     }
0676       else if( 0 < angle_abs && angle_abs <= 10 )
0677     {
0678       hist_ang0_->FillAll( adc, is_single_hit_cluster_adc7, is_double_hit_cluster_adc14 );
0679     }
0680       else
0681     cout << (*track_incoming_theta)[i] << endl;
0682 
0683       ///////////////////////////////////////////////////////////////
0684       // Fine selection                                            //
0685       ///////////////////////////////////////////////////////////////
0686       if( 10 < angle_abs && angle_abs <= 11 )
0687     {
0688       hist_ang10_11_->FillAll( adc, is_single_hit_cluster_adc7, is_double_hit_cluster_adc14 );
0689     }
0690       else if( 20 < angle_abs && angle_abs <= 21 )
0691     {
0692       hist_ang20_21_->FillAll( adc, is_single_hit_cluster_adc7, is_double_hit_cluster_adc14 );
0693     }
0694       else if( 30 < angle_abs && angle_abs <= 31 )
0695     {
0696       hist_ang30_31_->FillAll( adc, is_single_hit_cluster_adc7, is_double_hit_cluster_adc14 );
0697     }
0698 
0699     } // end of for( int i=0; i<adc->size(); i++ )
0700 
0701 }
0702 
0703 
0704 double Analysis::GetExpectedPeakPosition( double theta )
0705 {
0706   double thickness = 320; // um
0707   double path_length = thickness / TMath::Sin( TMath::Pi() * theta / 180 );
0708   //cout << "Path length: " << path_length << endl;
0709   
0710   double energy_loss_per_cm = 1.15; // MeV/g/cm/cm
0711   double density = 2.329; // g/cm^3
0712   double energy_loss = path_length * energy_loss_per_cm * density / 10; // keV
0713   //cout << "Energy loss: " << energy_loss << " keV" << endl;
0714   
0715   double energy_electron_hole = 3.62; // eV/pair
0716   double elementaly_charge = 1.602e-19; // C
0717   double charge = energy_loss * 1000 / energy_electron_hole * elementaly_charge / 1e-15; // fC
0718   //cout << "Charge: " << charge << " fC" << endl;
0719   
0720   double pulse_height = charge * 100 + 210;
0721   double dac = (pulse_height - 210 ) / 4;
0722 
0723   return dac;
0724 }
0725 
0726 ///////////////////////////////////////////////////////////////////////
0727 // public functions                                                  //
0728 ///////////////////////////////////////////////////////////////////////
0729 
0730 void Analysis::Loop()
0731 {
0732   //     This is the loop skeleton where:
0733   //    jentry is the global entry number in the chain
0734   //    ientry is the entry number in the current Tree
0735   //  Note that the argument to GetEntry must be:
0736   //    jentry for TChain::GetEntry
0737   //    ientry for TTree::GetEntry and TBranch::GetEntry
0738   
0739   if( fChain == 0 )
0740     return;
0741 
0742    Long64_t nentries = fChain->GetEntriesFast();
0743 
0744    bool does_z_cut = true;
0745    
0746    for (Long64_t jentry=0; jentry<nentries;jentry++)
0747      {
0748        Long64_t ientry = LoadTree(jentry);
0749 
0750        int hit_num_inner = std::count( is_associated_in->begin(), is_associated_in->end(), true );
0751        int hit_num_outer = std::count( is_associated_out->begin(), is_associated_out->end(), true );
0752        hist_correlation_->Fill( hit_num_inner, hit_num_outer );
0753        
0754        if (ientry < 0)
0755      break;
0756 
0757        fChain->GetEntry(jentry);
0758        
0759        double cluster_num_in = adc_in->size();
0760        double cluster_num_out = adc_out->size();
0761        double cluster_num = cluster_num_in + cluster_num_out;
0762        double cluster_num_asymmetry = fabs(cluster_num_in - cluster_num_out) / cluster_num;
0763        double associated_cluster_num_in = count( is_associated_in->begin(), is_associated_in->end(), true );
0764        double associated_cluster_num_out = count( is_associated_out->begin(), is_associated_out->end(), true );
0765 
0766        double associated_cluster_num = associated_cluster_num_in + associated_cluster_num_out;
0767 
0768        double associated_ratio = associated_cluster_num / cluster_num;
0769        double associated_ratio_in = associated_cluster_num_in / cluster_num_in;
0770        double associated_ratio_out = associated_cluster_num_out / cluster_num_out;
0771 
0772        hist_association_->Fill( associated_ratio );
0773        hist_association_in_out_->Fill( associated_ratio_in, associated_ratio_out );
0774        hist_association_in_->Fill( associated_ratio_in );
0775        hist_association_out_->Fill( associated_ratio_out );
0776 
0777        bool is_good_event = true;
0778        if( cluster_num_in <= 3
0779        || cluster_num_out <= 3
0780        || cluster_num <= 6
0781        //|| cluster_num_asymmetry > 0.1
0782        )
0783      is_good_event = false;
0784 
0785        if( associated_ratio > 1 || associated_ratio_in > 1 || associated_ratio_out > 1 )
0786      cout << associated_cluster_num << " / " << cluster_num << "=" << associated_ratio << "\t"
0787           << associated_ratio_in << "\t"
0788           << associated_ratio_out << "\t"
0789           << endl;
0790 
0791        double associated_threshold = 0.4;
0792        if( associated_ratio_in  < associated_threshold || 
0793        associated_ratio_out < associated_threshold || 
0794        associated_ratio     < associated_threshold )
0795      is_good_event = false;
0796               
0797        ////////////////////////////////////////////////////////////////////
0798        // Loop over inner hits
0799        if( is_good_event == false )
0800      continue;
0801 
0802        // cout << cluster_num << "\t"
0803        //       << associated_cluster_num << "\t"
0804        //       << associated_ratio
0805        //       << endl;
0806        
0807        this->FillClusterInfo( 0 );
0808        this->FillClusterInfo( 1 );
0809 
0810      }
0811 
0812    for( auto& hist : hists_ )
0813      {
0814        hist->SetAdc7Correction( does_adc7_correction_ );
0815        hist->SetAdc14Correction( does_adc14_correction_ );
0816        hist->ModifyAdc();
0817        //hist->Print();
0818      }
0819 
0820 }
0821 
0822 void Analysis::Draw()
0823 {
0824 
0825   SetsPhenixStyle();
0826   this->InitParameter();
0827   
0828   string output = "results/mip_" + to_string( run_ );
0829   if( is_preliminary_ == true )
0830     output += "_preliminary.pdf";
0831   else
0832     output += "_internal.pdf";
0833 
0834   TFile* tf_output = new TFile( (output.substr(0, output.size()-4) + ".root").c_str(), "RECREATE" );
0835   
0836   c_ = new TCanvas( output.c_str(), "title", 800, 800 );
0837   c_->Print( ((string)c_->GetName() + "[").c_str() );
0838 
0839   // drawing a single hist
0840   for( auto& mip_hist : hists_ )
0841     this->DrawSingle( mip_hist );
0842 
0843   //  if( false )
0844   { // same as Misaki
0845     vector < MipHist* > mip_hists;
0846     mip_hists.push_back( hist_ang85_ );
0847     mip_hists.push_back( hist_ang45_ );
0848     mip_hists.push_back( hist_ang35_ );
0849     this->DrawMultiple( mip_hists, "85-90_40-45_30-35" );
0850   }
0851 
0852   //  if( false )
0853   { // show cuts effect
0854     vector < MipHist* > mip_hists;
0855     mip_hists.push_back( hist_all_ );
0856     //mip_hists.push_back( hist_aso_ );
0857     mip_hists.push_back( hist_ang85_ );
0858     this->DrawMultiple( mip_hists, "raw_only_top_angle" );
0859   }
0860 
0861   //  if( false )
0862   { // tracklet association or not
0863     vector < MipHist* > mip_hists;
0864     mip_hists.push_back( hist_aso_ );
0865     mip_hists.push_back( hist_no_aso_ );
0866     this->DrawMultiple( mip_hists, "associated_not_associated" );
0867   }
0868 
0869   gStyle -> SetPalette( kLightTemperature );
0870   //  if( false )
0871   { // all my config
0872     vector < MipHist* > mip_hists;
0873     mip_hists.push_back( hist_ang80_ );
0874     mip_hists.push_back( hist_ang70_ );
0875     mip_hists.push_back( hist_ang60_ );
0876     mip_hists.push_back( hist_ang50_ );
0877     mip_hists.push_back( hist_ang40_ );
0878     mip_hists.push_back( hist_ang30_ );
0879     mip_hists.push_back( hist_ang20_ );
0880     mip_hists.push_back( hist_ang10_ );
0881     mip_hists.push_back( hist_ang0_ );
0882 
0883     this->DrawMultiple( mip_hists, "all", true );
0884   }
0885 
0886   { // all my config
0887     vector < MipHist* > mip_hists;
0888     mip_hists.push_back( hist_ang80_ );
0889     mip_hists.push_back( hist_ang50_ );
0890     mip_hists.push_back( hist_ang20_ );
0891     mip_hists.push_back( hist_ang0_ );
0892     this->DrawMultiple( mip_hists, "80-90_50-60_20-30_0-10", true );
0893   }
0894 
0895   { // all my config
0896     vector < MipHist* > mip_hists;
0897     mip_hists.push_back( hist_ang10_11_ );
0898     mip_hists.push_back( hist_ang20_21_ );
0899     mip_hists.push_back( hist_ang30_31_ );
0900     this->DrawMultiple( mip_hists, "10-11_20-21_30-31", true );
0901   }
0902 
0903   { // all my config
0904     vector < MipHist* > mip_hists;
0905     mip_hists.push_back( hist_ang80_ );
0906     mip_hists.push_back( hist_ang70_ );
0907     mip_hists.push_back( hist_ang60_ );
0908     mip_hists.push_back( hist_ang50_ );
0909     mip_hists.push_back( hist_ang40_ );
0910     mip_hists.push_back( hist_ang30_ );
0911     mip_hists.push_back( hist_ang20_ );
0912     this->DrawMultiPanel( mip_hists, "multi_panel" );
0913   }
0914 
0915   //////////////////////////////////////////////////////////////////////
0916   // #cluster correlation                                             //
0917   //////////////////////////////////////////////////////////////////////
0918   gStyle -> SetPalette( kBird );
0919 
0920   top_margin_ = right_margin_ = 0.13;
0921   c_->SetTopMargin( top_margin_ );
0922   c_->SetBottomMargin( 0.13 );
0923   c_->SetRightMargin( right_margin_ );
0924   c_->SetLeftMargin( 0.13 );
0925 
0926   hist_correlation_->GetYaxis()->SetTitleOffset( 1.2 );
0927   hist_correlation_->Draw( "colz" );
0928   gPad->SetLogz( true );
0929   //  DrawStats( hist_correlation_, 0.7, 0.7, 0.9, 0.9 );
0930   this->DrawWords();
0931   c_->Print( c_->GetName() );
0932   
0933   string output_correlation = output.substr( 0, output.find_last_of("_") ) + "_cluster_correlation";
0934   if( is_preliminary_ == true )
0935     output_correlation += "_preliminary.pdf";
0936   else
0937     output_correlation += "_internal.pdf";
0938     
0939   c_->Print( output_correlation.c_str() );
0940   
0941   c_->Print( ((string)c_->GetName() + "]").c_str() );
0942 
0943   for( auto& mip_hist : hists_ )
0944     tf_output->WriteTObject( mip_hist->GetHist(), mip_hist->GetHist()->GetName() );
0945 
0946   tf_output->WriteTObject( hist_association_, hist_association_->GetName() );
0947   tf_output->WriteTObject( hist_association_in_out_, hist_association_in_out_->GetName() );
0948   tf_output->WriteTObject( hist_association_in_, hist_association_in_->GetName() );
0949   tf_output->WriteTObject( hist_association_out_, hist_association_out_->GetName() );
0950   tf_output->WriteTObject( hist_correlation_, hist_correlation_->GetName() );
0951   tf_output->Close();
0952 }
0953 
0954 void Analysis::InitParameter()
0955 {
0956 
0957   adc7_ = 210;
0958   counter_too_large_cluster_ = 0;
0959   cluster_size_max_ = 5;
0960   hit_num_inner_ = 0;
0961   hit_num_outer_ = 0;
0962   
0963   top_margin_ = 0.05; // 0.1;
0964   right_margin_ = 0.05; // 0.15;
0965   adc7_modification_factor_ = 0.35;
0966 }
0967 
0968 void Analysis::Show(Long64_t entry)
0969 {
0970 // Print contents of entry.
0971 // If entry is not specified, print current entry
0972    if (!fChain) return;
0973    fChain->Show(entry);
0974 }
0975 
0976 void Analysis::Print()
0977 {
0978 
0979   for( auto& hist : hists_ )
0980     hist->Print();
0981 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!