Updated version from P. Cerello
[u/mrichter/AliRoot.git] / ITS / ITSsddanalysis.C
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();
23   }
24   
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");
29   
30   // Get AliRun object from file or create it if not on file
31   //   if (!gAlice) {
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");
35   //}
36   
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.);
65   
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.);
70
71   TH1F *p = new TH1F("p","Momentum ",100,0.,20.);
72   
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.);
75   
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);
84   
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.);
95   
96   TH1F *noihist = new TH1F("noisehist","noise",80,10.,30.);
97   
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);
100   
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.);
110
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];
116   /*
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_");
122     Char_t lad[2];
123     sprintf(lad,"%d",i+1);
124     hname->Append(lad);
125     hname->Append("_");
126     cname->Append(lad);
127     cname->Append("_");
128     Char_t det[2];
129     sprintf(det,"%d",j+1);
130     hname->Append(det);
131     hname->Append("_");
132     cname->Append(det);
133     cname->Append("_");
134     Char_t wng[2];
135     sprintf(wng,"%d",k+1);
136     hname->Append(wng);
137     cname->Append(wng);
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);
140     if(i<14 && j<6) {
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);
143     }
144     }
145     }
146     }
147   */
148  
149   AliITS *ITS  = (AliITS*) gAlice->GetModule("ITS");
150   if (!ITS) { cout << "no ITS" << endl; return; }
151     
152   Int_t nparticles = gAlice->GetEvent(0);
153
154   Int_t cp[8]={0,0,0,0,0,0,0,0};
155   
156   AliITSDetType *iDetType=ITS->DetType(1);
157   
158   AliITSresponseSDD *res1 = (AliITSresponseSDD*)iDetType->GetResponseModel();
159   if (!res1) {
160     res1=new AliITSresponseSDD();
161     ITS->SetResponseModel(1,res1);
162   }
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();
169   
170   AliITSsegmentationSDD *seg1=(AliITSsegmentationSDD*)iDetType->GetSegmentationModel();
171   AliITSgeom *aliitsgeo = ITS->GetITSgeom();
172   
173   Int_t cp[8]={0,0,0,0,0,0,0,0};
174     
175   Int_t dum = 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;
180   
181   Float_t n,b;
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());
185   Float_t dif0,dif1;
186   res1->DiffCoeff(dif0,dif1);
187   printf("SDD: dif0 %f dif1 %f\n",dif0,dif1);
188   
189   AliITSsimulationSDD *sim1=new AliITSsimulationSDD(seg1,res1);
190   ITS->SetSimulationModel(1,sim1);
191
192   //
193   //   Loop over events
194   //
195   
196   Int_t Nh=0;
197   Int_t Nh1=0;
198   for (int nev=0; nev<= evNumber2; nev++) {
199     if(nev>0) {
200       nparticles = gAlice->GetEvent(nev);
201       gAlice->SetEvent(nev);
202     }
203     cout << "nparticles  " <<nparticles<<endl;
204     if (nev < evNumber1) continue;
205     if (nparticles <= 0) return;
206
207     // Reset Pointers
208     AliITShit *itsHit;
209     AliITSRecPoint *itsPnt = 0;
210     AliITSRawClusterSDD *itsClu = 0;
211
212     // Reset Event Counters
213     Int_t nGoodTotalHits = 0;
214     Int_t nGoodTotalPnts = 0;
215         
216     // Get Hit, Cluster & Recpoints Tree Pointers
217     
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);
221     
222     ITS->GetTreeC(nev);
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);
226     
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);
230     
231     // Get Pointers to Clusters & Recpoints TClonesArrays
232     
233     Int_t iDet = 1;  // 1 = SDD
234
235     TClonesArray *ITSclu  = ITS->ClustersAddress(iDet);
236     TClonesArray *ITSrec  = ITS->RecPoints(); 
237     
238     // check recpoints
239     
240     Int_t nbytes = 0;
241     Int_t totpoints = 0;
242     Int_t totclust = 0;
243     
244     // check hits
245     
246     Int_t nmodules=0;    
247     ITS->InitModules(-1,nmodules); 
248     ITS->FillModules(nev,0,nmodules,"","");
249     
250     TObjArray *fITSmodules = ITS->GetModules();
251     
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);
255     
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(); 
261       
262       ITS->ResetClusters();
263       TC->GetEvent(mod);
264       ITS->ResetRecPoints();
265       nbytes += TR->GetEvent(mod);
266       
267       Int_t nrecp = ITSrec->GetEntries();
268       totpoints += nrecp;
269       //if (nrecp) printf("Found %d rec points for module %d\n",nrecp,mod);
270       if (!nrecp) continue;
271       
272       Int_t nrecc = ITSclu->GetEntries();
273       totclust += nrecc;
274       //if (nrecc) printf("Found %d clusters for module %d\n",nrecc,mod);
275       
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 ); 
280
281       //printf("point loop\n");
282
283       Int_t nGoodPnts = 0;
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();  
290         nGoodPnts++;
291         nGoodTotalPnts++;
292         
293         Int_t pntlayer;
294         Int_t pntladder;
295         Int_t pntdetector;
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.;
304
305         Float_t wy = 0.;
306         if(itsClu->Anodes() != 0.) {
307           wy = pntmult/((Float_t) pntands);
308         }
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);
314         
315         //  Detector  occupancy
316         
317         if(pntlayer == 3) {
318           occupancy3->Fill((Float_t) pntdetector,(Float_t) pntmult);
319           //mappa3pnt[pntladder-1][pntdetector-1][0]->Fill(pntanod,pnttime);
320         }
321         if(pntlayer == 4) {
322           occupancy4->Fill((Float_t) pntdetector,(Float_t) pntmult);
323           //mappa4pnt[pntladder-1][pntdetector-1][0]->Fill(pntanod,pnttime);
324         }
325         
326         //  Point  Efficiency  vs  time.
327         
328         effp->Fill(pntpath);
329         anodes->Fill(pnttime,pntands);
330         andtsm->Fill(wy,pntands);
331         chp->Fill(pnttime,pntchrg);
332       }
333       
334       //printf("hit loop\n");
335
336       Float_t sddLength = seg1->Dx();
337       Float_t sddWidth = seg1->Dz();
338       Float_t driftSpeed=res1->DriftSpeed();    
339       
340       Int_t nGoodHits = 0;
341
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);
346         
347         Float_t avx = 0.;
348         Float_t avy = 0.;
349         Float_t avz = 0.;
350         Int_t ifl = 0;
351         Float_t DepEnergy = 100000.*itsHit->GetIonization();
352         AliITShit *itsHit1 = 0;
353         if(DepEnergy == 0.) { 
354           hit++;
355           if(hit == nhits) break;
356           itsHit1 = (AliITShit*) Mod->GetHit(hit);
357           avx = itsHit1->GetXG();
358           avy = itsHit1->GetYG();
359           avz = itsHit1->GetZG();
360           ifl = 1;
361         }
362         avx += itsHit->GetXG();
363         avy += itsHit->GetYG();
364         avz += itsHit->GetZG();
365         if(DepEnergy == 0.) { 
366           avx /= 2.;
367           avy /= 2.;
368           avz /= 2.;
369         }
370         if(ifl == 0) continue;
371         
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);
375         p->Fill(ptot*100);
376         if(ptot < 0.05) continue;
377         
378         Int_t Layer = itsHit->GetLayer();
379         Int_t Ladder = itsHit->GetLadder();
380         Int_t Det = itsHit->GetDetector();
381         
382         Float_t And;
383         Float_t Tim;
384         Float_t x = itsHit->GetXL();
385         Float_t z = itsHit->GetZL();
386         xz->Fill(z,x);
387         seg1->GetPadTxz(x,z);
388         And = z;
389         Tim = x*tstep;
390         and_tim->Fill(And,Tim);
391
392         Float_t And1;
393         Float_t Tim1;
394         Float_t x1;
395         Float_t z1;
396         if(itsHit1) {
397           x1 = itsHit1->GetXL();
398           z1 = itsHit1->GetZL();
399           xz->Fill(z1,x1);
400           seg1->GetPadTxz(x1,z1);
401           And1 = z1;
402           Tim1 = x1*tstep;
403           and_tim->Fill(And1,Tim1);
404         }
405         Float_t DepEnergy = 100000.*itsHit->GetIonization();
406         if(DepEnergy == 0.) DepEnergy = 100000.*itsHit1->GetIonization(); 
407         if(DepEnergy < 5.) continue;    
408
409         if(itsHit1) {
410           Tim += Tim1;
411           Tim /= 2.;
412           And += And1;
413           And /= 2.;
414         }
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);
420         
421         nGoodHits++;
422         nGoodTotalHits++;
423         
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);
426         
427         effh->Fill(path);
428         Float_t ww = DepEnergy;
429         chh->Fill(TMath::Abs(Tim),ww);
430         Int_t inmatches = 0;
431         
432         Float_t diffmin = 100000.;
433         Int_t pntmin = -1;
434
435         //printf("point loop\n");
436         for (Int_t pnt=0;pnt<nrecp;pnt++) {
437           itsPnt  = (AliITSRecPoint*)ITSrec->At(pnt);
438           if(!itsPnt) continue;
439           itsClu  = (AliITSRawClusterSDD*)ITSclu->At(pnt);
440           if(!itsClu) continue;
441           
442           Int_t LayerP;
443           Int_t LadderP;
444           Int_t DetP;
445           aliitsgeo->GetModuleId(mod+first,LayerP,LadderP,DetP);
446           Int_t LayerH =  itsHit->GetLayer();
447           Int_t LadderH =  itsHit->GetLadder();
448           Int_t DetH =  itsHit->GetDetector();
449           if(LayerH != LayerP) continue;
450           if(LadderH != LadderP) continue;
451           if(DetH != DetP) continue;
452
453           Float_t Pand = (Float_t) itsClu->A();
454           if(Pand < 0 || Pand > maxand) { cout << "Pand: " << Pand << endl; continue; }
455           Float_t Ptim = (Float_t) itsClu->T();
456           Float_t Pwng = (Float_t) itsClu->W();
457           if(Pwng == 1) Ptim *= -1.;
458           pand_ptim->Fill(Pand,Ptim);
459           Float_t adiff = And-Pand;
460           Float_t tdiff = Tim-Ptim;  
461           if(And < 0) {
462             printf("tim %f\n",Tim);
463             printf("and %f\n",And);
464           }
465           if(Pwng == 1) tdiff *=-1.;
466           local1->Fill(adiff,tdiff);
467           
468           if(TMath::Abs(adiff) >= 1) continue;
469           if(TMath::Abs(tdiff) >= 100) continue;
470
471           Float_t apdiff = adiff*apitch;
472           Float_t tpdiff = tdiff*vdrift;
473          
474           Float_t diff = TMath::Sqrt( apdiff*apdiff+tpdiff*tpdiff ); 
475           if( diff < diffmin ){
476             diffmin = diff;
477             pntmin = pnt;
478           }
479                   
480           if(TMath::Abs(adiff) < 1. && TMath::Abs(tdiff) < 100.) {
481             inmatches++;
482           }
483         } // loop (points)
484
485         if( pntmin > -1 ) {
486           if( flagP[pntmin] == 1) continue;
487           flagP[pntmin] = 1;
488           itsClu  = (AliITSRawClusterSDD*)ITSclu->At( pntmin );
489           Float_t Pand = (Float_t) itsClu->A();
490           Float_t Ptim = (Float_t) itsClu->T();
491           Float_t Pwng = (Float_t) itsClu->W();
492           Float_t sigma = itsClu->Asigma();
493           Float_t tau = itsClu->Tsigma();
494           Int_t pntands = itsClu->Anodes();
495           if(Pwng == 1) Ptim *= -1.;
496           Float_t adiff = And-Pand;
497           Float_t tdiff = Tim-Ptim;
498           if(Pwng == 1) tdiff *=-1.;
499           local->Fill(adiff,tdiff);
500           
501           Float_t dpath = Ptim*vdrift/1000.;
502           Float_t dpathh = Tim*vdrift/1000.;
503           Float_t adpath = TMath::Abs(dpath);
504           Float_t adpathh = TMath::Abs(dpathh);
505           Float_t apdiff = adiff*apitch;
506           Float_t tpdiff = tdiff*vdrift;
507           aa->Fill(Pand,And);
508           
509           Int_t pntands = itsClu->Anodes();
510           if(pntands == 1) {
511             at1->Fill(adpath,apdiff);
512             tt1->Fill(adpath,tpdiff);
513           } 
514           if(pntands == 2) {
515             at2->Fill(adpath,apdiff);
516             tt2->Fill(adpath,tpdiff);
517           } 
518           
519           at->Fill(adpathh,apdiff);
520           tt->Fill(adpathh,tpdiff);
521           asigma->Fill(adpathh,sigma);
522           tsigma->Fill(adpathh,tau);
523           if(pntands == 2) asigma2->Fill(adpathh,sigma);
524           
525           Float_t *lP = new Float_t[3];
526           lP[0] = itsPnt->GetX();
527           lP[1] = 0.;
528           lP[2] = itsPnt->GetZ();
529           Float_t *gP = new Float_t[3];
530           aliitsgeo->LtoG(LayerH,LadderH,DetH,lP,gP);
531           Float_t dx = avx - gP[0];
532           Float_t dy = avy - gP[1];
533           Float_t dz = avz - gP[2];
534           delete lP;
535           delete gP;
536           
537           Float_t pntchrg = itsClu->Q();
538           Float_t dq = DepEnergy/0.122 - pntchrg;
539           Float_t rq = 0;
540           if(pntchrg != 0) rq = DepEnergy/0.122/((Float_t) pntchrg);
541           if(LayerH == 3) {
542             xy3->Fill(dx,dy);
543             xz3->Fill(dz,dx);
544             yz3->Fill(dz,dy);
545           } else if(LayerH == 4) {
546             xy4->Fill(dx,dy);
547             xz4->Fill(dz,dx);
548             yz4->Fill(dz,dy);
549           }
550           chd->Fill(dq);
551           if(rq != 0.) chr->Fill(rq); 
552           
553           rec_vs_time->Fill(adpathh);
554           if(Layer==3) rec_vs_time3->Fill(adpathh);
555           if(Layer==4) rec_vs_time4->Fill(adpathh);
556         }
557         
558         nmatch->Fill(inmatches);
559       }  // loop hits
560       
561       if(nGoodHits != nGoodPnts) {
562         printf("module: %d",mod+1);
563         printf(", nGoodHits: %d",nGoodHits);
564         printf(", nGoodPnts: %d\n",nGoodPnts);
565       }
566       Float_t nHP = (Float_t) nGoodHits-nGoodPnts;
567       hitpnt->Fill(nHP);
568       
569       Int_t www = 0.;
570       if(nGoodHits != 0) www = nGoodPnts/((Float_t) nGoodHits);
571       if(Layer == 3) {
572         pntmap3->Fill(Ladder,Det,(Float_t) nGoodPnts);
573         hitmap3->Fill(Ladder,Det,(Float_t) nGoodHits);
574         map3->Fill(Ladder,Det,www);
575       }
576       if(Layer == 4) {
577         pntmap4->Fill(Ladder,Det,(Float_t) nGoodPnts);
578         hitmap4->Fill(Ladder,Det,(Float_t) nGoodHits);
579         map4->Fill(Ladder,Det,www);
580       }
581       
582       //printf("double hit loop\n");
583       
584       Stat_t wh = 1.;
585       if(nGoodHits > 1) 
586         wh /= (((Float_t) nGoodHits)*((Float_t) nGoodHits)-1)/2.;
587       else
588         wh = 0.;
589       
590       Int_t *flag = new Int_t[nhits];
591       Int_t nGoodHitsOK = 0;
592       for (Int_t hit=0;hit<nhits;hit++) {
593         flag[hit] = 0;
594         itsHit   = (AliITShit*)Mod->GetHit(hit);
595         Float_t avx = 0.;
596         Float_t avy = 0.;
597         Float_t avz = 0.;
598         Int_t ifl = 0;
599         Float_t DepEnergy = 100000.*itsHit->GetIonization();
600         AliITShit *itsHit1 = 0;
601         if(DepEnergy == 0.) { 
602           hit++;
603           flag[hit] = 0;
604           if(hit == nhits) break;
605           itsHit1 = (AliITShit*) Mod->GetHit(hit);
606           avx = itsHit1->GetXG();
607           avy = itsHit1->GetYG();
608           avz = itsHit1->GetZG();
609           ifl = 1;
610         }
611         avx += itsHit->GetXG();
612         avy += itsHit->GetYG();
613         avz += itsHit->GetZG();
614         if(DepEnergy == 0.) { 
615           avx /= 2.;
616           avy /= 2.;
617           avz /= 2.;
618         }
619         if(DepEnergy < 5. && DepEnergy > 0.) continue;
620         if(ifl == 0) continue;
621
622         Float_t px; Float_t py; Float_t pz;
623         itsHit->GetMomentumG(px,py,pz);
624         Float_t ptot = TMath::Sqrt(px*px+py*py+pz*pz);
625         if(ptot < 0.05) continue;
626         
627         for (Int_t hit1=hit+1;hit1<nhits;hit1++) {
628           itsHit2   = (AliITShit*)Mod->GetHit(hit1);
629           
630           Float_t avx2 = 0.;
631           Float_t avy2 = 0.;
632           Float_t avz2 = 0.;
633           Int_t ifl2 = 0;
634           Float_t DepEnergy2 = 100000.*itsHit2->GetIonization();
635           AliITShit *itsHit3 = 0;
636           if(DepEnergy2 == 0.) { 
637             hit1++;
638             itsHit3 = (AliITShit*) Mod->GetHit(hit1);
639             avx2 = itsHit3->GetXG();
640             avy2 = itsHit3->GetYG();
641             avz2 = itsHit3->GetZG();
642             ifl2 = 1;
643           }
644           avx2 += itsHit2->GetXG();
645           avy2 += itsHit2->GetYG();
646           avz2 += itsHit2->GetZG();
647           if(DepEnergy2 == 0.) { 
648             avx2 /= 2.;
649             avy2 /= 2.;
650             avz2 /= 2.;
651           }
652           if(DepEnergy2 < 5. && DepEnergy2 > 0.) continue;
653           if(itsHit->GetLayer() != itsHit2->GetLayer()) continue;
654           if(itsHit->GetLadder() != itsHit2->GetLadder()) continue;
655           if(itsHit->GetDetector() != itsHit2->GetDetector()) continue;
656           if(ifl2 == 0) continue;
657           
658           Float_t px1; Float_t py1; Float_t pz1;
659           itsHit2->GetMomentumG(px1,py1,pz1);
660           Float_t ptot1 = TMath::Sqrt(px1*px1+py1*py1+pz1*pz1);
661           if(ptot1 < 0.05) continue;
662           
663           Float_t And;
664           Float_t Tim;
665           Float_t x = itsHit->GetXL();
666           Float_t z = itsHit->GetZL();
667           seg1->GetPadTxz(x,z);
668           And = z;
669           Tim = x*tstep;
670           if(And < 0 || And > maxand) continue;
671           Float_t And2;
672           Float_t Tim2;
673           Float_t x2 = itsHit2->GetXL();
674           Float_t z2 = itsHit2->GetZL();
675           seg1->GetPadTxz(x2,z2);
676           And2 = z2;
677           Tim2 = x2*tstep;
678           if(And2 < 0 || And2 > maxand) continue;
679           Float_t da = apitch*(And-And2);
680           Float_t dt = vdrift*(Tim-Tim2);
681           Float_t danh = And-And2;
682           Float_t dtmh = Tim-Tim2;
683           Float_t dist = TMath::Sqrt(da*da+dt*dt);
684           if(dt < 1000.) {
685             Float_t wx = dt*clock/(1000.*vdrift);
686             Float_t wy = da/apitch;
687             dtah->Fill(wx,wy);
688           }
689           if(dist<20.) { cout << "skip hit " << hit1 << endl; flag[hit] = 1; }
690          
691           if(ifl == 1 && ifl2 == 1) {
692             if(dist>10)dtrh->Fill(dist,wh);
693             dtrhall->Fill(dist/1000.,wh);
694             dadth->Fill(dtmh,danh);
695           }
696         }  // end cluster loop
697         
698         Float_t path = TMath::Abs(Tim)*vdrift/1000.;
699         if(flag[hit] == 0) { nGoodHitsOK++; hit_vs_time1->Fill(path); }
700       } // end hit loop
701       delete [] flag;
702       printf("nGoodHits: %d",nGoodHits);
703       printf(", nGoodHitsOK: %d\n",nGoodHitsOK);
704
705       //printf("cluster loop\n");
706  
707       AliITSRawClusterSDD *itsCluFake = 0;
708       Int_t nGoodPntsOK = 0;
709       for( int ip=0; ip<nrecp; ip++) {
710         itsClu  = (AliITSRawClusterSDD*)ITSclu->At(ip);
711         if(!itsClu) continue;
712         Float_t Ptim = (Float_t) itsClu->T();
713         Float_t dpath = Ptim*vdrift/1000.;
714         Float_t adpath = TMath::Abs(dpath);
715         if( flagP[ip] == 1) { nGoodPntsOK++; rec_vs_time1->Fill(dpath); }
716         else {
717           cout << "ip: " << ip << endl;
718           itsCluFake  = (AliITSRawClusterSDD*)ITSclu->At( ip );
719           if(!itsCluFake) continue;
720           cout << "pointer: " << itsCluFake << endl;
721           itsCluFake->PrintInfo();
722           Float_t Ptim = (Float_t) itsCluFake->T();
723           Float_t dpath = Ptim*vdrift/1000.;
724           fake_vs_time->Fill(dpath);  
725         }
726       }
727       
728       Stat_t wp = 1.;
729       if(nGoodPntsOK > 1) 
730         wp /= (((Float_t) nGoodPntsOK)*((Float_t) nGoodPntsOK)-1)/2.;
731       else
732         wp = 0.;
733       
734       //printf("double cluster loop\n");
735       for (Int_t pnt=0;pnt<nrecp;pnt++) {
736         if( flagP[pnt] == 0) continue;
737         itsPnt  = (AliITSRecPoint*)ITSrec->At(pnt);
738         if(!itsPnt) continue;
739         itsClu  = (AliITSRawClusterSDD*)ITSclu->At(pnt);
740         if(!itsClu) continue;
741         AliITSRecPoint *itsPnt1;
742         AliITSRawClusterSDD *itsClu1;
743         for (Int_t pnt1=pnt+1;pnt1<nrecp;pnt1++) {
744           if( flagP[pnt1] == 0) continue;
745           itsPnt1 = (AliITSRecPoint*)ITSrec->At(pnt1);
746           if(!itsPnt1) continue;
747           itsClu1 = (AliITSRawClusterSDD*)ITSclu->At(pnt1);
748           if(!itsClu1) continue;
749           
750           Float_t dan = itsClu->A()-itsClu1->A();
751           Float_t Pwng = (Float_t) itsClu->W();
752           Float_t dt1 = itsClu->T();
753           if(Pwng == 1) dt1 *= -1.;       
754           Float_t dt2 = itsClu1->T();
755           Float_t Pwng = (Float_t) itsClu1->W();
756           if(Pwng == 1) dt2 *= -1.;       
757           Float_t dtm = dt1-dt2;
758           dadt->Fill(dtm,dan);
759           Float_t dap = apitch*(itsClu->A()-itsClu1->A());
760           Float_t dtp = vdrift*(dt1-dt2);
761           Float_t distp = TMath::Sqrt(dap*dap+dtp*dtp);
762           if(TMath::Abs(dan) < 1.) dclutim->Fill(dtp);
763           if(TMath::Abs(dan) > 10.) dfkctim->Fill(dtp); 
764           if(dtp < 1000.) {
765             Float_t wx = dtp*clock;
766             Float_t wy = dap/apitch;
767             dtap->Fill(wx,wy);
768           }
769           dtrp->Fill(distp,wp);
770           dtrpall->Fill(distp/1000.,wp);
771         } // end loop points
772       } // end loop points
773       
774       delete [] flagP;
775     } // end loop modules
776
777     ITS->ClearModules();
778     gAlice->CleanDetectors();
779     
780   } // end loop events
781
782   cout << "open output file" << endl;
783   TFile fhistos("SDD_histos_test.root","RECREATE");
784   local->Write();
785   local1->Write();
786   
787   aa->Write();
788   at->Write();
789   tt->Write();
790   at1->Write();
791   tt1->Write();
792   at2->Write();
793   tt2->Write();
794   asigma->Write();
795   tsigma->Write();
796   asigma2->Write();
797   dadt->Write();
798   dadth->Write();
799   dclutim->Write();
800   dfkctim->Write();
801   xy3->Write();
802   xz3->Write();
803   yz3->Write();
804   xy4->Write();
805   xz4->Write();
806   yz4->Write();
807   chd->Write();
808   chr->Write();
809   chh->Write();
810   chp->Write();
811   dtrp->Write();
812   dtrpall->Write();
813   dtah->Write();
814   dtap->Write();
815   dtrh->Write();
816   dtrhall->Write();
817   effh->Write();
818   effp->Write();
819   rec_vs_time->Write();
820   hit_vs_time->Write();
821   hit_vs_time1->Write();
822   rec_vs_time1->Write();
823   rec_vs_time3->Write();
824   hit_vs_time3->Write();
825   rec_vs_time4->Write();
826   hit_vs_time4->Write();
827   fake_vs_time->Write();
828   p->Write();
829   nmatch->Write();
830   anodes->Write();
831   tsampl->Write();
832   ntotal->Write();
833   clmult->Write();
834   andtsm->Write();
835   amplit->Write();
836   amplit1->Write();
837   hitpnt->Write();
838   noihist->Write();
839   occupancy3->Write();
840   occupancy4->Write();
841   map3->Write();
842   hitmap3->Write();
843   pntmap3->Write();
844   map4->Write();
845   hitmap4->Write();
846   pntmap4->Write();
847   /*
848     for(Int_t i=0;i<22;i++) {
849     for(Int_t j=0;j<8;j++) {
850     for(Int_t k=0;k<2;k++) {
851     //mappa4hit[i][j][k]->Write();
852     //mappa4pnt[i][j][k]->Write();
853     if(i<14 && j<6) {
854     //mappa3hit[i][j][k]->Write();
855     //mappa3pnt[i][j][k]->Write();
856     }
857     }
858     }
859     }
860   */
861   xz->Write();
862   and_tim->Write();
863   pand_ptim->Write();
864   fhistos.Close();
865   file->Close();
866 }
867
868
869
870
871
872