Renamed houghtest.C to runhough.C
[u/mrichter/AliRoot.git] / HLT / exa / trigger_pp.C
1 void trigger_pp(char *outfile="results.root")
2 {
3   
4   TNtuple *ntuppel = new TNtuple("ntuppel","","pt:eta:xvert:yvert:zvert:nhits:px:py:pz:event");
5   Float_t meas[10];  
6   
7   for(int event=0; event<1; event++)
8     {
9       char fname[256];
10       sprintf(fname,"/prog/alice/data/Rawdata/1_patch/pp/pileups/recon_%d/tracks.raw",event);
11       //sprintf(fname,"/prog/alice/data/Rawdata/1_patch/pp/recon_%d/tracks.raw",event);
12
13       //Get the tracks:
14       AliL3TrackArray *tracks = new AliL3TrackArray();
15       AliL3FileHandler *file = new AliL3FileHandler();
16       file->SetBinaryInput(fname);
17       file->Binary2TrackArray(tracks);
18       file->CloseBinaryInput();
19       delete file;
20       
21       sprintf(fname,"/prog/alice/data/Rawdata/1_patch/pp/pileups/recon_%d/",event);
22       //sprintf(fname,"/prog/alice/data/Rawdata/1_patch/pp/recon_%d/",event);
23
24       Int_t ntracks=0;
25       Double_t xc,yc,zc;
26       Double_t impact;
27       AliL3Vertex vertex;
28       AliL3Fitter *fitter = new AliL3Fitter(&vertex);
29       fitter->LoadClusters(fname);
30       //fitter->NoVertex();
31       
32       AliL3TrackArray *ftracks = new AliL3TrackArray();
33       
34       for(int i=0; i<tracks->GetNTracks(); i++)
35         {
36           track = (AliL3Track*)tracks->GetCheckedTrack(i);
37           if(!track) continue;
38           track->CalculateHelix();
39           //cout<<"Pt "<<track->GetPt()<<" eta "<<track->GetPseudoRapidity()<<" Nhits "<<track->GetNHits()<<endl;
40           //cout<<"Before second fit; pt "<<track->GetPt()<<" firstpoint "<<track->GetFirstPointX()<<" "<<track->GetFirstPointY()<<" "<<track->GetFirstPointZ()<<endl;
41           fitter->FitHelix(track);//refit the tracks
42           track->CalculateHelix();
43           track->GetClosestPoint(&vertex,xc,yc,zc);
44           meas[0]=track->GetPt();
45           meas[1]=track->GetPseudoRapidity();
46           meas[2]=xc;
47           meas[3]=yc;
48           meas[4]=zc;
49           meas[5]=track->GetNHits();
50           meas[6]=track->GetPx();
51           meas[7]=track->GetPy();
52           meas[8]=track->GetPz();
53           meas[9]=event;
54           ntuppel->Fill(meas);
55           //continue;
56           if(fabs(track->GetPseudoRapidity())>0.9) continue;
57           if(track->GetNHits() < 100) continue;
58           if(track->GetPt()<0.2) continue;
59           impact = sqrt(xc*xc+yc*yc+zc*zc);
60           if(fabs(zc)>3) continue;
61           ftracks->AddLast(track);
62           cout<<"-------------------------------------"<<endl;
63           cout<<"Number of hits "<<track->GetNHits()<<endl;
64           cout<<"Transversal impact "<<sqrt(xc*xc+yc*yc)<<endl;
65           cout<<"Longitudinal impact "<<zc<<endl;
66           cout<<"Total impact "<<sqrt(xc*xc+yc*yc+zc*zc)<<endl;
67           cout<<"xc "<<xc<<" yc "<<yc<<" zc "<<zc<<endl;
68           
69           //cout<<"After second fit; pt "<<track->GetPt()<<" firstpoint "<<track->GetFirstPointX()<<" "<<track->GetFirstPointY()<<" "<<track->GetFirstPointZ()<<endl;
70           cout<<"pt "<<track->GetPt()<<" eta "<<track->GetPseudoRapidity()<<endl;
71           
72           ntracks++;
73         }
74       cout<<endl<<"There was "<<ntracks<<" accepted tracks, out of total "<<tracks->GetNTracks()<<" found tracks"<<endl;
75       
76       display(ftracks,fname);
77       delete tracks;
78       delete fitter;
79       delete ftracks;
80     }
81   
82   TFile *rfile = TFile::Open(outfile,"RECREATE");
83   ntuppel->Write();
84   rfile->Close();
85   delete ntuppel;
86
87 }
88
89 void display(AliL3TrackArray *tracks,char *path)
90 {
91   int slice[2]={0,35};
92   d = new AliL3Display(slice);
93   d->Setup("tracks.raw",path);
94   d->SetTracks(tracks);
95   //d->DisplayClusters();
96     d->DisplayAll();
97     //d->DisplayTracks();
98   
99 }
100
101 void ploteff()
102 {
103   gROOT->LoadMacro("XFunct.C");
104   //double x[6]={1,2,4,6,8,10};
105   //double y[6]={0.873684,1.01379,1.17751,1.28614,1.31638,1.32022};
106   
107   double x[6]={0.8,1.6,3.2,4.8,6.4,8};
108   double y[6]={0.497006,0.88024,1.19162,1.23952,1.39521,1.40719};//->sigmas
109
110   
111   TGraph *gr = new TGraph(6,x,y);
112   c1 = new TCanvas("c1","",2);
113   SetCanvasOptions(c1);
114   gr->SetMarkerStyle(20);
115   gr->SetTitle("");
116   gr->Draw("APL");
117   gr->GetHistogram()->SetXTitle("Z_{CUT} [x #sigma_{Z}]");
118   gr->GetHistogram()->SetYTitle("Efficiency");
119   SetTGraphOptions(gr);
120 }
121
122 double geteff(char *infile,char *singlefile,double cut)
123 {
124   gROOT->LoadMacro("XFunct.C");
125   gStyle->SetStatColor(10);
126   gStyle->SetOptFit(1);
127   Int_t pileups[25],single[25];
128   file = TFile::Open(infile,"READ");
129   
130   char name[256];
131   for(int i=0; i<25; i++)
132     {
133       sprintf(name,"zvert < %f && zvert > %f && pt>0.2 && nhits>100 && eta>-0.9 && eta<0.9 && event==%d",cut,-1.*cut,i);
134       ntuppel->Draw(">>eventlist",name);
135       int entries = eventlist->GetN();
136       pileups[i]=entries;
137     }
138   //eventlist->Print("all");
139   file->Close();
140   
141   file = TFile::Open(singlefile,"read");
142   for(int i=0; i<25; i++)
143     {
144       sprintf(name,"zvert < %f && zvert > %f && pt>0.2 && nhits>100 && eta>-0.9 && eta<0.9 && event==%d",3,-3,i);
145       ntuppel->Draw(">>eventlist",name);
146       int entries = eventlist->GetN();
147       single[i]=entries;
148     }
149   file->Close();
150   double totsingle=0,totpileup=0;
151   for(int i=0; i<25; i++)
152     {
153       totsingle += single[i];
154       totpileup += pileups[i];
155     }
156   double toteff = totpileup/totsingle;
157   cout<<"Total eff "<<toteff<<endl;
158   
159   return toteff;
160 }
161
162
163 void plot(char *infile)
164 {
165   gROOT->LoadMacro("XFunct.C");
166   gStyle->SetStatColor(10);
167   gStyle->SetOptFit(1);
168   file = TFile::Open(infile,"READ");
169   
170   histo = new TH1F("histo","histo",20,-3,3);
171   
172   vhist = new TH1F("vhist","",20,-3,3);
173   SetTH1Options(vhist);
174   
175   char fname[256];
176   char fname2[256];
177   for(int ev=0; ev<25; ev++)
178     {
179       sprintf(fname,"pt>0.2 && nhits>100 && eta>-0.9 && eta<0.9 && event==%d",ev);
180       ntuppel->Draw("zvert>>histo",fname,"goff");
181       float mean = histo->GetMean();
182       vhist->Fill(mean);
183       continue;
184       sprintf(fname2,"(zvert-(%f))>>+vhist",mean);
185       cout<<fname2<<endl;
186       ntuppel->Draw(fname2,fname,"goff");
187     }
188   
189   c1 = new TCanvas("c1","",2);
190   SetCanvasOptions(c1);
191   vhist->SetXTitle("Z* [cm]");
192   vhist->Draw();
193   return;
194   //ntuppel->Draw("zvert>>histo","pt>0.2");
195   TF1 *f1 = new TF1("f1","gaus",-3,3);
196   histo->Fit("f1","R");
197
198   //histo->Draw();
199   //file->Close();
200 }
201
202 enum tagprimary {kPrimaryCharged=0x4000};
203 void LoadEvent(Int_t event=0)
204 {
205   //Load the generated particles
206   
207   gROOT->LoadMacro("$(ALICE_ROOT)/macros/loadlibs.C");
208   loadlibs();
209   TFile *rootfile = TFile::Open("/prog/alice/data/pro/25event_pp.root","READ");
210   gAlice = (AliRun*)rootfile->Get("gAlice");
211   
212   //  TNtuple *ntup = new TNtuple(
213   
214   Int_t nparticles = gAlice->GetEvent(event);
215   Int_t nprim = FindPrimaries(nparticles,0.,0.,0.);
216   cout<<"Number of primaries "<<nprim<<endl;
217   int co=0;
218   for(Int_t i=0; i<nparticles; i++)
219     {
220       TParticle *part = gAlice->Particle(i);
221       if(!part->TestBit(kPrimaryCharged)) continue;
222       if(fabs(part->Eta())>0.9) continue;
223       cout<<part->GetName()<<" pt "<<part->Pt()<<" eta "<<part->Eta()<<" xvert "<<part->Vx()<<" yvert "<<part->Vy()<<" zvert "<<part->Vz()<<endl;
224       co++;
225     }
226   cout<<endl<<"Number of primary tracks in the detector: "<<co<<endl;
227   gAlice=0;
228   rootfile->Close();
229   
230 }
231
232 Int_t FindPrimaries(Int_t nparticles, Double_t xori, Double_t yori, Double_t zori)
233 {
234   //Define primary particles in a pp-event. Code taken from offline.
235
236
237   // cuts:
238   Double_t vertcut = 0.001;  // cut on the vertex position
239   Double_t decacut = 3.;     // cut if the part. decays close to the vert.
240   Double_t timecut = 0.;
241   
242   TList *listprim = new TList();
243   listprim->SetOwner(kFALSE);
244   
245   Int_t nprch1=0;
246   Int_t nprch2=0;
247   for(Int_t iprim = 0; iprim<nparticles; iprim++){   //loop on  tracks
248     
249     TParticle * part = gAlice->Particle(iprim);
250     char *xxx=strstr(part->GetName(),"XXX");
251     if(xxx)continue;
252     
253     TParticlePDG *ppdg = part->GetPDG();
254     if(TMath::Abs(ppdg->Charge())!=3)continue;  // only charged (no quarks)
255     
256     Double_t dist=TMath::Sqrt((part->Vx()-xori)*(part->Vx()-xori)+(part->Vy()-yori)*(part->Vy()-yori)+(part->Vz()-zori)*(part->Vz()-zori));
257     if(dist>vertcut)continue;  // cut on the vertex
258
259     if(part->T()>timecut)continue;
260
261     Double_t ptot=TMath::Sqrt(part->Px()*part->Px()+part->Py()*part->Py()+part->Pz()*part->Pz());
262     if(ptot==(TMath::Abs(part->Pz())))continue; // no beam particles
263
264     Bool_t prmch = kTRUE;   // candidate primary track
265     Int_t fidau=part->GetFirstDaughter();  // cut on daughters
266     Int_t lasdau=0;
267     Int_t ndau=0;
268     if(fidau>=0){
269       lasdau=part->GetLastDaughter();
270       ndau=lasdau-fidau+1;
271     }
272     if(ndau>0){
273       for(Int_t j=fidau;j<=lasdau;j++){
274         TParticle *dau=gAlice->Particle(j);
275         Double_t distd=TMath::Sqrt((dau->Vx()-xori)*(dau->Vx()-xori)+(dau->Vy()-yori)*(dau->Vy()-yori)+(dau->Vz()-zori)*(dau->Vz()-zori));
276         Double_t rad = TMath::Sqrt((dau->Vx()-xori)*(dau->Vx()-xori)+(dau->Vy()-yori)*(dau->Vy()-yori));
277         if(distd<decacut)prmch=kFALSE;  // eliminate if the decay is near the vertex
278         //if(rad < 20)prmch=kFALSE;
279       }
280     }
281
282     if(prmch){
283       nprch1++;
284       part->SetBit(kPrimaryCharged);
285       listprim->Add(part);    // list of primary particles (before cleanup)
286     }
287   }
288
289
290   nprch2=0;
291   for(Int_t iprim = 0; iprim<nparticles; iprim++){ // cleanup loop
292     TParticle * part = gAlice->Particle(iprim);
293     if(part->TestBit(kPrimaryCharged)){
294       Int_t mothind=part->GetFirstMother();
295       if(mothind<0)continue;
296       TParticle *moth=gAlice->Particle(mothind);
297       TParticle *mothb=(TParticle*)listprim->FindObject(moth);
298       if(mothb){
299         listprim->Remove(moth);
300         moth->ResetBit(kPrimaryCharged);
301         nprch2++;
302       }
303     }
304   }
305
306   listprim->Clear("nodelete");
307   delete listprim;
308   return nprch1-nprch2;
309   
310 }