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