HLT/src/AliL3FileHandler.cxx

   1 //$Id$
   2
   3 // Author: Uli Frankenfeld <mailto:franken@fi.uib.no>, Anders Vestbo <mailto:vestbo$fi.uib.no>
   4 //*-- Copyright &copy Uli
   5
   6 #include "AliL3StandardIncludes.h"
   7 #include <TClonesArray.h>
   8 #include <TSystem.h>
   9
  10 #include <AliTPCDigitsArray.h>
  11 #include <AliTPCClustersArray.h>
  12 #include <AliTPCcluster.h>
  13 #include <AliTPCClustersRow.h>
  14 #include <AliSimDigits.h>
  15
  16 #include "AliL3Logging.h"
  17 #include "AliL3Transform.h"
  18 #include "AliL3MemHandler.h"
  19 #include "AliL3FileHandler.h"
  20 #include "AliL3DigitData.h"
  21 #include "AliL3TrackSegmentData.h"
  22 #include "AliL3SpacePointData.h"
  23 #include "AliL3TrackArray.h"
  24
  25
  26 //_____________________________________________________________
  27 // AliL3FileHandler
  28 //
  29 // The HLT ROOT <-> binary files handling class
  30 //
  31 // This class provides the interface between AliROOT files,
  32 // and HLT binary files. It should be used for converting
  33 // TPC data stored in AliROOT format (outputfile from a simulation),
  34 // into the data format currently used by in the HLT framework.
  35 // This enables the possibility to always use the same data format,
  36 // whether you are using a binary file as an input, or a AliROOT file.
  37 //
  38 // For example on how to create binary files from a AliROOT simulation,
  39 // see example macro exa/Binary.C.
  40 //
  41 // For reading a AliROOT file into HLT format in memory, do the following:
  42 //
  43 // AliL3FileHandler file;
  44 // file.Init(slice,patch);
  45 // file.SetAliInput("galice.root");
  46 // AliL3DigitRowData *dataPt = (AliL3DigitRowData*)file.AliDigits2Memory(nrows,eventnr);
  47 //
  48 // All the data are then stored in memory and accessible via the pointer dataPt.
  49 // Accesing the data is then identical to the example 1) showed in AliL3MemHandler class.
  50 //
  51 // For converting the data back, and writing it to a new AliROOT file do:
  52 //
  53 // AliL3FileHandler file;
  54 // file.Init(slice,patch);
  55 // file.SetAliInput("galice.root");
  56 // file.Init(slice,patch,NumberOfRowsInPatch);
  57 // file.AliDigits2RootFile(dataPt,"new_galice.root");
  58 // file.CloseAliInput();
  59
  60 ClassImp(AliL3FileHandler)
  61
  62 AliL3FileHandler::AliL3FileHandler()
  63 {
  64   //Default constructor
  65   fInAli = 0;
  66   fParam = 0;
  67   fMC =0;
  68   fLastIndex=0;
  69   fDigits=0;
  70   fDigitsTree=0;
  71 }
  72
  73 AliL3FileHandler::~AliL3FileHandler()
  74 {
  75   //Destructor
  76   if(fMC) CloseMCOutput();
  77   FreeDigitsTree();
  78   if(fInAli) CloseAliInput();
  79
  80 }
  81
  82 void AliL3FileHandler::FreeDigitsTree()
  83 {
  84   if(!fDigitsTree)
  85     {
  86       LOG(AliL3Log::kWarning,"AliL3FileHandler::FreeDigitsTree()","Pointer")
  87         <<"Cannot free digitstree, it is not present"<<ENDLOG;
  88       return;
  89     }
  90   fDigits=0;
  91   fDigitsTree->Delete();
  92   fDigitsTree=0;
  93 }
  94
  95 Bool_t AliL3FileHandler::SetMCOutput(char *name)
  96 {
  97   fMC = fopen(name,"w");
  98   if(!fMC){
  99     LOG(AliL3Log::kWarning,"AliL3FileHandler::SetMCOutput","File Open")
 100       <<"Pointer to File = 0x0 "<<ENDLOG;
 101     return kFALSE;
 102   }
 103   return kTRUE;
 104 }
 105
 106 Bool_t AliL3FileHandler::SetMCOutput(FILE *file)
 107 {
 108   fMC = file;
 109   if(!fMC){
 110     LOG(AliL3Log::kWarning,"AliL3FileHandler::SetMCOutput","File Open")
 111       <<"Pointer to File = 0x0 "<<ENDLOG;
 112     return kFALSE;
 113   }
 114   return kTRUE;
 115 }
 116
 117 void AliL3FileHandler::CloseMCOutput()
 118 {
 119   if(!fMC){
 120     LOG(AliL3Log::kWarning,"AliL3FileHandler::CloseMCOutPut","File Close")
 121       <<"Nothing to Close"<<ENDLOG;
 122     return;
 123   }
 124   fclose(fMC);
 125   fMC =0;
 126 }
 127
 128 Bool_t AliL3FileHandler::SetAliInput()
 129 {
 130   if(!fInAli->IsOpen()){
 131     LOG(AliL3Log::kError,"AliL3FileHandler::SetAliInput","File Open")
 132       <<"Ali File "<<fInAli->GetName()<<" does not exist"<<ENDLOG;
 133     return kFALSE;
 134   }
 135   fParam = (AliTPCParam*)fInAli->Get(AliL3Transform::GetParamName());
 136   if(!fParam){
 137     LOG(AliL3Log::kError,"AliL3FileHandler::SetAliInput","File Open")
 138       <<"No AliTPCParam "<<AliL3Transform::GetParamName()<<" in File "<<fInAli->GetName()<<ENDLOG;
 139     return kFALSE;
 140   }
 141   return kTRUE;
 142 }
 143
 144 Bool_t AliL3FileHandler::SetAliInput(char *name)
 145 {
 146   //Open the AliROOT file with name.
 147
 148   fInAli= new TFile(name,"READ");
 149   if(!fInAli){
 150     LOG(AliL3Log::kWarning,"AliL3FileHandler::SetAliInput","File Open")
 151     <<"Pointer to TFile = 0x0 "<<ENDLOG;
 152     return kFALSE;
 153   }
 154   return SetAliInput();
 155 }
 156
 157 Bool_t AliL3FileHandler::SetAliInput(TFile *file)
 158 {
 159   //Specify already opened AliROOT file to use as an input.
 160
 161   fInAli=file;
 162   if(!fInAli){
 163     LOG(AliL3Log::kWarning,"AliL3FileHandler::SetAliInput","File Open")
 164     <<"Pointer to TFile = 0x0 "<<ENDLOG;
 165     return kFALSE;
 166   }
 167   return SetAliInput();
 168 }
 169
 170 void AliL3FileHandler::CloseAliInput()
 171 {
 172   if(!fInAli){
 173     LOG(AliL3Log::kWarning,"AliL3FileHandler::CloseAliInput","File Close")
 174       <<"Nothing to Close"<<ENDLOG;
 175     return;
 176   }
 177   if(fInAli->IsOpen()) fInAli->Close();
 178   delete fInAli;
 179   fInAli = 0;
 180
 181 }
 182
 183 Bool_t AliL3FileHandler::IsDigit(Int_t event)
 184 {
 185   //Check if there is a TPC digit tree in the current file.
 186   //Return kTRUE if tree was found, and kFALSE if not found.
 187
 188   if(!fInAli){
 189     LOG(AliL3Log::kWarning,"AliL3FileHandler::IsDigit","File")
 190     <<"Pointer to TFile = 0x0 "<<ENDLOG;
 191     return kTRUE;  //maybe you are using binary input which is Digits!!
 192   }
 193   Char_t name[1024];
 194   sprintf(name,"TreeD_%s_%d",AliL3Transform::GetParamName(),event);
 195   TTree *t=(TTree*)fInAli->Get(name);
 196   if(t){
 197     LOG(AliL3Log::kInformational,"AliL3FileHandler::IsDigit","File Type")
 198     <<"Found Digit Tree -> Use Fast Cluster Finder"<<ENDLOG;
 199     return kTRUE;
 200   }
 201   else{
 202     LOG(AliL3Log::kInformational,"AliL3FileHandler::IsDigit","File Type")
 203     <<"No Digit Tree -> Use Cluster Tree"<<ENDLOG;
 204     return kFALSE;
 205   }
 206 }
 207
 208 ///////////////////////////////////////// Digit IO
 209 Bool_t AliL3FileHandler::AliDigits2Binary(Int_t event)
 210 {
 211   Bool_t out = kTRUE;
 212   UInt_t nrow;
 213   AliL3DigitRowData* data = AliDigits2Memory(nrow,event);
 214   out = Memory2Binary(nrow,data);
 215   Free();
 216   return out;
 217 }
 218
 219 Bool_t AliL3FileHandler::AliDigits2CompBinary(Int_t event,Bool_t altro)
 220 {
 221   //Convert AliROOT TPC data, into HLT data format.
 222   //event specifies the event you want in the aliroot file.
 223
 224   Bool_t out = kTRUE;
 225   UInt_t ndigits=0;
 226   AliL3DigitRowData *digits=0;
 227   if(!altro)
 228     digits = AliDigits2Memory(ndigits,event);
 229   else
 230     digits = AliAltroDigits2Memory(ndigits,event);
 231   out = Memory2CompBinary(ndigits,digits);
 232   Free();
 233   return out;
 234 }
 235
 236 AliL3DigitRowData * AliL3FileHandler::AliDigits2Memory(UInt_t & nrow,Int_t event)
 237 {
 238   //Read data from AliROOT file into memory, and store it in the HLT data format.
 239   //Returns a pointer to the data.
 240
 241   AliL3DigitRowData *data = 0;
 242   nrow=0;
 243
 244   if(!fInAli){
 245     LOG(AliL3Log::kWarning,"AliL3FileHandler::AliDigits2Memory","File")
 246     <<"No Input avalible: no object TFile"<<ENDLOG;
 247     return 0;
 248   }
 249   if(!fInAli->IsOpen()){
 250     LOG(AliL3Log::kWarning,"AliL3FileHandler::AliDigits2Memory","File")
 251     <<"No Input avalible: TFile not opened"<<ENDLOG;
 252     return 0;
 253   }
 254
 255   if(!fDigitsTree)
 256     GetDigitsTree(event);
 257
 258   UShort_t dig;
 259   Int_t time,pad,sector,row;
 260   Int_t nrows=0;
 261   Int_t ndigitcount=0;
 262   Int_t entries = (Int_t)fDigitsTree->GetEntries();
 263   Int_t ndigits[entries];
 264   Int_t lslice,lrow;
 265   Float_t xyz[3];
 266
 267   for(Int_t n=fLastIndex; n<fDigitsTree->GetEntries(); n++)
 268     {
 269       fDigitsTree->GetEvent(n);
 270       fParam->AdjustSectorRow(fDigits->GetID(),sector,row);
 271       AliL3Transform::Sector2Slice(lslice,lrow,sector,row);
 272       //if(fSlice != lslice || lrow<fRowMin || lrow>fRowMax) continue;
 273       if(lslice < fSlice) continue;
 274       if(lslice != fSlice) break;
 275       if(lrow < fRowMin) continue;
 276       if(lrow > fRowMax) break;
 277
 278       ndigits[lrow] = 0;
 279       fDigits->First();
 280       do {
 281         time=fDigits->CurrentRow();
 282         pad=fDigits->CurrentColumn();
 283         dig = fDigits->GetDigit(time,pad);
 284         if(dig <= fParam->GetZeroSup()) continue;
 285         if(dig >= AliL3Transform::GetADCSat())
 286           dig = AliL3Transform::GetADCSat();
 287
 288         AliL3Transform::Raw2Local(xyz,sector,row,pad,time);
 289         if(fParam->GetPadRowRadii(sector,row)<230./250.*fabs(xyz[2]))
 290           continue;
 291
 292         ndigits[lrow]++; //for this row only
 293         ndigitcount++;   //total number of digits to be published
 294
 295       } while (fDigits->Next());
 296       //cout << lrow << " " << ndigits[lrow] << " - " << ndigitcount << endl;
 297       nrows++;
 298     }
 299
 300   Int_t size = sizeof(AliL3DigitData)*ndigitcount
 301     + nrows*sizeof(AliL3DigitRowData);
 302
 303   LOG(AliL3Log::kDebug,"AliL3FileHandler::AliDigits2Memory","Digits")
 304     <<AliL3Log::kDec<<"Found "<<ndigitcount<<" Digits"<<ENDLOG;
 305
 306   data=(AliL3DigitRowData*) Allocate(size);
 307   nrow = (UInt_t)nrows;
 308   AliL3DigitRowData *tempPt = data;
 309   for(Int_t n=fLastIndex; n<fDigitsTree->GetEntries(); n++)
 310     {
 311       fDigitsTree->GetEvent(n);
 312       fParam->AdjustSectorRow(fDigits->GetID(),sector,row);
 313       AliL3Transform::Sector2Slice(lslice,lrow,sector,row);
 314       //if(fSlice != lslice || lrow<fRowMin || lrow>fRowMax) continue;
 315       if(lslice < fSlice) continue;
 316       if(lslice != fSlice) break;
 317       if(lrow < fRowMin) continue;
 318       if(lrow > fRowMax) break;
 319
 320       tempPt->fRow = lrow;
 321       tempPt->fNDigit = ndigits[lrow];
 322
 323       Int_t localcount=0;
 324       fDigits->First();
 325       do {
 326         time=fDigits->CurrentRow();
 327         pad=fDigits->CurrentColumn();
 328         dig = fDigits->GetDigit(time,pad);
 329         if (dig <= fParam->GetZeroSup()) continue;
 330         if(dig >= AliL3Transform::GetADCSat())
 331           dig = AliL3Transform::GetADCSat();
 332
 333         //Exclude data outside cone:
 334         AliL3Transform::Raw2Local(xyz,sector,row,pad,time);
 335         if(fParam->GetPadRowRadii(sector,row)<230./250.*fabs(xyz[2]))
 336           continue;
 337
 338         if(localcount >= ndigits[lrow])
 339           LOG(AliL3Log::kFatal,"AliL3FileHandler::AliDigits2Binary","Memory")
 340             <<AliL3Log::kDec<<"Mismatch: localcount "<<localcount<<" ndigits "
 341             <<ndigits[lrow]<<ENDLOG;
 342
 343         tempPt->fDigitData[localcount].fCharge=dig;
 344         tempPt->fDigitData[localcount].fPad=pad;
 345         tempPt->fDigitData[localcount].fTime=time;
 346 #ifdef do_mc
 347         tempPt->fDigitData[localcount].fTrackID[0] = fDigits->GetTrackID(time,pad,0);
 348         tempPt->fDigitData[localcount].fTrackID[1] = fDigits->GetTrackID(time,pad,1);
 349         tempPt->fDigitData[localcount].fTrackID[2] = fDigits->GetTrackID(time,pad,2);
 350 #endif
 351         localcount++;
 352       } while (fDigits->Next());
 353
 354       Byte_t *tmp = (Byte_t*)tempPt;
 355       Int_t size = sizeof(AliL3DigitRowData)
 356                                       + ndigits[lrow]*sizeof(AliL3DigitData);
 357       tmp += size;
 358       tempPt = (AliL3DigitRowData*)tmp;
 359 #ifdef ASVVERSION
 360       fLastIndex=n;
 361 #endif
 362     }
 363 #ifdef ASVVERSION
 364   fLastIndex++;
 365 #endif
 366   return data;
 367 }
 368
 369 AliL3DigitRowData * AliL3FileHandler::AliAltroDigits2Memory(UInt_t & nrow,Int_t event)
 370 {
 371   //Read data from AliROOT file into memory, and store it in the HLT data format.
 372   //Returns a pointer to the data.
 373   //This functions filter out single timebins, which is noise. The timebins which
 374   //are removed are timebins which have the 4 zero neighbours;
 375   //(pad-1,time),(pad+1,time),(pad,time-1),(pad,time+1).
 376
 377   AliL3DigitRowData *data = 0;
 378   nrow=0;
 379
 380   if(!fInAli){
 381     LOG(AliL3Log::kWarning,"AliL3FileHandler::AliAltroDigits2Memory","File")
 382     <<"No Input avalible: no object TFile"<<ENDLOG;
 383     return 0;
 384   }
 385   if(!fInAli->IsOpen()){
 386     LOG(AliL3Log::kWarning,"AliL3FileHandler::AliAltroDigits2Memory","File")
 387     <<"No Input avalible: TFile not opened"<<ENDLOG;
 388     return 0;
 389   }
 390
 391   if(!fDigitsTree)
 392     GetDigitsTree(event);
 393
 394   UShort_t dig;
 395   Int_t time,pad,sector,row;
 396   Int_t nrows=0;
 397   Int_t ndigitcount=0;
 398   Int_t entries = (Int_t)fDigitsTree->GetEntries();
 399   Int_t ndigits[entries];
 400   Int_t lslice,lrow;
 401   Float_t xyz[3];
 402
 403   for(Int_t n=fLastIndex; n<fDigitsTree->GetEntries(); n++)
 404     {
 405       fDigitsTree->GetEvent(n);
 406       fParam->AdjustSectorRow(fDigits->GetID(),sector,row);
 407       AliL3Transform::Sector2Slice(lslice,lrow,sector,row);
 408       //if(fSlice != lslice || lrow<fRowMin || lrow>fRowMax) continue;
 409       if(lslice < fSlice) continue;
 410       if(lslice != fSlice) break;
 411       if(lrow < fRowMin) continue;
 412       if(lrow > fRowMax) break;
 413
 414       ndigits[lrow] = 0;
 415       fDigits->ExpandBuffer();
 416       fDigits->ExpandTrackBuffer();
 417       for(Int_t i=0; i<fDigits->GetNCols(); i++)
 418         {
 419           for(Int_t j=0; j<fDigits->GetNRows(); j++)
 420             {
 421               pad=i;
 422               time=j;
 423               dig = fDigits->GetDigitFast(time,pad);
 424               if(dig <= fParam->GetZeroSup()) continue;
 425               if(dig >= AliL3Transform::GetADCSat())
 426                 dig = AliL3Transform::GetADCSat();
 427
 428               //Check for single timebins, and remove them because they are noise for sure.
 429               if(i>0 && i<fDigits->GetNCols()-1 && j>0 && j<fDigits->GetNRows()-1)
 430                 if(fDigits->GetDigitFast(time,pad-1)<=fParam->GetZeroSup() &&
 431                    fDigits->GetDigitFast(time-1,pad)<=fParam->GetZeroSup() &&
 432                    fDigits->GetDigitFast(time+1,pad)<=fParam->GetZeroSup() &&
 433                    fDigits->GetDigitFast(time,pad+1)<=fParam->GetZeroSup())
 434                   continue;
 435
 436               //Boundaries:
 437               if(i==0) //pad ==0
 438                 {
 439                   if(j < fDigits->GetNRows()-1 && j > 0)
 440                     {
 441                       if(fDigits->GetDigitFast(time-1,pad)<=fParam->GetZeroSup() &&
 442                          fDigits->GetDigitFast(time+1,pad)<=fParam->GetZeroSup() &&
 443                          fDigits->GetDigitFast(time,pad+1)<=fParam->GetZeroSup())
 444                         continue;
 445                     }
 446                   else if(j > 0)
 447                     {
 448                       if(fDigits->GetDigitFast(time-1,pad)<=fParam->GetZeroSup() &&
 449                          fDigits->GetDigitFast(time,pad+1)<=fParam->GetZeroSup())
 450                         continue;
 451                     }
 452                 }
 453               if(j==0)
 454                 {
 455                   if(i < fDigits->GetNCols()-1 && i > 0)
 456                     {
 457                       if(fDigits->GetDigitFast(time,pad-1)<=fParam->GetZeroSup() &&
 458                          fDigits->GetDigitFast(time,pad+1)<=fParam->GetZeroSup() &&
 459                          fDigits->GetDigitFast(time+1,pad)<=fParam->GetZeroSup())
 460                         continue;
 461                     }
 462                   else if(i > 0)
 463                     {
 464                       if(fDigits->GetDigitFast(time,pad-1)<=fParam->GetZeroSup() &&
 465                          fDigits->GetDigitFast(time+1,pad)<=fParam->GetZeroSup())
 466                         continue;
 467                     }
 468                 }
 469
 470               if(i == fDigits->GetNCols()-1)
 471                 {
 472
 473                   if(j>0 && j<fDigits->GetNRows()-1)
 474                     {
 475                       if(fDigits->GetDigitFast(time-1,pad)<=fParam->GetZeroSup() &&
 476                          fDigits->GetDigitFast(time+1,pad)<=fParam->GetZeroSup() &&
 477                          fDigits->GetDigitFast(time,pad-1)<=fParam->GetZeroSup())
 478                         continue;
 479                     }
 480                   else if(j==0 && j<fDigits->GetNRows()-1)
 481                     {
 482                       if(fDigits->GetDigitFast(time+1,pad)<=fParam->GetZeroSup() &&
 483                          fDigits->GetDigitFast(time,pad-1)<=fParam->GetZeroSup())
 484                         continue;
 485                     }
 486                   else
 487                     {
 488                       if(fDigits->GetDigitFast(time-1,pad)<=fParam->GetZeroSup() &&
 489                          fDigits->GetDigitFast(time,pad-1)<=fParam->GetZeroSup())
 490                         continue;
 491                     }
 492                 }
 493
 494               if(j==fDigits->GetNRows()-1)
 495                 {
 496                   if(i>0 && i<fDigits->GetNCols()-1)
 497                     {
 498                       if(fDigits->GetDigitFast(time,pad-1)<=fParam->GetZeroSup() &&
 499                          fDigits->GetDigitFast(time,pad+1)<=fParam->GetZeroSup() &&
 500                          fDigits->GetDigitFast(time-1,pad)<=fParam->GetZeroSup())
 501                         continue;
 502                     }
 503                   else if(i==0 && fDigits->GetNCols()-1)
 504                     {
 505                       if(fDigits->GetDigitFast(time,pad+1)<=fParam->GetZeroSup() &&
 506                          fDigits->GetDigitFast(time-1,pad)<=fParam->GetZeroSup())
 507                         continue;
 508                     }
 509                   else
 510                     {
 511                       if(fDigits->GetDigitFast(time,pad-1)<=fParam->GetZeroSup() &&
 512                          fDigits->GetDigitFast(time-1,pad)<=fParam->GetZeroSup())
 513                         continue;
 514                     }
 515                 }
 516
 517               AliL3Transform::Raw2Local(xyz,sector,row,pad,time);
 518               if(fParam->GetPadRowRadii(sector,row)<230./250.*fabs(xyz[2]))
 519                 continue;
 520
 521               ndigits[lrow]++; //for this row only
 522               ndigitcount++;   //total number of digits to be published
 523             }
 524         }
 525       nrows++;
 526     }
 527   Int_t size = sizeof(AliL3DigitData)*ndigitcount
 528     + nrows*sizeof(AliL3DigitRowData);
 529
 530   LOG(AliL3Log::kDebug,"AliL3FileHandler::AliAltroDigits2Memory","Digits")
 531     <<AliL3Log::kDec<<"Found "<<ndigitcount<<" Digits"<<ENDLOG;
 532
 533   data=(AliL3DigitRowData*) Allocate(size);
 534   nrow = (UInt_t)nrows;
 535   AliL3DigitRowData *tempPt = data;
 536   for(Int_t n=fLastIndex; n<fDigitsTree->GetEntries(); n++)
 537     {
 538       fDigitsTree->GetEvent(n);
 539       fParam->AdjustSectorRow(fDigits->GetID(),sector,row);
 540       AliL3Transform::Sector2Slice(lslice,lrow,sector,row);
 541       //if(fSlice != lslice || lrow<fRowMin || lrow>fRowMax) continue;
 542       if(lslice < fSlice) continue;
 543       if(lslice != fSlice) break;
 544       if(lrow < fRowMin) continue;
 545       if(lrow > fRowMax) break;
 546
 547       tempPt->fRow = lrow;
 548       tempPt->fNDigit = ndigits[lrow];
 549
 550       Int_t localcount=0;
 551       fDigits->ExpandBuffer();
 552       fDigits->ExpandTrackBuffer();
 553       for(Int_t i=0; i<fDigits->GetNCols(); i++)
 554         {
 555           for(Int_t j=0; j<fDigits->GetNRows(); j++)
 556             {
 557               pad=i;
 558               time=j;
 559               dig = fDigits->GetDigitFast(time,pad);
 560               if(dig <= fParam->GetZeroSup()) continue;
 561               if(dig >= AliL3Transform::GetADCSat())
 562                 dig = AliL3Transform::GetADCSat();
 563
 564               //Check for single timebins, and remove them because they are noise for sure.
 565               if(i>0 && i<fDigits->GetNCols()-1 && j>0 && j<fDigits->GetNRows()-1)
 566                 if(fDigits->GetDigitFast(time,pad-1)<=fParam->GetZeroSup() &&
 567                    fDigits->GetDigitFast(time-1,pad)<=fParam->GetZeroSup() &&
 568                    fDigits->GetDigitFast(time+1,pad)<=fParam->GetZeroSup() &&
 569                    fDigits->GetDigitFast(time,pad+1)<=fParam->GetZeroSup())
 570                   continue;
 571
 572               //Boundaries:
 573               if(i==0) //pad ==0
 574                 {
 575                   if(j < fDigits->GetNRows()-1 && j > 0)
 576                     {
 577                       if(fDigits->GetDigitFast(time-1,pad)<=fParam->GetZeroSup() &&
 578                          fDigits->GetDigitFast(time+1,pad)<=fParam->GetZeroSup() &&
 579                          fDigits->GetDigitFast(time,pad+1)<=fParam->GetZeroSup())
 580                         continue;
 581                     }
 582                   else if(j > 0)
 583                     {
 584                       if(fDigits->GetDigitFast(time-1,pad)<=fParam->GetZeroSup() &&
 585                          fDigits->GetDigitFast(time,pad+1)<=fParam->GetZeroSup())
 586                         continue;
 587                     }
 588                 }
 589               if(j==0)
 590                 {
 591                   if(i < fDigits->GetNCols()-1 && i > 0)
 592                     {
 593                       if(fDigits->GetDigitFast(time,pad-1)<=fParam->GetZeroSup() &&
 594                          fDigits->GetDigitFast(time,pad+1)<=fParam->GetZeroSup() &&
 595                          fDigits->GetDigitFast(time+1,pad)<=fParam->GetZeroSup())
 596                         continue;
 597                     }
 598                   else if(i > 0)
 599                     {
 600                       if(fDigits->GetDigitFast(time,pad-1)<=fParam->GetZeroSup() &&
 601                          fDigits->GetDigitFast(time+1,pad)<=fParam->GetZeroSup())
 602                         continue;
 603                     }
 604                 }
 605
 606               if(i == fDigits->GetNCols()-1)
 607                 {
 608
 609                   if(j>0 && j<fDigits->GetNRows()-1)
 610                     {
 611                       if(fDigits->GetDigitFast(time-1,pad)<=fParam->GetZeroSup() &&
 612                          fDigits->GetDigitFast(time+1,pad)<=fParam->GetZeroSup() &&
 613                          fDigits->GetDigitFast(time,pad-1)<=fParam->GetZeroSup())
 614                         continue;
 615                     }
 616                   else if(j==0 && j<fDigits->GetNRows()-1)
 617                     {
 618                       if(fDigits->GetDigitFast(time+1,pad)<=fParam->GetZeroSup() &&
 619                          fDigits->GetDigitFast(time,pad-1)<=fParam->GetZeroSup())
 620                         continue;
 621                     }
 622                   else
 623                     {
 624                       if(fDigits->GetDigitFast(time-1,pad)<=fParam->GetZeroSup() &&
 625                          fDigits->GetDigitFast(time,pad-1)<=fParam->GetZeroSup())
 626                         continue;
 627                     }
 628                 }
 629
 630               if(j==fDigits->GetNRows()-1)
 631                 {
 632                   if(i>0 && i<fDigits->GetNCols()-1)
 633                     {
 634                       if(fDigits->GetDigitFast(time,pad-1)<=fParam->GetZeroSup() &&
 635                          fDigits->GetDigitFast(time,pad+1)<=fParam->GetZeroSup() &&
 636                          fDigits->GetDigitFast(time-1,pad)<=fParam->GetZeroSup())
 637                         continue;
 638                     }
 639                   else if(i==0 && fDigits->GetNCols()-1)
 640                     {
 641                       if(fDigits->GetDigitFast(time,pad+1)<=fParam->GetZeroSup() &&
 642                          fDigits->GetDigitFast(time-1,pad)<=fParam->GetZeroSup())
 643                         continue;
 644                     }
 645                   else
 646                     {
 647                       if(fDigits->GetDigitFast(time,pad-1)<=fParam->GetZeroSup() &&
 648                          fDigits->GetDigitFast(time-1,pad)<=fParam->GetZeroSup())
 649                         continue;
 650                     }
 651                 }
 652
 653
 654               AliL3Transform::Raw2Local(xyz,sector,row,pad,time);
 655               if(fParam->GetPadRowRadii(sector,row)<230./250.*fabs(xyz[2]))
 656                 continue;
 657
 658               if(localcount >= ndigits[lrow])
 659                 LOG(AliL3Log::kFatal,"AliL3FileHandler::AliDigits2Binary","Memory")
 660                   <<AliL3Log::kDec<<"Mismatch: localcount "<<localcount<<" ndigits "
 661                   <<ndigits[lrow]<<ENDLOG;
 662
 663               tempPt->fDigitData[localcount].fCharge=dig;
 664               tempPt->fDigitData[localcount].fPad=pad;
 665               tempPt->fDigitData[localcount].fTime=time;
 666 #ifdef do_mc
 667               tempPt->fDigitData[localcount].fTrackID[0] = fDigits->GetTrackIDFast(time,pad,0)-2;
 668               tempPt->fDigitData[localcount].fTrackID[1] = fDigits->GetTrackIDFast(time,pad,1)-2;
 669               tempPt->fDigitData[localcount].fTrackID[2] = fDigits->GetTrackIDFast(time,pad,2)-2;
 670 #endif
 671               localcount++;
 672             }
 673         }
 674
 675       Byte_t *tmp = (Byte_t*)tempPt;
 676       Int_t size = sizeof(AliL3DigitRowData)
 677         + ndigits[lrow]*sizeof(AliL3DigitData);
 678       tmp += size;
 679       tempPt = (AliL3DigitRowData*)tmp;
 680 #ifdef ASVVERSION
 681       fLastIndex=n;
 682 #endif
 683     }
 684 #ifdef ASVVERSION
 685   fLastIndex++;
 686 #endif
 687   return data;
 688 }
 689
 690 Bool_t AliL3FileHandler::GetDigitsTree(Int_t event)
 691 {
 692   //Connects to the TPC digit tree in the AliROOT file.
 693
 694   fInAli->cd();
 695   Char_t dname[100];
 696   sprintf(dname,"TreeD_%s_%d",AliL3Transform::GetParamName(),event);
 697   fDigitsTree = (TTree*)fInAli->Get(dname);
 698   if(!fDigitsTree)
 699     {
 700       LOG(AliL3Log::kError,"AliL3FileHandler::GetDigitsTree","Digits Tree")
 701         <<AliL3Log::kHex<<"Error getting digitstree "<<(Int_t)fDigitsTree<<ENDLOG;
 702       return kFALSE;
 703     }
 704   fDigitsTree->GetBranch("Segment")->SetAddress(&fDigits);
 705   return kTRUE;
 706 }
 707
 708 void AliL3FileHandler::AliDigits2RootFile(AliL3DigitRowData *rowPt,Char_t *new_digitsfile)
 709 {
 710   //Write the data stored in rowPt, into a new AliROOT file.
 711   //The data is stored in the AliROOT format
 712   //This is specially a nice thing if you have modified data, and wants to run it
 713   //through the offline reconstruction chain.
 714   //The arguments is a pointer to the data, and the name of the new AliROOT file.
 715   //Remember to pass the original AliROOT file (the one that contains the original
 716   //simulated data) to this object, in order to retrieve the MC id's of the digits.
 717
 718   if(!fInAli)
 719     {
 720       printf("AliL3FileHandler::AliDigits2RootFile : No rootfile\n");
 721       return;
 722     }
 723   if(!fParam)
 724     {
 725       printf("AliL3FileHandler::AliDigits2RootFile : No parameter object. Run on rootfile\n");
 726       return;
 727     }
 728
 729   //Get the original digitstree:
 730   Char_t dname[100];
 731   sprintf(dname,"TreeD_%s_0",AliL3Transform::GetParamName());
 732
 733   fInAli->cd();
 734   AliTPCDigitsArray *old_array = new AliTPCDigitsArray();
 735   old_array->Setup(fParam);
 736   old_array->SetClass("AliSimDigits");
 737
 738   Bool_t ok = old_array->ConnectTree(dname);
 739   if(!ok)
 740     {
 741       printf("AliL3FileHandler::AliDigits2RootFile : No digits tree object\n");
 742       return;
 743     }
 744
 745   Bool_t create=kFALSE;
 746   TFile *digFile;
 747
 748   if(gSystem->AccessPathName(new_digitsfile))
 749     {
 750       cout<<"AliL3FileHandler::AliDigits2RootFile : Creating new file :"<<new_digitsfile<<endl;
 751       create = kTRUE;
 752       digFile = TFile::Open(new_digitsfile,"RECREATE");
 753       fParam->Write(fParam->GetTitle());
 754     }
 755   else
 756     {
 757       create = kFALSE;
 758       digFile = TFile::Open(new_digitsfile,"UPDATE");
 759
 760     }
 761   if(!digFile->IsOpen())
 762     {
 763       LOG(AliL3Log::kError,"AliL3FileHandler::AliDigits2RootFile","Rootfile")
 764         <<"Error opening rootfile "<<new_digitsfile<<ENDLOG;
 765       return;
 766     }
 767
 768   digFile->cd();
 769
 770   //setup a new one, or connect it to the existing one:
 771   AliTPCDigitsArray *arr = new AliTPCDigitsArray();
 772   arr->SetClass("AliSimDigits");
 773   arr->Setup(fParam);
 774   if(create)
 775     arr->MakeTree();
 776   else
 777     {
 778       Bool_t ok = arr->ConnectTree(dname);
 779       if(!ok)
 780         {
 781           printf("AliL3FileHandler::AliDigits2RootFile : No digits tree object in existing file\n");
 782           return;
 783         }
 784     }
 785   Int_t digcounter=0,trackID[3];
 786
 787   for(Int_t i=fRowMin; i<=fRowMax; i++)
 788     {
 789
 790       if((Int_t)rowPt->fRow != i)
 791         cerr<<"AliL3FileHandler::AliDigits2RootFile : Mismatching row numbering!!!"<<endl;
 792
 793       Int_t sector,row;
 794       AliL3Transform::Slice2Sector(fSlice,i,sector,row);
 795
 796       AliSimDigits *old_dig = (AliSimDigits*)old_array->LoadRow(sector,row);
 797       AliSimDigits * dig = (AliSimDigits*)arr->CreateRow(sector,row);
 798       old_dig->ExpandBuffer();
 799       old_dig->ExpandTrackBuffer();
 800       dig->ExpandBuffer();
 801       dig->ExpandTrackBuffer();
 802
 803       if(!old_dig)
 804         printf("AliL3FileHandler::AliDigits2RootFile : No padrow %d %d\n",sector,row);
 805
 806       AliL3DigitData *digPt = rowPt->fDigitData;
 807       digcounter=0;
 808       for(UInt_t j=0; j<rowPt->fNDigit; j++)
 809         {
 810           Short_t charge = (Short_t)digPt[j].fCharge;
 811           Int_t pad = (Int_t)digPt[j].fPad;
 812           Int_t time = (Int_t)digPt[j].fTime;
 813
 814           if(charge == 0) //Only write the digits that has not been removed
 815             {
 816               cerr<<"AliL3FileHandler::AliDigits2RootFile : Zero charge!!! "<<endl;
 817               continue;
 818             }
 819
 820           digcounter++;
 821
 822           //Tricks to get and set the correct track id's.
 823           for(Int_t t=0; t<3; t++)
 824             {
 825               Int_t label = old_dig->GetTrackIDFast(time,pad,t);
 826               if(label > 1)
 827                 trackID[t] = label - 2;
 828               else if(label==0)
 829                 trackID[t] = -2;
 830               else
 831                 trackID[t] = -1;
 832             }
 833
 834           dig->SetDigitFast(charge,time,pad);
 835
 836           for(Int_t t=0; t<3; t++)
 837             ((AliSimDigits*)dig)->SetTrackIDFast(trackID[t],time,pad,t);
 838
 839         }
 840       //cout<<"Wrote "<<digcounter<<" on row "<<i<<endl;
 841       UpdateRowPointer(rowPt);
 842       arr->StoreRow(sector,row);
 843       arr->ClearRow(sector,row);
 844       old_array->ClearRow(sector,row);
 845     }
 846   digFile->cd();
 847   char treeName[100];
 848   sprintf(treeName,"TreeD_%s_0",fParam->GetTitle());
 849
 850   arr->GetTree()->SetName(treeName);
 851   arr->GetTree()->AutoSave();
 852   delete arr;
 853   delete old_array;
 854   digFile->Close();
 855
 856 }
 857
 858 ///////////////////////////////////////// Point IO
 859 Bool_t AliL3FileHandler::AliPoints2Binary(){
 860   Bool_t out = kTRUE;
 861   UInt_t npoint;
 862   AliL3SpacePointData *data = AliPoints2Memory(npoint);
 863   out = Memory2Binary(npoint,data);
 864   Free();
 865   return out;
 866 }
 867
 868 AliL3SpacePointData * AliL3FileHandler::AliPoints2Memory(UInt_t & npoint){
 869   AliL3SpacePointData *data = 0;
 870   npoint=0;
 871   if(!fInAli){
 872     LOG(AliL3Log::kWarning,"AliL3FileHandler::AliPoints2Memory","File")
 873     <<"No Input avalible: no object TFile"<<ENDLOG;
 874     return 0;
 875   }
 876   if(!fInAli->IsOpen()){
 877     LOG(AliL3Log::kWarning,"AliL3FileHandler::AliPoints2Memory","File")
 878     <<"No Input avalible: TFile not opend"<<ENDLOG;
 879     return 0;
 880   }
 881
 882   TDirectory *savedir = gDirectory;
 883   fInAli->cd();
 884
 885   Char_t cname[100];
 886   Int_t eventn = 0;
 887   sprintf(cname,"TreeC_TPC_%d",eventn);
 888   AliTPCClustersArray carray;
 889   carray.Setup(fParam);
 890   carray.SetClusterType("AliTPCcluster");
 891   Bool_t clusterok = carray.ConnectTree(cname);
 892   if(!clusterok) return 0;
 893
 894   AliTPCClustersRow ** clusterrow =
 895                new AliTPCClustersRow*[ (int)carray.GetTree()->GetEntries()];
 896   Int_t *rows = new int[ (int)carray.GetTree()->GetEntries()];
 897   Int_t *sects = new int[  (int)carray.GetTree()->GetEntries()];
 898   Int_t sum=0;
 899
 900   Int_t lslice,lrow;
 901   for(Int_t i=0; i<carray.GetTree()->GetEntries(); i++){
 902     AliSegmentID *s = carray.LoadEntry(i);
 903     Int_t sector,row;
 904     fParam->AdjustSectorRow(s->GetID(),sector,row);
 905     rows[i] = row;
 906     sects[i] = sector;
 907     clusterrow[i] = 0;
 908     AliL3Transform::Sector2Slice(lslice,lrow,sector,row);
 909     if(fSlice != lslice || lrow<fRowMin || lrow>fRowMax) continue;
 910     clusterrow[i] = carray.GetRow(sector,row);
 911     if(clusterrow[i])
 912       sum+=clusterrow[i]->GetArray()->GetEntriesFast();
 913   }
 914   UInt_t size = sum*sizeof(AliL3SpacePointData);
 915
 916   LOG(AliL3Log::kDebug,"AliL3FileHandler::AliPoints2Memory","File")
 917   <<AliL3Log::kDec<<"Found "<<sum<<" SpacePoints"<<ENDLOG;
 918
 919   data = (AliL3SpacePointData *) Allocate(size);
 920   npoint = sum;
 921   UInt_t n=0;
 922   for(Int_t i=0; i<carray.GetTree()->GetEntries(); i++){
 923     if(!clusterrow[i]) continue;
 924     Int_t row = rows[i];
 925     Int_t sector = sects[i];
 926     AliL3Transform::Sector2Slice(lslice,lrow,sector,row);
 927     Int_t entries_in_row = clusterrow[i]->GetArray()->GetEntriesFast();
 928     for(Int_t j = 0;j<entries_in_row;j++){
 929       AliTPCcluster *c = (AliTPCcluster*)(*clusterrow[i])[j];
 930       data[n].fZ = c->GetZ();
 931       data[n].fY = c->GetY();
 932       data[n].fX = fParam->GetPadRowRadii(sector,row);
 933       data[n].fCharge = (UInt_t)c->GetQ();
 934       data[n].fID = n+((fSlice&0x7f)<<25)+((fPatch&0x7)<<22);//uli
 935       data[n].fPadRow = lrow;
 936       data[n].fXYErr = sqrt(c->GetSigmaY2());
 937       data[n].fZErr = sqrt(c->GetSigmaZ2());
 938       if(fMC) fprintf(fMC,"%d %d\n",data[n].fID,c->GetLabel(0));
 939       n++;
 940     }
 941   }
 942   for(Int_t i=0;i<carray.GetTree()->GetEntries();i++){
 943     Int_t row = rows[i];
 944     Int_t sector = sects[i];
 945     if(carray.GetRow(sector,row))
 946       carray.ClearRow(sector,row);
 947   }
 948
 949   delete [] clusterrow;
 950   delete [] rows;
 951   delete [] sects;
 952   savedir->cd();
 953
 954   return data;
 955 }
 956