TPC/TPCsim/AliTPCDigitizer.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2000, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 /*
  19   Class for creating of the sumable digits and digits from MC data
  20   //
  21   The input :  ideal signals (Hits->Diffusion->Attachment -Ideal signal)
  22   The output:  "raw digits"
  23
  24   Effect implemented:
  25   1. Pad by pad gain map
  26   2. Noise map
  27   3. The dead channels identified  - zerro noise for corresponding pads
  28      In this case the outpu equal zerro
  29
  30 */
  31
  32
  33
  34
  35 #include <stdlib.h>
  36 #include <TTree.h>
  37 #include <TObjArray.h>
  38 #include <TFile.h>
  39 #include <TDirectory.h>
  40 #include <Riostream.h>
  41 #include <TParameter.h>
  42
  43 #include "AliTPCDigitizer.h"
  44
  45 #include "AliTPC.h"
  46 #include "AliTPCParam.h"
  47 #include "AliTPCParamSR.h"
  48 #include "AliRun.h"
  49 #include "AliLoader.h"
  50 #include "AliPDG.h"
  51 #include "AliDigitizationInput.h"
  52 #include "AliSimDigits.h"
  53 #include "AliLog.h"
  54
  55 #include "AliTPCcalibDB.h"
  56 #include "AliTPCCalPad.h"
  57 #include "AliTPCCalROC.h"
  58 #include "TTreeStream.h"
  59 #include "AliTPCReconstructor.h"
  60 #include <TGraphErrors.h>
  61
  62 using std::cout;
  63 using std::cerr;
  64 using std::endl;
  65 ClassImp(AliTPCDigitizer)
  66
  67 //___________________________________________
  68   AliTPCDigitizer::AliTPCDigitizer() :AliDigitizer(),fDebug(0), fDebugStreamer(0)
  69 {
  70   //
  71 // Default ctor - don't use it
  72 //
  73
  74 }
  75
  76 //___________________________________________
  77 AliTPCDigitizer::AliTPCDigitizer(AliDigitizationInput* digInput)
  78   :AliDigitizer(digInput),fDebug(0), fDebugStreamer(0)
  79 {
  80   //
  81 // ctor which should be used
  82 //
  83   AliDebug(2,"(AliDigitizationInput* digInput) was processed");
  84   if (AliTPCReconstructor::StreamLevel()>0)  fDebugStreamer = new TTreeSRedirector("TPCDigitDebug.root","recreate");
  85
  86 }
  87
  88 //------------------------------------------------------------------------
  89 AliTPCDigitizer::~AliTPCDigitizer()
  90 {
  91 // Destructor
  92   if (fDebugStreamer) delete fDebugStreamer;
  93 }
  94
  95
  96
  97 //------------------------------------------------------------------------
  98 Bool_t AliTPCDigitizer::Init()
  99 {
 100 // Initialization
 101
 102  return kTRUE;
 103 }
 104
 105
 106 //------------------------------------------------------------------------
 107 void AliTPCDigitizer::Digitize(Option_t* option)
 108 {
 109   //DigitizeFast(option);
 110   DigitizeWithTailAndCrossTalk(option);
 111
 112 }
 113 //------------------------------------------------------------------------
 114 void AliTPCDigitizer::DigitizeFast(Option_t* option)
 115 {
 116
 117   // merge input tree's with summable digits
 118   //output stored in TreeTPCD
 119   char s[100];
 120   char ss[100];
 121   TString optionString = option;
 122   if (!strcmp(optionString.Data(),"deb")) {
 123     cout<<"AliTPCDigitizer:::DigitizeFast called with option deb "<<endl;
 124     fDebug = 3;
 125   }
 126   //get detector and geometry
 127
 128
 129   AliRunLoader *rl, *orl;
 130   AliLoader *gime, *ogime;
 131
 132   if (gAlice == 0x0)
 133    {
 134      Warning("DigitizeFast","gAlice is NULL. Loading from input 0");
 135      rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(0));
 136      if (rl == 0x0)
 137       {
 138         Error("DigitizeFast","Can not find Run Loader for input 0. Can not proceed.");
 139         return;
 140       }
 141      rl->LoadgAlice();
 142      rl->GetAliRun();
 143    }
 144   AliTPC *pTPC  = (AliTPC *) gAlice->GetModule("TPC");
 145   AliTPCParam * param = pTPC->GetParam();
 146
 147   //sprintf(s,param->GetTitle());
 148   snprintf(s,100,"%s",param->GetTitle());
 149   //sprintf(ss,"75x40_100x60");
 150   snprintf(ss,100,"75x40_100x60");
 151   if(strcmp(s,ss)==0){
 152     printf("2 pad-length geom hits with 3 pad-lenght geom digits...\n");
 153     delete param;
 154     param=new AliTPCParamSR();
 155   }
 156   else{
 157     //sprintf(ss,"75x40_100x60_150x60");
 158     snprintf(ss,100,"75x40_100x60_150x60");
 159    if(strcmp(s,ss)!=0) {
 160      printf("No TPC parameters found...\n");
 161      exit(2);
 162    }
 163   }
 164
 165   pTPC->GenerNoise(500000); //create table with noise
 166   //
 167   Int_t nInputs = fDigInput->GetNinputs();
 168   Int_t * masks = new Int_t[nInputs];
 169   for (Int_t i=0; i<nInputs;i++)
 170     masks[i]= fDigInput->GetMask(i);
 171   Short_t **pdig= new Short_t*[nInputs];   //pointers to the expanded digits array
 172   Int_t **ptr=  new Int_t*[nInputs];       //pointers to the expanded tracks array
 173   Bool_t *active=  new Bool_t[nInputs];    //flag for active input segments
 174   Char_t phname[100];
 175
 176   //create digits array for given sectors
 177   // make indexes
 178   //
 179   //create branch's in TPC treeD
 180   orl = AliRunLoader::GetRunLoader(fDigInput->GetOutputFolderName());
 181   ogime = orl->GetLoader("TPCLoader");
 182   TTree * tree  = ogime->TreeD();
 183   AliSimDigits * digrow = new AliSimDigits;
 184
 185   if (tree == 0x0)
 186    {
 187      ogime->MakeTree("D");
 188      tree  = ogime->TreeD();
 189    }
 190   tree->Branch("Segment","AliSimDigits",&digrow);
 191   //
 192   AliSimDigits ** digarr = new AliSimDigits*[nInputs];
 193   for (Int_t i1=0;i1<nInputs; i1++)
 194     {
 195       digarr[i1]=0;
 196      //    intree[i1]
 197       rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i1));
 198       gime = rl->GetLoader("TPCLoader");
 199       gime->LoadSDigits("read");
 200       TTree * treear =  gime->TreeS();
 201
 202       if (!treear)
 203        {
 204         cerr<<"AliTPCDigitizer: Input tree with SDigits not found in"
 205             <<" input "<< i1<<endl;
 206         for (Int_t i2=0;i2<i1+1; i2++){
 207
 208           if(digarr[i2])  delete digarr[i2];
 209         }
 210         delete [] digarr;
 211         delete [] active;
 212         delete []masks;
 213         delete []pdig;
 214         delete []ptr;
 215         return;
 216        }
 217
 218       //sprintf(phname,"lhcphase%d",i1);
 219       snprintf(phname,100,"lhcphase%d",i1);
 220       TParameter<float> *ph = (TParameter<float>*)treear->GetUserInfo()
 221                                ->FindObject("lhcphase0");
 222       if(!ph){
 223         cerr<<"AliTPCDigitizer: LHC phase  not found in"
 224             <<" input "<< i1<<endl;
 225         for (Int_t i2=0;i2<i1+1; i2++){
 226           if(digarr[i2])  delete digarr[i2];
 227         }
 228         delete [] digarr;
 229         delete [] active;
 230         delete []masks;
 231         delete []pdig;
 232         delete []ptr;
 233         return;
 234       }
 235       tree->GetUserInfo()->Add(new TParameter<float>(phname,ph->GetVal()));
 236               //
 237       if (treear->GetIndex()==0)
 238         treear->BuildIndex("fSegmentID","fSegmentID");
 239       treear->GetBranch("Segment")->SetAddress(&digarr[i1]);
 240     }
 241
 242
 243
 244
 245   //
 246
 247   param->SetZeroSup(2);
 248
 249   Int_t zerosup = param->GetZeroSup();
 250   AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetDedxGainFactor();
 251   AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
 252   //
 253   //Loop over segments of the TPC
 254
 255   for (Int_t segmentID=0; segmentID<param->GetNRowsTotal(); segmentID++)
 256    {
 257     Int_t sector, padRow;
 258     if (!param->AdjustSectorRow(segmentID,sector,padRow))
 259      {
 260       cerr<<"AliTPC warning: invalid segment ID ! "<<segmentID<<endl;
 261       continue;
 262      }
 263     AliTPCCalROC * gainROC = gainTPC->GetCalROC(sector);  // pad gains per given sector
 264     AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(sector);  // noise per given sector
 265     digrow->SetID(segmentID);
 266
 267     Int_t nTimeBins = 0;
 268     Int_t nPads = 0;
 269
 270     Bool_t digitize = kFALSE;
 271     for (Int_t i=0;i<nInputs; i++)
 272      {
 273
 274       rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i));
 275       gime = rl->GetLoader("TPCLoader");
 276
 277       if (gime->TreeS()->GetEntryWithIndex(segmentID,segmentID) >= 0) {
 278         digarr[i]->ExpandBuffer();
 279         digarr[i]->ExpandTrackBuffer();
 280         nTimeBins = digarr[i]->GetNRows();
 281         nPads = digarr[i]->GetNCols();
 282         active[i] = kTRUE;
 283         if (!GetRegionOfInterest() || (i == 0)) digitize = kTRUE;
 284       } else {
 285         active[i] = kFALSE;
 286       }
 287       if (GetRegionOfInterest() && !digitize) break;
 288      }
 289     if (!digitize) continue;
 290
 291     digrow->Allocate(nTimeBins,nPads);
 292     digrow->AllocateTrack(3);
 293
 294     Float_t q=0;
 295     Int_t label[1000]; //stack for 300 events
 296     Int_t labptr = 0;
 297
 298     Int_t nElems = nTimeBins*nPads;
 299
 300     for (Int_t i=0;i<nInputs; i++)
 301      if (active[i]) {
 302        pdig[i] = digarr[i]->GetDigits();
 303        ptr[i]  = digarr[i]->GetTracks();
 304       }
 305
 306     Short_t *pdig1= digrow->GetDigits();
 307     Int_t   *ptr1= digrow->GetTracks() ;
 308
 309
 310
 311     for (Int_t elem=0;elem<nElems; elem++)
 312      {
 313
 314        q=0;
 315        labptr=0;
 316        // looop over digits
 317         for (Int_t i=0;i<nInputs; i++) if (active[i])
 318          {
 319           //          q  += digarr[i]->GetDigitFast(rows,col);
 320             q  += *(pdig[i]);
 321
 322            for (Int_t tr=0;tr<3;tr++)
 323             {
 324              //             Int_t lab = digarr[i]->GetTrackIDFast(rows,col,tr);
 325              Int_t lab = ptr[i][tr*nElems];
 326              if ( (lab > 1) && *(pdig[i])>zerosup)
 327               {
 328                 label[labptr]=lab+masks[i];
 329                 labptr++;
 330               }
 331             }
 332            pdig[i]++;
 333            ptr[i]++;
 334          }
 335         q/=16.;  //conversion factor
 336         Float_t gain = gainROC->GetValue(padRow,elem/nTimeBins);  // get gain for given - pad-row pad
 337         //if (gain<0.5){
 338           //printf("problem\n");
 339         //}
 340         q*= gain;
 341         Float_t noisePad = noiseROC->GetValue(padRow,elem/nTimeBins);
 342         //       Float_t noise  = gRandom->Gaus(0,param->GetNoise()*param->GetNoiseNormFac());
 343         Float_t noise  = pTPC->GetNoise();
 344         q+=noise*noisePad;
 345         q=TMath::Nint(q);
 346         if (q > zerosup)
 347          {
 348           if(q >= param->GetADCSat()) q = (Short_t)(param->GetADCSat() - 1);
 349           //digrow->SetDigitFast((Short_t)q,rows,col);
 350           *pdig1 =Short_t(q);
 351           for (Int_t tr=0;tr<3;tr++)
 352            {
 353             if (tr<labptr)
 354              ptr1[tr*nElems] = label[tr];
 355            }
 356           }
 357         pdig1++;
 358         ptr1++;
 359      }
 360     //
 361     //  glitch filter
 362     //
 363     digrow->GlitchFilter();
 364     //
 365     digrow->CompresBuffer(1,zerosup);
 366     digrow->CompresTrackBuffer(1);
 367     tree->Fill();
 368     if (fDebug>0) cerr<<sector<<"\t"<<padRow<<"\n";
 369    } //for (Int_t n=0; n<param->GetNRowsTotal(); n++)
 370
 371
 372   orl = AliRunLoader::GetRunLoader(fDigInput->GetOutputFolderName());
 373   ogime = orl->GetLoader("TPCLoader");
 374   ogime->WriteDigits("OVERWRITE");
 375
 376   //fDigInput->GetTreeDTPC()->Write(0,TObject::kOverwrite);
 377
 378   delete digrow;
 379   for (Int_t i1=0;i1<nInputs; i1++) delete digarr[i1];
 380   delete []masks;
 381   delete []pdig;
 382   delete []ptr;
 383   delete []active;
 384   delete []digarr;
 385 }
 386
 387
 388
 389 //------------------------------------------------------------------------
 390 void AliTPCDigitizer::DigitizeSave(Option_t* option)
 391 {
 392   //
 393   // Merge input tree's with summable digits
 394   // Output digits stored in TreeTPCD
 395   //
 396   // Not active for long time.
 397   // Before adding modification (for ion tail calucation and for the crorsstalk) it should be
 398   //  checked one by one with currenlty used AliTPCDigitizer::DigitizeFast
 399   //
 400   TString optionString = option;
 401   if (!strcmp(optionString.Data(),"deb")) {
 402     cout<<"AliTPCDigitizer::Digitize: called with option deb "<<endl;
 403     fDebug = 3;
 404   }
 405   //get detector and geometry
 406   AliRunLoader *rl, *orl;
 407   AliLoader *gime, *ogime;
 408
 409
 410   orl = AliRunLoader::GetRunLoader(fDigInput->GetOutputFolderName());
 411   ogime = orl->GetLoader("TPCLoader");
 412
 413   rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(0));
 414   //gime = rl->GetLoader("TPCLoader");
 415   rl->GetLoader("TPCLoader");
 416   rl->LoadgAlice();
 417   AliRun* alirun = rl->GetAliRun();
 418
 419   AliTPC *pTPC  = (AliTPC *) alirun->GetModule("TPC");
 420   AliTPCParam * param = pTPC->GetParam();
 421   pTPC->GenerNoise(500000); //create teble with noise
 422   printf("noise %f \n",  param->GetNoise()*param->GetNoiseNormFac());
 423   //
 424   Int_t nInputs = fDigInput->GetNinputs();
 425   // stupid protection...
 426   if (nInputs <= 0) return;
 427   //
 428   Int_t * masks = new Int_t[nInputs];
 429   for (Int_t i=0; i<nInputs;i++)
 430     masks[i]= fDigInput->GetMask(i);
 431
 432   AliSimDigits ** digarr = new AliSimDigits*[nInputs];
 433   for(Int_t ii=0;ii<nInputs;ii++) digarr[ii]=0;
 434
 435   for (Int_t i1=0;i1<nInputs; i1++)
 436    {
 437      //digarr[i1]=0;
 438     //    intree[i1]
 439     rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i1));
 440     gime = rl->GetLoader("TPCLoader");
 441
 442     TTree * treear =  gime->TreeS();
 443     //
 444     if (!treear) {
 445       cerr<<" TPC -  not existing input = \n"<<i1<<" ";
 446       delete [] masks;
 447       for(Int_t i=0; i<nInputs; i++) delete digarr[i];
 448       delete [] digarr;
 449       return;
 450     }
 451     //
 452     TBranch * br = treear->GetBranch("fSegmentID");
 453     if (br) br->GetFile()->cd();
 454     treear->GetBranch("Segment")->SetAddress(&digarr[i1]);
 455    }
 456
 457   rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(0));
 458   gime = rl->GetLoader("TPCLoader");
 459   Stat_t nentries = gime->TreeS()->GetEntries();
 460
 461
 462   //create branch's in TPC treeD
 463   AliSimDigits * digrow = new AliSimDigits;
 464   TTree * tree  = ogime->TreeD();
 465
 466   tree->Branch("Segment","AliSimDigits",&digrow);
 467   param->SetZeroSup(2);
 468
 469   Int_t zerosup = param->GetZeroSup();
 470   //Loop over segments of the TPC
 471
 472   AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetDedxGainFactor();
 473   AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
 474   for (Int_t n=0; n<nentries; n++) {
 475     rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(0));
 476     gime = rl->GetLoader("TPCLoader");
 477     gime->TreeS()->GetEvent(n);
 478
 479     digarr[0]->ExpandBuffer();
 480     digarr[0]->ExpandTrackBuffer();
 481
 482
 483     for (Int_t i=1;i<nInputs; i++){
 484 //      fDigInput->GetInputTreeTPCS(i)->GetEntryWithIndex(digarr[0]->GetID(),digarr[0]->GetID());
 485       rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i));
 486       gime = rl->GetLoader("TPCLoader");
 487       gime->TreeS()->GetEntryWithIndex(digarr[0]->GetID(),digarr[0]->GetID());
 488       digarr[i]->ExpandBuffer();
 489       digarr[i]->ExpandTrackBuffer();
 490       if ((digarr[0]->GetID()-digarr[i]->GetID())>0)
 491        printf("problem\n");
 492
 493     }
 494
 495     Int_t sector, padRow;
 496     if (!param->AdjustSectorRow(digarr[0]->GetID(),sector,padRow)) {
 497       cerr<<"AliTPC warning: invalid segment ID ! "<<digarr[0]->GetID()<<endl;
 498       continue;
 499     }
 500
 501     AliTPCCalROC * gainROC = gainTPC->GetCalROC(sector);  // pad gains per given sector
 502     AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(sector);  // noise per given sector
 503     digrow->SetID(digarr[0]->GetID());
 504
 505     Int_t nTimeBins = digarr[0]->GetNRows();
 506     Int_t nPads = digarr[0]->GetNCols();
 507     digrow->Allocate(nTimeBins,nPads);
 508     digrow->AllocateTrack(3);
 509
 510     Float_t q=0;
 511     Int_t label[1000]; //stack for 300 events
 512     Int_t labptr = 0;
 513
 514
 515
 516     for (Int_t iTimeBin=0;iTimeBin<nTimeBins; iTimeBin++){   // iTimeBin
 517       for (Int_t iPad=0;iPad<nPads; iPad++){    // pad
 518
 519        q=0;
 520        labptr=0;
 521        // looop over digits
 522         for (Int_t i=0;i<nInputs; i++){
 523          q  += digarr[i]->GetDigitFast(iTimeBin,iPad);
 524           //q  += *(pdig[i]);
 525
 526           for (Int_t tr=0;tr<3;tr++) {
 527            Int_t lab = digarr[i]->GetTrackIDFast(iTimeBin,iPad,tr);
 528            //Int_t lab = ptr[i][tr*nElems];
 529             if ( (lab > 1) ) {
 530               label[labptr]=lab+masks[i];
 531               labptr++;
 532             }
 533           }
 534          // pdig[i]++;
 535          //ptr[i]++;
 536
 537         }
 538        q/=16.;  //conversion factor
 539        //       Float_t noise  = gRandom->Gaus(0,param->GetNoise()*param->GetNoiseNormFac());
 540        Float_t gain = gainROC->GetValue(padRow,iPad);
 541        q*= gain;
 542        Float_t noisePad = noiseROC->GetValue(padRow, iPad);
 543
 544        Float_t noise  = pTPC->GetNoise();
 545        q+=noise*noisePad;
 546        //
 547        // here we can get digits from past and add signal
 548        //
 549        //
 550        //for (Int_t jTimeBin=0; jTimeBin<iTimeBin; jTimeBin++)
 551        //  q+=ionTail
 552        //
 553
 554         q=TMath::Nint(q);
 555         if (q > zerosup){
 556
 557          if(q >= param->GetADCSat()) q = (Short_t)(param->GetADCSat() - 1);
 558          digrow->SetDigitFast((Short_t)q,iTimeBin,iPad);
 559          // *pdig1 =Short_t(q);
 560          for (Int_t tr=0;tr<3;tr++){
 561            if (tr<labptr)
 562              ((AliSimDigits*)digrow)->SetTrackIDFast(label[tr],iTimeBin,iPad,tr);
 563            //ptr1[tr*nElems] = label[tr];
 564            //else
 565              //           ((AliSimDigits*)digrow)->SetTrackIDFast(-1,iTimeBin,iPad,tr);
 566            //  ptr1[tr*nElems] = 1;
 567          }
 568        }
 569        //pdig1++;
 570        //ptr1++;
 571     }
 572     }
 573
 574     digrow->CompresBuffer(1,zerosup);
 575     digrow->CompresTrackBuffer(1);
 576     tree->Fill();
 577     if (fDebug>0) cerr<<sector<<"\t"<<padRow<<"\n";
 578   }
 579 //  printf("end TPC merging - end -Tree %s\t%p\n",fDigInput->GetInputTreeH(0)->GetName(),fDigInput->GetInputTreeH(0)->GetListOfBranches()->At(3));
 580   //fDigInput->GetTreeDTPC()->Write(0,TObject::kOverwrite);
 581     ogime->WriteDigits("OVERWRITE");
 582
 583     for (Int_t i=1;i<nInputs; i++)
 584      {
 585       rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i));
 586       gime = rl->GetLoader("TPCLoader");
 587       gime->UnloadSDigits();
 588      }
 589     ogime->UnloadDigits();
 590
 591   delete digrow;
 592   for (Int_t i1=0;i1<nInputs; i1++) delete digarr[i1];
 593   delete [] masks;
 594   delete [] digarr;
 595 }
 596
 597
 598
 599
 600
 601
 602
 603 //------------------------------------------------------------------------
 604 void AliTPCDigitizer::DigitizeWithTailAndCrossTalk(Option_t* option)
 605 {
 606   // Modified version of the digitization function
 607   // Modification: adding the ion tail and crosstalk:
 608   //
 609   // pcstream used in order to visually inspect data
 610   //
 611   //
 612   // Crosstalk simulation:
 613   //   1.) Calculate per time bin mean charge (per pad)  within anode wire segment
 614   //   2.) Subsract for the clusters at given time bin fraction of (mean charge) normalized by add hoc constant
 615   //       AliTPCRecoParam::GetCrosstalkCorrection() (0 if not crosstalk, 1 if ideal crosstalk)
 616   //       for simplicity we are assuming that wire segents are related to pad-rows
 617   //       Wire segmentationn is obtatined from the
 618   //       AliTPCParam::GetWireSegment(Int_t sector, Int_t row); // to be implemented
 619   //       AliTPCParam::GetNPadsPerSegment(Int_t segmentID); // to be implemented
 620   //   3.) Substract form the signal contribution from the previous tracks - Ion tail in case speified in the AliTPCRecoParam
 621   //        AliTPCRecoParam::GetUseIonTailCorrection()
 622   //
 623   // Ion tail simulation:
 624   //    1.) Needs signal from pad+-1, taking signal from history
 625   // merge input tree's with summable digits
 626   // output stored in TreeTPCD
 627   //
 628   AliTPCcalibDB* const calib=AliTPCcalibDB::Instance();
 629   AliTPCRecoParam *recoParam = calib->GetRecoParam(0);
 630   AliDebug(1, Form(" recoParam->GetCrosstalkCorrection()  =   %f", recoParam->GetCrosstalkCorrection()));
 631   AliDebug(1,Form(" recoParam->GetUseIonTailCorrection() =  %d ", recoParam->GetUseIonTailCorrection()));
 632   Int_t nROCs = 72;
 633   char s[100];
 634   char ss[100];
 635   TString optionString = option;
 636   if (!strcmp(optionString.Data(),"deb")) {
 637     cout<<"AliTPCDigitizer:::DigitizeFast called with option deb "<<endl;
 638     fDebug = 3;
 639   }
 640
 641   // ======== get detector and geometry =======
 642
 643   AliRunLoader *rl, *orl;
 644   AliLoader *gime, *ogime;
 645
 646   if (gAlice == 0x0)
 647   {
 648     Warning("DigitizeFast","gAlice is NULL. Loading from input 0");
 649     rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(0));
 650     if (rl == 0x0)
 651     {
 652       Error("DigitizeFast","Can not find Run Loader for input 0. Can not proceed.");
 653       return;
 654     }
 655     rl->LoadgAlice();
 656     rl->GetAliRun();
 657   }
 658   AliTPC *pTPC  = (AliTPC *) gAlice->GetModule("TPC");
 659   AliTPCParam * param = pTPC->GetParam();
 660
 661   //sprintf(s,param->GetTitle());
 662   snprintf(s,100,"%s",param->GetTitle());
 663   //sprintf(ss,"75x40_100x60");
 664   snprintf(ss,100,"75x40_100x60");
 665   if(strcmp(s,ss)==0){
 666     printf("2 pad-length geom hits with 3 pad-lenght geom digits...\n");
 667     delete param;
 668     param=new AliTPCParamSR();
 669   } else {
 670     //sprintf(ss,"75x40_100x60_150x60");
 671     snprintf(ss,100,"75x40_100x60_150x60");
 672     if(strcmp(s,ss)!=0) {
 673       printf("No TPC parameters found...\n");
 674       exit(2);
 675     }
 676   }
 677
 678   pTPC->GenerNoise(500000); //create table with noise
 679   //
 680   Int_t nInputs = fDigInput->GetNinputs();
 681   Int_t * masks = new Int_t[nInputs];
 682   for (Int_t i=0; i<nInputs;i++)
 683     masks[i]= fDigInput->GetMask(i);
 684   Short_t **pdig= new Short_t*[nInputs];   //pointers to the expanded digits array
 685   Int_t **ptr=  new Int_t*[nInputs];       //pointers to the expanded tracks array
 686   Bool_t *active=  new Bool_t[nInputs];    //flag for active input segments
 687   Char_t phname[100];
 688
 689   //create digits array for given sectors
 690   // make indexes
 691   //
 692   //create branch's in TPC treeD
 693   orl = AliRunLoader::GetRunLoader(fDigInput->GetOutputFolderName());
 694   ogime = orl->GetLoader("TPCLoader");
 695   TTree * tree  = ogime->TreeD();
 696   AliSimDigits * digrow = new AliSimDigits;
 697
 698   if (tree == 0x0)
 699   {
 700     ogime->MakeTree("D");
 701     tree  = ogime->TreeD();
 702   }
 703   tree->Branch("Segment","AliSimDigits",&digrow);
 704   //
 705   AliSimDigits ** digarr = new AliSimDigits*[nInputs];
 706   for (Int_t i1=0;i1<nInputs; i1++)
 707   {
 708     digarr[i1]=0;
 709     //    intree[i1]
 710     rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i1));
 711     gime = rl->GetLoader("TPCLoader");
 712     gime->LoadSDigits("read");
 713     TTree * treear =  gime->TreeS();
 714
 715     if (!treear)
 716     {
 717       cerr<<"AliTPCDigitizer: Input tree with SDigits not found in"
 718         <<" input "<< i1<<endl;
 719       for (Int_t i2=0;i2<i1+1; i2++){
 720
 721         if(digarr[i2])  delete digarr[i2];
 722       }
 723       delete [] digarr;
 724       delete [] active;
 725       delete []masks;
 726       delete []pdig;
 727       delete []ptr;
 728       return;
 729     }
 730
 731     //sprintf(phname,"lhcphase%d",i1);
 732     snprintf(phname,100,"lhcphase%d",i1);
 733     TParameter<float> *ph = (TParameter<float>*)treear->GetUserInfo()
 734       ->FindObject("lhcphase0");
 735     if(!ph)
 736     {
 737       cerr<<"AliTPCDigitizer: LHC phase  not found in"
 738         <<" input "<< i1<<endl;
 739       for (Int_t i2=0;i2<i1+1; i2++){
 740         if(digarr[i2])  delete digarr[i2];
 741       }
 742       delete [] digarr;
 743       delete [] active;
 744       delete []masks;
 745       delete []pdig;
 746       delete []ptr;
 747       return;
 748     }
 749     tree->GetUserInfo()->Add(new TParameter<float>(phname,ph->GetVal()));
 750     //
 751     if (treear->GetIndex()==0)
 752       treear->BuildIndex("fSegmentID","fSegmentID");
 753     treear->GetBranch("Segment")->SetAddress(&digarr[i1]);
 754   }
 755
 756
 757
 758
 759   //
 760   // zero supp, take gain and noise map of TPC from OCDB
 761   param->SetZeroSup(2);
 762   Int_t zerosup = param->GetZeroSup();
 763   AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetDedxGainFactor();
 764   AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
 765
 766
 767   //
 768   // Cache the ion tail objects form OCDB
 769   //
 770
 771   TObjArray *ionTailArr = (TObjArray*)AliTPCcalibDB::Instance()->GetIonTailArray();
 772   if (!ionTailArr) {AliFatal("TPC - Missing IonTail OCDB object");}
 773  //  TObject *rocFactorIROC  = ionTailArr->FindObject("factorIROC");
 774 //   TObject *rocFactorOROC  = ionTailArr->FindObject("factorOROC");
 775 //   Float_t factorIROC      = (atof(rocFactorIROC->GetTitle()));
 776 //   Float_t factorOROC      = (atof(rocFactorOROC->GetTitle()));
 777   Int_t nIonTailBins =0;
 778   TObjArray timeResFunc(nROCs);
 779   for (Int_t isec = 0;isec<nROCs;isec++){        //loop overs sectors
 780     // Array of TGraphErrors for a given sector
 781     TGraphErrors ** graphRes   = new TGraphErrors *[20];
 782     Float_t * trfIndexArr    = new Float_t[20];
 783     for (Int_t icache=0; icache<20; icache++)
 784     {
 785       graphRes[icache]       = NULL;
 786       trfIndexArr[icache] = 0;
 787     }
 788     if (!AliTPCcalibDB::Instance()->GetTailcancelationGraphs(isec,graphRes,trfIndexArr)) continue;
 789
 790     // fill all TGraphErrors of trfs (time response function) of a given sector to a TObjArray
 791     TObjArray *timeResArr = new TObjArray(20);  timeResArr -> SetOwner(kTRUE);
 792     for (Int_t ires = 0;ires<20;ires++) timeResArr->AddAt(graphRes[ires],ires);
 793     timeResFunc.AddAt(timeResArr,isec); // Fill all trfs into a single TObjArray
 794     nIonTailBins = graphRes[3]->GetN();
 795   }
 796
 797   //
 798   // 1.) Make first loop to calculate mean amplitude per pad per segment for cross talk
 799   //
 800
 801   TObjArray   crossTalkSignalArray(nROCs);  // for 36 sectors
 802   TVectorD  * qTotSectorOld  = new TVectorD(nROCs);
 803   Float_t qTotTPC = 0.;
 804   Float_t qTotPerSector = 0.;
 805   Int_t nTimeBinsAll = 1100;
 806   Int_t nWireSegments=11;
 807   // 1.a) crorstalk matrix initialization
 808   for (Int_t sector=0; sector<nROCs; sector++){
 809     TMatrixD *pcrossTalkSignal = new TMatrixD(nWireSegments,nTimeBinsAll);
 810     for (Int_t imatrix = 0; imatrix<11; imatrix++)
 811       for (Int_t jmatrix = 0; jmatrix<nTimeBinsAll; jmatrix++){
 812         (*pcrossTalkSignal)[imatrix][jmatrix]=0.;
 813       }
 814     crossTalkSignalArray.AddAt(pcrossTalkSignal,sector);
 815   }
 816   //
 817   // main loop over rows of whole TPC
 818   for (Int_t globalRowID=0; globalRowID<param->GetNRowsTotal(); globalRowID++) {
 819     Int_t sector, padRow;
 820     if (!param->AdjustSectorRow(globalRowID,sector,padRow)) {
 821       cerr<<"AliTPC warning: invalid segment ID ! "<<globalRowID<<endl;
 822       continue;
 823     }
 824     // Calculate number of pads in a anode wire segment for normalization
 825     Int_t wireSegmentID    = param->GetWireSegment(sector,padRow);
 826     Float_t nPadsPerSegment = (Float_t)(param->GetNPadsPerSegment(wireSegmentID));
 827     // structure with mean signal per pad to be filled for each timebin in first loop (11 anodeWireSegment and 1100 timebin)
 828     TMatrixD &crossTalkSignal =  *((TMatrixD*)crossTalkSignalArray.At(sector));
 829     AliTPCCalROC * gainROC = gainTPC->GetCalROC(sector);  // pad gains per given sector
 830     //  digrow->SetID(globalRowID);
 831     Int_t nTimeBins = 0;
 832     Int_t nPads = 0;
 833     Bool_t digitize = kFALSE;
 834     for (Int_t i=0;i<nInputs; i++){   //here we can have more than one input  - merging of separate events, signal1+signal2+background
 835       rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i));
 836       gime = rl->GetLoader("TPCLoader");
 837       if (gime->TreeS()->GetEntryWithIndex(globalRowID,globalRowID) >= 0) {
 838         digarr[i]->ExpandBuffer();
 839         digarr[i]->ExpandTrackBuffer();
 840         nTimeBins = digarr[i]->GetNRows();
 841         nPads = digarr[i]->GetNCols();
 842         active[i] = kTRUE;
 843         if (!GetRegionOfInterest() || (i == 0)) digitize = kTRUE;
 844       } else {
 845         active[i] = kFALSE;
 846       }
 847       if (GetRegionOfInterest() && !digitize) break;
 848     }
 849     if (!digitize) continue;
 850     //digrow->Allocate(nTimeBins,nPads);
 851     Float_t q    = 0;
 852     Int_t labptr = 0;
 853     Int_t nElems = nTimeBins*nPads; // element is a unit of a given row's pad-timebin space
 854     for (Int_t i=0;i<nInputs; i++)
 855       if (active[i]) {
 856         pdig[i] = digarr[i]->GetDigits();
 857       }
 858     //
 859     // loop over elements i.e pad-timebin space of a row, "padNumber=elem/nTimeBins", "timeBin=elem%nTimeBins"
 860     //
 861     for (Int_t elem=0;elem<nElems; elem++) {
 862       q=0;
 863       labptr=0;
 864       // looop over digits
 865       for (Int_t i=0;i<nInputs; i++) if (active[i])       {
 866         q  += *(pdig[i]);
 867         pdig[i]++;
 868       }
 869       if (q<=0) continue;
 870       Int_t padNumber = elem/nTimeBins;
 871       Int_t timeBin   = elem%nTimeBins;
 872       Float_t gain = gainROC->GetValue(padRow,padNumber);  // get gain for given - pad-row pad
 873       q*= gain;
 874       crossTalkSignal[wireSegmentID][timeBin]+= q/nPadsPerSegment;        // Qtot per segment for a given timebin
 875       qTotSectorOld -> GetMatrixArray()[sector] += q;                      // Qtot for each sector
 876       qTotTPC += q;                                                        // Qtot for whole TPC
 877     } // end of q loop
 878   } // end of global row loop
 879
 880   //
 881   // 1.b) Dump the content of the crossTalk signal to the debug stremer - to be corealted later with the same crosstalk correction
 882   //      assumed during reconstruction
 883   //
 884   if (AliTPCReconstructor::StreamLevel()==1) {
 885     //
 886     // dump the crosstalk matrices to tree for further investigation
 887     //     a.) to estimate fluctuation of pedestal in indiviula wire segments
 888     //     b.) to check correlation between regions
 889     //     c.) to check relative conribution of signal below threshold to crosstalk
 890     for (Int_t isector=0; isector<nROCs; isector++){  //set all ellemts of crosstalk matrix to 0
 891       TMatrixD * crossTalkMatrix = (TMatrixD*)crossTalkSignalArray.At(isector);
 892       //TMatrixD * crossTalkMatrixBelow = (TMatrixD*)fCrossTalkSignalArray->At(isector+nROCs);
 893       TVectorD vecAll(crossTalkMatrix->GetNrows());
 894       //TVectorD vecBelow(crossTalkMatrix->GetNrows());
 895       //
 896       for (Int_t itime=0; itime<crossTalkMatrix->GetNcols(); itime++){
 897         for (Int_t iwire=0; iwire<crossTalkMatrix->GetNrows(); iwire++){
 898           vecAll[iwire]=(*crossTalkMatrix)(iwire,itime);
 899           //vecBelow[iwire]=(*crossTalkMatrixBelow)(iwire,itime);
 900         }
 901         (*fDebugStreamer)<<"crosstalkMatrix"<<
 902           "sector="<<isector<<
 903           "itime="<<itime<<
 904           "vecAll.="<<&vecAll<<                // crosstalk charge + charge below threshold
 905           //"vecBelow.="<<&vecBelow<<            // crosstalk contribution from signal below threshold
 906           "\n";
 907       }
 908     }
 909   }
 910
 911
 912
 913
 914   //
 915   // 2.) Loop over segments (padrows) of the TPC
 916   //
 917   for (Int_t globalRowID=0; globalRowID<param->GetNRowsTotal(); globalRowID++) {
 918     Int_t sector, padRow;
 919     if (!param->AdjustSectorRow(globalRowID,sector,padRow)) {
 920       cerr<<"AliTPC warning: invalid segment ID ! "<<globalRowID<<endl;
 921       continue;
 922     }
 923     TObjArray *arrTRF = (TObjArray*)timeResFunc.At(sector);
 924     //    TGraphErrors  *graphTRF = (TGraphErrors*)arrTRF->At(1);
 925     Int_t wireSegmentID    = param->GetWireSegment(sector,padRow);
 926     Float_t nPadsPerSegment = (Float_t)(param->GetNPadsPerSegment(wireSegmentID));
 927     //    const Float_t ampfactor = (sector<36)?factorIROC:factorOROC;      // factor for the iontail which is ROC type dependent
 928     AliTPCCalROC * gainROC  = gainTPC->GetCalROC(sector);  // pad gains per given sector
 929     AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(sector);  // noise per given sector
 930     digrow->SetID(globalRowID);
 931     Int_t nTimeBins = 0;
 932     Int_t nPads = 0;
 933     Bool_t digitize = kFALSE;
 934     for (Int_t i=0;i<nInputs; i++) {   //here we can have more than one input  - merging of separate events, signal1+signal2+background
 935       rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i));
 936       gime = rl->GetLoader("TPCLoader");
 937       if (gime->TreeS()->GetEntryWithIndex(globalRowID,globalRowID) >= 0) {
 938         digarr[i]->ExpandBuffer();
 939         digarr[i]->ExpandTrackBuffer();
 940         nTimeBins = digarr[i]->GetNRows();
 941         nPads = digarr[i]->GetNCols();
 942         active[i] = kTRUE;
 943         if (!GetRegionOfInterest() || (i == 0)) digitize = kTRUE;
 944       } else {
 945         active[i] = kFALSE;
 946       }
 947       if (GetRegionOfInterest() && !digitize) break;
 948     }
 949     if (!digitize) continue;
 950
 951     digrow->Allocate(nTimeBins,nPads);
 952     digrow->AllocateTrack(3);
 953
 954     Int_t localPad = 0;
 955     Float_t q    = 0.;
 956     Float_t qXtalk   = 0.;
 957     Float_t qIonTail = 0.;
 958     Float_t qOrig = 0.;
 959     Int_t label[1000]; //stack for 300 events
 960     Int_t labptr = 0;
 961     Int_t nElems = nTimeBins*nPads; // element is a unit of a given row's pad-timebin space
 962     for (Int_t i=0;i<nInputs; i++)
 963       if (active[i]) {
 964         pdig[i] = digarr[i]->GetDigits();
 965         ptr[i]  = digarr[i]->GetTracks();
 966       }
 967     Short_t *pdig1= digrow->GetDigits();
 968     Int_t   *ptr1= digrow->GetTracks() ;
 969     // loop over elements i.e pad-timebin space of a row
 970     for (Int_t elem=0;elem<nElems; elem++)     {
 971       q=0;
 972       labptr=0;
 973       // looop over digits
 974       for (Int_t i=0;i<nInputs; i++) if (active[i]){
 975         q  += *(pdig[i]);
 976         for (Int_t tr=0;tr<3;tr++)         {
 977           Int_t lab = ptr[i][tr*nElems];
 978           if ( (lab > 1) && *(pdig[i])>zerosup) {
 979             label[labptr]=lab+masks[i];
 980             labptr++;
 981           }
 982         }
 983         pdig[i]++;
 984         ptr[i]++;
 985       }
 986       Int_t padNumber = elem/nTimeBins;
 987       Int_t timeBin   = elem%nTimeBins;
 988       localPad = padNumber-nPads/2;
 989
 990       Float_t gain = gainROC->GetValue(padRow,padNumber);  // get gain for given - pad-row pad
 991       //if (gain<0.5){
 992       //printf("problem\n");
 993       //}
 994       q*= gain;
 995       qOrig = q;
 996       Float_t noisePad = noiseROC->GetValue(padRow,padNumber);
 997       Float_t noise  = pTPC->GetNoise()*noisePad;
 998       if ( (q/16.+noise)> zerosup){
 999         // Crosstalk correction
1000         qXtalk = (*(TMatrixD*)crossTalkSignalArray.At(sector))[wireSegmentID][timeBin];
1001         qTotPerSector = qTotSectorOld -> GetMatrixArray()[sector];
1002
1003         // Ion tail correction: being elem=padNumber*nTimeBins+timeBin;
1004         Int_t lowerElem=elem-nIonTailBins;
1005         Int_t zeroElem =(elem/nTimeBins)*nTimeBins;
1006         if (zeroElem<0) zeroElem=0;
1007         if (lowerElem<zeroElem) lowerElem=zeroElem;
1008         //
1009         qIonTail=0;
1010         if (q>0 && recoParam->GetUseIonTailCorrection()){
1011           for (Int_t i=0;i<nInputs; i++) if (active[i]){
1012             Short_t *pdigC= digarr[i]->GetDigits();
1013             if (padNumber==0) continue;
1014             if (padNumber>=nPads-1) continue;
1015             for (Int_t dpad=-1; dpad<=1; dpad++){          // calculate iontail due signals from neigborhood pads
1016               for (Int_t celem=elem-1; celem>lowerElem; celem--){
1017                 Int_t celemPad=celem+dpad*nTimeBins;
1018                 Double_t qCElem=pdigC[celemPad];
1019                 if ( qCElem<=0) continue;
1020                 //here we substract ion tail
1021                 Double_t COG=0;
1022                 if (celemPad-nTimeBins>nTimeBins && celemPad+nTimeBins<nElems){   // COG calculation in respect to current pad in pad units
1023                   Double_t sumAmp=pdigC[celemPad-nTimeBins]+pdigC[celemPad]+pdigC[celemPad+nTimeBins];
1024                   COG=(-1.0*pdigC[celemPad-nTimeBins]+pdigC[celemPad+nTimeBins])/sumAmp;
1025                 }
1026                 Int_t indexTRFPRF = (TMath::Nint(TMath::Abs(COG*10.))%20);
1027                 TGraphErrors  *graphTRFPRF = (TGraphErrors*)arrTRF->At(indexTRFPRF);
1028                 if (graphTRFPRF==NULL) continue;
1029                 // here we should get index and point of TRF corresponding to given COG position
1030                 if (graphTRFPRF->GetY()[elem-celem]<0)qIonTail+=qCElem*graphTRFPRF->GetY()[elem-celem];
1031               }
1032             }
1033           }
1034         }
1035       }
1036       //
1037       q -= qXtalk*recoParam->GetCrosstalkCorrection();
1038       q+=qIonTail;
1039       q/=16.;                                              //conversion factor
1040       q+=noise;
1041       q=TMath::Nint(q);  // round to the nearest integer
1042
1043
1044       // fill info for degugging
1045       if (AliTPCReconstructor::StreamLevel()==1 && qOrig > zerosup ) {
1046         TTreeSRedirector &cstream = *fDebugStreamer;
1047         UInt_t uid = AliTPCROC::GetTPCUniqueID(sector, padRow, padNumber);
1048         cstream <<"ionTailXtalk"<<
1049           "uid="<<uid<<                        // globla unique identifier
1050           "sector="<< sector<<
1051           "globalRowID="<<globalRowID<<
1052           "padRow="<< padRow<<                 //pad row
1053           "wireSegmentID="<< wireSegmentID<<   //wire segment 0-11, 0-3 in IROC 4-10 in OROC
1054           "localPad="<<localPad<<              // pad number -npads/2 .. npads/2
1055           "padNumber="<<padNumber<<            // pad number 0..npads
1056           "timeBin="<< timeBin<<               // time bin
1057           "nPadsPerSegment="<<nPadsPerSegment<<// number of pads per wire segment
1058           "qTotPerSector="<<qTotPerSector<<    // total charge in sector
1059           //
1060           "qTotTPC="<<qTotTPC<<                // acumulated charge without crosstalk and ion tail in full TPC
1061           "qOrig="<< qOrig<<                   // charge in given pad-row,pad,time-bin
1062           "q="<<q<<                            // q=qOrig-qXtalk-qIonTail - to check sign of the effects
1063           "qXtalk="<<qXtalk<<                  // crosstal contribtion at given position
1064           "qIonTail="<<qIonTail<<              // ion tail cotribution from past signal
1065           "\n";
1066       }// dump the results to the debug streamer if in debug mode
1067
1068       if (q > zerosup){
1069         if(q >= param->GetADCSat()) q = (Short_t)(param->GetADCSat() - 1);
1070         //digrow->SetDigitFast((Short_t)q,rows,col);
1071         *pdig1 =Short_t(q);
1072         for (Int_t tr=0;tr<3;tr++)
1073           {
1074           if (tr<labptr)
1075             ptr1[tr*nElems] = label[tr];
1076         }
1077       }
1078       pdig1++;
1079       ptr1++;
1080     }
1081
1082    if (AliTPCReconstructor::StreamLevel()==1 && qOrig > zerosup ) {
1083       cout << " sector = " << sector  << " row = " << padRow << " localPad = " << localPad << " SegmentID = " << wireSegmentID << " nPadsPerSegment = " <<  nPadsPerSegment << endl;
1084       cout << " qXtalk =   " <<  qXtalk ;
1085       cout << " qOrig = " <<  qOrig ;
1086       cout << " q = "    <<  q ;
1087       cout << " qsec = " <<  qTotPerSector ;
1088       cout << " qTotTPC "<<  qTotTPC << endl;
1089    }
1090     //
1091     //  glitch filter
1092     //
1093     digrow->GlitchFilter();
1094     //
1095     digrow->CompresBuffer(1,zerosup);
1096     digrow->CompresTrackBuffer(1);
1097     tree->Fill();
1098     if (fDebug>0) cerr<<sector<<"\t"<<padRow<<"\n";
1099   } //for (Int_t n=0; n<param->GetNRowsTotal(); n++)
1100
1101
1102   orl = AliRunLoader::GetRunLoader(fDigInput->GetOutputFolderName());
1103   ogime = orl->GetLoader("TPCLoader");
1104   ogime->WriteDigits("OVERWRITE");
1105
1106   //fDigInput->GetTreeDTPC()->Write(0,TObject::kOverwrite);
1107
1108   delete digrow;
1109   for (Int_t i1=0;i1<nInputs; i1++) delete digarr[i1];
1110   delete []masks;
1111   delete []pdig;
1112   delete []ptr;
1113   delete []active;
1114   delete []digarr;
1115 }