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   Int_t digcounter=0;
 389
 390   for(Int_t i=fRowMin; i<=fRowMax; i++)
 391     {
 392
 393       if((Int_t)rowPt->fRow != i) printf("AliL3FileHandler::AliDigits2RootFile : Mismatching row numbering!!!\n");
 394
 395       Int_t sector,row;
 396       fTransformer->Slice2Sector(fSlice,i,sector,row);
 397       AliSimDigits * dig = (AliSimDigits*)arr->CreateRow(sector,row);
 398       AliSimDigits *old_dig = (AliSimDigits*)old_array->LoadRow(sector,row);
 399       if(!old_dig)
 400         printf("AliL3FileHandler::AliDigits2RootFile : No padrow %d %d\n",sector,row);
 401
 402       AliL3DigitData *digPt = rowPt->fDigitData;
 403       digcounter=0;
 404       for(UInt_t j=0; j<rowPt->fNDigit; j++)
 405         {
 406           Int_t charge = (Int_t)digPt[j].fCharge;
 407           Int_t pad = (Int_t)digPt[j].fPad;
 408           Int_t time = (Int_t)digPt[j].fTime;
 409
 410           if(charge == 0) //Only write the digits that has not been removed
 411             continue;
 412           digcounter++;
 413           dig->SetDigitFast(charge,time,pad);
 414
 415           Int_t trackID[3] = {old_dig->GetTrackID(time,pad,0),old_dig->GetTrackID(time,pad,1),old_dig->GetTrackID(time,pad,2)};
 416           Int_t s_pad = pad;
 417           Int_t s_time = time - 1;
 418           while(trackID[0] < 0)
 419             {
 420               if(s_time >= 0 && s_time < fTransformer->GetNTimeBins() && s_pad >= 0 && s_pad < fTransformer->GetNPads(i))
 421                 {
 422                   if(old_dig->GetTrackID(s_time,s_pad,0) > 0)
 423                     {
 424                       trackID[0]=old_dig->GetTrackID(s_time,s_pad,0);
 425                       trackID[1]=old_dig->GetTrackID(s_time,s_pad,1);
 426                       trackID[2]=old_dig->GetTrackID(s_time,s_pad,2);
 427                     }
 428                 }
 429               if(s_pad == pad && s_time == time - 1)
 430                 s_time = time + 1;
 431               else if(s_pad == pad && s_time == time + 1)
 432                 {s_pad = pad - 1; s_time = time;}
 433               else if(s_pad == pad - 1 && s_time == time)
 434                 s_time = time - 1;
 435               else if(s_pad == pad - 1 && s_time == time - 1)
 436                 s_time = time + 1;
 437               else if(s_pad == pad - 1 && s_time == time + 1)
 438                 {s_pad = pad + 1; s_time = time;}
 439               else if(s_pad == pad + 1 && s_time == time)
 440                 s_time = time - 1;
 441               else if(s_pad == pad + 1 && s_time == time - 1)
 442                 s_time = time + 1;
 443               else
 444                 break;
 445             }
 446
 447           dig->SetTrackIDFast(trackID[0],time,pad,0);
 448           dig->SetTrackIDFast(trackID[1],time,pad,1);
 449           dig->SetTrackIDFast(trackID[2],time,pad,2);
 450
 451         }
 452       //cout<<"Wrote "<<digcounter<<" on row "<<i<<endl;
 453       UpdateRowPointer(rowPt);
 454       arr->StoreRow(sector,row);
 455       arr->ClearRow(sector,row);
 456       old_array->ClearRow(sector,row);
 457     }
 458   digFile->cd();
 459   char treeName[100];
 460   sprintf(treeName,"TreeD_%s_0",fParam->GetTitle());
 461   printf("Writing tree to file.....");
 462   arr->GetTree()->Write(treeName,TObject::kOverwrite);
 463   printf("done\n");
 464   digFile->Close();
 465   //arr->GetTree()->Delete();
 466   //delete arr;
 467 }
 468
 469 ///////////////////////////////////////// Point IO
 470 Bool_t AliL3FileHandler::AliPoints2Binary(){
 471   Bool_t out = kTRUE;
 472   UInt_t npoint;
 473   AliL3SpacePointData *data = AliPoints2Memory(npoint);
 474   out = Memory2Binary(npoint,data);
 475   Free();
 476   return out;
 477 }
 478
 479 AliL3SpacePointData * AliL3FileHandler::AliPoints2Memory(UInt_t & npoint){
 480   AliL3SpacePointData *data = 0;
 481   npoint=0;
 482   if(!fInAli){
 483     LOG(AliL3Log::kWarning,"AliL3FileHandler::AliPoints2Memory","File")
 484     <<"No Input avalible: no object TFile"<<ENDLOG;
 485     return 0;
 486   }
 487   if(!fInAli->IsOpen()){
 488     LOG(AliL3Log::kWarning,"AliL3FileHandler::AliPoints2Memory","File")
 489     <<"No Input avalible: TFile not opend"<<ENDLOG;
 490     return 0;
 491   }
 492   TDirectory *savedir = gDirectory;
 493   fInAli->cd();
 494
 495   Char_t cname[100];
 496   Int_t eventn = 0;
 497   sprintf(cname,"TreeC_TPC_%d",eventn);
 498   AliTPCClustersArray carray;
 499   carray.Setup(fParam);
 500   carray.SetClusterType("AliTPCcluster");
 501   Bool_t clusterok = carray.ConnectTree(cname);
 502   if(!clusterok) return 0;
 503
 504   AliTPCClustersRow ** clusterrow =
 505                new AliTPCClustersRow*[ (int)carray.GetTree()->GetEntries()];
 506   Int_t *rows = new int[ (int)carray.GetTree()->GetEntries()];
 507   Int_t *sects = new int[  (int)carray.GetTree()->GetEntries()];
 508   Int_t sum=0;
 509
 510   Int_t lslice,lrow;
 511   for(Int_t i=0; i<carray.GetTree()->GetEntries(); i++){
 512     AliSegmentID *s = carray.LoadEntry(i);
 513     Int_t sector,row;
 514     fParam->AdjustSectorRow(s->GetID(),sector,row);
 515     rows[i] = row;
 516     sects[i] = sector;
 517     clusterrow[i] = 0;
 518     fTransformer->Sector2Slice(lslice,lrow,sector,row);
 519     if(fSlice != lslice || lrow<fRowMin || lrow>fRowMax) continue;
 520     clusterrow[i] = carray.GetRow(sector,row);
 521     if(clusterrow[i])
 522       sum+=clusterrow[i]->GetArray()->GetEntriesFast();
 523   }
 524   UInt_t size = sum*sizeof(AliL3SpacePointData);
 525
 526   LOG(AliL3Log::kDebug,"AliL3FileHandler::AliPoints2Memory","File")
 527   <<AliL3Log::kDec<<"Found "<<sum<<" SpacePoints"<<ENDLOG;
 528
 529   data = (AliL3SpacePointData *) Allocate(size);
 530   npoint = sum;
 531   UInt_t n=0;
 532   for(Int_t i=0; i<carray.GetTree()->GetEntries(); i++){
 533     if(!clusterrow[i]) continue;
 534     Int_t row = rows[i];
 535     Int_t sector = sects[i];
 536     fTransformer->Sector2Slice(lslice,lrow,sector,row);
 537     Int_t entries_in_row = clusterrow[i]->GetArray()->GetEntriesFast();
 538     for(Int_t j = 0;j<entries_in_row;j++){
 539       AliTPCcluster *c = (AliTPCcluster*)(*clusterrow[i])[j];
 540       data[n].fZ = c->GetZ();
 541       data[n].fY = c->GetY();
 542       data[n].fX = fParam->GetPadRowRadii(sector,row);
 543       data[n].fID = n+((fSlice&0x7f)<<25)+((fPatch&0x7)<<22);//uli
 544       data[n].fPadRow = lrow;
 545       data[n].fXYErr = c->GetSigmaY2();
 546       data[n].fZErr = c->GetSigmaZ2();
 547       if(fMC) fprintf(fMC,"%d %d\n",data[n].fID,c->GetLabel(0));
 548       n++;
 549     }
 550   }
 551   for(Int_t i=0;i<carray.GetTree()->GetEntries();i++){
 552     Int_t row = rows[i];
 553     Int_t sector = sects[i];
 554     if(carray.GetRow(sector,row))
 555       carray.ClearRow(sector,row);
 556   }
 557
 558   delete [] clusterrow;
 559   delete [] rows;
 560   delete [] sects;
 561   savedir->cd();
 562
 563   return data;
 564 }
 565