The operator[] is replaced by At() or AddAt() in case of TObjArray.
[u/mrichter/AliRoot.git] / TPC / TPCDigits2Clusters.C
1
2 /*
3 #include "alles.h"
4 extern AliTPCParam *gTPCParam;
5 extern AliTPCClustersArray * gCalcClusters;
6 extern AliTPCClustersArray * gDifClusters;
7 AliTPCDigitsArray * GetDigitsArray();
8   AliTPCClustersArray * GetExactClustersArray();
9 AliTPCClustersArray *  GetCalcClustersArray(Bool_t newtree=kFALSE);
10 AliTPCClustersArray *  GetDifClustersArray(Bool_t newtree=kFALSE);
11 TClonesArray * CompareClusters(TClonesArray *calcclusters,TClonesArray *exactclusters, 
12                                 TClonesArray *diffclusters, Float_t shx, Float_t shy);
13 */
14
15 void   TPCDigits2Clusters(AliTPCDigitsArray * digarr, AliTPCClustersArray *calcarr,
16                           AliTPCClustersArray *exarr=0, AliTPCClustersArray *difarr=0)
17 {
18   //
19   //calculate clusters position from digits information
20   //
21
22   
23   //initialisation of cluster finder
24   AliClusterFinder cf;
25   cf.SetThreshold(gTPCParam->GetZeroSup());
26   cf.SetDetectorParam(gTPCParam);
27   cf.SetBFit(kFALSE);
28   cf.GetMinuit()->SetPrintLevel(-1);
29   cf.GetMinuit()->SetMaxIterations(20);
30   
31   if (digarr ==0) digarr = GetDigitsArray();
32   if ( (digarr == 0) || (digarr->GetTree()==0)) return;  
33   if (exarr==0) exarr =GetExactClustersArray();
34   if ((exarr!=0)&&(difarr==0)) difarr= GetDifClustersArray();
35
36   if (calcarr==0) calcarr = GetCalcClustersArray(kTRUE);
37   //loop over all writen digits segments
38   Int_t nsegment = (Int_t)digarr->GetTree()->GetEntries();
39   Int_t segment;
40   for (segment = 0; segment<nsegment;segment++){
41     //load segment with index (TTree internal)  number segment
42     AliSimDigits *digrow= (AliSimDigits*)digarr->LoadEntry(segment);
43     Int_t sector,row;
44     ((AliTPCParam*)digarr->GetParam())->AdjustSectorRow(digrow->GetID(),sector,row);
45     
46     AliTPCClustersRow * clrow = calcarr->CreateRow(sector,row);
47     AliTPCClustersRow * difrow= difarr->CreateRow(sector,row);    
48     AliTPCClustersRow * exrow =  exarr->GetRow(sector,row);
49     if (exrow==0) exrow = exarr->LoadRow(sector,row);
50     TH2F *his = (TH2F*)digrow->GenerHisto();
51     digrow->ExpandTrackBuffer();
52     if (his==0) return;
53
54     //set current index for cluster finder
55     Int_t  index[3]= {0,sector,row};
56     cf.GetHisto(his);
57     cf.SetDetectorIndex(index);    
58     cf.FindPeaks3(clrow->GetArray());
59
60     //get cluster tracks ID
61     Int_t ncl = clrow->GetArray()->GetEntriesFast();
62     for (Int_t icl = 0 ;icl<ncl; icl++)
63       {
64         AliDigitCluster *  cl =(AliDigitCluster*)clrow->GetArray()->UncheckedAt(icl);
65         Int_t i = (Int_t)cl->fMaxX;
66         Int_t j = (Int_t)cl->fMaxY;             
67         cl->fTracks[0]=digrow->GetTrackIDFast(i,j,0);
68         cl->fTracks[1]=digrow->GetTrackIDFast(i,j,1);
69         cl->fTracks[2]=digrow->GetTrackIDFast(i,j,2);   
70       }
71     Float_t shx = gTPCParam->GetZOffset()/gTPCParam->GetZWidth()+0.5;
72     Float_t shy= 0.5;
73     CompareClusters(clrow->GetArray(),exrow->GetArray(),difrow->GetArray(),shx,shy);        
74     calcarr->StoreRow(sector,row);
75     difarr->StoreRow(sector,row);
76     digarr->ClearRow(sector,row);
77   }    
78   //STORING RESULTS
79  char treeName[100];
80  sprintf(treeName,"TreeCCalc_%s",gTPCParam->GetTitle());
81  calcarr->GetTree()->Write(treeName);
82  if (difarr!=0){
83    char treeName[100];
84    sprintf(treeName,"TreeCDif_%s",gTPCParam->GetTitle());
85    difarr->GetTree()->Write(treeName);
86  }
87
88 }
89
90
91 AliTPCClustersArray *  GetCalcClustersArray(Bool_t newtree=kFALSE, const char* name=0) 
92 {
93   //
94   //construct AliTPCClusters Array object
95   //
96   AliTPCClustersArray * arr=0;
97   //  if ( (gAlice!=0) && (gAlice->GetDetector("TPC")!=0) ) {    
98   //  arr = ((AliTPC*)gAlice->GetDetector("TPC"))->GetClustersArray();
99   //  if (arr!=0) arr->Update();
100   //}
101   if (arr==0) {    
102     arr = new AliTPCClustersArray;
103     TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
104     if (file==0) file = TFile::Open("rfio:galice.root","update");   
105     arr->SetClusterType("AliDigitCluster");         
106     arr->Setup(gTPCParam);
107     cout<<"Update status : "<<arr->Update()<<"\n";
108     char treeName[100];
109     if (name==0)
110       sprintf(treeName,"TreeCExact_%s",gTPCParam->GetTitle());
111     else
112       sprintf(treeName,"TreeCCalc_%s",gTPCParam->GetTitle());    
113     if (newtree!=kTRUE){
114       cout<<"Connect tree status : "<<arr->ConnectTree(treeName)<<"\n";
115     }
116     else {
117       arr->MakeTree();
118     }
119   }
120   gCalcClusters = arr;
121   return arr;
122 }
123
124
125 AliTPCClustersArray *  GetDifClustersArray(Bool_t newtree=kFALSE, const char* name=0) 
126 {
127   //
128   //construct AliTPCClusters Array object
129   //
130   AliTPCClustersArray * arr=0;
131   //  if ( (gAlice!=0) && (gAlice->GetDetector("TPC")!=0) ) {    
132   //  arr = ((AliTPC*)gAlice->GetDetector("TPC"))->GetClustersArray();
133   //  if (arr!=0) arr->Update();
134   //}
135   if (arr==0) {    
136     arr = new AliTPCClustersArray;
137     TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
138     if (file==0) file = TFile::Open("rfio:galice.root","update");   
139     arr->SetClusterType("AliDifCluster");         
140     arr->Setup(gTPCParam);
141     cout<<"Update status : "<<arr->Update()<<"\n";
142     char treeName[100];
143     if (name==0) 
144       sprintf(treeName,"TreeCExact_%s",gTPCParam->GetTitle());
145     else
146       sprintf(treeName,"TreeCDif_%s",gTPCParam->GetTitle());
147     if (newtree!=kTRUE){
148       cout<<"Connect tree status : "<<arr->ConnectTree(treeName)<<"\n";
149     }
150     else {
151       arr->MakeTree();
152     }
153   }
154   gDifClusters  = arr; 
155   return arr;
156 }
157
158
159
160 TClonesArray * CompareClusters(TClonesArray *calcclusters,TClonesArray *exactclusters, 
161                                 TClonesArray *diffclusters, Float_t shx, Float_t shy)
162 {
163   if (calcclusters==0) return 0;
164   if (exactclusters==0) return 0;
165   if (diffclusters==0) return 0;
166   TClonesArray & diff=*diffclusters;
167    
168    
169   Int_t nclusters2 = calcclusters->GetEntriesFast();
170   Int_t nclusters1 = exactclusters->GetEntriesFast();
171   
172   for(Int_t i=0; i<nclusters1; i++){
173     AliCluster * cl1 = (AliCluster*)exactclusters->UncheckedAt(i);
174     AliDigitCluster *cl2;
175     AliDifCluster diffc;
176     Int_t index=-1;
177     Float_t dx=10000;
178     Float_t dy=10000;
179     for (Int_t j=0; j<nclusters2;j++){
180       cl2 = (AliDigitCluster*)calcclusters->UncheckedAt(j);
181       if ( (cl2!=0)&& (cl1!=0)){ 
182         Float_t ddx = (cl2->fX-shx)-cl1->fX;
183         Float_t ddy = (cl2->fY-shy)-cl1->fY;
184         if ((ddx*ddx+ddy*ddy)<(dx*dx+dy*dy)){
185           dx=ddx;
186           dy=ddy;
187           index=j;
188         }      
189       }
190     }
191     
192     cl2 = (AliDigitCluster*)calcclusters->UncheckedAt(index);
193     if (cl2!=0){
194       diffc.fDx      =dx;
195       diffc.fDy      =dy;
196       
197       diffc.fAngleX  = cl1->fSigmaX2;
198       diffc.fAngleY  = cl1->fSigmaY2; 
199       diffc.fTracks[0] = cl2->fTracks[0];
200       diffc.fTracks[1] = cl2->fTracks[1];
201       diffc.fTracks[2] = cl2->fTracks[2];
202       diffc.fGenerTrack= cl1->fTracks[0];
203
204       diffc.fX       = cl2->fX;
205       diffc.fY       = cl2->fY;
206       diffc.fNx      = cl2->fNx;
207       diffc.fNy      = cl2->fNy;
208       diffc.fQ       = cl2->fQ;    
209       diffc.fOrigQ   = cl1->fQ;
210       diffc.fSigmaX2 = cl2->fSigmaX2;
211       diffc.fSigmaY2 = cl2->fSigmaY2;
212       diffc.fArea    = cl2->fArea;
213       diffc.fMax     = cl2->fMax;
214     }
215     else cout<<"pici[ici/n";
216     
217     new(diff[i]) AliDifCluster(diffc);
218   }
219   return diffclusters;
220