Latest changes

diff --git a/HLT/exa/trigger_pp.C b/HLT/exa/trigger_pp.C
index 9d87415a799b43767ad52637c93c3137a2bdd88e..ddec0eea9fb11d676ab23532a3bfadc749a3c19a 100644 (file)
--- a/HLT/exa/trigger_pp.C
+++ b/HLT/exa/trigger_pp.C
@@ -1,63 +1,216 @@
-void trigger()
+void trigger_pp(char *outfile="results.root")
 {
   
-  //Get the tracks:
-  AliL3TrackArray *tracks = new AliL3TrackArray();
-  AliL3FileHandler *file = new AliL3FileHandler();
-  file->SetBinaryInput("tracks.raw");
-  file->Binary2TrackArray(tracks);
-  file->CloseBinaryInput();
-  delete file;
-  
-  Int_t ntracks=0;
-  Double_t xc,yc,zc;
-  Double_t impact;
-  AliL3Vertex vertex;
-  AliL3Fitter *fitter = new AliL3Fitter(&vertex);
-  fitter->LoadClusters("./");
-  //fitter->NoVertex();
-  for(int i=0; i<tracks->GetNTracks(); i++)
+  TNtuple *ntuppel = new TNtuple("ntuppel","","pt:eta:xvert:yvert:zvert:nhits:px:py:pz:event");
+  Float_t meas[10];  
+  
+  for(int event=0; event<1; event++)
     {
-      track = (AliL3Track*)tracks->GetCheckedTrack(i);
-      if(!track) continue;
-      if(fabs(track->GetPseudoRapidity())>0.9) continue;
-      if(track->GetNHits() < 50) continue;
-      if(track->GetPt()<0.1) continue;
-      track->CalculateHelix();
-      //cout<<"Pt "<<track->GetPt()<<" eta "<<track->GetPseudoRapidity()<<" Nhits "<<track->GetNHits()<<endl;
-      //cout<<"Before second fit; pt "<<track->GetPt()<<" firstpoint "<<track->GetFirstPointX()<<" "<<track->GetFirstPointY()<<" "<<track->GetFirstPointZ()<<endl;
-      fitter->FitHelix(track);
-      track->CalculateHelix();
+      char fname[256];
+      sprintf(fname,"/prog/alice/data/Rawdata/1_patch/pp/pileups/recon_%d/tracks.raw",event);
+      //sprintf(fname,"/prog/alice/data/Rawdata/1_patch/pp/recon_%d/tracks.raw",event);
+
+      //Get the tracks:
+      AliL3TrackArray *tracks = new AliL3TrackArray();
+      AliL3FileHandler *file = new AliL3FileHandler();
+      file->SetBinaryInput(fname);
+      file->Binary2TrackArray(tracks);
+      file->CloseBinaryInput();
+      delete file;
+      
+      sprintf(fname,"/prog/alice/data/Rawdata/1_patch/pp/pileups/recon_%d/",event);
+      //sprintf(fname,"/prog/alice/data/Rawdata/1_patch/pp/recon_%d/",event);
+
+      Int_t ntracks=0;
+      Double_t xc,yc,zc;
+      Double_t impact;
+      AliL3Vertex vertex;
+      AliL3Fitter *fitter = new AliL3Fitter(&vertex);
+      fitter->LoadClusters(fname);
+      //fitter->NoVertex();
       
-      track->GetClosestPoint(&vertex,xc,yc,zc);
-      impact = sqrt(xc*xc+yc*yc+zc*zc);
-      if(impact>3) continue;
-      cout<<"Number of hits "<<track->GetNHits()<<endl;
-      cout<<"Transversal impact "<<sqrt(xc*xc+yc*yc)<<endl;
-      cout<<"Longitudinal impact "<<zc<<endl;
-      cout<<"Total impact "<<sqrt(xc*xc+yc*yc+zc*zc)<<endl;
-      cout<<"xc "<<xc<<" yc "<<yc<<" zc "<<zc<<endl;
+      AliL3TrackArray *ftracks = new AliL3TrackArray();
       
-      //cout<<"After second fit; pt "<<track->GetPt()<<" firstpoint "<<track->GetFirstPointX()<<" "<<track->GetFirstPointY()<<" "<<track->GetFirstPointZ()<<endl;
-      cout<<"pt "<<track->GetPt()<<" eta "<<track->GetPseudoRapidity()<<endl;
+      for(int i=0; i<tracks->GetNTracks(); i++)
+       {
+         track = (AliL3Track*)tracks->GetCheckedTrack(i);
+         if(!track) continue;
+         track->CalculateHelix();
+         //cout<<"Pt "<<track->GetPt()<<" eta "<<track->GetPseudoRapidity()<<" Nhits "<<track->GetNHits()<<endl;
+         //cout<<"Before second fit; pt "<<track->GetPt()<<" firstpoint "<<track->GetFirstPointX()<<" "<<track->GetFirstPointY()<<" "<<track->GetFirstPointZ()<<endl;
+         fitter->FitHelix(track);//refit the tracks
+         track->CalculateHelix();
+         track->GetClosestPoint(&vertex,xc,yc,zc);
+         meas[0]=track->GetPt();
+         meas[1]=track->GetPseudoRapidity();
+         meas[2]=xc;
+         meas[3]=yc;
+         meas[4]=zc;
+         meas[5]=track->GetNHits();
+         meas[6]=track->GetPx();
+         meas[7]=track->GetPy();
+         meas[8]=track->GetPz();
+         meas[9]=event;
+         ntuppel->Fill(meas);
+         //continue;
+         if(fabs(track->GetPseudoRapidity())>0.9) continue;
+         if(track->GetNHits() < 100) continue;
+         if(track->GetPt()<0.2) continue;
+         impact = sqrt(xc*xc+yc*yc+zc*zc);
+         if(fabs(zc)>3) continue;
+         ftracks->AddLast(track);
+         cout<<"-------------------------------------"<<endl;
+         cout<<"Number of hits "<<track->GetNHits()<<endl;
+         cout<<"Transversal impact "<<sqrt(xc*xc+yc*yc)<<endl;
+         cout<<"Longitudinal impact "<<zc<<endl;
+         cout<<"Total impact "<<sqrt(xc*xc+yc*yc+zc*zc)<<endl;
+         cout<<"xc "<<xc<<" yc "<<yc<<" zc "<<zc<<endl;
+         
+         //cout<<"After second fit; pt "<<track->GetPt()<<" firstpoint "<<track->GetFirstPointX()<<" "<<track->GetFirstPointY()<<" "<<track->GetFirstPointZ()<<endl;
+         cout<<"pt "<<track->GetPt()<<" eta "<<track->GetPseudoRapidity()<<endl;
+         
+         ntracks++;
+       }
+      cout<<endl<<"There was "<<ntracks<<" accepted tracks, out of total "<<tracks->GetNTracks()<<" found tracks"<<endl;
       
-      cout<<"-------------------------------------"<<endl;
-      ntracks++;
+      display(ftracks,fname);
+      delete tracks;
+      delete fitter;
+      delete ftracks;
+    }
+  
+  TFile *rfile = TFile::Open(outfile,"RECREATE");
+  ntuppel->Write();
+  rfile->Close();
+  delete ntuppel;
+
+}
+
+void display(AliL3TrackArray *tracks,char *path)
+{
+  int slice[2]={0,35};
+  d = new AliL3Display(slice);
+  d->Setup("tracks.raw",path);
+  d->SetTracks(tracks);
+  //d->DisplayClusters();
+    d->DisplayAll();
+    //d->DisplayTracks();
+  
+}
+
+void ploteff()
+{
+  gROOT->LoadMacro("XFunct.C");
+  //double x[6]={1,2,4,6,8,10};
+  //double y[6]={0.873684,1.01379,1.17751,1.28614,1.31638,1.32022};
+  
+  double x[6]={0.8,1.6,3.2,4.8,6.4,8};
+  double y[6]={0.497006,0.88024,1.19162,1.23952,1.39521,1.40719};//->sigmas
+
+  
+  TGraph *gr = new TGraph(6,x,y);
+  c1 = new TCanvas("c1","",2);
+  SetCanvasOptions(c1);
+  gr->SetMarkerStyle(20);
+  gr->SetTitle("");
+  gr->Draw("APL");
+  gr->GetHistogram()->SetXTitle("Z_{CUT} [x #sigma_{Z}]");
+  gr->GetHistogram()->SetYTitle("Efficiency");
+  SetTGraphOptions(gr);
+}
+
+double geteff(char *infile,char *singlefile,double cut)
+{
+  gROOT->LoadMacro("XFunct.C");
+  gStyle->SetStatColor(10);
+  gStyle->SetOptFit(1);
+  Int_t pileups[25],single[25];
+  file = TFile::Open(infile,"READ");
+  
+  char name[256];
+  for(int i=0; i<25; i++)
+    {
+      sprintf(name,"zvert < %f && zvert > %f && pt>0.2 && nhits>100 && eta>-0.9 && eta<0.9 && event==%d",cut,-1.*cut,i);
+      ntuppel->Draw(">>eventlist",name);
+      int entries = eventlist->GetN();
+      pileups[i]=entries;
+    }
+  //eventlist->Print("all");
+  file->Close();
+  
+  file = TFile::Open(singlefile,"read");
+  for(int i=0; i<25; i++)
+    {
+      sprintf(name,"zvert < %f && zvert > %f && pt>0.2 && nhits>100 && eta>-0.9 && eta<0.9 && event==%d",3,-3,i);
+      ntuppel->Draw(">>eventlist",name);
+      int entries = eventlist->GetN();
+      single[i]=entries;
     }
-  cout<<endl<<"There was "<<ntracks<<" found tracks"<<endl;
-  delete tracks;
-  delete fitter;
+  file->Close();
+  double totsingle=0,totpileup=0;
+  for(int i=0; i<25; i++)
+    {
+      totsingle += single[i];
+      totpileup += pileups[i];
+    }
+  double toteff = totpileup/totsingle;
+  cout<<"Total eff "<<toteff<<endl;
+  
+  return toteff;
+}
+
+
+void plot(char *infile)
+{
+  gROOT->LoadMacro("XFunct.C");
+  gStyle->SetStatColor(10);
+  gStyle->SetOptFit(1);
+  file = TFile::Open(infile,"READ");
+  
+  histo = new TH1F("histo","histo",20,-3,3);
+  
+  vhist = new TH1F("vhist","",20,-3,3);
+  SetTH1Options(vhist);
+  
+  char fname[256];
+  char fname2[256];
+  for(int ev=0; ev<25; ev++)
+    {
+      sprintf(fname,"pt>0.2 && nhits>100 && eta>-0.9 && eta<0.9 && event==%d",ev);
+      ntuppel->Draw("zvert>>histo",fname,"goff");
+      float mean = histo->GetMean();
+      vhist->Fill(mean);
+      continue;
+      sprintf(fname2,"(zvert-(%f))>>+vhist",mean);
+      cout<<fname2<<endl;
+      ntuppel->Draw(fname2,fname,"goff");
+    }
+  
+  c1 = new TCanvas("c1","",2);
+  SetCanvasOptions(c1);
+  vhist->SetXTitle("Z* [cm]");
+  vhist->Draw();
+  return;
+  //ntuppel->Draw("zvert>>histo","pt>0.2");
+  TF1 *f1 = new TF1("f1","gaus",-3,3);
+  histo->Fit("f1","R");
+
+  //histo->Draw();
+  //file->Close();
 }
 
 enum tagprimary {kPrimaryCharged=0x4000};
 void LoadEvent(Int_t event=0)
 {
   //Load the generated particles
-
+  
   gROOT->LoadMacro("$(ALICE_ROOT)/macros/loadlibs.C");
   loadlibs();
   TFile *rootfile = TFile::Open("/prog/alice/data/pro/25event_pp.root","READ");
   gAlice = (AliRun*)rootfile->Get("gAlice");
+  
+  //  TNtuple *ntup = new TNtuple(
+  
   Int_t nparticles = gAlice->GetEvent(event);
   Int_t nprim = FindPrimaries(nparticles,0.,0.,0.);
   cout<<"Number of primaries "<<nprim<<endl;
@@ -85,21 +238,21 @@ Int_t FindPrimaries(Int_t nparticles, Double_t xori, Double_t yori, Double_t zor
   Double_t vertcut = 0.001;  // cut on the vertex position
   Double_t decacut = 3.;     // cut if the part. decays close to the vert.
   Double_t timecut = 0.;
-
+  
   TList *listprim = new TList();
   listprim->SetOwner(kFALSE);
-
+  
   Int_t nprch1=0;
   Int_t nprch2=0;
   for(Int_t iprim = 0; iprim<nparticles; iprim++){   //loop on  tracks
-
+    
     TParticle * part = gAlice->Particle(iprim);
     char *xxx=strstr(part->GetName(),"XXX");
     if(xxx)continue;
-
+    
     TParticlePDG *ppdg = part->GetPDG();
     if(TMath::Abs(ppdg->Charge())!=3)continue;  // only charged (no quarks)
-
+    
     Double_t dist=TMath::Sqrt((part->Vx()-xori)*(part->Vx()-xori)+(part->Vy()-yori)*(part->Vy()-yori)+(part->Vz()-zori)*(part->Vz()-zori));
     if(dist>vertcut)continue;  // cut on the vertex
 
@@ -120,7 +273,9 @@ Int_t FindPrimaries(Int_t nparticles, Double_t xori, Double_t yori, Double_t zor
       for(Int_t j=fidau;j<=lasdau;j++){
         TParticle *dau=gAlice->Particle(j);
         Double_t distd=TMath::Sqrt((dau->Vx()-xori)*(dau->Vx()-xori)+(dau->Vy()-yori)*(dau->Vy()-yori)+(dau->Vz()-zori)*(dau->Vz()-zori));
-        if(distd<decacut)prmch=kFALSE;  // eliminate if the decay is near the vertex
+        Double_t rad = TMath::Sqrt((dau->Vx()-xori)*(dau->Vx()-xori)+(dau->Vy()-yori)*(dau->Vy()-yori));
+       if(distd<decacut)prmch=kFALSE;  // eliminate if the decay is near the vertex
+       //if(rad < 20)prmch=kFALSE;
       }
     }