Added function AliDigits2RootFile
[u/mrichter/AliRoot.git] / HLT / src / AliL3FileHandler.cxx
1
2 //Author:        Uli Frankenfeld
3 //Last Modified: 17.12.2000
4
5 #include <math.h>
6 #include <iostream.h>
7 #include <TObject.h>
8 #include <TFile.h>
9 #include <TTree.h>
10 #include <stdio.h>
11
12 #include "AliL3Transform.h"
13 #include "AliL3Logging.h"
14 #include "AliL3MemHandler.h"
15 #include "AliL3FileHandler.h"
16
17 #include "AliTPCDigitsArray.h"
18 #include "AliTPCClustersArray.h"
19 #include "AliTPCcluster.h"
20 #include "AliTPCClustersRow.h"
21 #include "AliTPCParam.h"
22 #include "AliSimDigits.h"
23
24 #include "AliL3DigitData.h"
25 #include "AliL3TrackSegmentData.h"
26 #include "AliL3SpacePointData.h"
27 #include "AliL3TrackArray.h"
28 //_____________________________________________________________
29 //
30 // The L3 Binary File handler 
31 //
32
33 ClassImp(AliL3FileHandler)
34
35 AliL3FileHandler::AliL3FileHandler(){
36   //Default constructor
37   fInAli = 0;
38   fParam = 0;
39   fTransformer = 0;
40   fMC =0;
41 }
42
43
44 AliL3FileHandler::~AliL3FileHandler(){
45   //Destructor
46   if(fTransformer) delete fTransformer;
47   if(fMC) CloseMCOutput();
48 }
49
50 Bool_t AliL3FileHandler::SetMCOutput(char *name){
51   fMC = fopen(name,"w");
52   if(!fMC){
53     LOG(AliL3Log::kWarning,"AliL3FileHandler::SetMCOutput","File Open")
54     <<"Pointer to File = 0x0 "<<ENDLOG;
55     return kFALSE;
56   }
57   return kTRUE;
58 }
59
60 Bool_t AliL3FileHandler::SetMCOutput(FILE *file){
61   fMC = file;
62   if(!fMC){
63     LOG(AliL3Log::kWarning,"AliL3FileHandler::SetMCOutput","File Open")
64     <<"Pointer to File = 0x0 "<<ENDLOG;
65     return kFALSE;
66   }
67   return kTRUE;
68 }
69
70 void AliL3FileHandler::CloseMCOutput(){
71   if(!fMC){
72     LOG(AliL3Log::kWarning,"AliL3FileHandler::CloseMCOutPut","File Close")
73     <<"Nothing to Close"<<ENDLOG;
74     return;
75   }
76   fclose(fMC);
77   fMC =0;
78 }
79
80 Bool_t AliL3FileHandler::SetAliInput(){
81   if(!fInAli->IsOpen()){
82     LOG(AliL3Log::kError,"AliL3FileHandler::SetAliInput","File Open")
83     <<"Ali File "<<fInAli->GetName()<<" does not exist"<<ENDLOG;
84     return kFALSE;
85   }
86   fParam = (AliTPCParam*)fInAli->Get("75x40_100x60");
87   if(!fParam){ 
88     LOG(AliL3Log::kError,"AliL3FileHandler::SetAliInput","File Open")
89     <<"No AliTPCParam 75x40_100x60 in File "<<fInAli->GetName()<<ENDLOG;
90      return kFALSE;
91   }
92   fTransformer = new AliL3Transform();
93   return kTRUE;
94 }
95
96 Bool_t AliL3FileHandler::SetAliInput(char *name){
97   fInAli= new TFile(name,"READ");
98   if(!fInAli){
99     LOG(AliL3Log::kWarning,"AliL3FileHandler::SetAliInput","File Open")
100     <<"Pointer to TFile = 0x0 "<<ENDLOG;
101     return kFALSE;
102   }
103   return SetAliInput();
104 }
105
106 Bool_t AliL3FileHandler::SetAliInput(TFile *file){
107   fInAli=file;
108   if(!fInAli){
109     LOG(AliL3Log::kWarning,"AliL3FileHandler::SetAliInput","File Open")
110     <<"Pointer to TFile = 0x0 "<<ENDLOG;
111     return kFALSE;
112   }
113   return SetAliInput();
114 }
115
116 void AliL3FileHandler::CloseAliInput(){
117   if(!fInAli){
118     LOG(AliL3Log::kWarning,"AliL3FileHandler::CloseAliInput","File Close")
119     <<"Nothing to Close"<<ENDLOG;
120      return;
121   }
122   if(fInAli->IsOpen()) fInAli->Close();
123   delete fInAli;
124   fInAli = 0;
125 }
126
127 Bool_t AliL3FileHandler::IsDigit(){
128   if(!fInAli){
129     LOG(AliL3Log::kWarning,"AliL3FileHandler::IsDigit","File")
130     <<"Pointer to TFile = 0x0 "<<ENDLOG;
131     return kTRUE;  //may you are use binary input which is Digits!!
132   }
133   TTree *t=(TTree*)fInAli->Get("TreeD_75x40_100x60");
134   if(t){
135     LOG(AliL3Log::kInformational,"AliL3FileHandler::IsDigit","File Type")
136     <<"Found Digit Tree -> Use Fast Cluster Finder"<<ENDLOG;
137     return kTRUE;
138   }
139   else{
140     LOG(AliL3Log::kInformational,"AliL3FileHandler::IsDigit","File Type")
141     <<"No Digit Tree -> Use Cluster Tree"<<ENDLOG;
142     return kFALSE;
143   }
144 }
145
146 ///////////////////////////////////////// Digit IO  
147 Bool_t AliL3FileHandler::AliDigits2Binary(){
148   Bool_t out = kTRUE;
149   UInt_t nrow;
150   AliL3DigitRowData* data = AliDigits2Memory(nrow);
151   out = Memory2Binary(nrow,data);
152   Free();
153   return out;
154 }
155
156
157 Bool_t AliL3FileHandler::AliDigits2CompBinary(){
158   Bool_t out = kTRUE;
159   UInt_t ndigits=0;
160   AliL3DigitRowData* digits = AliDigits2Memory(ndigits);
161   out = Memory2CompBinary(ndigits,digits);
162   Free();
163   return out;
164 }
165
166
167 AliL3DigitRowData * AliL3FileHandler::AliDigits2Memory(UInt_t & nrow){
168   AliL3DigitRowData *data = 0;
169   nrow=0;
170   if(!fInAli){
171     LOG(AliL3Log::kWarning,"AliL3FileHandler::AliDigits2Memory","File")
172     <<"No Input avalible: no object TFile"<<ENDLOG;
173     return 0; 
174   }
175   if(!fInAli->IsOpen()){
176     LOG(AliL3Log::kWarning,"AliL3FileHandler::AliDigits2Memory","File")
177     <<"No Input avalible: TFile not opend"<<ENDLOG;
178     return 0;
179   }
180
181   TDirectory *savedir = gDirectory;
182   fInAli->cd();
183   TTree *t=(TTree*)fInAli->Get("TreeD_75x40_100x60");
184   if(!t){
185     LOG(AliL3Log::kWarning,"AliL3FileHandler::AliDigits2Binary","AliRoot")
186     <<"No Digit Tree inside!"<<ENDLOG;
187     return 0;
188   }
189   AliSimDigits digarr, *dummy=&digarr;
190   t->GetBranch("Segment")->SetAddress(&dummy);
191   UShort_t dig;
192   Int_t time,pad,sector,row;
193   Int_t nrows=0;
194   Int_t ndigitcount=0;
195   Int_t entries = (Int_t)t->GetEntries();
196   Int_t ndigits[entries];
197   Int_t lslice,lrow;
198   for(Int_t n=0; n<t->GetEntries(); n++)
199     {
200       t->GetEvent(n);
201       fParam->AdjustSectorRow(digarr.GetID(),sector,row);
202       fTransformer->Sector2Slice(lslice,lrow,sector,row);
203       if(fSlice != lslice || lrow<fRowMin || lrow>fRowMax) continue;
204 //      if(fSlice != lslice) continue;
205
206       Float_t xyz[3];
207       ndigits[lrow] = 0;
208       digarr.First();
209       do {
210         time=digarr.CurrentRow();
211         pad=digarr.CurrentColumn();
212         dig = digarr.GetDigit(time,pad);
213         if(dig<=fParam->GetZeroSup()) continue;
214         if(time < fParam->GetMaxTBin()-1 && time > 0)
215           if(digarr.GetDigit(time+1,pad) <= fParam->GetZeroSup()
216              && digarr.GetDigit(time-1,pad) <= fParam->GetZeroSup())
217             continue;
218
219         fTransformer->Raw2Local(xyz,sector,row,pad,time);
220         if(fParam->GetPadRowRadii(sector,row)<230./250.*fabs(xyz[2]))
221           continue;
222
223         ndigits[lrow]++; //for this row only
224         ndigitcount++;  //total number of digits to be published
225
226       } while (digarr.Next());
227
228       nrows++;
229     }
230   Int_t size = sizeof(AliL3DigitData)*ndigitcount
231                                       + nrows*sizeof(AliL3DigitRowData);
232
233   LOG(AliL3Log::kDebug,"AliL3FileHandler::AliDigits2Memory","Digits")
234   <<AliL3Log::kDec<<"Found "<<ndigitcount<<" Digits"<<ENDLOG;
235
236   data=(AliL3DigitRowData*) Allocate(size);
237   nrow = (UInt_t)nrows;
238   AliL3DigitRowData *tempPt = data;
239   for(Int_t n=0; n<t->GetEntries(); n++)
240     {
241       t->GetEvent(n);
242
243       Float_t xyz[3];
244       fParam->AdjustSectorRow(digarr.GetID(),sector,row);
245       fTransformer->Sector2Slice(lslice,lrow,sector,row);
246 //      if(fSlice != lslice) continue;
247       if(fSlice != lslice || lrow<fRowMin || lrow>fRowMax) continue;
248       tempPt->fRow = lrow;
249       tempPt->fNDigit = ndigits[lrow];
250
251       Int_t localcount=0;
252       digarr.First();
253       do {
254         dig=digarr.CurrentDigit();
255         if (dig<=fParam->GetZeroSup()) continue;
256         time=digarr.CurrentRow();
257         pad=digarr.CurrentColumn();
258         if(time < fParam->GetMaxTBin()-1 && time > 0)
259           if(digarr.GetDigit(time-1,pad) <= fParam->GetZeroSup() &&
260              digarr.GetDigit(time+1,pad) <= fParam->GetZeroSup()) continue;
261
262         //Exclude data outside cone:
263         fTransformer->Raw2Local(xyz,sector,row,pad,time);
264         if(fParam->GetPadRowRadii(sector,row)<230./250.*fabs(xyz[2]))
265           continue;
266
267         if(localcount >= ndigits[lrow])
268           LOG(AliL3Log::kFatal,"AliL3FileHandler::AliDigits2Binary","Memory")
269           <<AliL3Log::kDec<<"Mismatch: localcount "<<localcount<<" ndigits "
270           <<ndigits[lrow]<<ENDLOG;
271
272         tempPt->fDigitData[localcount].fCharge=dig;
273         tempPt->fDigitData[localcount].fPad=pad;
274         tempPt->fDigitData[localcount].fTime=time;
275         localcount++;
276       } while (digarr.Next());
277
278       Byte_t *tmp = (Byte_t*)tempPt;
279       Int_t size = sizeof(AliL3DigitRowData)
280                                       + ndigits[lrow]*sizeof(AliL3DigitData);
281       tmp += size;
282       tempPt = (AliL3DigitRowData*)tmp;
283     }
284   savedir->cd(); 
285   return data;
286 }
287
288 void AliL3FileHandler::AliDigits2RootFile(AliL3DigitRowData *rowPt,Char_t *new_digitsfile)
289 {
290   //Write digits to a new alirootfile.
291
292   if(!fInAli)
293     {
294       printf("AliL3FileHandler::AliDigits2RootFile : No input alirootfile\n");
295       return;
296     }
297   if(!fParam)
298     {
299       printf("AliL3FileHandler::AliDigits2RootFile : No parameter object. Run on rootfile\n");
300       return;
301     }
302   if(!fTransformer)
303     {
304       printf("AliL3FileHandler::AliDigits2RootFile : No transform object\n");
305       return;
306     }
307   
308   //Get the original digitstree:
309   fInAli->cd();
310   AliTPCDigitsArray *old_array = new AliTPCDigitsArray();
311   old_array->Setup(fParam);
312   old_array->SetClass("AliSimDigits");
313   Bool_t ok = old_array->ConnectTree("TreeD_75x40_100x60");
314   if(!ok)
315     {
316       printf("AliL3FileHandler::AliDigits2RootFile : No digits tree object\n");
317       return;
318     }
319     
320   TFile *digFile = new TFile(new_digitsfile,"RECREATE");
321   digFile->cd();
322   
323   //setup a new one:
324   AliTPCDigitsArray *arr = new AliTPCDigitsArray; 
325   arr->SetClass("AliSimDigits");
326   arr->Setup(fParam);
327   arr->MakeTree();
328   if(fRowMin !=0 || fRowMax != 175)
329     printf("\n AliL3FileHandler::AliDigits2RootFile : Rather stupid row numbers...%d %d\n\n",fRowMin,fRowMax);
330   for(Int_t i=fRowMin; i<=fRowMax; i++)
331     {
332       Int_t sector,row;
333       fTransformer->Slice2Sector(fSlice,i,sector,row);
334       AliDigits * dig = arr->CreateRow(sector,row);
335       AliDigits *old_dig = old_array->LoadRow(sector,row);
336       if(!old_dig)
337         printf("AliL3FileHandler::AliDigits2RootFile : No padrow %d %d\n",sector,row);
338       
339       AliL3DigitData *digPt = rowPt->fDigitData;
340       for(UInt_t j=0; j<rowPt->fNDigit; j++)
341         {
342           UShort_t charge = digPt[j].fCharge;
343           UChar_t pad = digPt[j].fPad;
344           UShort_t time = digPt[j].fTime;
345           dig->SetDigitFast(charge,time,pad);
346           ((AliSimDigits*)dig)->SetTrackIDFast(((AliSimDigits*)old_dig)->GetTrackID((Int_t)time,(Int_t)pad,0),time,pad,0);
347           ((AliSimDigits*)dig)->SetTrackIDFast(((AliSimDigits*)old_dig)->GetTrackID((Int_t)time,(Int_t)pad,1),time,pad,1);
348           ((AliSimDigits*)dig)->SetTrackIDFast(((AliSimDigits*)old_dig)->GetTrackID((Int_t)time,(Int_t)pad,2),time,pad,2);
349           
350         }
351       UpdateRowPointer(rowPt);
352       arr->StoreRow(sector,row);
353       arr->ClearRow(sector,row);  
354       old_array->ClearRow(sector,row);
355     }
356   digFile->cd();
357   char treeName[100];
358   sprintf(treeName,"TreeD_%s",fParam->GetTitle());
359   arr->GetTree()->Write(treeName,TObject::kOverwrite);
360   fParam->Write(fParam->GetTitle());
361   digFile->Close();
362   delete digFile;
363 }
364
365 ///////////////////////////////////////// Point IO  
366 Bool_t AliL3FileHandler::AliPoints2Binary(){
367   Bool_t out = kTRUE;
368   UInt_t npoint;
369   AliL3SpacePointData *data = AliPoints2Memory(npoint);
370   out = Memory2Binary(npoint,data);
371   Free();
372   return out;
373 }
374
375 AliL3SpacePointData * AliL3FileHandler::AliPoints2Memory(UInt_t & npoint){
376   AliL3SpacePointData *data = 0;
377   npoint=0;
378   if(!fInAli){
379     LOG(AliL3Log::kWarning,"AliL3FileHandler::AliPoints2Memory","File")
380     <<"No Input avalible: no object TFile"<<ENDLOG;
381     return 0;
382   }
383   if(!fInAli->IsOpen()){
384     LOG(AliL3Log::kWarning,"AliL3FileHandler::AliPoints2Memory","File")
385     <<"No Input avalible: TFile not opend"<<ENDLOG;
386     return 0;
387   }
388   TDirectory *savedir = gDirectory;
389   fInAli->cd();
390
391   AliTPCClustersArray carray;
392   carray.Setup(fParam);
393   carray.SetClusterType("AliTPCcluster");
394   Bool_t clusterok = carray.ConnectTree("Segment Tree");
395   if(!clusterok) return 0;
396
397   AliTPCClustersRow ** clusterrow = 
398                new AliTPCClustersRow*[ (int)carray.GetTree()->GetEntries()];
399   Int_t *rows = new int[ (int)carray.GetTree()->GetEntries()];
400   Int_t *sects = new int[  (int)carray.GetTree()->GetEntries()];
401   Int_t sum=0;
402
403   Int_t lslice,lrow;
404   for(Int_t i=0; i<carray.GetTree()->GetEntries(); i++){
405     AliSegmentID *s = carray.LoadEntry(i);
406     Int_t sector,row;
407     fParam->AdjustSectorRow(s->GetID(),sector,row);
408     rows[i] = row;
409     sects[i] = sector;
410     clusterrow[i] = 0;
411     fTransformer->Sector2Slice(lslice,lrow,sector,row);
412     if(fSlice != lslice || lrow<fRowMin || lrow>fRowMax) continue;
413     clusterrow[i] = carray.GetRow(sector,row);
414     if(clusterrow[i])
415       sum+=clusterrow[i]->GetArray()->GetEntriesFast();
416   }
417   UInt_t size = sum*sizeof(AliL3SpacePointData);
418
419   LOG(AliL3Log::kDebug,"AliL3FileHandler::AliPoints2Memory","File")
420   <<AliL3Log::kDec<<"Found "<<sum<<" SpacePoints"<<ENDLOG;
421
422   data = (AliL3SpacePointData *) Allocate(size);
423   npoint = sum;
424   UInt_t n=0; 
425   for(Int_t i=0; i<carray.GetTree()->GetEntries(); i++){
426     if(!clusterrow[i]) continue;
427     Int_t row = rows[i];
428     Int_t sector = sects[i];
429     fTransformer->Sector2Slice(lslice,lrow,sector,row);
430     Int_t entries_in_row = clusterrow[i]->GetArray()->GetEntriesFast();
431     for(Int_t j = 0;j<entries_in_row;j++){
432       AliTPCcluster *c = (AliTPCcluster*)(*clusterrow[i])[j];
433       data[n].fZ = c->GetZ();
434       data[n].fY = c->GetY();
435       data[n].fX = fParam->GetPadRowRadii(sector,row);
436       data[n].fID = n+((fSlice&0x7f)<<25)+((fPatch&0x7)<<22);//uli
437       data[n].fPadRow = lrow;
438       data[n].fXYErr = c->GetSigmaY2();
439       data[n].fZErr = c->GetSigmaZ2();
440       if(fMC) fprintf(fMC,"%d %d\n",data[n].fID,c->GetLabel(0));
441       n++;
442     }
443   }
444   for(Int_t i=0;i<carray.GetTree()->GetEntries();i++){
445     Int_t row = rows[i];
446     Int_t sector = sects[i];
447     if(carray.GetRow(sector,row))
448       carray.ClearRow(sector,row);
449   }
450
451   delete [] clusterrow;
452   delete [] rows;
453   delete [] sects;
454   savedir->cd();   
455
456   return data;
457 }
458