//macro to read mini root data
{

// Clear ROOT environment
        gROOT.Reset();
// Switch on ROOT timer
        gROOT.Time();
        gROOT->SetStyle("Plain");

        Int_t MaxEvents = 100000000;
// Create a T49Run object
        T49Run *run = new T49Run();
        run->SetVerbose(kFALSE);


        Char_t *RunType="20G+-20GeV-central-03A";   
        Int_t StartAtRun = 4883; //4881    
        Int_t EndAtRun = 4899;
        Int_t eventCount = 0;

// Define a pointer to an T49EventRoot object
        T49EventRoot    *event;
// Define a pointer to an T49ParticleRoot object
        T49ParticleRoot *track;


 T49CutEvent   *cutevent  = new T49CutEvent("eventcut","Event cuts");
 T49CutTrack   *cuttrack = new T49CutTrack("trackcut","Track cuts");
 TClonesArray  *AllTrackList;
 TObjArray     *CutTrackList;
     


// Define my own variables (number of charged part.)
	Int_t plus1;
	Int_t minus1;
        Int_t all;
        double plus_pt;
        double minus_pt;
        double plus_ptkw;
        double minus_ptkw;
        double all_pt;
        double all_ptkw; 



        ofstream *out = new ofstream;
        out->open("20_ana_pions.txt");

//        ofstream *online = new ofstream;
//        online->open("online_20_protons.text");


//Set event cuts
cutevent->SetVertexIflag(0);
cutevent->SetVertexX(-0.25, 0.2);
cutevent->SetVertexY(-0.2, 0.2);
cutevent->SetVertexZ(-581.3, -580.8);
//cutevent->SetEVeto(0,9250);
cutevent->SetVerbose(1);

//Set track cuts
cuttrack->SetBx(-2, 2);
cuttrack->SetBy(-1, 1);
cuttrack->SetNMaxPoint(3,30,1000000);
//cuttrack->SetNPointToNMaxPoint(3, 0.5);   //cos nie tak z tym cieciem
cuttrack->SetPt(0.005,1.5);
cuttrack->SetZFirst(-1000000.0, 200.0);
////SetDedxSigma(Int_t data, Int_t pid, Int_t tpc, Double_t min, Double_t max)
//(data type):0->pp,1->PbPb
//(pid):0->el,1->pi,2->ka,3->pr (tpc):2->mtpc

cuttrack->SetNDedxPoint(2,60,150);
//cuttrack->SetDedxSigma(1,3,2,0.9,1.04);  //protons - Panos cuts
cuttrack->SetDedxSigma(1,1,2,0.9,1.1);  //pions - Panos

// Loop over events in the data file
        while ((event = (T49EventRoot *) 
run->GetNextEvent(RunType,StartAtRun,EndAtRun)) && (eventCount < 
MaxEvents)) {

// find primary vertex
     T49Vertex *vertex = event->GetPrimaryVertex();
     T49Vertex  *bpd = event->GetBPDVertex();

     double ntf=vertex->GetNTrkFit();
     double nto=vertex->GetNTrkOut();
     double ratio=ntf/nto;


//Check events with given event cuts
if((cutevent->CheckEvent(event) == 1)
   &&(vertex->GetPchi2()>0)
   &&(nto!=0)
   &&(ratio>0.25))

{

//*online<<event->GetVertexZ()<<endl;


// Get a list of all primary particles in the event
           TClonesArray *AllTrackList = event->GetPrimaryParticles();
//Get accepted track list
           CutTrackList = cuttrack->GetAcceptedParticles(AllTrackList);

// start to count negative and positive particles:
	   plus1=0;
	   minus1=0;
           all=0; 
             plus_pt=0;
             minus_pt=0;
             plus_ptkw=0;
             minus_ptkw=0;
             all_pt=0;
             all_ptkw=0;
    

// Loop over all particles
           for(Int_t i=0;i<CutTrackList->GetEntries();i++)
           {
               track = (T49ParticleRoot *) CutTrackList->At(i);

 double np=track->GetNPoint();
 double nmp=track->GetNMaxPoint();

 double angle=(atan2(track->GetPy(), track->GetPx()))*(180.0/3.14);
 double rapidity=track->GetRap(0.1396);
 double ycms=rapidity-1.88;

 if(track->GetCharge()<0)
 {
   if(angle<0)
   angle=angle+360;
 angle=angle-180;
 }

double curve;

if(angle==0)
curve=1.5;
else
{
 if(ycms>1.0 && ycms<1.2)
 curve=1/((angle*angle)/600)-0.2;

 if(ycms>1.2 && ycms<1.4)
 curve=1/((angle*angle)/700)-0.2;

 if(ycms>1.4 && ycms<1.6)
 curve=1/((angle*angle)/1000)-0.2;

 if(ycms>1.6 && ycms<1.8)   
 curve=1/((angle*angle)/2600)-0.5;

 if(ycms>1.8 && ycms<2.0)
 curve=1/((angle*angle)/3000)-0.4;

 if(ycms>2.0 && ycms<2.2)
 curve=1/((angle*angle)/2500)-0.3;

 if(ycms>2.2 && ycms<2.4)
 curve=1/((angle*angle)/1800)-0.3;

 if(ycms>2.4 && ycms<2.6)
 curve=1/((angle*angle)/1000)-0.1;

}

//additional track cuts:
//not included in NA49 analysis classes

    if(((track->GetIflag()&0xFF000000) == 0)
    &&(np/nmp>0.5)
    &&(track->GetPt()<curve)
    &&(ycms>1.1)
    &&(ycms<2.6))


   {

//*online<<track->GetTmeanCharge()<<'\t'<<track->GetP()<<endl;

     all=all+1;
     all_pt=all_pt+(track->GetPt());
     all_ptkw=all_ptkw+pow((track->GetPt()),2.0);

     if(track->GetCharge()>0)
     {
     plus1=plus1+1;
     plus_pt=plus_pt+(track->GetPt());
     plus_ptkw=plus_ptkw+pow((track->GetPt()),2.0);
     }

     if(track->GetCharge()<0)
     {
     minus1=minus1+1;
     minus_pt=minus_pt+(track->GetPt());
     minus_ptkw=minus_ptkw+pow((track->GetPt()),2.0);
     }


    }
           }  
               eventCount++;
  *out << event->GetEveto() << "\t" << "\t" << all << "\t" << all_pt << 
"\t" << all_ptkw << "\t" << "\t" << minus1 << "\t" << minus_pt << "\t" << 
minus_ptkw << "\t" << "\t" << plus1 << "\t" << plus_pt << "\t" << plus_ptkw << "\t" << endl;
        } 
 
} 
        out->close();
//        online->close();
}