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