Updated version of the Bari code to work with the HEAD. A new test macros has also...
[u/mrichter/AliRoot.git] / ITS / ITSsddanalysis.C
CommitLineData
13787886 1void 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