source/LQGENEP/Distill_pi.C

0001 int Distill_pi(string seed)
0002 {
0003
0004     std::string inputFile = "TestOut_hadronic_"+seed+"seed.txt";
0005     std::ifstream fin(inputFile.c_str());
0006     std::string str;
0007
0008     std::string foutname = "LQ_3pion_"+seed+"seed.txt";
0009     std::ofstream fout;
0010     fout.open(foutname.c_str());
0011
0012     int I, KS, KF, ORIG;
0013     const int tau_KF = 15;
0014     const int pi_KF = 211;
0015     int I_prev = 0;
0016     int event_count = 0;
0017     int pi_num = 0;
0018
0019     vector<int> is_three_pi;
0020     vector<int> tau_orig;
0021
0022     fin.clear();
0023         fin.seekg(0, ios::beg);
0024
0025
0026     while (std::getline(fin, str))
0027       {
0028         if(str.length() == 127)
0029           {
0030         //                    cout << str << endl;
0031         std::istringstream iss(str);
0032
0033         //Read in variables in event
0034         if(iss >> I >> KS >> KF >> ORIG)
0035           {
0036             //Reset tau daughters if new event
0037             if(I < I_prev) tau_orig.clear();
0038             I_prev = I;
0039
0040             if(TMath::Abs(KF) == tau_KF){
0041             //Put tau particle in array
0042             tau_orig.push_back(I);
0043             }
0044             else
0045               {
0046             for(int i = 0; i < tau_orig.size(); i++)
0047               {
0048                 //If daughter from tau
0049                 if(ORIG == tau_orig[i])
0050                   {
0051
0052                 //add up number of pions
0053                 if(TMath::Abs(KF) == pi_KF){
0054                   pi_num++;
0055                 }
0056                 //Add particle to daughters from tau
0057                 tau_orig.push_back(I);
0058                   }
0059               }
0060               }
0061           }
0062           }
0063         else
0064           {
0065         //End of the event
0066         if(str.find("Event finished") != std::string::npos) {
0067           //Store if this event was a three pion event
0068
0069           if(pi_num == 3) is_three_pi.push_back(1);
0070           else {
0071            is_three_pi.push_back(0);
0072           }
0073           event_count++;
0074           if(event_count % 100==0) cout<<event_count<<endl;
0075
0076           //Reset pion counter
0077           pi_num = 0;
0078         }
0079           }
0080       }
0081
0082
0083     int events_left = 0;
0084     int line_count = 0;
0085     fin.clear();
0086     fin.seekg(0, ios::beg);
0087
0088     event_count = 0;
0089     while (std::getline(fin, str))
0090           {
0091         line_count++;
0092         if(line_count<7 && is_three_pi[event_count]==0) fout << str << endl;
0093
0094             if(str.length() == 127)
0095               {
0096         if(is_three_pi[event_count] == 1) fout << str << endl;
0097         /*
0098         std::istringstream iss(str);
0099
0100                 if(iss >> I >> KS >> KF >> ORIG)
0101                   {
0102                     if(I < I_prev) tau_orig.clear();
0103                     I_prev = I;
0104
0105                     if(I < 3)
0106                       {
0107             if(is_three_pi[event_count] == 1) fout << str << endl;
0108                       }
0109                     else if(KS == 21)
0110                       { if(is_three_pi[event_count] == 1) fout << str << endl; }
0111                     else if(TMath::Abs(KF) == tau_KF)
0112                       {
0113             if(is_three_pi[event_count] == 1) fout << str << endl;
0114             tau_orig.push_back(I);
0115                       }
0116                     else
0117                       {
0118                         for(int i = 0; i < tau_orig.size(); i++)
0119                           {
0120                             if(ORIG == tau_orig[i])
0121                               {
0122                 if(is_three_pi[event_count] == 1)   fout << str << endl;
0123                 tau_orig.push_back(I);
0124                               }
0125                           }
0126                       }
0127
0128                   }
0129         */
0130           }
0131             else
0132               {
0133
0134         if(is_three_pi[event_count]) fout<<str<<endl;
0135
0136         if(str.find("Event finished") != std::string::npos) {
0137           if(is_three_pi[event_count] == 1) {
0138             events_left++;
0139           }
0140           event_count++;
0141           if(event_count % 100==0) cout<<event_count<<endl;
0142
0143           if(events_left == 1000) break;
0144         }
0145               }
0146       }
0147
0148
0149     gSystem->Load("libeicsmear");
0150         BuildTree(foutname.c_str(),"./",events_left);
0151
0152     return 0;
0153 }