1 void ITSsddanalysis (Int_t evNumber1=0,Int_t evNumber2=0)
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.
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
13 <img src="gif/anal.gif">
16 /////////////////////////////////////////////////////////////////////////
18 // Dynamically link some shared libs
20 if (gClassTable->GetID("AliRun") < 0) {
21 gROOT->LoadMacro("loadlibs.C");
25 // Connect the Root Galice file containing Geometry, Kine and Hits
26 TString *str = new TString("galice.root");
27 TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(str->Data());
28 if (!file) file = new TFile(str->Data(),"UPDATE");
30 // Get AliRun object from file or create it if not on file
32 gAlice = (AliRun*)file->Get("gAlice");
33 if (gAlice) printf("AliRun object found on file\n");
34 if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
37 TH2F *local = new TH2F("local","time vs anode difference",101,-1.01,1.01,210,-210.,210.);
38 TH2F *local1 = new TH2F("local1","time vs anode difference",101,-5.01,5.01,210,-2100.,2100.);
39 TH2F *dtah = new TH2F("dtah","anode difference vs drift time (hits)",210,-21000.,21000.,101,-10000.01,10000.01);
40 TH2F *dtap = new TH2F("dtap","anode difference vs drift time (points)",21000,-21000.,210.,101,-10000.01,10000.01);
41 TH1F *dclutim = new TH1F("dclutim","cluster time difference (dan < 1)",200,-2000.,2000.);
42 TH1F *dfkctim = new TH1F("dfkctim","cluster time difference (dan > 5)",200,-2000.,2000.);
43 TH2F *xy3 = new TH2F("xy3","y vs x",100,-200.02,200.02,100,-200.02,200.02);
44 TH2F *xz3 = new TH2F("xz3","x vs z",100,-200.02,200.02,100,-200.02,200.02);
45 TH2F *yz3 = new TH2F("yz3","y vs z",100,-200.02,200.02,100,-200.02,200.02);
46 TH2F *xy4 = new TH2F("xy4","y vs x",100,-200.02,200.02,100,-200.02,200.02);
47 TH2F *xz4 = new TH2F("xz4","x vs z",100,-200.02,200.02,100,-200.02,200.02);
48 TH2F *yz4 = new TH2F("yz4","y vs z",100,-200.02,200.02,100,-200.02,200.02);
49 TH1F *chd = new TH1F("chargediff","Charge difference (Gen-Rec)",100,-1000.,1000.);
50 TH1F *chr = new TH1F("chargeratio","Charge ratio (Gen/Rec)",100,0.,0.1);
51 TH2F *chp = new TH2F("chp","Point Charge vs time",28,0.,7000.,300,0.,3000.);
52 TH2F *chh = new TH2F("chh","Hit Charge vs time",28,0.,7000.,80,0.,40.);
53 TH2F *dadth = new TH2F("dadth","danode vs dtime (hits)",560,-14000.,14000.,601,-300.5,300.5);
54 TH2F *dadt = new TH2F("dadt","danode vs dtime (points)",280,-7000.,7000.,601,-300.5,300.5);
55 TH2F *aa = new TH2F("aa","anode hit vs point",204,-0.5,50.5.,204,-.5,50.5);
56 TH2F *at = new TH2F("at","anode difference vs drift path (mm)",18,0.,36.,100,-200.,200.);
57 TH2F *tt = new TH2F("tt","time coord. difference (um) vs drift path (mm)",18,0.,36.,100,-200.,200.);
58 TH2F *at1 = new TH2F("at1","(1a) anode difference vs drift path (mm)",18,0.,36.,100,-200.,200.);
59 TH2F *tt1 = new TH2F("tt1","(1a) time coord. difference (um) vs drift path (mm)",18,0.,36.,100,-200.,200.);
60 TH2F *at2 = new TH2F("at2","(2a) anode difference vs drift path (mm)",18,0.,36.,100,-200.,200.);
61 TH2F *tt2 = new TH2F("tt2","(2a) time coord. difference (um) vs drift path (mm)",18,0.,36.,100,-200.,200.);
62 TH2F *asigma = new TH2F("asigma","anode sigma vs drift path (mm)",18,0.,36.,300,0.,300.);
63 TH2F *tsigma = new TH2F("tsigma","tau. vs drift path (mm)",18,0.,36.,200,0.,100.);
64 TH2F *asigma2 = new TH2F("asigma2","2 anode sigma vs drift path (mm)",18,0.,36.,150,0.,300.);
66 TH1F *dtrp = new TH1F("dtrp","Double track separation (microns)",40,0.,2000.);
67 TH1F *dtrpall = new TH1F("dtrpall","Double track separation (mm)",100,0.,100.);
68 TH1F *dtrh = new TH1F("dtrh","Double track separation (microns)",40,0.,2000.);
69 TH1F *dtrhall = new TH1F("dtrhall","Double track separation (mm)",100,0.,100.);
71 TH1F *p = new TH1F("p","Momentum ",100,0.,20.);
73 TH1F *effh = new TH1F("effh","Hit Multiplicity vs drift path (mm)",18,0.,36.);
74 TH1F *effp = new TH1F("effp","Point Multiplicity vs drift path (mm)",18,0.,36.);
76 TH2F *anodes = new TH2F("nanodes","Anode Multiplicity vs drift time",28,0.,7000.,5,0.5,5.5);
77 TH2F *andtsm = new TH2F("nand_ntsm","Anode Mult vs Time Mult",15,0.5,15.5,5,0.5,5.5);
78 TH2F *tsampl = new TH2F("nsampls","Sample Multiplicity vs drift time",28,0.,7000.,15,0.5,15.5);
79 TH2F *ntotal = new TH2F("ntotal","Cluster Multiplicity vs drift time",28,0.,7000.,60,0.5,60.5);
80 TH2F *clmult = new TH2F("clmult","Anode Multiplicity vs Total Multiplicity",50,0.5,50.5,5,0.5,5.5);
81 TH2F *amplit1 = new TH2F("amplit1","Point Amplitude vs drift path (mm)",28,0.,7000.,64,0.5,1024.5);
82 TH2F *amplit = new TH2F("amplit","Point Amplitude vs drift path (mm)",28,0.,7000.,60,0.5,600.5);
83 TH1F *hitpnt = new TH1F("hitpnt","Hit-Point Multiplicity",21,-10.5,10.5);
85 TH1F *nmatch = new TH1F("nmatch","Number of matched points",5,-0.5.,4.5);
86 TH1F *rec_vs_time = new TH1F("rec_vs_time","Point Rec. vs drift path",36,0.,36.);
87 TH1F *hit_vs_time = new TH1F("hit_vs_time","Hit vs drift path",36,0.,36.);
88 TH1F *rec_vs_time3 = new TH1F("rec_vs_time3","Point Rec. vs drift path",36,0.,36.);
89 TH1F *hit_vs_time3 = new TH1F("hit_vs_time3","Hit vs drift path",36,0.,36.);
90 TH1F *rec_vs_time4 = new TH1F("rec_vs_time4","Point Rec. vs drift path",36,0.,36.);
91 TH1F *hit_vs_time4 = new TH1F("hit_vs_time4","Hit vs drift path",36,0.,36.);
92 TH1F *rec_vs_time1 = new TH1F("rec_vs_time1","Point Rec. vs drift path",36,0.,36.);
93 TH1F *hit_vs_time1 = new TH1F("hit_vs_time1","Hit vs drift path",36,0.,36.);
94 TH1F *fake_vs_time = new TH1F("fake_vs_time","fake points vs drift path",36,0.,36.);
96 TH1F *noihist = new TH1F("noisehist","noise",80,10.,30.);
98 TH1F *occupancy3 = new TH1F("occupancy3","Occupancy vs Detector Number, Layer 3",20,0.5,20.5);
99 TH1F *occupancy4 = new TH1F("occupancy4","Occupancy vs Detector Number, Layer 4",20,0.5,20.5);
101 TH2F *pntmap3 = new TH2F("pntmap3","Point map Layer 3",20,0.5,20.5,10,0.5,10.5);
102 TH2F *hitmap3 = new TH2F("hitmap3","Hit map Layer 3",20,0.5,20.5,10,0.5,10.5);
103 TH2F *map3 = new TH2F("map3","Hit/Point map Layer 3",20,0.5,20.5,10,0.5,10.5);
104 TH2F *pntmap4 = new TH2F("pntmap4","Point map Layer 4",30,0.5,30.5,10,0.5,10.5);
105 TH2F *hitmap4 = new TH2F("hitmap4","Hit map Layer 4",30,0.5,30.5,10,0.5,10.5);
106 TH2F *map4 = new TH2F("map4","Hit/Point map Layer 4",30,0.5,30.5,10,0.5,10.5);
107 TH2F *xz = new TH2F("xz","X vs Z",50,-5,5.,50,-5.,5.);
108 TH2F *and_tim = new TH2F("and_tim","Tim vs Anode",30,-100,356.,30,-8000.,8000.);
109 TH2F *pand_ptim = new TH2F("pand_ptim","Tim vs Anode",30,-100,356.,30,-8000.,8000.);
111 //Int_t nanodes = 256;
112 //TH2F *mappa3hit[14][6][2];
113 //TH2F *mappa4hit[22][8][2];
114 //TH2F *mappa3pnt[14][6][2];
115 //TH2F *mappa4pnt[22][8][2];
117 for(Int_t i=0;i<22;i++) {
118 for(Int_t j=0;j<8;j++) {
119 for(Int_t k=0;k<2;k++) {
120 TString *hname = new TString("hitmap_");
121 TString *cname = new TString("pntmap_");
123 sprintf(lad,"%d",i+1);
129 sprintf(det,"%d",j+1);
135 sprintf(wng,"%d",k+1);
138 //mappa4hit[i][j][k] = new TH2F(hname->Data(),hname->Data(),nanodes,0.5,nanodes+0.5,256,0.5,256.5);
139 //mappa4pnt[i][j][k] = new TH2F(cname->Data(),cname->Data(),nanodes,0.5,nanodes+0.5,256,0.5,256.5);
141 mappa3hit[i][j][k] = new TH2F(hname->Data(),hname->Data(),nanodes,0.5,nanodes+0.5,256,0.5,256.5);
142 mappa3pnt[i][j][k] = new TH2F(cname->Data(),cname->Data(),nanodes,0.5,nanodes+0.5,256,0.5,256.5);
149 AliITS *ITS = (AliITS*) gAlice->GetModule("ITS");
150 if (!ITS) { cout << "no ITS" << endl; return; }
152 Int_t nparticles = gAlice->GetEvent(0);
154 Int_t cp[8]={0,0,0,0,0,0,0,0};
156 AliITSDetType *iDetType=ITS->DetType(1);
158 AliITSresponseSDD *res1 = (AliITSresponseSDD*)iDetType->GetResponseModel();
160 res1=new AliITSresponseSDD();
161 ITS->SetResponseModel(1,res1);
163 res1->SetZeroSupp("2D"); // 1D
164 res1->SetParamOptions("same","same");
165 //res1->SetFilenames(" ","$(ALICE_ROOT)/ITS/base.dat","$(ALICE_ROOT)/ITS/2D.dat ");
166 res1->SetCompressParam(cp);
167 res1->SetDriftSpeed(7.3);
168 Float_t vdrift = res1->DriftSpeed();
170 AliITSsegmentationSDD *seg1=(AliITSsegmentationSDD*)iDetType->GetSegmentationModel();
171 AliITSgeom *aliitsgeo = ITS->GetITSgeom();
173 Int_t cp[8]={0,0,0,0,0,0,0,0};
176 Float_t apitch = seg1->Dpz(dum);
177 Float_t tstep = seg1->Dpx(dum);
178 Float_t maxand = seg1->Npz()/2.;
179 //cout << "anodes: " << maxand << ", tstep: " << tstep << endl;
182 res1->GetNoiseParam(n,b);
183 printf("SDD: noise baseline %f %f zs option %s data type %s\n",n,b,res1->ZeroSuppOption(),res1->DataType());
184 printf("SDD: DriftSpeed %f TopValue %f\n",res1->DriftSpeed(),res1->DynamicRange());
186 res1->DiffCoeff(dif0,dif1);
187 printf("SDD: dif0 %f dif1 %f\n",dif0,dif1);
189 AliITSsimulationSDD *sim1=new AliITSsimulationSDD(seg1,res1);
190 ITS->SetSimulationModel(1,sim1);
198 for (int nev=0; nev<= evNumber2; nev++) {
200 nparticles = gAlice->GetEvent(nev);
201 gAlice->SetEvent(nev);
203 cout << "nparticles " <<nparticles<<endl;
204 if (nev < evNumber1) continue;
205 if (nparticles <= 0) return;
209 AliITSRecPoint *itsPnt = 0;
210 AliITSRawClusterSDD *itsClu = 0;
212 // Reset Event Counters
213 Int_t nGoodTotalHits = 0;
214 Int_t nGoodTotalPnts = 0;
216 // Get Hit, Cluster & Recpoints Tree Pointers
218 TTree *TH = gAlice->TreeH();
219 Int_t nenthit=TH->GetEntries();
220 printf("Found %d entries in the Hit tree (must be one per track per event!)\n",nenthit);
223 TTree *TC=ITS->TreeC();
224 Int_t nentclu=TC->GetEntries();
225 printf("Found %d entries in the Hit tree (must be one per module per event!)\n",nentclu);
227 TTree *TR = gAlice->TreeR();
228 Int_t nentrec=TR->GetEntries();
229 printf("Found %d entries in the Rec tree (must be one per module per event!)\n",nentrec);
231 // Get Pointers to Clusters & Recpoints TClonesArrays
233 Int_t iDet = 1; // 1 = SDD
235 TClonesArray *ITSclu = ITS->ClustersAddress(iDet);
236 TClonesArray *ITSrec = ITS->RecPoints();
247 ITS->InitModules(-1,nmodules);
248 ITS->FillModules(nev,0,nmodules,"","");
250 TObjArray *fITSmodules = ITS->GetModules();
252 Int_t first = aliitsgeo->GetStartDet(iDet);
253 Int_t last = aliitsgeo->GetLastDet(iDet);
254 printf("det type %d first, last %d %d \n",iDet,first,last);
256 for (Int_t mod=0; mod<last-first+1; mod++) {
257 cout << "Module: " << mod+1 << endl;
258 TTree *TR = gAlice->TreeR();
259 Int_t nentrec=TR->GetEntries();
260 TClonesArray *ITSrec = ITS->RecPoints();
262 ITS->ResetClusters();
264 ITS->ResetRecPoints();
265 nbytes += TR->GetEvent(mod);
267 Int_t nrecp = ITSrec->GetEntries();
269 //if (nrecp) printf("Found %d rec points for module %d\n",nrecp,mod);
270 if (!nrecp) continue;
272 Int_t nrecc = ITSclu->GetEntries();
274 //if (nrecc) printf("Found %d clusters for module %d\n",nrecc,mod);
276 Int_t nrecp = ITSrec->GetEntries();
277 Int_t startSDD = aliitsgeo->GetStartSDD();
278 Int_t *flagP = new Int_t [nrecp];
279 memset( flagP, 0, sizeof(Int_t)*nrecp );
281 //printf("point loop\n");
284 for (Int_t pnt=0;pnt<nrecp;pnt++) {
285 itsPnt = (AliITSRecPoint*)ITSrec->At(pnt);
286 if(!itsPnt) continue;
287 itsClu = (AliITSRawClusterSDD*)ITSclu->At(pnt);
288 if(!itsClu) continue;
289 //itsClu->PrintInfo();
296 aliitsgeo->GetModuleId(mod+first,pntlayer,pntladder,pntdetector);
297 Int_t pntmult = itsClu->Samples();
298 Int_t pntands = itsClu->Anodes();
299 Float_t pnttime = itsClu->T();
300 Float_t pntanod = itsClu->A();
301 Float_t pntchrg = itsClu->Q();
302 Float_t pntampl = itsClu->PeakAmpl();
303 Float_t pntpath = pnttime*vdrift/1000.;
306 if(itsClu->Anodes() != 0.) {
307 wy = pntmult/((Float_t) pntands);
309 clmult->Fill((Float_t)pntmult,(Float_t) pntands);
310 ntotal->Fill(pnttime,(Float_t)pntmult);
311 tsampl->Fill(pnttime,wy);
312 amplit->Fill(pnttime,pntampl);
313 amplit1->Fill(pnttime,pntampl);
315 // Detector occupancy
318 occupancy3->Fill((Float_t) pntdetector,(Float_t) pntmult);
319 //mappa3pnt[pntladder-1][pntdetector-1][0]->Fill(pntanod,pnttime);
322 occupancy4->Fill((Float_t) pntdetector,(Float_t) pntmult);
323 //mappa4pnt[pntladder-1][pntdetector-1][0]->Fill(pntanod,pnttime);
326 // Point Efficiency vs time.
329 anodes->Fill(pnttime,pntands);
330 andtsm->Fill(wy,pntands);
331 chp->Fill(pnttime,pntchrg);
334 //printf("hit loop\n");
336 Float_t sddLength = seg1->Dx();
337 Float_t sddWidth = seg1->Dz();
338 Float_t driftSpeed=res1->DriftSpeed();
342 AliITSmodule *Mod = (AliITSmodule *)fITSmodules->At(mod+first);
343 Int_t nhits = Mod->GetNhits();
344 for (Int_t hit=0;hit<nhits;hit++) {
345 itsHit = (AliITShit*)Mod->GetHit(hit);
351 Float_t DepEnergy = 100000.*itsHit->GetIonization();
352 AliITShit *itsHit1 = 0;
353 if(DepEnergy == 0.) {
355 if(hit == nhits) break;
356 itsHit1 = (AliITShit*) Mod->GetHit(hit);
357 avx = itsHit1->GetXG();
358 avy = itsHit1->GetYG();
359 avz = itsHit1->GetZG();
362 avx += itsHit->GetXG();
363 avy += itsHit->GetYG();
364 avz += itsHit->GetZG();
365 if(DepEnergy == 0.) {
370 if(ifl == 0) continue;
372 Float_t px; Float_t py; Float_t pz;
373 itsHit->GetMomentumG(px,py,pz);
374 Float_t ptot = TMath::Sqrt(px*px+py*py+pz*pz);
376 if(ptot < 0.05) continue;
378 Int_t Layer = itsHit->GetLayer();
379 Int_t Ladder = itsHit->GetLadder();
380 Int_t Det = itsHit->GetDetector();
384 Float_t x = itsHit->GetXL();
385 Float_t z = itsHit->GetZL();
387 seg1->GetPadTxz(x,z);
390 and_tim->Fill(And,Tim);
397 x1 = itsHit1->GetXL();
398 z1 = itsHit1->GetZL();
400 seg1->GetPadTxz(x1,z1);
403 and_tim->Fill(And1,Tim1);
405 Float_t DepEnergy = 100000.*itsHit->GetIonization();
406 if(DepEnergy == 0.) DepEnergy = 100000.*itsHit1->GetIonization();
407 if(DepEnergy < 5.) continue;
415 if(And < 0. || And > maxand) { cout << "And: " << And << endl; continue; }
416 Float_t path = TMath::Abs(Tim)*vdrift/1000.;
417 hit_vs_time->Fill(path);
418 if(Layer==3) hit_vs_time3->Fill(path);
419 if(Layer==4) hit_vs_time4->Fill(path);
424 //if(Layer == 3) mappa3hit[Ladder-1][Det-1][0]->Fill(And,Tim);
425 //if(Layer == 4) mappa4hit[Ladder-1][Det-1][0]->Fill(And,Tim);
428 Float_t ww = DepEnergy;
429 chh->Fill(TMath::Abs(Tim),ww);
432 Float_t diffmin = 100000.;
435 //printf("point loop\n");
437 for (Int_t pnt=0;pnt<nrecp;pnt++) {
438 itsPnt = (AliITSRecPoint*)ITSrec->At(pnt);
439 if(!itsPnt) continue;
440 itsClu = (AliITSRawClusterSDD*)ITSclu->At(pnt);
441 if(!itsClu) continue;
446 aliitsgeo->GetModuleId(mod+first,LayerP,LadderP,DetP);
447 Int_t LayerH = itsHit->GetLayer();
448 Int_t LadderH = itsHit->GetLadder();
449 Int_t DetH = itsHit->GetDetector();
450 if(LayerH != LayerP) continue;
451 if(LadderH != LadderP) continue;
452 if(DetH != DetP) continue;
454 Float_t Pand = (Float_t) itsClu->A();
455 if(Pand < 0 || Pand > maxand) { cout << "Pand: " << Pand << endl; continue; }
456 Float_t Ptim = (Float_t) itsClu->T();
457 Float_t Pwng = (Float_t) itsClu->W();
458 if(Pwng == 1) Ptim *= -1.;
459 pand_ptim->Fill(Pand,Ptim);
460 Float_t adiff = And-Pand;
461 Float_t tdiff = Tim-Ptim;
463 printf("tim %f\n",Tim);
464 printf("and %f\n",And);
466 if(Pwng == 1) tdiff *=-1.;
467 local1->Fill(adiff,tdiff);
469 if(TMath::Abs(adiff) >= 1) continue;
470 if(TMath::Abs(tdiff) >= 100) continue;
472 Float_t apdiff = adiff*apitch;
473 Float_t tpdiff = tdiff*vdrift;
475 Float_t diff = TMath::Sqrt( apdiff*apdiff+tpdiff*tpdiff );
476 if( diff < diffmin ){
481 if(TMath::Abs(adiff) < 1. && TMath::Abs(tdiff) < 100.) {
485 if( flagP[pntmin] == 1) continue;
487 itsClu = (AliITSRawClusterSDD*)ITSclu->At( pntmin );
488 Float_t Pand = (Float_t) itsClu->A();
489 Float_t Ptim = (Float_t) itsClu->T();
490 Float_t Pwng = (Float_t) itsClu->W();
491 Float_t sigma = itsClu->Asigma();
492 Float_t tau = itsClu->Tsigma();
493 Int_t pntands = itsClu->Anodes();
494 if(Pwng == 1) Ptim *= -1.;
495 Float_t adiff = And-Pand;
496 Float_t tdiff = Tim-Ptim;
497 if(Pwng == 1) tdiff *=-1.;
498 local->Fill(adiff,tdiff);
500 Float_t dpath = Ptim*vdrift/1000.;
501 Float_t dpathh = Tim*vdrift/1000.;
502 Float_t adpath = TMath::Abs(dpath);
503 Float_t adpathh = TMath::Abs(dpathh);
504 Float_t apdiff = adiff*apitch;
505 Float_t tpdiff = tdiff*vdrift;
508 Int_t pntands = itsClu->Anodes();
510 at1->Fill(adpath,apdiff);
511 tt1->Fill(adpath,tpdiff);
514 at2->Fill(adpath,apdiff);
515 tt2->Fill(adpath,tpdiff);
518 at->Fill(adpathh,apdiff);
519 tt->Fill(adpathh,tpdiff);
520 asigma->Fill(adpathh,sigma);
521 tsigma->Fill(adpathh,tau);
522 if(pntands == 2) asigma2->Fill(adpathh,sigma);
524 Float_t *lP = new Float_t[3];
525 lP[0] = itsPnt->GetX();
527 lP[2] = itsPnt->GetZ();
528 Float_t *gP = new Float_t[3];
529 aliitsgeo->LtoG(LayerH,LadderH,DetH,lP,gP);
530 Float_t dx = avx - gP[0];
531 Float_t dy = avy - gP[1];
532 Float_t dz = avz - gP[2];
536 Float_t pntchrg = itsClu->Q();
537 Float_t dq = DepEnergy/0.122 - pntchrg;
539 if(pntchrg != 0) rq = DepEnergy/0.122/((Float_t) pntchrg);
544 } else if(LayerH == 4) {
550 if(rq != 0.) chr->Fill(rq);
552 rec_vs_time->Fill(adpathh);
553 if(Layer==3) rec_vs_time3->Fill(adpathh);
554 if(Layer==4) rec_vs_time4->Fill(adpathh);
557 nmatch->Fill(inmatches);
560 if(nGoodHits != nGoodPnts) {
561 printf("module: %d",mod+1);
562 printf(", nGoodHits: %d",nGoodHits);
563 printf(", nGoodPnts: %d\n",nGoodPnts);
565 Float_t nHP = (Float_t) nGoodHits-nGoodPnts;
569 if(nGoodHits != 0) www = nGoodPnts/((Float_t) nGoodHits);
571 pntmap3->Fill(Ladder,Det,(Float_t) nGoodPnts);
572 hitmap3->Fill(Ladder,Det,(Float_t) nGoodHits);
573 map3->Fill(Ladder,Det,www);
576 pntmap4->Fill(Ladder,Det,(Float_t) nGoodPnts);
577 hitmap4->Fill(Ladder,Det,(Float_t) nGoodHits);
578 map4->Fill(Ladder,Det,www);
581 //printf("double hit loop\n");
585 wh /= (((Float_t) nGoodHits)*((Float_t) nGoodHits)-1)/2.;
589 Int_t *flag = new Int_t[nhits];
590 Int_t nGoodHitsOK = 0;
591 for (Int_t hit=0;hit<nhits;hit++) {
593 itsHit = (AliITShit*)Mod->GetHit(hit);
598 Float_t DepEnergy = 100000.*itsHit->GetIonization();
599 AliITShit *itsHit1 = 0;
600 if(DepEnergy == 0.) {
603 if(hit == nhits) break;
604 itsHit1 = (AliITShit*) Mod->GetHit(hit);
605 avx = itsHit1->GetXG();
606 avy = itsHit1->GetYG();
607 avz = itsHit1->GetZG();
610 avx += itsHit->GetXG();
611 avy += itsHit->GetYG();
612 avz += itsHit->GetZG();
613 if(DepEnergy == 0.) {
618 if(DepEnergy < 5. && DepEnergy > 0.) continue;
619 if(ifl == 0) continue;
621 Float_t px; Float_t py; Float_t pz;
622 itsHit->GetMomentumG(px,py,pz);
623 Float_t ptot = TMath::Sqrt(px*px+py*py+pz*pz);
624 if(ptot < 0.05) continue;
626 for (Int_t hit1=hit+1;hit1<nhits;hit1++) {
627 itsHit2 = (AliITShit*)Mod->GetHit(hit1);
633 Float_t DepEnergy2 = 100000.*itsHit2->GetIonization();
634 AliITShit *itsHit3 = 0;
635 if(DepEnergy2 == 0.) {
637 itsHit3 = (AliITShit*) Mod->GetHit(hit1);
638 avx2 = itsHit3->GetXG();
639 avy2 = itsHit3->GetYG();
640 avz2 = itsHit3->GetZG();
643 avx2 += itsHit2->GetXG();
644 avy2 += itsHit2->GetYG();
645 avz2 += itsHit2->GetZG();
646 if(DepEnergy2 == 0.) {
651 if(DepEnergy2 < 5. && DepEnergy2 > 0.) continue;
652 if(itsHit->GetLayer() != itsHit2->GetLayer()) continue;
653 if(itsHit->GetLadder() != itsHit2->GetLadder()) continue;
654 if(itsHit->GetDetector() != itsHit2->GetDetector()) continue;
655 if(ifl2 == 0) continue;
657 Float_t px1; Float_t py1; Float_t pz1;
658 itsHit2->GetMomentumG(px1,py1,pz1);
659 Float_t ptot1 = TMath::Sqrt(px1*px1+py1*py1+pz1*pz1);
660 if(ptot1 < 0.05) continue;
664 Float_t x = itsHit->GetXL();
665 Float_t z = itsHit->GetZL();
666 seg1->GetPadTxz(x,z);
669 if(And < 0 || And > maxand) continue;
672 Float_t x2 = itsHit2->GetXL();
673 Float_t z2 = itsHit2->GetZL();
674 seg1->GetPadTxz(x2,z2);
677 if(And2 < 0 || And2 > maxand) continue;
678 Float_t da = apitch*(And-And2);
679 Float_t dt = vdrift*(Tim-Tim2);
680 Float_t danh = And-And2;
681 Float_t dtmh = Tim-Tim2;
682 Float_t dist = TMath::Sqrt(da*da+dt*dt);
684 Float_t wx = dt*clock/(1000.*vdrift);
685 Float_t wy = da/apitch;
688 if(dist<20.) { cout << "skip hit " << hit1 << endl; flag[hit] = 1; }
690 if(ifl == 1 && ifl2 == 1) {
691 if(dist>10)dtrh->Fill(dist,wh);
692 dtrhall->Fill(dist/1000.,wh);
693 dadth->Fill(dtmh,danh);
695 } // end cluster loop
697 Float_t path = TMath::Abs(Tim)*vdrift/1000.;
698 if(flag[hit] == 0) { nGoodHitsOK++; hit_vs_time1->Fill(path); }
701 printf("nGoodHits: %d",nGoodHits);
702 printf(", nGoodHitsOK: %d\n",nGoodHitsOK);
704 //printf("cluster loop\n");
706 AliITSRawClusterSDD *itsCluFake = 0;
707 Int_t nGoodPntsOK = 0;
708 for( int ip=0; ip<nrecp; ip++) {
709 itsClu = (AliITSRawClusterSDD*)ITSclu->At(ip);
710 if(!itsClu) continue;
711 Float_t Ptim = (Float_t) itsClu->T();
712 Float_t dpath = Ptim*vdrift/1000.;
713 Float_t adpath = TMath::Abs(dpath);
714 if( flagP[ip] == 1) { nGoodPntsOK++; rec_vs_time1->Fill(dpath); }
716 cout << "ip: " << ip << endl;
717 itsCluFake = (AliITSRawClusterSDD*)ITSclu->At( ip );
718 if(!itsCluFake) continue;
719 cout << "pointer: " << itsCluFake << endl;
720 itsCluFake->PrintInfo();
721 Float_t Ptim = (Float_t) itsCluFake->T();
722 Float_t dpath = Ptim*vdrift/1000.;
723 fake_vs_time->Fill(dpath);
729 wp /= (((Float_t) nGoodPntsOK)*((Float_t) nGoodPntsOK)-1)/2.;
733 //printf("double cluster loop\n");
734 for (Int_t pnt=0;pnt<nrecp;pnt++) {
735 if( flagP[pnt] == 0) continue;
736 itsPnt = (AliITSRecPoint*)ITSrec->At(pnt);
737 if(!itsPnt) continue;
738 itsClu = (AliITSRawClusterSDD*)ITSclu->At(pnt);
739 if(!itsClu) continue;
740 AliITSRecPoint *itsPnt1;
741 AliITSRawClusterSDD *itsClu1;
742 for (Int_t pnt1=pnt+1;pnt1<nrecp;pnt1++) {
743 if( flagP[pnt1] == 0) continue;
744 itsPnt1 = (AliITSRecPoint*)ITSrec->At(pnt1);
745 if(!itsPnt1) continue;
746 itsClu1 = (AliITSRawClusterSDD*)ITSclu->At(pnt1);
747 if(!itsClu1) continue;
749 Float_t dan = itsClu->A()-itsClu1->A();
750 Float_t Pwng = (Float_t) itsClu->W();
751 Float_t dt1 = itsClu->T();
752 if(Pwng == 1) dt1 *= -1.;
753 Float_t dt2 = itsClu1->T();
754 Float_t Pwng = (Float_t) itsClu1->W();
755 if(Pwng == 1) dt2 *= -1.;
756 Float_t dtm = dt1-dt2;
758 Float_t dap = apitch*(itsClu->A()-itsClu1->A());
759 Float_t dtp = vdrift*(dt1-dt2);
760 Float_t distp = TMath::Sqrt(dap*dap+dtp*dtp);
761 if(TMath::Abs(dan) < 1.) dclutim->Fill(dtp);
762 if(TMath::Abs(dan) > 10.) dfkctim->Fill(dtp);
764 Float_t wx = dtp*clock;
765 Float_t wy = dap/apitch;
768 dtrp->Fill(distp,wp);
769 dtrpall->Fill(distp/1000.,wp);
774 } // end loop modules
777 gAlice->CleanDetectors();
781 cout << "open output file" << endl;
782 TFile fhistos("SDD_histos_test.root","RECREATE");
818 rec_vs_time->Write();
819 hit_vs_time->Write();
820 hit_vs_time1->Write();
821 rec_vs_time1->Write();
822 rec_vs_time3->Write();
823 hit_vs_time3->Write();
824 rec_vs_time4->Write();
825 hit_vs_time4->Write();
826 fake_vs_time->Write();
847 for(Int_t i=0;i<22;i++) {
848 for(Int_t j=0;j<8;j++) {
849 for(Int_t k=0;k<2;k++) {
850 //mappa4hit[i][j][k]->Write();
851 //mappa4pnt[i][j][k]->Write();
853 //mappa3hit[i][j][k]->Write();
854 //mappa3pnt[i][j][k]->Write();