New class for ITS coordiante transformations used by AliITSgeom nearly
[u/mrichter/AliRoot.git] / ITS / AliITSpoints.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16
17 ///////////////////////////////////////////////////////////////////////////////
18 //                                                                           //
19 //  This class contains the points for the ALICE event display               //
20 //                                                                           //
21 //Begin_Html
22 /*
23 <img src="gif/AliITSpointsClass.gif">
24 */
25 //End_Html
26 //                                                                           //
27 //                                                                           //
28 ///////////////////////////////////////////////////////////////////////////////
29
30 #include <TPad.h>
31 #include <TVirtualPad.h>
32 #include <TPolyLine3D.h>
33 #include <TPaveText.h>
34 #include <TView.h>
35 #include <TMath.h>
36 #include <TH1.h>
37 #include <TF1.h>
38 #include <TGraph.h>
39 #include <TCanvas.h>
40 #include <TSpectrum.h>
41
42 #include "AliRun.h"
43 #include "AliITSpoints.h"
44 #include "AliITSdisplay.h"
45 #include "AliITSMap.h"
46 #include "AliITSRecPointSDDnew.h"
47
48
49 //const Int_t MAXNipx=1026, MAXNipy=1026;
50 const Int_t MAXNipx=752, MAXNipy=256;
51
52 static AliITSMapA2 *sMap = 0; 
53 static Int_t  sModule=0; 
54  
55 ClassImp(AliITSpoints)
56
57 //_____________________________________________________________________________
58 AliITSpoints::AliITSpoints()
59 {
60   //
61   // Default constructor
62   //
63   fModuleIndex = -1;
64   fHitIndex = -1;
65   fTrackIndex = -1;
66   fDigitIndex = -1;
67   fMarker[0] = fMarker[1] = fMarker[2]=0;
68   fMatrix = 0;
69
70   //fConnect=kFALSE;
71 }
72
73 //_____________________________________________________________________________
74 AliITSpoints::AliITSpoints(Int_t npoints)
75   :AliPoints(npoints)
76 {
77   //
78   // Standard constructor
79   //
80   fModuleIndex = -1;
81   fHitIndex = -1;
82   fTrackIndex = -1;
83   fDigitIndex = -1;
84   fMarker[0] = fMarker[1] = fMarker[2]=0;
85   fMatrix = 0;
86
87   //fConnect=kFALSE;
88
89 }
90          
91 //_____________________________________________________________________________
92 AliITSpoints::~AliITSpoints()
93 {
94   //
95   // Default destructor
96   //
97   fHitIndex = 0;
98   fTrackIndex = 0;
99   fDigitIndex = 0;
100   fModuleIndex = 0;
101   for (Int_t i=0;i<3;i++){
102       if (fMarker[i]) delete fMarker[i];
103   }
104   fMatrix = 0;
105 }
106
107 //_____________________________________________________________________________
108 void AliITSpoints::ExecuteEvent(Int_t event, Int_t px, Int_t py)
109 {
110   //
111   //*-*-*-*-*-*-*-*-*-*Execute action corresponding to one event*-*-*-*-*-*-*-*
112   //*-*                =========================================
113   //*-*
114   //*-*  This member function must be implemented to realize the action
115   //*-*  corresponding to the mouse click on the object in the window
116   //*-*
117   //*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
118
119   if (!gPad) return;
120
121   gPad->SetCursor(kCross);
122
123   TObject *select = gPad->GetSelected();
124   if(!select) {printf("no select \n"); return;}
125   if (!select->InheritsFrom("AliITSpoints")) {gPad->SetUniqueID(0); return;}
126
127    //erase old position and draw a line at current position
128    switch (event) {
129
130
131    case kButton1Down:
132         return;
133    case kButton1Motion:
134        AnodeProjection(px,py);
135        return;
136    case kButton1Up:
137        return;
138
139
140    case kButton1Double:
141        TimeProjection(px,py);
142        return;
143
144        /* why button2 does not work ?
145         it does not know kButton2Down, kButton2Motion, kButton2Up
146    case kButton2Down:
147        printf("Button2Down \n");
148        return;
149    case kButton2Motion:
150        printf("Button2Motion \n");
151        TimeProjection(px,py);
152        return;
153    case kButton2Up:
154        return;
155    */
156    }
157    
158 }
159
160 //_____________________________________________________________________________
161 void AliITSpoints::AnodeProjection(Int_t px, Int_t py) 
162 {
163   //if(!fDrawHist) return;
164
165     gPad->SetCursor(kMove);
166     gPad->GetCanvas()->FeedbackMode(kTRUE);
167
168    Float_t xmin = gPad->GetX1();
169    Float_t xmax = gPad->GetX2();
170    Float_t ymin = gPad->GetY1();
171    Float_t ymax = gPad->GetY2();
172    //printf("xmin,xmax,ymin,ymax %f %f %f %f\n",xmin,xmax,ymin,ymax);
173     Float_t x,y;
174    // x,y in the -1,1 scale
175     x = gPad->AbsPixeltoX(px);
176     y = gPad->AbsPixeltoY(py);
177
178     Float_t scale[3],center[3];
179     Int_t irep;
180     gPad->GetView()->FindScope(scale,center,irep);
181
182     AliITSdisplay* display=((AliITSdisplay*)(gAlice->Display()));
183      Int_t anode, timebin;
184     display->GetPadIxy(anode,timebin,x*scale[0],y*scale[0]);
185     
186     Int_t amin, amax, tmin, tmax;
187     display->GetPadIxy(amin,tmin,xmin*scale[0],ymin*scale[0]);
188     display->GetPadIxy(amax,tmax,xmax*scale[0],ymax*scale[0]);
189
190     if (xmin < 0 && xmax > 0) {
191        if (0-xmin > xmax) amax-=374;
192        else {amin+=374; tmin=tmax;}
193        tmax=256;
194     } 
195     Int_t tminold=tmin;
196     Int_t tmaxold=tmax;
197     tmin=TMath::Min(tminold,tmaxold);
198     tmax=TMath::Max(tminold,tmaxold);
199
200     Int_t nbiny = tmax-tmin;
201
202     //printf("amin amax tmin tmax anode, timebin %d %d %d %d %d %d\n",amin,amax,tmin,tmax,anode,timebin);
203
204     //create or set the new canvas c2
205    TVirtualPad *padsav = gPad;
206    //printf("AnodeProj: gPad %p\n",gPad);
207    TCanvas *c2 = (TCanvas*)gROOT->GetListOfCanvases()->FindObject("c2");
208    if(c2) delete c2->GetPrimitive("Anode Projection");
209    if(c2) delete c2->GetPrimitive("Time Projection");
210    else     c2=new TCanvas("c2");
211    c2->cd();
212    c2->SetFillColor(0);
213
214    char title[30];
215    sprintf(title,"anode=%d",anode);
216
217    TH1F *h1=new TH1F("Anode Projection","",nbiny,(float)tmin,(float)tmax);
218
219    h1->SetTitle(title);
220    h1->SetTitleSize(2.);
221    h1->SetTitleOffset(2.);
222    h1->SetLabelSize(0.05);
223    h1->SetYTitle(title);
224    //h1->SetStats(0);
225
226    AliITSMapA2 *sMap=display->GetMap();
227
228    for (int j=tmin;j<tmax;j++) {
229      //printf("an tbin signal %d %d %f\n",anode,j,(float)sMap->GetSignal(anode-1,j-1));
230       h1->Fill((float)j,(float)sMap->GetSignal(anode-1,j-1));
231    }
232    h1->Smooth(1);
233    h1->Fit("gaus","ql");
234    h1->Draw();
235
236    c2->Update();
237    padsav->cd();
238    //printf("AnodeProj: gPad padsav %p %p\n",gPad,padsav);
239
240 }
241
242 //_____________________________________________________________________________
243 void AliITSpoints::TimeProjection(Int_t px, Int_t py)
244 {
245   //if(!fDrawHist) return;
246   //printf("TimeProjection \n");
247
248     gPad->SetCursor(kRightSide);
249     gPad->GetCanvas()->FeedbackMode(kTRUE);
250
251    Float_t xmin = gPad->GetX1();
252    Float_t xmax = gPad->GetX2();
253    Float_t ymin = gPad->GetY1();
254    Float_t ymax = gPad->GetY2();
255    //printf("xmin,xmax,ymin,ymax %f %f %f %f\n",xmin,xmax,ymin,ymax);
256     Float_t x,y;
257    // x,y in the -1,1 scale
258     x = gPad->AbsPixeltoX(px);
259     y = gPad->AbsPixeltoY(py);
260
261     Float_t scale[3],center[3];
262     Int_t irep;
263     gPad->GetView()->FindScope(scale,center,irep);
264
265     AliITSdisplay* display=((AliITSdisplay*)(gAlice->Display()));
266      Int_t anode, timebin;
267     display->GetPadIxy(anode,timebin,x*scale[0],y*scale[2]);
268     
269     Int_t amin, amax, tmin, tmax;
270     display->GetPadIxy(amin,tmin,xmin*scale[0],ymin*scale[2]);
271     display->GetPadIxy(amax,tmax,xmax*scale[0],ymax*scale[2]);
272
273     if (xmin < 0 && xmax > 0) {
274        if (0-xmin > xmax) amax-=374;
275        else {amin+=374;tmin=tmax;}
276        tmax=256;
277     } 
278
279     Int_t aminold=amin;
280     Int_t amaxold=amax;
281     amin=TMath::Min(aminold,amaxold);
282     amax=TMath::Max(aminold,amaxold);
283
284     Int_t nbinx = amax-amin;
285
286     //printf("amin amax tmin tmax anode, timebin %d %d %d %d %d %d\n",amin,amax,tmin,tmax,anode,timebin);
287
288     //create or set the new canvas c2
289    TVirtualPad *padsav = gPad;
290    TCanvas *c2 = (TCanvas*)gROOT->GetListOfCanvases()->FindObject("c2");
291    if(c2) delete c2->GetPrimitive("Anode Projection");
292    if(c2) delete c2->GetPrimitive("Time Projection");
293    else     c2=new TCanvas("c2");
294    c2->cd();
295    c2->SetFillColor(0);
296
297    char title[30];
298    //printf("title: anode, timebin %d %d\n",anode,timebin);
299    sprintf(title,"time sample=%d",timebin);
300
301    TH1F *h2=new TH1F("Time Projection","",nbinx,(float)amin,(float)amax);
302
303    h2->SetTitle(title);
304    h2->SetTitleSize(2.);
305    h2->SetTitleOffset(2.);
306    h2->SetLabelSize(0.05);
307    h2->SetYTitle(title);
308    //h2->SetStats(0);
309
310    AliITSMapA2 *sMap=display->GetMap();
311
312    for (int i=amin;i<amax;i++) {
313      //printf("an tbin signal %d %d %f\n",i,timebin,(float)sMap->GetSignal(i-1,timebin-1));
314       h2->Fill((float)i,(float)sMap->GetSignal(i-1,timebin-1));
315    }
316    h2->Smooth(1);
317    h2->Fit("gaus","ql");
318    h2->Draw();
319
320    c2->Update();
321    padsav->cd();
322
323 }
324 //_____________________________________________________________________________
325 void AliITSpoints::DisplayModule()
326 {
327   /*
328  
329    if (!((AliITSDisplay*)(gAlice->Display()))->GetDrawTracksOpt()) return;
330    // ????
331    ((AliITSDisplay*)(gAlice->Display()))->SetModuleNumber(fModuleIndex);
332    Int_t event=((AliITSDisplay*)(gAlice->Display()))->GetEvent();
333    ((AliITSDisplay*)(gAlice->Display()))->SetEvent(event);
334    ((AliITSDisplay*)(gAlice->Display()))->SetRange(4,4);
335    //((AliITSDisplay*)(gAlice->Display()))->DrawModule(fModuleIndex);
336   */
337 }
338
339 //_____________________________________________________________________________
340 const Text_t *AliITSpoints::GetName() const
341 {
342   //
343   // Return name of the Geant3 particle corresponding to this point
344   //
345   TParticle *particle = GetParticle();
346   if (particle) return particle->GetName();
347   else  return IsA()->GetName();
348   //if (!particle) return "Particle";
349 }
350
351 //_____________________________________________________________________________
352 Text_t *AliITSpoints::GetObjectInfo(Int_t px, Int_t py)
353 {
354   //
355   //   Redefines TObject::GetObjectInfo.
356   //   Displays the info (particle,etc
357   //   corresponding to cursor position px,py
358   //
359    if (!gPad) return (char*)"";
360   static char info[64];
361   char an[6], tbin[9], track[6];
362   an="Anode"; tbin="Time bin"; track="Track"; 
363   if(strcmp(GetName(),"AliITSpoints ")) {
364         AliITSdigitSDD *dig=GetDigit();
365         if (!dig) {sprintf(info,"%s %s %d",GetName(),track,fIndex); return info;} 
366         int anode=dig->fCellX;
367         int time=dig->fCellY;
368         sprintf(info,"%s %d %s %d %s %d",an,anode+1,tbin,time+1,track,fIndex);
369   } else {
370         if(strcmp(GetName(),"TView ")) {
371             Float_t x = gPad->AbsPixeltoX(px);
372             Float_t y = gPad->AbsPixeltoY(py);
373             sprintf(info,"%s x=%.3g, y=%.3g",GetName(),gPad->PadtoX(x),gPad->PadtoY(y));
374         } else {
375             sprintf(info,"%s %s %d",GetName(),track,fIndex); 
376         }
377
378   }
379
380   return info;
381 }
382
383 //_____________________________________________________________________________
384 void AliITSpoints::DumpHit()
385 {
386   //
387   //   Dump hit corresponding to this point
388   //
389   AliITShit *hit = GetHit();
390   if (hit) hit->Dump();
391
392 }
393
394 //_____________________________________________________________________________
395 void AliITSpoints::DumpDigit()
396 {
397   //
398   //   Dump digit corresponding to this point
399   //
400   AliITSdigitSDD *digit = GetDigit();
401   if (digit) digit->Dump();
402 }
403
404 //_____________________________________________________________________________
405 void AliITSpoints::InspectHit()
406 {
407   //
408   //   Inspect hit corresponding to this point
409   //
410
411   if (fHitIndex < 0 ) return;
412   TVirtualPad *padsav = gPad;
413   AliITShit *hit = GetHit();
414   if (hit) hit->Inspect();
415   TVirtualPad *padinspect = (TVirtualPad*)(gROOT->GetListOfCanvases())->FindObject("inspect");
416    padinspect->cd();
417    Float_t xmin = gPad->GetX1();
418    Float_t xmax = gPad->GetX2();
419    Float_t ymin = gPad->GetY1();
420    Float_t ymax = gPad->GetY2();
421    Float_t dy   = ymax-ymin;
422
423       TPaveText *pad = new TPaveText(xmin, ymin+0.1*dy, xmax, ymin+0.15*dy);
424       pad->SetBit(kCanDelete);
425       pad->SetFillColor(42);
426       pad->Draw();
427       char ptitle[100];
428       sprintf(ptitle," %s , fTrack: %d  fTrackIndex: %d ",GetName(),fIndex,fTrackIndex);
429       pad->AddText(ptitle);
430       padinspect->cd();
431       padinspect->Update();
432   if (padsav) padsav->cd();
433
434 }
435
436 //_____________________________________________________________________________
437 void AliITSpoints::InspectDigit()
438 {
439   //
440   //   Inspect digit corresponding to this point
441   //
442
443   if (fDigitIndex < 0) return;
444   TVirtualPad *padsav = gPad;
445   AliITSdigitSDD *digit = GetDigit();
446   if (digit) digit->Inspect();
447   TVirtualPad *padinspect = (TVirtualPad*)(gROOT->GetListOfCanvases())->FindObject("inspect");
448    padinspect->cd();
449    Float_t xmin = gPad->GetX1();
450    Float_t xmax = gPad->GetX2();
451    Float_t ymin = gPad->GetY1();
452    Float_t ymax = gPad->GetY2();
453    Float_t dy   = ymax-ymin;
454
455       TPaveText *pad = new TPaveText(xmin, ymin+0.1*dy, xmax, ymin+0.25*dy);
456       pad->SetBit(kCanDelete);
457       pad->SetFillColor(42);
458       pad->Draw();
459       char ptitle[11][100];
460       //      sprintf(ptitle[11],"Tracks making this digit");
461       //      pad->AddText(ptitle[11]);
462   for (int i=0;i<3;i++) {
463       if (digit->fTracks[i] == 0) continue;  
464       sprintf(ptitle[i],"fTrackIndex: %d  Charge: %f",digit->fTracks[i],digit->fTcharges[i]);
465       pad->AddText(ptitle[i]);
466   }
467       padinspect->cd();
468       padinspect->Update();
469   if (padsav) padsav->cd();
470     
471 }
472
473 //_____________________________________________________________________________
474 Int_t AliITSpoints::GetTrackIndex()
475 {
476   //
477   //   Dump digit corresponding to this point
478   //
479
480   this->Inspect();
481   /*
482   if (fDigitIndex != 0) {
483     Int_t ncol=this->fMatrix->GetNcols();
484     for (int i=0;i<ncol;i++) {
485         printf(" track charge %f %f \n",(*(this->fMatrix))(0,i),(*(this->fMatrix))(1,i));
486     }
487   }
488   */
489   return fTrackIndex;
490 }
491
492 //_____________________________________________________________________________
493 AliITShit *AliITSpoints::GetHit() const
494 {
495   //
496   //   Returns pointer to hit index in AliRun::fParticles
497   //
498   AliITS *ITS  = (AliITS*)gAlice->GetModule("ITS");
499   gAlice->TreeH()->GetEvent(fTrackIndex);
500   TClonesArray *ITShits  = ITS->Hits();
501   Int_t nhits = ITShits->GetEntriesFast();
502   if (fHitIndex < 0 || fHitIndex >= nhits) return 0;
503   return (AliITShit*)ITShits->UncheckedAt(fHitIndex);
504 }
505
506 //_____________________________________________________________________________
507 AliITSdigitSDD *AliITSpoints::GetDigit() const
508 {
509   //
510   //   Returns pointer to digit index in AliRun::fParticles
511   //
512
513   AliITSdisplay *display=(AliITSdisplay*)gAlice->Display();
514   Int_t module=display->GetModule();
515   Int_t layer=display->GetLayer();
516   Int_t id=(Int_t)((layer-1)/2);
517    
518   AliITS *ITS  = (AliITS*)gAlice->GetModule("ITS");
519   TClonesArray *ITSdigits  = ITS->DigitsAddress(id);
520   gAlice->TreeD()->GetEvent(module);
521   Int_t ndigits = ITSdigits->GetEntriesFast();
522   if (fDigitIndex < 0 || fDigitIndex >= ndigits) return 0;
523   return (AliITSdigitSDD*)ITSdigits->UncheckedAt(fDigitIndex);
524
525   /*
526   // have smth like ??
527   AliITSgeom *geom  = ITS->GetITSgeom;
528   Int_t lastSPD=geom->GetLastSPD();
529   Int_t lastSDD=geom->GetLastSDD();
530
531   Int_t nent=(Int_t)gAlice->TreeD()->GetEntries();
532   if (nent < lastSDD) {
533        // take the appropriate action ...using lastSPD, lastSDD ...
534        // or better introduce an Option_t in the display : All, SPD, SDD, SSD
535        // or combination
536        printf("Digitisation wasn't done for all det types !\n");
537   }
538   */
539   // have smth like 
540   // if (module < lastSPD) return (AliITSdigitSPD*)ITSdigits->UncheckedAt(fDigitIndex)
541   // else  (module < lastSDD) return (AliITSdigitSDD*)ITSdigits->UncheckedAt(fDigitIndex)  .... ??
542
543   //return 0;
544 }
545
546 //_____________________________________________________________________________
547
548
549 struct Bin {
550    AliITSdigitSDD *dig;
551    int idx;
552    int digidx;
553    Bin() {dig=0; idx=-1; digidx=-1;}
554 };
555
556
557 struct PreCluster : public AliITSRecPointSDDnew {
558    AliITSdigitSDD* summit;
559    int idx;   
560    int npeaks;
561    int npoly;
562    float xpoly[300];
563    float ypoly[300];
564    float zpoly[300];
565    PreCluster() : AliITSRecPointSDDnew() {npeaks=npoly=0; 
566                                        for (int k=0;k<300;k++) {
567                                          xpoly[k]=ypoly[k]=zpoly[k]=0;
568                                        }
569    }
570                               
571 };
572
573
574 //_____________________________________________________________________________
575
576 //static void FindCluster(int i, int j, Bin bins[][MAXNipy], PreCluster &c) 
577 static void FindCluster(int i, int j, Bin *bins, PreCluster &c, int thresh) 
578
579 {
580
581   //
582   // Find clusters
583   //
584
585   //Bin& b=bins[i][j];
586   Bin& b=bins[i*MAXNipy+j];
587   Int_t q=b.dig->fSignal - thresh;
588
589   printf("i j q %d %d %d\n",i,j,q);
590
591   //  if (b.idx >= 0 && b.idx != c.idx) {
592   if (b.idx >= 0 && b.idx > c.idx) {
593     c.idx=b.idx;
594     c.npeaks++;
595   }
596   
597   if (q > TMath::Abs(c.summit->fSignal)) c.summit=b.dig;
598
599   // get pad coordinates and prepare the up and down steps   
600
601   Float_t xl[3], dx[3];
602   AliITSdisplay *display=(AliITSdisplay*)gAlice->Display();
603   display->GetPadCxy(i,j,xl,dx);
604
605   // calculate center of gravity
606   c.npoly++;
607   c.fMulDigits++;
608   c.fDigitsList[c.npoly-1]=b.digidx;
609
610   //printf("npoly c.fDigitsList[c.npoly-1] %d %d \n",c.npoly,c.fDigitsList[c.npoly-1]);
611
612   if (c.npoly > 300 ) {
613     printf("FindCluster - npoly >300,  npoly %d \n",c.npoly);
614     c.npoly=300;
615   }
616   c.xpoly[c.npoly-1]=xl[0];
617   c.ypoly[c.npoly-1]=xl[1];
618   c.zpoly[c.npoly-1]=xl[2];
619
620   c.fX += q*xl[0];
621   c.fZ += q*xl[2];
622   c.fY =  xl[1];
623   c.fQ += q;
624   
625   b.dig = 0;  b.idx = c.idx;
626
627   /*
628   // left and right  
629   if (i && bins[i-1][j].dig) FindCluster(i-1,j,bins,c);
630   if (i < MAXNipx && bins[i+1][j].dig) FindCluster(i+1,j,bins,c);
631   // up and down
632   if (j+1 < MAXNipy && bins[i][j+1].dig) FindCluster(i,j+1,bins,c);
633   if (j && bins[i][j-1].dig) FindCluster(i,j-1,bins,c);
634   */
635   // left and right  
636   if (i && bins[(i-1)*MAXNipy+j].dig) FindCluster(i-1,j,bins,c,thresh);
637   if (i < MAXNipx && bins[(i+1)*MAXNipy+j].dig) FindCluster(i+1,j,bins,c,thresh);
638   // up and down
639   if (j+1 < MAXNipy && bins[i*MAXNipy+(j+1)].dig) FindCluster(i,j+1,bins,c,thresh);
640   if (j && bins[i*MAXNipy+(j-1)].dig) FindCluster(i,j-1,bins,c,thresh);
641
642 }
643
644 //_____________________________________________________________________________
645 void AliITSpoints::GetCenterOfGravity()
646 {
647   //
648   // simple ITS cluster finder from digits -- finds neighbours and 
649   // calculates center of gravity for the cluster
650   //
651
652   const int THRESHOLD=20;
653
654   Bin *sBin=new Bin[MAXNipx*MAXNipy];
655   //else memset(sBin,0,sizeof(Bin)*MAXNipx*MAXNipy);
656   //printf("sBin %p\n",sBin);
657
658
659   //Bin bins[MAXNipx][MAXNipy]; 
660
661   AliITSdisplay *display=(AliITSdisplay*)gAlice->Display();
662   Int_t module=display->GetModule(); // or module=fModuleIndex;
663   Int_t layer=display->GetLayer();
664   Int_t id=(Int_t)((layer-1)/2);
665    
666   AliITS *ITS  = (AliITS*)gAlice->GetModule("ITS");
667
668   TClonesArray *ITSdigits  = ITS->DigitsAddress(id);
669   Int_t nent=(Int_t)gAlice->TreeD()->GetEntries();
670   //gAlice->TreeD()->GetEvent(nent-nofmodules+module-1);
671   gAlice->TreeD()->GetEvent(module);
672   Int_t ndigits = ITSdigits->GetEntriesFast();
673   if (fDigitIndex < 0 || fDigitIndex >= ndigits) return;
674
675   AliITSdigitSDD  *dig;
676   dig=(AliITSdigitSDD*)ITSdigits->UncheckedAt(fDigitIndex);
677   Int_t ipx=dig->fCellX;
678   Int_t ipy=dig->fCellY;
679   //bins[ipx][ipy].dig=dig;
680   sBin[ipx*MAXNipy+ipy].dig=dig;
681   //printf("ipx, ipy, Sbin.dig %d %d %p\n",ipx,ipy,sBin[ipx*MAXNipy+ipy].dig);
682     
683   int ndig;
684   int ncls=0;
685   for (ndig=0;ndig<ndigits;ndig++) {
686     dig = (AliITSdigitSDD*)ITSdigits->UncheckedAt(ndig);
687     int i=dig->fCellX; int j=dig->fCellY;
688     //printf("ndig i j %d %d %d \n",ndig,i,j);
689     //bins[i][j].dig=dig;
690     sBin[i*MAXNipy+j].dig=dig;
691     sBin[i*MAXNipy+j].digidx=ndig;
692   }
693   //PreCluster c; c.summit=bins[ipx][ipy].dig; c.idx=ncls;
694   // FindCluster(ipx, ipy, bins, c);
695    PreCluster c; c.summit=sBin[ipx*MAXNipy+ipy].dig; c.idx=ncls;
696    FindCluster(ipx, ipy, sBin, c, 0);
697     
698    if (c.npeaks>1) {
699       printf("GetCenterOfGravity -- more than one peak");
700    }
701    c.fX /= c.fQ;
702    c.fZ /= c.fQ;
703    ncls++;
704
705    AliITSpoints *points = 0;
706    points = new AliITSpoints(1);
707    points->SetMarkerColor(kYellow);
708    points->SetMarkerStyle(5);
709    points->SetMarkerSize(1.);
710    points->SetPoint(0,c.fX,c.fY,c.fZ);
711    points->SetParticle(-1);
712    points->Draw();
713
714    TPolyLine3D *pline=0; 
715
716    /*
717    pline=new TPolyLine3D(c.npoly);
718    Int_t np=c.npoly;
719    TVector *xp=new TVector(c.npoly);
720    TVector *yp=new TVector(c.npoly);
721    TVector *zp=new TVector(c.npoly);
722    for (int i=0;i<np;i++) {
723      (*xp)(i)=c.xpoly[i];
724      (*yp)(i)=c.ypoly[i];
725      (*zp)(i)=c.zpoly[i];
726      pline->SetPoint(i,(*xp)(i),(*yp)(i),(*zp)(i));
727      //printf("np, i, xp, yp, zp %d %d %f %f %f \n",np,i,(*xp)(i),(*yp)(i),(*zp)(i));
728    }
729    */
730    /*
731    delete xp;
732    delete yp;
733    delete zp;
734    */ 
735
736    //if (fConnect) {
737       pline=new TPolyLine3D(c.npoly,c.xpoly,c.ypoly,c.zpoly);
738       pline->SetLineColor(kWhite);
739       pline->Draw();
740       //}
741
742
743
744    for (int k=0;k<c.npoly;k++) {
745      c.xpoly[k]=c.ypoly[k]=c.zpoly[k]=0;
746    }
747    c.npoly=0;
748
749 }
750
751 //_____________________________________________________________________________
752 void AliITSpoints::FindLocalMaxima()
753 {
754   //
755
756   Bin *lBin=new Bin[MAXNipx*MAXNipy];
757
758   AliITSdisplay *display=(AliITSdisplay*)gAlice->Display();
759   Int_t module=display->GetModule(); // or module=fModuleIndex;
760   Int_t layer=display->GetLayer();
761   Int_t id=(Int_t)((layer-1)/2);
762    
763   AliITS *ITS  = (AliITS*)gAlice->GetModule("ITS");
764
765   TClonesArray *ITSdigits  = ITS->DigitsAddress(id);
766   Int_t nent=(Int_t)gAlice->TreeD()->GetEntries();
767   //gAlice->TreeD()->GetEvent(nent-nofmodules+module-1);
768   gAlice->TreeD()->GetEvent(module);
769   Int_t ndigits = ITSdigits->GetEntriesFast();
770   if (fDigitIndex < 0 || fDigitIndex >= ndigits) return;
771
772   AliITSdigitSDD  *dig;
773   dig=(AliITSdigitSDD*)ITSdigits->UncheckedAt(fDigitIndex);
774   Int_t ipx=dig->fCellX;
775   Int_t ipy=dig->fCellY;
776   lBin[ipx*MAXNipy+ipy].dig=dig;
777     
778   int ndig;
779   int ncls=0;
780
781   const int kMaxNeighbours=4;
782   
783   for (ndig=0;ndig<ndigits;ndig++) {
784     dig = (AliITSdigitSDD*)ITSdigits->UncheckedAt(ndig);
785     int i=dig->fCellX; int j=dig->fCellY;
786     //printf("ndig i j %d %d %d \n",ndig,i,j);
787     lBin[i*MAXNipy+j].dig=dig;
788     lBin[i*MAXNipy+j].digidx=ndig;
789   }
790    PreCluster c; c.summit=lBin[ipx*MAXNipy+ipy].dig; c.idx=ncls;
791    FindCluster(ipx, ipy, lBin, c, 0);
792     
793    printf("finished FindCluster !!!\n");
794
795    Int_t fMul=c.npoly;
796    printf("fMul %d \n",fMul);
797    AliITSdigitSDD* digt;
798
799   int fIx[fMul];
800   int fIy[fMul];
801   int fSortIx[fMul];
802   int fSortIy[fMul];
803   int fQ[fMul];
804   int fQini[fMul];
805   int fIndLocal[fMul];
806   int fIndSaddle[fMul];
807   AliITSdigitSDD*    fDig[fMul];        // current list of digits 
808
809   for (int i=0; i<fMul; i++)
810   {
811         int idx=c.fDigitsList[i];
812         digt = (AliITSdigitSDD*)ITSdigits->UncheckedAt(idx);
813         fDig[i]=digt;
814         fIx[i]= digt->fCellX;
815         fIy[i]= digt->fCellY;
816         fSortIx[i]= digt->fCellX;
817         fSortIy[i]= digt->fCellY;
818         fQini[i] = digt->fSignal;
819   }
820
821     Int_t idx_sort[fMul];
822
823     TMath::Sort(fMul,fSortIy,idx_sort,0);
824     Int_t tmin=fIy[idx_sort[0]];
825     Int_t tmax=fIy[idx_sort[fMul-1]];
826
827     TMath::Sort(fMul,fSortIx,idx_sort,0);
828     Int_t amin=fIx[idx_sort[0]];
829     Int_t amax=fIx[idx_sort[fMul-1]];
830     //printf("tmin tmax amin amax %d %d %d %d\n",tmin,tmax,amin,amax);
831
832
833    AliITSMapA2 *sMap=display->GetMap();
834
835    Int_t nt=tmax-tmin;
836    Int_t nbinx=amax-amin;
837
838    float **source = new float* [nbinx+1];
839    for(int i=0;i<nbinx+1;i++) source[i] = new float[nt+1];
840
841    Double_t Qsmooth[nbinx+1][nt+1];
842    Double_t temp[nt+1];
843    TH1F *h1=new TH1F("Anode Projection","",nbinx+1,(float)amin,(float)amax+1);
844
845    for (int i=amin;i<=amax;i++) {
846      for (int j=tmin;j<=tmax;j++) {
847          source[i-amin][j-tmin]=0;
848          Qsmooth[i-amin][j-tmin]=0.;
849          source[i-amin][j-tmin]=(float)sMap->GetSignal(i,j);
850          Qsmooth[i-amin][j-tmin]=sMap->GetSignal(i,j);
851          //printf("i j Qsmooth signal map %d %d %f %f\n",i,j,Qsmooth[i-amin][j-tmin],sMap->GetSignal(i,j));
852          temp[j-tmin]=sMap->GetSignal(i,j);
853      }
854      h1->SmoothArray(nt+1,temp,1);
855      for (int j=tmin;j<=tmax;j++) {
856          Qsmooth[i-amin][j-tmin]=temp[j-tmin];
857          //printf("i j Qsmooth temp %d %d %f %f\n",i,j,Qsmooth[i-amin][j-tmin],temp[j-tmin]);
858      }
859
860    }
861
862
863 //
864 //  Find local maxima//
865     int fNLocal=0;
866     int fNSaddle=0;
867     Bool_t IsLocal[fMul];
868     Bool_t IsSaddle[fMul];
869     Int_t AssocPeak[fMul];
870     Int_t nn;
871     Int_t Qneighb[2*kMaxNeighbours];
872     Int_t X[2*kMaxNeighbours], Y[2*kMaxNeighbours];
873     for (int i=0; i<fMul; i++) {
874         fQ[i]=(Int_t)Qsmooth[fIx[i]-amin][fIy[i]-tmin];
875         Neighbours(fIx[i], fIy[i], &nn, X, Y);
876         IsLocal[i]=kTRUE;
877         IsSaddle[i]=kTRUE;
878         for (int j=0; j<nn; j++) {
879             Qneighb[j]=0;
880             if (X[j]<amin || X[j]>amax || Y[j]<tmin || Y[j]>tmax) continue;
881             Qneighb[j]=(Int_t)Qsmooth[X[j]-amin][Y[j]-tmin];
882         }
883         int nlzero=0;
884         for (int j=0; j<nn; j++) {
885           //if (!fQ[i] || !Qneighb[j] || !fQini[i]) {nlzero++;continue;}
886             if (!fQ[i] || Qneighb[j]<4 || !fQini[i]) {nlzero++;continue;}
887             //if (Qneighb[j] > fQ[i] || nlzero>=2) { //???
888             if (Qneighb[j] > fQ[i] ) {
889               //printf("i j %d %d  local kFALSE Qneighb fQ %d %d\n",i,j,Qneighb[j],fQ[i]);
890                 IsLocal[i]=kFALSE;
891                 break;
892 //
893 // handle special case of neighbouring pads with equal signal
894              } else if (Qneighb[j] == fQ[i]) {
895                 if (fNLocal >0) {
896                     for (Int_t k=0; k<fNLocal; k++) {
897                         if (X[j]==fIx[fIndLocal[k]] && (Y[j]==fIy[fIndLocal[k]] || Y[j]==fIy[fIndLocal[k]+1] || Y[j]==fIy[fIndLocal[k]-1]))
898                         {
899                             IsLocal[i]=kFALSE;
900                         } 
901                         if(j-2 >= tmin && j+2 <= tmax && (fQ[i] == Qsmooth[X[j]-amin][Y[j-2]-tmin]) || (fQ[i] == Qsmooth[X[j]-amin][Y[j+2]-tmin]))
902                         {
903                             IsLocal[i]=kFALSE;
904                         } 
905                     } // loop over local maxima
906                 } // are there are already local maxima
907             } // IsLocal
908         } // loop over neighb
909
910         // find saddle points
911         int nzero=0; int nlower=0;
912         for (int j=0; j<nn; j++) {
913             if (X[j]<amin || X[j]>amax || Y[j]<tmin || Y[j]>tmax) continue;
914             Qneighb[j]=(int)source[X[j]-amin][Y[j]-tmin];
915             if (!Qneighb[j] || !fQ[i] || !fQini[i]) {nzero++;continue;}
916             //if (!fQ[i] || Qneighb[j]<4 || !fQini[i]) {nzero++;continue;}
917             if (Qneighb[j] < fQ[i] ) {
918               //printf("fIx fIy X[j] Y[j] %d %d %d %d  saddle kFALSE Qneighb fQ %d %d\n",fIx[i],fIy[i],X[j],Y[j],Qneighb[j],fQ[i]);
919               nlower++;
920               
921               IsSaddle[i]=kFALSE;
922               break;
923 //
924 // handle special case of neighbouring pads with equal signal
925               // this is dangerous for saddle points ! better out !!
926             } else if (Qneighb[j] == fQ[i]) {
927                 if (fNSaddle >0) {
928                     for (Int_t k=0; k<fNSaddle; k++) {
929                         if ((X[j]==fIx[fIndSaddle[k]] && (Y[j]==fIy[fIndSaddle[k]] || Y[j]==fIy[fIndSaddle[k]-1] || Y[j]==fIy[fIndSaddle[k]+1]))) 
930
931                           //|| (X[j]==fIx[fIndSaddle[k]-1] && (Y[j]==fIy[fIndSaddle[k]] || Y[j]==fIy[fIndSaddle[k]-1] ||Y[j]==fIy[fIndSaddle[k]+1])) || (X[j]==fIx[fIndSaddle[k]+1] && (Y[j]==fIy[fIndSaddle[k]] || Y[j]==fIy[fIndSaddle[k]-1] ||Y[j]==fIy[fIndSaddle[k]+1])))
932                         {
933                             IsSaddle[i]=kFALSE;
934                         } 
935
936                     } // loop over saddle points
937                 } // are there are already saddle points
938             } // IsSaddle
939         } // loop over next neighbours
940
941         // Maxima should not be on the edge - therefore cond on nlzero ?
942         if (IsLocal[i] && fQ[i]>0) {
943           //printf(" i fNLocal nlzero %d %d %d\n",i,fNLocal,nlzero);
944             fIndLocal[fNLocal]=i;
945             fNLocal++;
946         } else fQ[i]=0;
947         // Maxima should not be on the edge
948         //if (IsSaddle[i] && fQ[i]>0 && !nzero ) {
949         if (IsSaddle[i] && fQini[i]>0) {
950           //printf(" i fNSaddle nzero %d %d %d\n",i,fNSaddle,nzero);
951           if (fIy[i]-1 < tmin || fIy[i]+1 > tmax) continue;
952           if (source[fIx[i]-amin][fIy[i]-1-tmin] && source[fIx[i]-amin][fIy[i]+1-tmin]) {  
953             fIndSaddle[fNSaddle]=i;
954             fNSaddle++;
955           }
956         }
957     } // loop over all digits
958     printf("fNLocal %d  fNSaddle %d fMul %d \n",fNLocal,fNSaddle,fMul);
959
960
961     // take the highest local maxima and spread the charge, i.e. define the
962     // response 
963
964     Int_t idx=TMath::LocMax(fMul,fQ);
965     int thresh=0;
966     int qsaddle[fNSaddle];
967     if (fNSaddle) {
968        for (int i=0;i<fNSaddle;i++) {
969                qsaddle[i]=fQini[fIndSaddle[i]];
970                printf("isaddle qsaddle %d %d \n",i,qsaddle[i]);
971        }
972        Int_t idmin=TMath::LocMin(fNSaddle,qsaddle);
973        thresh=qsaddle[idmin];
974     }
975     printf("idx fQ[idx] thresh %d %d %d\n",idx,fQ[idx],thresh);
976
977     // this part should be recursive for fNLocal>2 and should apply only to
978     // clusters with fNLocal==2
979
980     Float_t xc[fNLocal], yc[fNLocal], zc[fNLocal];
981
982     if (fNLocal > 1 && fNSaddle) {
983       for (int i=0; i<fMul; i++) {
984           printf("before i %d\n",i);
985           if (source[fIx[i]-amin][fIy[i]-tmin]  > thresh+3) lBin[fIx[i]*MAXNipy+fIy[i]].dig=fDig[i];
986           else lBin[fIx[i]*MAXNipy+fIy[i]].dig=0;
987           printf("fDig %p \n",fDig[i]);
988
989           printf("i j source %d %d %f lBin %p \n",fIx[i],fIy[i],source[fIx[i]-amin][fIy[i]-tmin],lBin[fIx[i]*MAXNipy+fIy[i]].dig);
990       }
991
992       for (int i=0; i<fNLocal; i++) {
993         int ipxl=fIx[fIndLocal[i]]; int ipyl=fIy[fIndLocal[i]];
994         printf("ipxl ipyl %d %d\n",ipxl,ipyl);
995         //printf("q %p %d \n",fDig[fIndLocal[i]],fDig[fIndLocal[i]]->fSignal);
996         printf("%p  \n",lBin[ipxl*MAXNipy+ipyl].dig);
997
998         PreCluster cnew; cnew.summit=lBin[ipxl*MAXNipy+ipyl].dig; cnew.idx=ncls;
999         printf("q %p %d \n",fDig[fIndLocal[i]],lBin[ipx*MAXNipy+ipy].dig->fSignal);
1000         //FindCluster(ipx, ipy, lBin, cnew, 2*thresh+1); // it's further away
1001         FindCluster(ipxl, ipyl, lBin, cnew, 0);
1002         printf("finish FindCluster i %d\n",i);
1003         cnew.fX /= cnew.fQ;
1004         cnew.fZ /= cnew.fQ;
1005         ncls++;
1006         printf("i cnew.fX cnew.fY cnew.fZ %d %f %f %f\n",i,cnew.fX, cnew.fY, cnew.fZ);
1007
1008         xc[i]=cnew.fX;
1009         yc[i]=cnew.fY;
1010         zc[i]=cnew.fZ;
1011
1012         
1013         Int_t fMulnew=cnew.npoly;
1014         printf("fMulnew %d \n",fMulnew);
1015
1016         for (int j=0; j<fMulnew; j++)
1017           {
1018             int idx=cnew.fDigitsList[j];
1019             //printf("fMul idx %d  %d \n",i,idx);
1020             digt = (AliITSdigitSDD*)ITSdigits->UncheckedAt(idx);
1021             fDig[j]=digt;
1022             //printf("digt %p \n",digt);
1023             fIx[j]= digt->fCellX;
1024             fIy[j]= digt->fCellY;
1025             fSortIx[j]= digt->fCellX;
1026             fSortIy[j]= digt->fCellY;
1027             fQini[j] = digt->fSignal;
1028             // printf("i fIx fIy fQini %d %d %d %d\n",j,fIx[j],fIy[j],fQini[j]);
1029           }
1030
1031         Int_t idx_sort[fMulnew];
1032
1033         TMath::Sort(fMulnew,fSortIy,idx_sort,0);
1034         Int_t tmin=fIy[idx_sort[0]];
1035         Int_t tmax=fIy[idx_sort[fMulnew-1]];
1036     
1037         if (ipy-tmin != tmax-ipy) {
1038              printf("non-symetric cluster in time direction\n");
1039              // do something - take Local maxima coord
1040         } else {
1041             // take cog coord
1042         }
1043         printf("ipy-tmin tmax-ipy %d %d\n",ipy-tmin, tmax-ipy);
1044
1045         TMath::Sort(fMulnew,fSortIx,idx_sort,0);
1046         Int_t amin=fIx[idx_sort[0]];
1047         Int_t amax=fIx[idx_sort[fMulnew-1]];
1048
1049         if (ipxl-amin != amax-ipxl) {
1050              printf("non-symetric cluster in anode direction\n");
1051              // do something - take Local maxima coord
1052         } else {
1053             // take cog coord
1054         }
1055         printf("ipxl-amin amax-ipxl %d %d\n",ipxl-amin, amax-ipxl);
1056
1057         printf("inside loop over local maxima: i ipxl ipyl %d %d %d\n",i,ipxl,ipyl);
1058       } // loop over local maxima
1059
1060     } // if fNLocal
1061
1062
1063     if (fNLocal == 1) {
1064         int ipx1=fIx[fIndLocal[0]]; int ipy1=fIy[fIndLocal[0]];
1065         if (ipy1-tmin != tmax-ipy1) {
1066              printf("non-symetric cluster in time direction\n");
1067              // do something - take Local maxima coord or fit
1068         } else {
1069             // take cog coord or fit
1070         }
1071         printf("ipy1-tmin tmax-ipy1 %d %d\n",ipy1-tmin, tmax-ipy1);
1072
1073         if (ipx1-amin != amax-ipx1) {
1074              printf("non-symetric cluster in anode direction\n");
1075              // do something - take Local maxima coord or fit
1076         } else {
1077             // take cog coord or fit
1078         }
1079         printf("ipx1-amin amax-ipx1 %d %d\n",ipx1-amin, amax-ipx1);
1080     }
1081
1082
1083     printf("Here! fNLocal fNSaddle %d %d\n",fNLocal,fNSaddle);
1084    
1085    AliITSpoints *points = 0;
1086    for (int i=0; i<fNLocal; i++) {
1087         points = new AliITSpoints();
1088         points->SetMarkerColor(kYellow);
1089         points->SetMarkerStyle(5);
1090         points->SetMarkerSize(1.);
1091         points->SetPoint(0,xc[i],yc[i],zc[i]);
1092         points->SetParticle(-1);
1093         points->Draw();
1094    }
1095
1096
1097    //points = new AliITSpoints(fNLocal);
1098    for (int i=0; i<fNLocal; i++) {
1099        points = new AliITSpoints();
1100        int idx=fIndLocal[i];
1101        points->SetMarkerColor(kGreen);
1102        points->SetMarkerStyle(5);
1103        points->SetMarkerSize(1.5);
1104        points->SetPoint(0,c.xpoly[idx],c.ypoly[idx],c.zpoly[idx]);
1105        points->SetParticle(-1);
1106        points->Draw();
1107    }
1108    //points = new AliITSpoints(fNSaddle);
1109    for (int i=0; i<fNSaddle; i++) {
1110        points = new AliITSpoints();
1111        int idx=fIndSaddle[i];
1112        points->SetMarkerColor(kWhite);
1113        points->SetMarkerStyle(5);
1114        points->SetMarkerSize(1.5);
1115        points->SetPoint(0,c.xpoly[idx],c.ypoly[idx],c.zpoly[idx]);
1116        points->SetParticle(-1);
1117        points->Draw();
1118    }
1119
1120
1121    for (int k=0;k<c.npoly;k++) {
1122      c.xpoly[k]=c.ypoly[k]=c.zpoly[k]=0;
1123    }
1124    c.npoly=0;
1125
1126    delete [] lBin;
1127 }
1128
1129 //_____________________________________________________________________________
1130 static Double_t sinoid(Double_t *x, Double_t *par)
1131 {
1132     Double_t arg = -2*TMath::Pi()*x[0];
1133     Double_t fitval= par[0]*TMath::Sin(arg)+
1134         par[1]*TMath::Sin(2*arg)+
1135         par[2]*TMath::Sin(3*arg)+
1136         par[3]*TMath::Sin(4*arg)+
1137         par[4]*TMath::Sin(5*arg);
1138     return fitval;
1139  }
1140
1141
1142
1143 //_____________________________________________________________________________
1144 void AliITSpoints::
1145 Neighbours(Int_t iX, Int_t iY, Int_t* Nlist, Int_t Xlist[8], Int_t Ylist[8])
1146 {
1147
1148     *Nlist=8;
1149     Xlist[0]=Xlist[1]=iX;
1150     if(iX) Xlist[2]=iX-1;
1151     else Xlist[2]=iX;
1152     if (iX < MAXNipx) Xlist[3]=iX+1;
1153     else Xlist[3]=iX;
1154     if(iY) Ylist[0]=iY-1;
1155     else Ylist[0]=iY;
1156     if (iY < MAXNipy) Ylist[1]=iY+1;
1157     else Ylist[1]=iY;
1158     Ylist[2]=Ylist[3]=iY;
1159
1160     // diagonal elements
1161
1162     Xlist[4]=Xlist[5]=Xlist[2];
1163     Xlist[6]=Xlist[7]=Xlist[3];
1164
1165     Ylist[4]=Ylist[7]=Ylist[1];
1166     Ylist[5]=Ylist[6]=Ylist[0];
1167     
1168
1169 }
1170