]>
Commit | Line | Data |
---|---|---|
13787886 | 1 | void ITSsddanalysis (Int_t evNumber1=0,Int_t evNumber2=0) |
2 | { | |
3 | ///////////////////////////////////////////////////////////////////////// | |
4 | // This macro is a small example of a ROOT macro | |
5 | // illustrating how to read the output of GALICE | |
6 | // and fill some histograms. | |
7 | // | |
8 | // Root > .L anal.C //this loads the macro in memory | |
9 | // Root > anal(); //by default process first event | |
10 | // Root > anal(2); //process third event | |
11 | //Begin_Html | |
12 | /* | |
13 | <img src="gif/anal.gif"> | |
14 | */ | |
15 | //End_Html | |
16 | ///////////////////////////////////////////////////////////////////////// | |
17 | ||
18 | // Dynamically link some shared libs | |
19 | ||
20 | if (gClassTable->GetID("AliRun") < 0) { | |
21 | gROOT->LoadMacro("loadlibs.C"); | |
22 | loadlibs(); | |
1f8e4927 | 23 | } else { |
24 | delete gAlice; | |
25 | gAlice=0; | |
13787886 | 26 | } |
27 | ||
28 | // Connect the Root Galice file containing Geometry, Kine and Hits | |
29 | TString *str = new TString("galice.root"); | |
30 | TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(str->Data()); | |
31 | if (!file) file = new TFile(str->Data(),"UPDATE"); | |
32 | ||
33 | // Get AliRun object from file or create it if not on file | |
34 | // if (!gAlice) { | |
35 | gAlice = (AliRun*)file->Get("gAlice"); | |
36 | if (gAlice) printf("AliRun object found on file\n"); | |
37 | if (!gAlice) gAlice = new AliRun("gAlice","Alice test program"); | |
38 | //} | |
39 | ||
40 | TH2F *local = new TH2F("local","time vs anode difference",101,-1.01,1.01,210,-210.,210.); | |
41 | TH2F *local1 = new TH2F("local1","time vs anode difference",101,-5.01,5.01,210,-2100.,2100.); | |
42 | TH2F *dtah = new TH2F("dtah","anode difference vs drift time (hits)",210,-21000.,21000.,101,-10000.01,10000.01); | |
43 | TH2F *dtap = new TH2F("dtap","anode difference vs drift time (points)",21000,-21000.,210.,101,-10000.01,10000.01); | |
44 | TH1F *dclutim = new TH1F("dclutim","cluster time difference (dan < 1)",200,-2000.,2000.); | |
45 | TH1F *dfkctim = new TH1F("dfkctim","cluster time difference (dan > 5)",200,-2000.,2000.); | |
46 | TH2F *xy3 = new TH2F("xy3","y vs x",100,-200.02,200.02,100,-200.02,200.02); | |
47 | TH2F *xz3 = new TH2F("xz3","x vs z",100,-200.02,200.02,100,-200.02,200.02); | |
48 | TH2F *yz3 = new TH2F("yz3","y vs z",100,-200.02,200.02,100,-200.02,200.02); | |
49 | TH2F *xy4 = new TH2F("xy4","y vs x",100,-200.02,200.02,100,-200.02,200.02); | |
50 | TH2F *xz4 = new TH2F("xz4","x vs z",100,-200.02,200.02,100,-200.02,200.02); | |
51 | TH2F *yz4 = new TH2F("yz4","y vs z",100,-200.02,200.02,100,-200.02,200.02); | |
52 | TH1F *chd = new TH1F("chargediff","Charge difference (Gen-Rec)",100,-1000.,1000.); | |
53 | TH1F *chr = new TH1F("chargeratio","Charge ratio (Gen/Rec)",100,0.,0.1); | |
54 | TH2F *chp = new TH2F("chp","Point Charge vs time",28,0.,7000.,300,0.,3000.); | |
55 | TH2F *chh = new TH2F("chh","Hit Charge vs time",28,0.,7000.,80,0.,40.); | |
56 | TH2F *dadth = new TH2F("dadth","danode vs dtime (hits)",560,-14000.,14000.,601,-300.5,300.5); | |
57 | TH2F *dadt = new TH2F("dadt","danode vs dtime (points)",280,-7000.,7000.,601,-300.5,300.5); | |
58 | TH2F *aa = new TH2F("aa","anode hit vs point",204,-0.5,50.5.,204,-.5,50.5); | |
59 | TH2F *at = new TH2F("at","anode difference vs drift path (mm)",18,0.,36.,100,-200.,200.); | |
60 | TH2F *tt = new TH2F("tt","time coord. difference (um) vs drift path (mm)",18,0.,36.,100,-200.,200.); | |
61 | TH2F *at1 = new TH2F("at1","(1a) anode difference vs drift path (mm)",18,0.,36.,100,-200.,200.); | |
62 | TH2F *tt1 = new TH2F("tt1","(1a) time coord. difference (um) vs drift path (mm)",18,0.,36.,100,-200.,200.); | |
63 | TH2F *at2 = new TH2F("at2","(2a) anode difference vs drift path (mm)",18,0.,36.,100,-200.,200.); | |
64 | TH2F *tt2 = new TH2F("tt2","(2a) time coord. difference (um) vs drift path (mm)",18,0.,36.,100,-200.,200.); | |
65 | TH2F *asigma = new TH2F("asigma","anode sigma vs drift path (mm)",18,0.,36.,300,0.,300.); | |
66 | TH2F *tsigma = new TH2F("tsigma","tau. vs drift path (mm)",18,0.,36.,200,0.,100.); | |
67 | TH2F *asigma2 = new TH2F("asigma2","2 anode sigma vs drift path (mm)",18,0.,36.,150,0.,300.); | |
68 | ||
69 | TH1F *dtrp = new TH1F("dtrp","Double track separation (microns)",40,0.,2000.); | |
70 | TH1F *dtrpall = new TH1F("dtrpall","Double track separation (mm)",100,0.,100.); | |
71 | TH1F *dtrh = new TH1F("dtrh","Double track separation (microns)",40,0.,2000.); | |
72 | TH1F *dtrhall = new TH1F("dtrhall","Double track separation (mm)",100,0.,100.); | |
73 | ||
74 | TH1F *p = new TH1F("p","Momentum ",100,0.,20.); | |
75 | ||
76 | TH1F *effh = new TH1F("effh","Hit Multiplicity vs drift path (mm)",18,0.,36.); | |
77 | TH1F *effp = new TH1F("effp","Point Multiplicity vs drift path (mm)",18,0.,36.); | |
78 | ||
79 | TH2F *anodes = new TH2F("nanodes","Anode Multiplicity vs drift time",28,0.,7000.,5,0.5,5.5); | |
80 | TH2F *andtsm = new TH2F("nand_ntsm","Anode Mult vs Time Mult",15,0.5,15.5,5,0.5,5.5); | |
81 | TH2F *tsampl = new TH2F("nsampls","Sample Multiplicity vs drift time",28,0.,7000.,15,0.5,15.5); | |
82 | TH2F *ntotal = new TH2F("ntotal","Cluster Multiplicity vs drift time",28,0.,7000.,60,0.5,60.5); | |
83 | TH2F *clmult = new TH2F("clmult","Anode Multiplicity vs Total Multiplicity",50,0.5,50.5,5,0.5,5.5); | |
84 | TH2F *amplit1 = new TH2F("amplit1","Point Amplitude vs drift path (mm)",28,0.,7000.,64,0.5,1024.5); | |
85 | TH2F *amplit = new TH2F("amplit","Point Amplitude vs drift path (mm)",28,0.,7000.,60,0.5,600.5); | |
86 | TH1F *hitpnt = new TH1F("hitpnt","Hit-Point Multiplicity",21,-10.5,10.5); | |
87 | ||
88 | TH1F *nmatch = new TH1F("nmatch","Number of matched points",5,-0.5.,4.5); | |
89 | TH1F *rec_vs_time = new TH1F("rec_vs_time","Point Rec. vs drift path",36,0.,36.); | |
90 | TH1F *hit_vs_time = new TH1F("hit_vs_time","Hit vs drift path",36,0.,36.); | |
91 | TH1F *rec_vs_time3 = new TH1F("rec_vs_time3","Point Rec. vs drift path",36,0.,36.); | |
92 | TH1F *hit_vs_time3 = new TH1F("hit_vs_time3","Hit vs drift path",36,0.,36.); | |
93 | TH1F *rec_vs_time4 = new TH1F("rec_vs_time4","Point Rec. vs drift path",36,0.,36.); | |
94 | TH1F *hit_vs_time4 = new TH1F("hit_vs_time4","Hit vs drift path",36,0.,36.); | |
95 | TH1F *rec_vs_time1 = new TH1F("rec_vs_time1","Point Rec. vs drift path",36,0.,36.); | |
96 | TH1F *hit_vs_time1 = new TH1F("hit_vs_time1","Hit vs drift path",36,0.,36.); | |
97 | TH1F *fake_vs_time = new TH1F("fake_vs_time","fake points vs drift path",36,0.,36.); | |
98 | ||
99 | TH1F *noihist = new TH1F("noisehist","noise",80,10.,30.); | |
100 | ||
101 | TH1F *occupancy3 = new TH1F("occupancy3","Occupancy vs Detector Number, Layer 3",20,0.5,20.5); | |
102 | TH1F *occupancy4 = new TH1F("occupancy4","Occupancy vs Detector Number, Layer 4",20,0.5,20.5); | |
103 | ||
104 | TH2F *pntmap3 = new TH2F("pntmap3","Point map Layer 3",20,0.5,20.5,10,0.5,10.5); | |
105 | TH2F *hitmap3 = new TH2F("hitmap3","Hit map Layer 3",20,0.5,20.5,10,0.5,10.5); | |
106 | TH2F *map3 = new TH2F("map3","Hit/Point map Layer 3",20,0.5,20.5,10,0.5,10.5); | |
107 | TH2F *pntmap4 = new TH2F("pntmap4","Point map Layer 4",30,0.5,30.5,10,0.5,10.5); | |
108 | TH2F *hitmap4 = new TH2F("hitmap4","Hit map Layer 4",30,0.5,30.5,10,0.5,10.5); | |
109 | TH2F *map4 = new TH2F("map4","Hit/Point map Layer 4",30,0.5,30.5,10,0.5,10.5); | |
110 | TH2F *xz = new TH2F("xz","X vs Z",50,-5,5.,50,-5.,5.); | |
111 | TH2F *and_tim = new TH2F("and_tim","Tim vs Anode",30,-100,356.,30,-8000.,8000.); | |
112 | TH2F *pand_ptim = new TH2F("pand_ptim","Tim vs Anode",30,-100,356.,30,-8000.,8000.); | |
113 | ||
114 | //Int_t nanodes = 256; | |
115 | //TH2F *mappa3hit[14][6][2]; | |
116 | //TH2F *mappa4hit[22][8][2]; | |
117 | //TH2F *mappa3pnt[14][6][2]; | |
118 | //TH2F *mappa4pnt[22][8][2]; | |
119 | /* | |
120 | for(Int_t i=0;i<22;i++) { | |
121 | for(Int_t j=0;j<8;j++) { | |
122 | for(Int_t k=0;k<2;k++) { | |
123 | TString *hname = new TString("hitmap_"); | |
124 | TString *cname = new TString("pntmap_"); | |
125 | Char_t lad[2]; | |
126 | sprintf(lad,"%d",i+1); | |
127 | hname->Append(lad); | |
128 | hname->Append("_"); | |
129 | cname->Append(lad); | |
130 | cname->Append("_"); | |
131 | Char_t det[2]; | |
132 | sprintf(det,"%d",j+1); | |
133 | hname->Append(det); | |
134 | hname->Append("_"); | |
135 | cname->Append(det); | |
136 | cname->Append("_"); | |
137 | Char_t wng[2]; | |
138 | sprintf(wng,"%d",k+1); | |
139 | hname->Append(wng); | |
140 | cname->Append(wng); | |
141 | //mappa4hit[i][j][k] = new TH2F(hname->Data(),hname->Data(),nanodes,0.5,nanodes+0.5,256,0.5,256.5); | |
142 | //mappa4pnt[i][j][k] = new TH2F(cname->Data(),cname->Data(),nanodes,0.5,nanodes+0.5,256,0.5,256.5); | |
143 | if(i<14 && j<6) { | |
144 | mappa3hit[i][j][k] = new TH2F(hname->Data(),hname->Data(),nanodes,0.5,nanodes+0.5,256,0.5,256.5); | |
145 | mappa3pnt[i][j][k] = new TH2F(cname->Data(),cname->Data(),nanodes,0.5,nanodes+0.5,256,0.5,256.5); | |
146 | } | |
147 | } | |
148 | } | |
149 | } | |
150 | */ | |
151 | ||
152 | AliITS *ITS = (AliITS*) gAlice->GetModule("ITS"); | |
153 | if (!ITS) { cout << "no ITS" << endl; return; } | |
154 | ||
155 | Int_t nparticles = gAlice->GetEvent(0); | |
156 | ||
157 | Int_t cp[8]={0,0,0,0,0,0,0,0}; | |
158 | ||
159 | AliITSDetType *iDetType=ITS->DetType(1); | |
160 | ||
161 | AliITSresponseSDD *res1 = (AliITSresponseSDD*)iDetType->GetResponseModel(); | |
162 | if (!res1) { | |
163 | res1=new AliITSresponseSDD(); | |
164 | ITS->SetResponseModel(1,res1); | |
165 | } | |
166 | res1->SetZeroSupp("2D"); // 1D | |
167 | res1->SetParamOptions("same","same"); | |
168 | //res1->SetFilenames(" ","$(ALICE_ROOT)/ITS/base.dat","$(ALICE_ROOT)/ITS/2D.dat "); | |
169 | res1->SetCompressParam(cp); | |
170 | res1->SetDriftSpeed(7.3); | |
171 | Float_t vdrift = res1->DriftSpeed(); | |
172 | ||
173 | AliITSsegmentationSDD *seg1=(AliITSsegmentationSDD*)iDetType->GetSegmentationModel(); | |
174 | AliITSgeom *aliitsgeo = ITS->GetITSgeom(); | |
175 | ||
176 | Int_t cp[8]={0,0,0,0,0,0,0,0}; | |
177 | ||
178 | Int_t dum = 0; | |
179 | Float_t apitch = seg1->Dpz(dum); | |
180 | Float_t tstep = seg1->Dpx(dum); | |
181 | Float_t maxand = seg1->Npz()/2.; | |
182 | //cout << "anodes: " << maxand << ", tstep: " << tstep << endl; | |
183 | ||
184 | Float_t n,b; | |
185 | res1->GetNoiseParam(n,b); | |
186 | printf("SDD: noise baseline %f %f zs option %s data type %s\n",n,b,res1->ZeroSuppOption(),res1->DataType()); | |
187 | printf("SDD: DriftSpeed %f TopValue %f\n",res1->DriftSpeed(),res1->DynamicRange()); | |
188 | Float_t dif0,dif1; | |
189 | res1->DiffCoeff(dif0,dif1); | |
190 | printf("SDD: dif0 %f dif1 %f\n",dif0,dif1); | |
191 | ||
192 | AliITSsimulationSDD *sim1=new AliITSsimulationSDD(seg1,res1); | |
193 | ITS->SetSimulationModel(1,sim1); | |
194 | ||
195 | // | |
196 | // Loop over events | |
197 | // | |
198 | ||
199 | Int_t Nh=0; | |
200 | Int_t Nh1=0; | |
201 | for (int nev=0; nev<= evNumber2; nev++) { | |
202 | if(nev>0) { | |
203 | nparticles = gAlice->GetEvent(nev); | |
204 | gAlice->SetEvent(nev); | |
205 | } | |
206 | cout << "nparticles " <<nparticles<<endl; | |
207 | if (nev < evNumber1) continue; | |
208 | if (nparticles <= 0) return; | |
209 | ||
210 | // Reset Pointers | |
211 | AliITShit *itsHit; | |
212 | AliITSRecPoint *itsPnt = 0; | |
213 | AliITSRawClusterSDD *itsClu = 0; | |
214 | ||
215 | // Reset Event Counters | |
216 | Int_t nGoodTotalHits = 0; | |
217 | Int_t nGoodTotalPnts = 0; | |
218 | ||
219 | // Get Hit, Cluster & Recpoints Tree Pointers | |
220 | ||
221 | TTree *TH = gAlice->TreeH(); | |
222 | Int_t nenthit=TH->GetEntries(); | |
223 | printf("Found %d entries in the Hit tree (must be one per track per event!)\n",nenthit); | |
224 | ||
225 | ITS->GetTreeC(nev); | |
226 | TTree *TC=ITS->TreeC(); | |
227 | Int_t nentclu=TC->GetEntries(); | |
43c39b72 | 228 | printf("Found %d entries in the Cluster tree (must be one per module per event!)\n",nentclu); |
13787886 | 229 | |
230 | TTree *TR = gAlice->TreeR(); | |
231 | Int_t nentrec=TR->GetEntries(); | |
232 | printf("Found %d entries in the Rec tree (must be one per module per event!)\n",nentrec); | |
233 | ||
13787886 | 234 | // check recpoints |
235 | ||
236 | Int_t nbytes = 0; | |
237 | Int_t totpoints = 0; | |
238 | Int_t totclust = 0; | |
239 | ||
240 | // check hits | |
241 | ||
242 | Int_t nmodules=0; | |
243 | ITS->InitModules(-1,nmodules); | |
244 | ITS->FillModules(nev,0,nmodules,"",""); | |
245 | ||
246 | TObjArray *fITSmodules = ITS->GetModules(); | |
247 | ||
43c39b72 | 248 | Int_t iDet = 1; // 1 = SDD |
249 | ||
250 | Int_t first = aliitsgeo->GetStartDet(0); | |
251 | Int_t last = aliitsgeo->GetLastDet(2); | |
252 | Int_t first_ok = aliitsgeo->GetStartDet(iDet); | |
253 | Int_t last_ok = aliitsgeo->GetLastDet(iDet); | |
13787886 | 254 | printf("det type %d first, last %d %d \n",iDet,first,last); |
255 | ||
43c39b72 | 256 | TClonesArray *ITSclu = 0; |
257 | TClonesArray *ITSrec = 0; | |
258 | for (Int_t mod=first_ok; mod<=last_ok; mod++) { | |
259 | cout << "Module: " << mod << endl; | |
260 | // if(mod < first_ok) continue; | |
261 | // if(mod > last_ok) continue; | |
262 | // Get Pointers to Clusters & Recpoints TClonesArrays | |
263 | ||
264 | ITSclu = ITS->ClustersAddress(iDet); | |
265 | ITSrec = ITS->RecPoints(); | |
266 | ||
13787886 | 267 | TTree *TR = gAlice->TreeR(); |
268 | Int_t nentrec=TR->GetEntries(); | |
269 | TClonesArray *ITSrec = ITS->RecPoints(); | |
270 | ||
271 | ITS->ResetClusters(); | |
272 | TC->GetEvent(mod); | |
273 | ITS->ResetRecPoints(); | |
274 | nbytes += TR->GetEvent(mod); | |
275 | ||
276 | Int_t nrecp = ITSrec->GetEntries(); | |
277 | totpoints += nrecp; | |
43c39b72 | 278 | if (nrecp) printf("Found %d rec points for module %d\n",nrecp,mod); |
13787886 | 279 | if (!nrecp) continue; |
280 | ||
281 | Int_t nrecc = ITSclu->GetEntries(); | |
282 | totclust += nrecc; | |
43c39b72 | 283 | if (nrecc) printf("Found %d clusters for module %d\n",nrecc,mod); |
13787886 | 284 | |
285 | Int_t nrecp = ITSrec->GetEntries(); | |
286 | Int_t startSDD = aliitsgeo->GetStartSDD(); | |
287 | Int_t *flagP = new Int_t [nrecp]; | |
288 | memset( flagP, 0, sizeof(Int_t)*nrecp ); | |
289 | ||
290 | //printf("point loop\n"); | |
291 | ||
292 | Int_t nGoodPnts = 0; | |
293 | for (Int_t pnt=0;pnt<nrecp;pnt++) { | |
294 | itsPnt = (AliITSRecPoint*)ITSrec->At(pnt); | |
295 | if(!itsPnt) continue; | |
296 | itsClu = (AliITSRawClusterSDD*)ITSclu->At(pnt); | |
297 | if(!itsClu) continue; | |
298 | //itsClu->PrintInfo(); | |
299 | nGoodPnts++; | |
300 | nGoodTotalPnts++; | |
301 | ||
302 | Int_t pntlayer; | |
303 | Int_t pntladder; | |
304 | Int_t pntdetector; | |
305 | aliitsgeo->GetModuleId(mod+first,pntlayer,pntladder,pntdetector); | |
306 | Int_t pntmult = itsClu->Samples(); | |
307 | Int_t pntands = itsClu->Anodes(); | |
308 | Float_t pnttime = itsClu->T(); | |
309 | Float_t pntanod = itsClu->A(); | |
310 | Float_t pntchrg = itsClu->Q(); | |
311 | Float_t pntampl = itsClu->PeakAmpl(); | |
312 | Float_t pntpath = pnttime*vdrift/1000.; | |
313 | ||
314 | Float_t wy = 0.; | |
315 | if(itsClu->Anodes() != 0.) { | |
316 | wy = pntmult/((Float_t) pntands); | |
317 | } | |
318 | clmult->Fill((Float_t)pntmult,(Float_t) pntands); | |
319 | ntotal->Fill(pnttime,(Float_t)pntmult); | |
320 | tsampl->Fill(pnttime,wy); | |
321 | amplit->Fill(pnttime,pntampl); | |
322 | amplit1->Fill(pnttime,pntampl); | |
323 | ||
324 | // Detector occupancy | |
325 | ||
326 | if(pntlayer == 3) { | |
327 | occupancy3->Fill((Float_t) pntdetector,(Float_t) pntmult); | |
328 | //mappa3pnt[pntladder-1][pntdetector-1][0]->Fill(pntanod,pnttime); | |
329 | } | |
330 | if(pntlayer == 4) { | |
331 | occupancy4->Fill((Float_t) pntdetector,(Float_t) pntmult); | |
332 | //mappa4pnt[pntladder-1][pntdetector-1][0]->Fill(pntanod,pnttime); | |
333 | } | |
334 | ||
335 | // Point Efficiency vs time. | |
336 | ||
337 | effp->Fill(pntpath); | |
338 | anodes->Fill(pnttime,pntands); | |
339 | andtsm->Fill(wy,pntands); | |
340 | chp->Fill(pnttime,pntchrg); | |
341 | } | |
342 | ||
343 | //printf("hit loop\n"); | |
344 | ||
345 | Float_t sddLength = seg1->Dx(); | |
346 | Float_t sddWidth = seg1->Dz(); | |
347 | Float_t driftSpeed=res1->DriftSpeed(); | |
348 | ||
349 | Int_t nGoodHits = 0; | |
350 | ||
351 | AliITSmodule *Mod = (AliITSmodule *)fITSmodules->At(mod+first); | |
352 | Int_t nhits = Mod->GetNhits(); | |
353 | for (Int_t hit=0;hit<nhits;hit++) { | |
354 | itsHit = (AliITShit*)Mod->GetHit(hit); | |
355 | ||
356 | Float_t avx = 0.; | |
357 | Float_t avy = 0.; | |
358 | Float_t avz = 0.; | |
359 | Int_t ifl = 0; | |
360 | Float_t DepEnergy = 100000.*itsHit->GetIonization(); | |
361 | AliITShit *itsHit1 = 0; | |
362 | if(DepEnergy == 0.) { | |
363 | hit++; | |
364 | if(hit == nhits) break; | |
365 | itsHit1 = (AliITShit*) Mod->GetHit(hit); | |
366 | avx = itsHit1->GetXG(); | |
367 | avy = itsHit1->GetYG(); | |
368 | avz = itsHit1->GetZG(); | |
369 | ifl = 1; | |
370 | } | |
371 | avx += itsHit->GetXG(); | |
372 | avy += itsHit->GetYG(); | |
373 | avz += itsHit->GetZG(); | |
374 | if(DepEnergy == 0.) { | |
375 | avx /= 2.; | |
376 | avy /= 2.; | |
377 | avz /= 2.; | |
378 | } | |
379 | if(ifl == 0) continue; | |
380 | ||
381 | Float_t px; Float_t py; Float_t pz; | |
382 | itsHit->GetMomentumG(px,py,pz); | |
383 | Float_t ptot = TMath::Sqrt(px*px+py*py+pz*pz); | |
384 | p->Fill(ptot*100); | |
385 | if(ptot < 0.05) continue; | |
386 | ||
387 | Int_t Layer = itsHit->GetLayer(); | |
388 | Int_t Ladder = itsHit->GetLadder(); | |
389 | Int_t Det = itsHit->GetDetector(); | |
390 | ||
391 | Float_t And; | |
392 | Float_t Tim; | |
393 | Float_t x = itsHit->GetXL(); | |
394 | Float_t z = itsHit->GetZL(); | |
395 | xz->Fill(z,x); | |
396 | seg1->GetPadTxz(x,z); | |
397 | And = z; | |
398 | Tim = x*tstep; | |
399 | and_tim->Fill(And,Tim); | |
400 | ||
401 | Float_t And1; | |
402 | Float_t Tim1; | |
403 | Float_t x1; | |
404 | Float_t z1; | |
405 | if(itsHit1) { | |
406 | x1 = itsHit1->GetXL(); | |
407 | z1 = itsHit1->GetZL(); | |
408 | xz->Fill(z1,x1); | |
409 | seg1->GetPadTxz(x1,z1); | |
410 | And1 = z1; | |
411 | Tim1 = x1*tstep; | |
412 | and_tim->Fill(And1,Tim1); | |
413 | } | |
414 | Float_t DepEnergy = 100000.*itsHit->GetIonization(); | |
415 | if(DepEnergy == 0.) DepEnergy = 100000.*itsHit1->GetIonization(); | |
416 | if(DepEnergy < 5.) continue; | |
417 | ||
418 | if(itsHit1) { | |
419 | Tim += Tim1; | |
420 | Tim /= 2.; | |
421 | And += And1; | |
422 | And /= 2.; | |
423 | } | |
424 | if(And < 0. || And > maxand) { cout << "And: " << And << endl; continue; } | |
425 | Float_t path = TMath::Abs(Tim)*vdrift/1000.; | |
426 | hit_vs_time->Fill(path); | |
427 | if(Layer==3) hit_vs_time3->Fill(path); | |
428 | if(Layer==4) hit_vs_time4->Fill(path); | |
429 | ||
430 | nGoodHits++; | |
431 | nGoodTotalHits++; | |
432 | ||
433 | //if(Layer == 3) mappa3hit[Ladder-1][Det-1][0]->Fill(And,Tim); | |
434 | //if(Layer == 4) mappa4hit[Ladder-1][Det-1][0]->Fill(And,Tim); | |
435 | ||
436 | effh->Fill(path); | |
437 | Float_t ww = DepEnergy; | |
438 | chh->Fill(TMath::Abs(Tim),ww); | |
439 | Int_t inmatches = 0; | |
440 | ||
441 | Float_t diffmin = 100000.; | |
442 | Int_t pntmin = -1; | |
443 | ||
444 | //printf("point loop\n"); | |
13787886 | 445 | for (Int_t pnt=0;pnt<nrecp;pnt++) { |
446 | itsPnt = (AliITSRecPoint*)ITSrec->At(pnt); | |
447 | if(!itsPnt) continue; | |
448 | itsClu = (AliITSRawClusterSDD*)ITSclu->At(pnt); | |
449 | if(!itsClu) continue; | |
450 | ||
451 | Int_t LayerP; | |
452 | Int_t LadderP; | |
453 | Int_t DetP; | |
454 | aliitsgeo->GetModuleId(mod+first,LayerP,LadderP,DetP); | |
455 | Int_t LayerH = itsHit->GetLayer(); | |
456 | Int_t LadderH = itsHit->GetLadder(); | |
457 | Int_t DetH = itsHit->GetDetector(); | |
458 | if(LayerH != LayerP) continue; | |
459 | if(LadderH != LadderP) continue; | |
460 | if(DetH != DetP) continue; | |
461 | ||
462 | Float_t Pand = (Float_t) itsClu->A(); | |
463 | if(Pand < 0 || Pand > maxand) { cout << "Pand: " << Pand << endl; continue; } | |
464 | Float_t Ptim = (Float_t) itsClu->T(); | |
465 | Float_t Pwng = (Float_t) itsClu->W(); | |
466 | if(Pwng == 1) Ptim *= -1.; | |
467 | pand_ptim->Fill(Pand,Ptim); | |
468 | Float_t adiff = And-Pand; | |
469 | Float_t tdiff = Tim-Ptim; | |
470 | if(And < 0) { | |
471 | printf("tim %f\n",Tim); | |
472 | printf("and %f\n",And); | |
473 | } | |
474 | if(Pwng == 1) tdiff *=-1.; | |
475 | local1->Fill(adiff,tdiff); | |
476 | ||
477 | if(TMath::Abs(adiff) >= 1) continue; | |
478 | if(TMath::Abs(tdiff) >= 100) continue; | |
479 | ||
480 | Float_t apdiff = adiff*apitch; | |
481 | Float_t tpdiff = tdiff*vdrift; | |
482 | ||
483 | Float_t diff = TMath::Sqrt( apdiff*apdiff+tpdiff*tpdiff ); | |
484 | if( diff < diffmin ){ | |
485 | diffmin = diff; | |
486 | pntmin = pnt; | |
487 | } | |
488 | ||
489 | if(TMath::Abs(adiff) < 1. && TMath::Abs(tdiff) < 100.) { | |
490 | inmatches++; | |
491 | } | |
cc736781 | 492 | } // loop (points) |
493 | ||
494 | if( pntmin > -1 ) { | |
495 | if( flagP[pntmin] == 1) continue; | |
496 | flagP[pntmin] = 1; | |
497 | itsClu = (AliITSRawClusterSDD*)ITSclu->At( pntmin ); | |
498 | Float_t Pand = (Float_t) itsClu->A(); | |
499 | Float_t Ptim = (Float_t) itsClu->T(); | |
500 | Float_t Pwng = (Float_t) itsClu->W(); | |
501 | Float_t sigma = itsClu->Asigma(); | |
502 | Float_t tau = itsClu->Tsigma(); | |
503 | Int_t pntands = itsClu->Anodes(); | |
504 | if(Pwng == 1) Ptim *= -1.; | |
505 | Float_t adiff = And-Pand; | |
506 | Float_t tdiff = Tim-Ptim; | |
507 | if(Pwng == 1) tdiff *=-1.; | |
508 | local->Fill(adiff,tdiff); | |
509 | ||
510 | Float_t dpath = Ptim*vdrift/1000.; | |
511 | Float_t dpathh = Tim*vdrift/1000.; | |
512 | Float_t adpath = TMath::Abs(dpath); | |
513 | Float_t adpathh = TMath::Abs(dpathh); | |
514 | Float_t apdiff = adiff*apitch; | |
515 | Float_t tpdiff = tdiff*vdrift; | |
516 | aa->Fill(Pand,And); | |
517 | ||
518 | Int_t pntands = itsClu->Anodes(); | |
519 | if(pntands == 1) { | |
520 | at1->Fill(adpath,apdiff); | |
521 | tt1->Fill(adpath,tpdiff); | |
522 | } | |
523 | if(pntands == 2) { | |
524 | at2->Fill(adpath,apdiff); | |
525 | tt2->Fill(adpath,tpdiff); | |
526 | } | |
527 | ||
528 | at->Fill(adpathh,apdiff); | |
529 | tt->Fill(adpathh,tpdiff); | |
530 | asigma->Fill(adpathh,sigma); | |
531 | tsigma->Fill(adpathh,tau); | |
532 | if(pntands == 2) asigma2->Fill(adpathh,sigma); | |
533 | ||
534 | Float_t *lP = new Float_t[3]; | |
535 | lP[0] = itsPnt->GetX(); | |
536 | lP[1] = 0.; | |
537 | lP[2] = itsPnt->GetZ(); | |
538 | Float_t *gP = new Float_t[3]; | |
539 | aliitsgeo->LtoG(LayerH,LadderH,DetH,lP,gP); | |
540 | Float_t dx = avx - gP[0]; | |
541 | Float_t dy = avy - gP[1]; | |
542 | Float_t dz = avz - gP[2]; | |
543 | delete lP; | |
544 | delete gP; | |
545 | ||
546 | Float_t pntchrg = itsClu->Q(); | |
547 | Float_t dq = DepEnergy/0.122 - pntchrg; | |
548 | Float_t rq = 0; | |
549 | if(pntchrg != 0) rq = DepEnergy/0.122/((Float_t) pntchrg); | |
550 | if(LayerH == 3) { | |
551 | xy3->Fill(dx,dy); | |
552 | xz3->Fill(dz,dx); | |
553 | yz3->Fill(dz,dy); | |
554 | } else if(LayerH == 4) { | |
555 | xy4->Fill(dx,dy); | |
556 | xz4->Fill(dz,dx); | |
557 | yz4->Fill(dz,dy); | |
13787886 | 558 | } |
cc736781 | 559 | chd->Fill(dq); |
560 | if(rq != 0.) chr->Fill(rq); | |
561 | ||
562 | rec_vs_time->Fill(adpathh); | |
563 | if(Layer==3) rec_vs_time3->Fill(adpathh); | |
564 | if(Layer==4) rec_vs_time4->Fill(adpathh); | |
13787886 | 565 | } |
cc736781 | 566 | |
13787886 | 567 | nmatch->Fill(inmatches); |
568 | } // loop hits | |
569 | ||
570 | if(nGoodHits != nGoodPnts) { | |
571 | printf("module: %d",mod+1); | |
572 | printf(", nGoodHits: %d",nGoodHits); | |
573 | printf(", nGoodPnts: %d\n",nGoodPnts); | |
574 | } | |
575 | Float_t nHP = (Float_t) nGoodHits-nGoodPnts; | |
576 | hitpnt->Fill(nHP); | |
577 | ||
578 | Int_t www = 0.; | |
579 | if(nGoodHits != 0) www = nGoodPnts/((Float_t) nGoodHits); | |
580 | if(Layer == 3) { | |
581 | pntmap3->Fill(Ladder,Det,(Float_t) nGoodPnts); | |
582 | hitmap3->Fill(Ladder,Det,(Float_t) nGoodHits); | |
583 | map3->Fill(Ladder,Det,www); | |
584 | } | |
585 | if(Layer == 4) { | |
586 | pntmap4->Fill(Ladder,Det,(Float_t) nGoodPnts); | |
587 | hitmap4->Fill(Ladder,Det,(Float_t) nGoodHits); | |
588 | map4->Fill(Ladder,Det,www); | |
589 | } | |
590 | ||
591 | //printf("double hit loop\n"); | |
592 | ||
593 | Stat_t wh = 1.; | |
594 | if(nGoodHits > 1) | |
595 | wh /= (((Float_t) nGoodHits)*((Float_t) nGoodHits)-1)/2.; | |
596 | else | |
597 | wh = 0.; | |
598 | ||
599 | Int_t *flag = new Int_t[nhits]; | |
600 | Int_t nGoodHitsOK = 0; | |
601 | for (Int_t hit=0;hit<nhits;hit++) { | |
602 | flag[hit] = 0; | |
603 | itsHit = (AliITShit*)Mod->GetHit(hit); | |
604 | Float_t avx = 0.; | |
605 | Float_t avy = 0.; | |
606 | Float_t avz = 0.; | |
607 | Int_t ifl = 0; | |
608 | Float_t DepEnergy = 100000.*itsHit->GetIonization(); | |
609 | AliITShit *itsHit1 = 0; | |
610 | if(DepEnergy == 0.) { | |
611 | hit++; | |
612 | flag[hit] = 0; | |
613 | if(hit == nhits) break; | |
614 | itsHit1 = (AliITShit*) Mod->GetHit(hit); | |
615 | avx = itsHit1->GetXG(); | |
616 | avy = itsHit1->GetYG(); | |
617 | avz = itsHit1->GetZG(); | |
618 | ifl = 1; | |
619 | } | |
620 | avx += itsHit->GetXG(); | |
621 | avy += itsHit->GetYG(); | |
622 | avz += itsHit->GetZG(); | |
623 | if(DepEnergy == 0.) { | |
624 | avx /= 2.; | |
625 | avy /= 2.; | |
626 | avz /= 2.; | |
627 | } | |
628 | if(DepEnergy < 5. && DepEnergy > 0.) continue; | |
629 | if(ifl == 0) continue; | |
630 | ||
631 | Float_t px; Float_t py; Float_t pz; | |
632 | itsHit->GetMomentumG(px,py,pz); | |
633 | Float_t ptot = TMath::Sqrt(px*px+py*py+pz*pz); | |
634 | if(ptot < 0.05) continue; | |
635 | ||
636 | for (Int_t hit1=hit+1;hit1<nhits;hit1++) { | |
637 | itsHit2 = (AliITShit*)Mod->GetHit(hit1); | |
638 | ||
639 | Float_t avx2 = 0.; | |
640 | Float_t avy2 = 0.; | |
641 | Float_t avz2 = 0.; | |
642 | Int_t ifl2 = 0; | |
643 | Float_t DepEnergy2 = 100000.*itsHit2->GetIonization(); | |
644 | AliITShit *itsHit3 = 0; | |
645 | if(DepEnergy2 == 0.) { | |
646 | hit1++; | |
647 | itsHit3 = (AliITShit*) Mod->GetHit(hit1); | |
648 | avx2 = itsHit3->GetXG(); | |
649 | avy2 = itsHit3->GetYG(); | |
650 | avz2 = itsHit3->GetZG(); | |
651 | ifl2 = 1; | |
652 | } | |
653 | avx2 += itsHit2->GetXG(); | |
654 | avy2 += itsHit2->GetYG(); | |
655 | avz2 += itsHit2->GetZG(); | |
656 | if(DepEnergy2 == 0.) { | |
657 | avx2 /= 2.; | |
658 | avy2 /= 2.; | |
659 | avz2 /= 2.; | |
660 | } | |
661 | if(DepEnergy2 < 5. && DepEnergy2 > 0.) continue; | |
662 | if(itsHit->GetLayer() != itsHit2->GetLayer()) continue; | |
663 | if(itsHit->GetLadder() != itsHit2->GetLadder()) continue; | |
664 | if(itsHit->GetDetector() != itsHit2->GetDetector()) continue; | |
665 | if(ifl2 == 0) continue; | |
666 | ||
667 | Float_t px1; Float_t py1; Float_t pz1; | |
668 | itsHit2->GetMomentumG(px1,py1,pz1); | |
669 | Float_t ptot1 = TMath::Sqrt(px1*px1+py1*py1+pz1*pz1); | |
670 | if(ptot1 < 0.05) continue; | |
671 | ||
672 | Float_t And; | |
673 | Float_t Tim; | |
674 | Float_t x = itsHit->GetXL(); | |
675 | Float_t z = itsHit->GetZL(); | |
676 | seg1->GetPadTxz(x,z); | |
677 | And = z; | |
678 | Tim = x*tstep; | |
679 | if(And < 0 || And > maxand) continue; | |
680 | Float_t And2; | |
681 | Float_t Tim2; | |
682 | Float_t x2 = itsHit2->GetXL(); | |
683 | Float_t z2 = itsHit2->GetZL(); | |
684 | seg1->GetPadTxz(x2,z2); | |
685 | And2 = z2; | |
686 | Tim2 = x2*tstep; | |
687 | if(And2 < 0 || And2 > maxand) continue; | |
688 | Float_t da = apitch*(And-And2); | |
689 | Float_t dt = vdrift*(Tim-Tim2); | |
690 | Float_t danh = And-And2; | |
691 | Float_t dtmh = Tim-Tim2; | |
692 | Float_t dist = TMath::Sqrt(da*da+dt*dt); | |
693 | if(dt < 1000.) { | |
694 | Float_t wx = dt*clock/(1000.*vdrift); | |
695 | Float_t wy = da/apitch; | |
696 | dtah->Fill(wx,wy); | |
697 | } | |
698 | if(dist<20.) { cout << "skip hit " << hit1 << endl; flag[hit] = 1; } | |
699 | ||
700 | if(ifl == 1 && ifl2 == 1) { | |
701 | if(dist>10)dtrh->Fill(dist,wh); | |
702 | dtrhall->Fill(dist/1000.,wh); | |
703 | dadth->Fill(dtmh,danh); | |
704 | } | |
705 | } // end cluster loop | |
706 | ||
707 | Float_t path = TMath::Abs(Tim)*vdrift/1000.; | |
708 | if(flag[hit] == 0) { nGoodHitsOK++; hit_vs_time1->Fill(path); } | |
709 | } // end hit loop | |
710 | delete [] flag; | |
711 | printf("nGoodHits: %d",nGoodHits); | |
712 | printf(", nGoodHitsOK: %d\n",nGoodHitsOK); | |
713 | ||
714 | //printf("cluster loop\n"); | |
715 | ||
716 | AliITSRawClusterSDD *itsCluFake = 0; | |
717 | Int_t nGoodPntsOK = 0; | |
718 | for( int ip=0; ip<nrecp; ip++) { | |
719 | itsClu = (AliITSRawClusterSDD*)ITSclu->At(ip); | |
720 | if(!itsClu) continue; | |
721 | Float_t Ptim = (Float_t) itsClu->T(); | |
722 | Float_t dpath = Ptim*vdrift/1000.; | |
723 | Float_t adpath = TMath::Abs(dpath); | |
724 | if( flagP[ip] == 1) { nGoodPntsOK++; rec_vs_time1->Fill(dpath); } | |
725 | else { | |
726 | cout << "ip: " << ip << endl; | |
727 | itsCluFake = (AliITSRawClusterSDD*)ITSclu->At( ip ); | |
728 | if(!itsCluFake) continue; | |
729 | cout << "pointer: " << itsCluFake << endl; | |
730 | itsCluFake->PrintInfo(); | |
731 | Float_t Ptim = (Float_t) itsCluFake->T(); | |
732 | Float_t dpath = Ptim*vdrift/1000.; | |
733 | fake_vs_time->Fill(dpath); | |
734 | } | |
735 | } | |
736 | ||
737 | Stat_t wp = 1.; | |
738 | if(nGoodPntsOK > 1) | |
739 | wp /= (((Float_t) nGoodPntsOK)*((Float_t) nGoodPntsOK)-1)/2.; | |
740 | else | |
741 | wp = 0.; | |
742 | ||
743 | //printf("double cluster loop\n"); | |
744 | for (Int_t pnt=0;pnt<nrecp;pnt++) { | |
745 | if( flagP[pnt] == 0) continue; | |
746 | itsPnt = (AliITSRecPoint*)ITSrec->At(pnt); | |
747 | if(!itsPnt) continue; | |
748 | itsClu = (AliITSRawClusterSDD*)ITSclu->At(pnt); | |
749 | if(!itsClu) continue; | |
750 | AliITSRecPoint *itsPnt1; | |
751 | AliITSRawClusterSDD *itsClu1; | |
752 | for (Int_t pnt1=pnt+1;pnt1<nrecp;pnt1++) { | |
753 | if( flagP[pnt1] == 0) continue; | |
754 | itsPnt1 = (AliITSRecPoint*)ITSrec->At(pnt1); | |
755 | if(!itsPnt1) continue; | |
756 | itsClu1 = (AliITSRawClusterSDD*)ITSclu->At(pnt1); | |
757 | if(!itsClu1) continue; | |
758 | ||
759 | Float_t dan = itsClu->A()-itsClu1->A(); | |
760 | Float_t Pwng = (Float_t) itsClu->W(); | |
761 | Float_t dt1 = itsClu->T(); | |
762 | if(Pwng == 1) dt1 *= -1.; | |
763 | Float_t dt2 = itsClu1->T(); | |
764 | Float_t Pwng = (Float_t) itsClu1->W(); | |
765 | if(Pwng == 1) dt2 *= -1.; | |
766 | Float_t dtm = dt1-dt2; | |
767 | dadt->Fill(dtm,dan); | |
768 | Float_t dap = apitch*(itsClu->A()-itsClu1->A()); | |
769 | Float_t dtp = vdrift*(dt1-dt2); | |
770 | Float_t distp = TMath::Sqrt(dap*dap+dtp*dtp); | |
771 | if(TMath::Abs(dan) < 1.) dclutim->Fill(dtp); | |
772 | if(TMath::Abs(dan) > 10.) dfkctim->Fill(dtp); | |
773 | if(dtp < 1000.) { | |
774 | Float_t wx = dtp*clock; | |
775 | Float_t wy = dap/apitch; | |
776 | dtap->Fill(wx,wy); | |
777 | } | |
778 | dtrp->Fill(distp,wp); | |
779 | dtrpall->Fill(distp/1000.,wp); | |
780 | } // end loop points | |
781 | } // end loop points | |
782 | ||
783 | delete [] flagP; | |
784 | } // end loop modules | |
785 | ||
786 | ITS->ClearModules(); | |
787 | gAlice->CleanDetectors(); | |
788 | ||
789 | } // end loop events | |
790 | ||
791 | cout << "open output file" << endl; | |
792 | TFile fhistos("SDD_histos_test.root","RECREATE"); | |
793 | local->Write(); | |
794 | local1->Write(); | |
795 | ||
796 | aa->Write(); | |
797 | at->Write(); | |
798 | tt->Write(); | |
799 | at1->Write(); | |
800 | tt1->Write(); | |
801 | at2->Write(); | |
802 | tt2->Write(); | |
803 | asigma->Write(); | |
804 | tsigma->Write(); | |
805 | asigma2->Write(); | |
806 | dadt->Write(); | |
807 | dadth->Write(); | |
808 | dclutim->Write(); | |
809 | dfkctim->Write(); | |
810 | xy3->Write(); | |
811 | xz3->Write(); | |
812 | yz3->Write(); | |
813 | xy4->Write(); | |
814 | xz4->Write(); | |
815 | yz4->Write(); | |
816 | chd->Write(); | |
817 | chr->Write(); | |
818 | chh->Write(); | |
819 | chp->Write(); | |
820 | dtrp->Write(); | |
821 | dtrpall->Write(); | |
822 | dtah->Write(); | |
823 | dtap->Write(); | |
824 | dtrh->Write(); | |
825 | dtrhall->Write(); | |
826 | effh->Write(); | |
827 | effp->Write(); | |
828 | rec_vs_time->Write(); | |
829 | hit_vs_time->Write(); | |
830 | hit_vs_time1->Write(); | |
831 | rec_vs_time1->Write(); | |
832 | rec_vs_time3->Write(); | |
833 | hit_vs_time3->Write(); | |
834 | rec_vs_time4->Write(); | |
835 | hit_vs_time4->Write(); | |
836 | fake_vs_time->Write(); | |
837 | p->Write(); | |
838 | nmatch->Write(); | |
839 | anodes->Write(); | |
840 | tsampl->Write(); | |
841 | ntotal->Write(); | |
842 | clmult->Write(); | |
843 | andtsm->Write(); | |
844 | amplit->Write(); | |
845 | amplit1->Write(); | |
846 | hitpnt->Write(); | |
847 | noihist->Write(); | |
848 | occupancy3->Write(); | |
849 | occupancy4->Write(); | |
850 | map3->Write(); | |
851 | hitmap3->Write(); | |
852 | pntmap3->Write(); | |
853 | map4->Write(); | |
854 | hitmap4->Write(); | |
855 | pntmap4->Write(); | |
856 | /* | |
857 | for(Int_t i=0;i<22;i++) { | |
858 | for(Int_t j=0;j<8;j++) { | |
859 | for(Int_t k=0;k<2;k++) { | |
860 | //mappa4hit[i][j][k]->Write(); | |
861 | //mappa4pnt[i][j][k]->Write(); | |
862 | if(i<14 && j<6) { | |
863 | //mappa3hit[i][j][k]->Write(); | |
864 | //mappa3pnt[i][j][k]->Write(); | |
865 | } | |
866 | } | |
867 | } | |
868 | } | |
869 | */ | |
870 | xz->Write(); | |
871 | and_tim->Write(); | |
872 | pand_ptim->Write(); | |
873 | fhistos.Close(); | |
874 | file->Close(); | |
875 | } | |
876 | ||
877 | ||
878 | ||
879 | ||
880 | ||
881 |