1 void trigger_pp(char *outfile="results.root")
4 TNtuple *ntuppel = new TNtuple("ntuppel","","pt:eta:xvert:yvert:zvert:nhits:px:py:pz:event");
7 for(int event=0; event<1; event++)
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);
14 AliL3TrackArray *tracks = new AliL3TrackArray();
15 AliL3FileHandler *file = new AliL3FileHandler();
16 file->SetBinaryInput(fname);
17 file->Binary2TrackArray(tracks);
18 file->CloseBinaryInput();
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);
28 AliL3Fitter *fitter = new AliL3Fitter(&vertex);
29 fitter->LoadClusters(fname);
32 AliL3TrackArray *ftracks = new AliL3TrackArray();
34 for(int i=0; i<tracks->GetNTracks(); i++)
36 track = (AliL3Track*)tracks->GetCheckedTrack(i);
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();
49 meas[5]=track->GetNHits();
50 meas[6]=track->GetPx();
51 meas[7]=track->GetPy();
52 meas[8]=track->GetPz();
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;
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;
74 cout<<endl<<"There was "<<ntracks<<" accepted tracks, out of total "<<tracks->GetNTracks()<<" found tracks"<<endl;
76 display(ftracks,fname);
82 TFile *rfile = TFile::Open(outfile,"RECREATE");
89 void display(AliL3TrackArray *tracks,char *path)
92 d = new AliL3Display(slice);
93 d->Setup("tracks.raw",path);
95 //d->DisplayClusters();
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};
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
111 TGraph *gr = new TGraph(6,x,y);
112 c1 = new TCanvas("c1","",2);
113 SetCanvasOptions(c1);
114 gr->SetMarkerStyle(20);
117 gr->GetHistogram()->SetXTitle("Z_{CUT} [x #sigma_{Z}]");
118 gr->GetHistogram()->SetYTitle("Efficiency");
119 SetTGraphOptions(gr);
122 double geteff(char *infile,char *singlefile,double cut)
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");
131 for(int i=0; i<25; i++)
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();
138 //eventlist->Print("all");
141 file = TFile::Open(singlefile,"read");
142 for(int i=0; i<25; i++)
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();
150 double totsingle=0,totpileup=0;
151 for(int i=0; i<25; i++)
153 totsingle += single[i];
154 totpileup += pileups[i];
156 double toteff = totpileup/totsingle;
157 cout<<"Total eff "<<toteff<<endl;
163 void plot(char *infile)
165 gROOT->LoadMacro("XFunct.C");
166 gStyle->SetStatColor(10);
167 gStyle->SetOptFit(1);
168 file = TFile::Open(infile,"READ");
170 histo = new TH1F("histo","histo",20,-3,3);
172 vhist = new TH1F("vhist","",20,-3,3);
173 SetTH1Options(vhist);
177 for(int ev=0; ev<25; ev++)
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();
184 sprintf(fname2,"(zvert-(%f))>>+vhist",mean);
186 ntuppel->Draw(fname2,fname,"goff");
189 c1 = new TCanvas("c1","",2);
190 SetCanvasOptions(c1);
191 vhist->SetXTitle("Z* [cm]");
194 //ntuppel->Draw("zvert>>histo","pt>0.2");
195 TF1 *f1 = new TF1("f1","gaus",-3,3);
196 histo->Fit("f1","R");
202 enum tagprimary {kPrimaryCharged=0x4000};
203 void LoadEvent(Int_t event=0)
205 //Load the generated particles
207 gROOT->LoadMacro("$(ALICE_ROOT)/macros/loadlibs.C");
209 TFile *rootfile = TFile::Open("/prog/alice/data/pro/25event_pp.root","READ");
210 gAlice = (AliRun*)rootfile->Get("gAlice");
212 // TNtuple *ntup = new TNtuple(
214 Int_t nparticles = gAlice->GetEvent(event);
215 Int_t nprim = FindPrimaries(nparticles,0.,0.,0.);
216 cout<<"Number of primaries "<<nprim<<endl;
218 for(Int_t i=0; i<nparticles; i++)
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;
226 cout<<endl<<"Number of primary tracks in the detector: "<<co<<endl;
232 Int_t FindPrimaries(Int_t nparticles, Double_t xori, Double_t yori, Double_t zori)
234 //Define primary particles in a pp-event. Code taken from offline.
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.;
242 TList *listprim = new TList();
243 listprim->SetOwner(kFALSE);
247 for(Int_t iprim = 0; iprim<nparticles; iprim++){ //loop on tracks
249 TParticle * part = gAlice->Particle(iprim);
250 char *xxx=strstr(part->GetName(),"XXX");
253 TParticlePDG *ppdg = part->GetPDG();
254 if(TMath::Abs(ppdg->Charge())!=3)continue; // only charged (no quarks)
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
259 if(part->T()>timecut)continue;
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
264 Bool_t prmch = kTRUE; // candidate primary track
265 Int_t fidau=part->GetFirstDaughter(); // cut on daughters
269 lasdau=part->GetLastDaughter();
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;
284 part->SetBit(kPrimaryCharged);
285 listprim->Add(part); // list of primary particles (before cleanup)
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);
299 listprim->Remove(moth);
300 moth->ResetBit(kPrimaryCharged);
306 listprim->Clear("nodelete");
308 return nprch1-nprch2;