TPC/Sim/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
  60 using std::cout;
  61 using std::cerr;
  62 using std::endl;
  63 ClassImp(AliTPCDigitizer)
  64
  65 //___________________________________________
  66   AliTPCDigitizer::AliTPCDigitizer() :AliDigitizer(),fDebug(0)
  67 {
  68   //
  69 // Default ctor - don't use it
  70 //
  71
  72 }
  73
  74 //___________________________________________
  75 AliTPCDigitizer::AliTPCDigitizer(AliDigitizationInput* digInput)
  76   :AliDigitizer(digInput),fDebug(0)
  77 {
  78   //
  79 // ctor which should be used
  80 //
  81   AliDebug(2,"(AliDigitizationInput* digInput) was processed");
  82 }
  83
  84 //------------------------------------------------------------------------
  85 AliTPCDigitizer::~AliTPCDigitizer()
  86 {
  87 // Destructor
  88 }
  89
  90
  91
  92 //------------------------------------------------------------------------
  93 Bool_t AliTPCDigitizer::Init()
  94 {
  95 // Initialization
  96
  97  return kTRUE;
  98 }
  99
 100
 101 //------------------------------------------------------------------------
 102 void AliTPCDigitizer::Digitize(Option_t* option)
 103 {
 104   DigitizeFast(option);
 105 }
 106 //------------------------------------------------------------------------
 107 void AliTPCDigitizer::DigitizeFast(Option_t* option)
 108 {
 109
 110   // merge input tree's with summable digits
 111   //output stored in TreeTPCD
 112   char s[100];
 113   char ss[100];
 114   TString optionString = option;
 115   if (!strcmp(optionString.Data(),"deb")) {
 116     cout<<"AliTPCDigitizer:::DigitizeFast called with option deb "<<endl;
 117     fDebug = 3;
 118   }
 119   //get detector and geometry
 120
 121
 122   AliRunLoader *rl, *orl;
 123   AliLoader *gime, *ogime;
 124
 125   if (gAlice == 0x0)
 126    {
 127      Warning("DigitizeFast","gAlice is NULL. Loading from input 0");
 128      rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(0));
 129      if (rl == 0x0)
 130       {
 131         Error("DigitizeFast","Can not find Run Loader for input 0. Can not proceed.");
 132         return;
 133       }
 134      rl->LoadgAlice();
 135      rl->GetAliRun();
 136    }
 137   AliTPC *pTPC  = (AliTPC *) gAlice->GetModule("TPC");
 138   AliTPCParam * param = pTPC->GetParam();
 139
 140   //sprintf(s,param->GetTitle());
 141   snprintf(s,100,"%s",param->GetTitle());
 142   //sprintf(ss,"75x40_100x60");
 143   snprintf(ss,100,"75x40_100x60");
 144   if(strcmp(s,ss)==0){
 145     printf("2 pad-length geom hits with 3 pad-lenght geom digits...\n");
 146     delete param;
 147     param=new AliTPCParamSR();
 148   }
 149   else{
 150     //sprintf(ss,"75x40_100x60_150x60");
 151     snprintf(ss,100,"75x40_100x60_150x60");
 152    if(strcmp(s,ss)!=0) {
 153      printf("No TPC parameters found...\n");
 154      exit(2);
 155    }
 156   }
 157
 158   pTPC->GenerNoise(500000); //create table with noise
 159   //
 160   Int_t nInputs = fDigInput->GetNinputs();
 161   Int_t * masks = new Int_t[nInputs];
 162   for (Int_t i=0; i<nInputs;i++)
 163     masks[i]= fDigInput->GetMask(i);
 164   Short_t **pdig= new Short_t*[nInputs];   //pointers to the expanded digits array
 165   Int_t **ptr=  new Int_t*[nInputs];       //pointers to the expanded tracks array
 166   Bool_t *active=  new Bool_t[nInputs];    //flag for active input segments
 167   Char_t phname[100];
 168
 169   //create digits array for given sectors
 170   // make indexes
 171   //
 172   //create branch's in TPC treeD
 173   orl = AliRunLoader::GetRunLoader(fDigInput->GetOutputFolderName());
 174   ogime = orl->GetLoader("TPCLoader");
 175   TTree * tree  = ogime->TreeD();
 176   AliSimDigits * digrow = new AliSimDigits;
 177
 178   if (tree == 0x0)
 179    {
 180      ogime->MakeTree("D");
 181      tree  = ogime->TreeD();
 182    }
 183   tree->Branch("Segment","AliSimDigits",&digrow);
 184   //
 185   AliSimDigits ** digarr = new AliSimDigits*[nInputs];
 186   for (Int_t i1=0;i1<nInputs; i1++)
 187     {
 188       digarr[i1]=0;
 189      //    intree[i1]
 190       rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i1));
 191       gime = rl->GetLoader("TPCLoader");
 192       gime->LoadSDigits("read");
 193       TTree * treear =  gime->TreeS();
 194
 195       if (!treear)
 196        {
 197         cerr<<"AliTPCDigitizer: Input tree with SDigits not found in"
 198             <<" input "<< i1<<endl;
 199         for (Int_t i2=0;i2<i1+1; i2++){
 200
 201           if(digarr[i2])  delete digarr[i2];
 202         }
 203         delete [] digarr;
 204         delete [] active;
 205         delete []masks;
 206         delete []pdig;
 207         delete []ptr;
 208         return;
 209        }
 210
 211       //sprintf(phname,"lhcphase%d",i1);
 212       snprintf(phname,100,"lhcphase%d",i1);
 213       TParameter<float> *ph = (TParameter<float>*)treear->GetUserInfo()
 214                                ->FindObject("lhcphase0");
 215       if(!ph){
 216         cerr<<"AliTPCDigitizer: LHC phase  not found in"
 217             <<" input "<< i1<<endl;
 218         for (Int_t i2=0;i2<i1+1; i2++){
 219           if(digarr[i2])  delete digarr[i2];
 220         }
 221         delete [] digarr;
 222         delete [] active;
 223         delete []masks;
 224         delete []pdig;
 225         delete []ptr;
 226         return;
 227       }
 228       tree->GetUserInfo()->Add(new TParameter<float>(phname,ph->GetVal()));
 229               //
 230       if (treear->GetIndex()==0)
 231         treear->BuildIndex("fSegmentID","fSegmentID");
 232       treear->GetBranch("Segment")->SetAddress(&digarr[i1]);
 233     }
 234
 235
 236
 237
 238   //
 239
 240   param->SetZeroSup(2);
 241
 242   Int_t zerosup = param->GetZeroSup();
 243   AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetDedxGainFactor();
 244   AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
 245   //
 246   //Loop over segments of the TPC
 247
 248   for (Int_t segmentID=0; segmentID<param->GetNRowsTotal(); segmentID++)
 249    {
 250     Int_t sector, padRow;
 251     if (!param->AdjustSectorRow(segmentID,sector,padRow))
 252      {
 253       cerr<<"AliTPC warning: invalid segment ID ! "<<segmentID<<endl;
 254       continue;
 255      }
 256     AliTPCCalROC * gainROC = gainTPC->GetCalROC(sector);  // pad gains per given sector
 257     AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(sector);  // noise per given sector
 258     digrow->SetID(segmentID);
 259
 260     Int_t nTimeBins = 0;
 261     Int_t nPads = 0;
 262
 263     Bool_t digitize = kFALSE;
 264     for (Int_t i=0;i<nInputs; i++)
 265      {
 266
 267       rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i));
 268       gime = rl->GetLoader("TPCLoader");
 269
 270       if (gime->TreeS()->GetEntryWithIndex(segmentID,segmentID) >= 0) {
 271         digarr[i]->ExpandBuffer();
 272         digarr[i]->ExpandTrackBuffer();
 273         nTimeBins = digarr[i]->GetNRows();
 274         nPads = digarr[i]->GetNCols();
 275         active[i] = kTRUE;
 276         if (!GetRegionOfInterest() || (i == 0)) digitize = kTRUE;
 277       } else {
 278         active[i] = kFALSE;
 279       }
 280       if (GetRegionOfInterest() && !digitize) break;
 281      }
 282     if (!digitize) continue;
 283
 284     digrow->Allocate(nTimeBins,nPads);
 285     digrow->AllocateTrack(3);
 286
 287     Float_t q=0;
 288     Int_t label[1000]; //stack for 300 events
 289     Int_t labptr = 0;
 290
 291     Int_t nElems = nTimeBins*nPads;
 292
 293     for (Int_t i=0;i<nInputs; i++)
 294      if (active[i]) {
 295        pdig[i] = digarr[i]->GetDigits();
 296        ptr[i]  = digarr[i]->GetTracks();
 297       }
 298
 299     Short_t *pdig1= digrow->GetDigits();
 300     Int_t   *ptr1= digrow->GetTracks() ;
 301
 302
 303
 304     for (Int_t elem=0;elem<nElems; elem++)
 305      {
 306
 307        q=0;
 308        labptr=0;
 309        // looop over digits
 310         for (Int_t i=0;i<nInputs; i++) if (active[i])
 311          {
 312           //          q  += digarr[i]->GetDigitFast(rows,col);
 313             q  += *(pdig[i]);
 314
 315            for (Int_t tr=0;tr<3;tr++)
 316             {
 317              //             Int_t lab = digarr[i]->GetTrackIDFast(rows,col,tr);
 318              Int_t lab = ptr[i][tr*nElems];
 319              if ( (lab > 1) && *(pdig[i])>zerosup)
 320               {
 321                 label[labptr]=lab+masks[i];
 322                 labptr++;
 323               }
 324             }
 325            pdig[i]++;
 326            ptr[i]++;
 327          }
 328         q/=16.;  //conversion factor
 329         Float_t gain = gainROC->GetValue(padRow,elem/nTimeBins);  // get gain for given - pad-row pad
 330         //if (gain<0.5){
 331           //printf("problem\n");
 332         //}
 333         q*= gain;
 334         Float_t noisePad = noiseROC->GetValue(padRow,elem/nTimeBins);
 335         //       Float_t noise  = gRandom->Gaus(0,param->GetNoise()*param->GetNoiseNormFac());
 336         Float_t noise  = pTPC->GetNoise();
 337         q+=noise*noisePad;
 338         q=TMath::Nint(q);
 339         if (q > zerosup)
 340          {
 341           if(q >= param->GetADCSat()) q = (Short_t)(param->GetADCSat() - 1);
 342           //digrow->SetDigitFast((Short_t)q,rows,col);
 343           *pdig1 =Short_t(q);
 344           for (Int_t tr=0;tr<3;tr++)
 345            {
 346             if (tr<labptr)
 347              ptr1[tr*nElems] = label[tr];
 348            }
 349           }
 350         pdig1++;
 351         ptr1++;
 352      }
 353     //
 354     //  glitch filter
 355     //
 356     digrow->GlitchFilter();
 357     //
 358     digrow->CompresBuffer(1,zerosup);
 359     digrow->CompresTrackBuffer(1);
 360     tree->Fill();
 361     if (fDebug>0) cerr<<sector<<"\t"<<padRow<<"\n";
 362    } //for (Int_t n=0; n<param->GetNRowsTotal(); n++)
 363
 364
 365   orl = AliRunLoader::GetRunLoader(fDigInput->GetOutputFolderName());
 366   ogime = orl->GetLoader("TPCLoader");
 367   ogime->WriteDigits("OVERWRITE");
 368
 369   //fDigInput->GetTreeDTPC()->Write(0,TObject::kOverwrite);
 370
 371   delete digrow;
 372   for (Int_t i1=0;i1<nInputs; i1++) delete digarr[i1];
 373   delete []masks;
 374   delete []pdig;
 375   delete []ptr;
 376   delete []active;
 377   delete []digarr;
 378 }
 379
 380
 381
 382 //------------------------------------------------------------------------
 383 void AliTPCDigitizer::DigitizeSave(Option_t* option)
 384 {
 385   //
 386   // Merge input tree's with summable digits
 387   // Output digits stored in TreeTPCD
 388   //
 389   // Not active for long time.
 390   // Before adding modification (for ion tail calucation and for the crorsstalk) it should be
 391   //  checked one by one with currenlty used AliTPCDigitizer::DigitizeFast
 392   //
 393   TString optionString = option;
 394   if (!strcmp(optionString.Data(),"deb")) {
 395     cout<<"AliTPCDigitizer::Digitize: called with option deb "<<endl;
 396     fDebug = 3;
 397   }
 398   //get detector and geometry
 399   AliRunLoader *rl, *orl;
 400   AliLoader *gime, *ogime;
 401
 402
 403   orl = AliRunLoader::GetRunLoader(fDigInput->GetOutputFolderName());
 404   ogime = orl->GetLoader("TPCLoader");
 405
 406   rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(0));
 407   //gime = rl->GetLoader("TPCLoader");
 408   rl->GetLoader("TPCLoader");
 409   rl->LoadgAlice();
 410   AliRun* alirun = rl->GetAliRun();
 411
 412   AliTPC *pTPC  = (AliTPC *) alirun->GetModule("TPC");
 413   AliTPCParam * param = pTPC->GetParam();
 414   pTPC->GenerNoise(500000); //create teble with noise
 415   printf("noise %f \n",  param->GetNoise()*param->GetNoiseNormFac());
 416   //
 417   Int_t nInputs = fDigInput->GetNinputs();
 418   // stupid protection...
 419   if (nInputs <= 0) return;
 420   //
 421   Int_t * masks = new Int_t[nInputs];
 422   for (Int_t i=0; i<nInputs;i++)
 423     masks[i]= fDigInput->GetMask(i);
 424
 425   AliSimDigits ** digarr = new AliSimDigits*[nInputs];
 426   for(Int_t ii=0;ii<nInputs;ii++) digarr[ii]=0;
 427
 428   for (Int_t i1=0;i1<nInputs; i1++)
 429    {
 430      //digarr[i1]=0;
 431     //    intree[i1]
 432     rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i1));
 433     gime = rl->GetLoader("TPCLoader");
 434
 435     TTree * treear =  gime->TreeS();
 436     //
 437     if (!treear) {
 438       cerr<<" TPC -  not existing input = \n"<<i1<<" ";
 439       delete [] masks;
 440       for(Int_t i=0; i<nInputs; i++) delete digarr[i];
 441       delete [] digarr;
 442       return;
 443     }
 444     //
 445     TBranch * br = treear->GetBranch("fSegmentID");
 446     if (br) br->GetFile()->cd();
 447     treear->GetBranch("Segment")->SetAddress(&digarr[i1]);
 448    }
 449
 450   rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(0));
 451   gime = rl->GetLoader("TPCLoader");
 452   Stat_t nentries = gime->TreeS()->GetEntries();
 453
 454
 455   //create branch's in TPC treeD
 456   AliSimDigits * digrow = new AliSimDigits;
 457   TTree * tree  = ogime->TreeD();
 458
 459   tree->Branch("Segment","AliSimDigits",&digrow);
 460   param->SetZeroSup(2);
 461
 462   Int_t zerosup = param->GetZeroSup();
 463   //Loop over segments of the TPC
 464
 465   AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetDedxGainFactor();
 466   AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
 467   for (Int_t n=0; n<nentries; n++) {
 468     rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(0));
 469     gime = rl->GetLoader("TPCLoader");
 470     gime->TreeS()->GetEvent(n);
 471
 472     digarr[0]->ExpandBuffer();
 473     digarr[0]->ExpandTrackBuffer();
 474
 475
 476     for (Int_t i=1;i<nInputs; i++){
 477 //      fDigInput->GetInputTreeTPCS(i)->GetEntryWithIndex(digarr[0]->GetID(),digarr[0]->GetID());
 478       rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i));
 479       gime = rl->GetLoader("TPCLoader");
 480       gime->TreeS()->GetEntryWithIndex(digarr[0]->GetID(),digarr[0]->GetID());
 481       digarr[i]->ExpandBuffer();
 482       digarr[i]->ExpandTrackBuffer();
 483       if ((digarr[0]->GetID()-digarr[i]->GetID())>0)
 484        printf("problem\n");
 485
 486     }
 487
 488     Int_t sector, padRow;
 489     if (!param->AdjustSectorRow(digarr[0]->GetID(),sector,padRow)) {
 490       cerr<<"AliTPC warning: invalid segment ID ! "<<digarr[0]->GetID()<<endl;
 491       continue;
 492     }
 493
 494     AliTPCCalROC * gainROC = gainTPC->GetCalROC(sector);  // pad gains per given sector
 495     AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(sector);  // noise per given sector
 496     digrow->SetID(digarr[0]->GetID());
 497
 498     Int_t nTimeBins = digarr[0]->GetNRows();
 499     Int_t nPads = digarr[0]->GetNCols();
 500     digrow->Allocate(nTimeBins,nPads);
 501     digrow->AllocateTrack(3);
 502
 503     Float_t q=0;
 504     Int_t label[1000]; //stack for 300 events
 505     Int_t labptr = 0;
 506
 507
 508
 509     for (Int_t iTimeBin=0;iTimeBin<nTimeBins; iTimeBin++){   // iTimeBin
 510       for (Int_t iPad=0;iPad<nPads; iPad++){    // pad
 511
 512        q=0;
 513        labptr=0;
 514        // looop over digits
 515         for (Int_t i=0;i<nInputs; i++){
 516          q  += digarr[i]->GetDigitFast(iTimeBin,iPad);
 517           //q  += *(pdig[i]);
 518
 519           for (Int_t tr=0;tr<3;tr++) {
 520            Int_t lab = digarr[i]->GetTrackIDFast(iTimeBin,iPad,tr);
 521            //Int_t lab = ptr[i][tr*nElems];
 522             if ( (lab > 1) ) {
 523               label[labptr]=lab+masks[i];
 524               labptr++;
 525             }
 526           }
 527          // pdig[i]++;
 528          //ptr[i]++;
 529
 530         }
 531        q/=16.;  //conversion factor
 532        //       Float_t noise  = gRandom->Gaus(0,param->GetNoise()*param->GetNoiseNormFac());
 533        Float_t gain = gainROC->GetValue(padRow,iPad);
 534        q*= gain;
 535        Float_t noisePad = noiseROC->GetValue(padRow, iPad);
 536
 537        Float_t noise  = pTPC->GetNoise();
 538        q+=noise*noisePad;
 539        //
 540        // here we can get digits from past and add signal
 541        //
 542        //
 543        //for (Int_t jTimeBin=0; jTimeBin<iTimeBin; jTimeBin++)
 544        //  q+=ionTail
 545        //
 546
 547         q=TMath::Nint(q);
 548         if (q > zerosup){
 549
 550          if(q >= param->GetADCSat()) q = (Short_t)(param->GetADCSat() - 1);
 551          digrow->SetDigitFast((Short_t)q,iTimeBin,iPad);
 552          // *pdig1 =Short_t(q);
 553          for (Int_t tr=0;tr<3;tr++){
 554            if (tr<labptr)
 555              ((AliSimDigits*)digrow)->SetTrackIDFast(label[tr],iTimeBin,iPad,tr);
 556            //ptr1[tr*nElems] = label[tr];
 557            //else
 558              //           ((AliSimDigits*)digrow)->SetTrackIDFast(-1,iTimeBin,iPad,tr);
 559            //  ptr1[tr*nElems] = 1;
 560          }
 561        }
 562        //pdig1++;
 563        //ptr1++;
 564     }
 565     }
 566
 567     digrow->CompresBuffer(1,zerosup);
 568     digrow->CompresTrackBuffer(1);
 569     tree->Fill();
 570     if (fDebug>0) cerr<<sector<<"\t"<<padRow<<"\n";
 571   }
 572 //  printf("end TPC merging - end -Tree %s\t%p\n",fDigInput->GetInputTreeH(0)->GetName(),fDigInput->GetInputTreeH(0)->GetListOfBranches()->At(3));
 573   //fDigInput->GetTreeDTPC()->Write(0,TObject::kOverwrite);
 574     ogime->WriteDigits("OVERWRITE");
 575
 576     for (Int_t i=1;i<nInputs; i++)
 577      {
 578       rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i));
 579       gime = rl->GetLoader("TPCLoader");
 580       gime->UnloadSDigits();
 581      }
 582     ogime->UnloadDigits();
 583
 584   delete digrow;
 585   for (Int_t i1=0;i1<nInputs; i1++) delete digarr[i1];
 586   delete [] masks;
 587   delete [] digarr;
 588 }
 589
 590
 591
 592 //------------------------------------------------------------------------
 593 void AliTPCDigitizer::DigitizeWithTailAndCrossTalk(Option_t* option, TTreeSredirector *pcstream)
 594 {
 595   // Modified version of the digitization function
 596   // Modification: adding the ion tail and crosstalk:
 597   //
 598   // pcstream used in order to visually inspect data
 599   //
 600   //
 601   // Crosstalk simulation:
 602   //   1.) Calculate per time bin mean charge (per pad)  within anode wire segment
 603   //   2.) Subsract for the clusters at given time bin fraction of (mean charge) normalized by add hoc constant
 604   //       AliTPCRecoParam::GetCrosstalkCorrection() (0 if not crosstalk, 1 if ideal crosstalk)
 605   //       for simplicity we are assuming that wire segents are related to pad-rows
 606   //       Wire segmentationn is obtatined from the
 607   //       AliTPCParam::GetWireSegment(Int_t sector, Int_t row); // to be implemented
 608   //       AliTPCParam::GetNPadsPerSegment(Int_t segmentID); // to be implemented
 609   //
 610   // Ion tail simulation:
 611   //    1.) Needs signal from pad+-1, taking signal from history
 612   // merge input tree's with summable digits
 613   // output stored in TreeTPCD
 614   //
 615   char s[100];
 616   char ss[100];
 617   TString optionString = option;
 618   if (!strcmp(optionString.Data(),"deb")) {
 619     cout<<"AliTPCDigitizer:::DigitizeFast called with option deb "<<endl;
 620     fDebug = 3;
 621   }
 622   //get detector and geometry
 623
 624
 625   AliRunLoader *rl, *orl;
 626   AliLoader *gime, *ogime;
 627
 628   if (gAlice == 0x0)
 629    {
 630      Warning("DigitizeFast","gAlice is NULL. Loading from input 0");
 631      rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(0));
 632      if (rl == 0x0)
 633       {
 634         Error("DigitizeFast","Can not find Run Loader for input 0. Can not proceed.");
 635         return;
 636       }
 637      rl->LoadgAlice();
 638      rl->GetAliRun();
 639    }
 640   AliTPC *pTPC  = (AliTPC *) gAlice->GetModule("TPC");
 641   AliTPCParam * param = pTPC->GetParam();
 642
 643   //sprintf(s,param->GetTitle());
 644   snprintf(s,100,"%s",param->GetTitle());
 645   //sprintf(ss,"75x40_100x60");
 646   snprintf(ss,100,"75x40_100x60");
 647   if(strcmp(s,ss)==0){
 648     printf("2 pad-length geom hits with 3 pad-lenght geom digits...\n");
 649     delete param;
 650     param=new AliTPCParamSR();
 651   }
 652   else{
 653     //sprintf(ss,"75x40_100x60_150x60");
 654     snprintf(ss,100,"75x40_100x60_150x60");
 655    if(strcmp(s,ss)!=0) {
 656      printf("No TPC parameters found...\n");
 657      exit(2);
 658    }
 659   }
 660
 661   pTPC->GenerNoise(500000); //create table with noise
 662   //
 663   Int_t nInputs = fDigInput->GetNinputs();
 664   Int_t * masks = new Int_t[nInputs];
 665   for (Int_t i=0; i<nInputs;i++)
 666     masks[i]= fDigInput->GetMask(i);
 667   Short_t **pdig= new Short_t*[nInputs];   //pointers to the expanded digits array
 668   Int_t **ptr=  new Int_t*[nInputs];       //pointers to the expanded tracks array
 669   Bool_t *active=  new Bool_t[nInputs];    //flag for active input segments
 670   Char_t phname[100];
 671
 672   //create digits array for given sectors
 673   // make indexes
 674   //
 675   //create branch's in TPC treeD
 676   orl = AliRunLoader::GetRunLoader(fDigInput->GetOutputFolderName());
 677   ogime = orl->GetLoader("TPCLoader");
 678   TTree * tree  = ogime->TreeD();
 679   AliSimDigits * digrow = new AliSimDigits;
 680
 681   if (tree == 0x0)
 682    {
 683      ogime->MakeTree("D");
 684      tree  = ogime->TreeD();
 685    }
 686   tree->Branch("Segment","AliSimDigits",&digrow);
 687   //
 688   AliSimDigits ** digarr = new AliSimDigits*[nInputs];
 689   for (Int_t i1=0;i1<nInputs; i1++)
 690     {
 691       digarr[i1]=0;
 692      //    intree[i1]
 693       rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i1));
 694       gime = rl->GetLoader("TPCLoader");
 695       gime->LoadSDigits("read");
 696       TTree * treear =  gime->TreeS();
 697
 698       if (!treear)
 699        {
 700         cerr<<"AliTPCDigitizer: Input tree with SDigits not found in"
 701             <<" input "<< i1<<endl;
 702         for (Int_t i2=0;i2<i1+1; i2++){
 703
 704           if(digarr[i2])  delete digarr[i2];
 705         }
 706         delete [] digarr;
 707         delete [] active;
 708         delete []masks;
 709         delete []pdig;
 710         delete []ptr;
 711         return;
 712        }
 713
 714       //sprintf(phname,"lhcphase%d",i1);
 715       snprintf(phname,100,"lhcphase%d",i1);
 716       TParameter<float> *ph = (TParameter<float>*)treear->GetUserInfo()
 717                                ->FindObject("lhcphase0");
 718       if(!ph){
 719         cerr<<"AliTPCDigitizer: LHC phase  not found in"
 720             <<" input "<< i1<<endl;
 721         for (Int_t i2=0;i2<i1+1; i2++){
 722           if(digarr[i2])  delete digarr[i2];
 723         }
 724         delete [] digarr;
 725         delete [] active;
 726         delete []masks;
 727         delete []pdig;
 728         delete []ptr;
 729         return;
 730       }
 731       tree->GetUserInfo()->Add(new TParameter<float>(phname,ph->GetVal()));
 732               //
 733       if (treear->GetIndex()==0)
 734         treear->BuildIndex("fSegmentID","fSegmentID");
 735       treear->GetBranch("Segment")->SetAddress(&digarr[i1]);
 736     }
 737
 738
 739
 740
 741   //
 742
 743   param->SetZeroSup(2);
 744   Int_t zerosup = param->GetZeroSup();
 745   AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetDedxGainFactor();
 746   AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
 747   //
 748   // 1.) Make first loop to calculate mean amplitude per pad per segment
 749   //
 750   // TObjArray * crossTalkSignalArray(nsectors);
 751   // TMatrixD  crossTalkSignal(nWireSegments, timeBin);  // structure with mean signal per pad to be filled in first loop
 752   //
 753   //
 754   //2.) Loop over segments (padrows) of the TPC
 755   //
 756   for (Int_t segmentID=0; segmentID<param->GetNRowsTotal(); segmentID++)
 757    {
 758     Int_t sector, padRow;
 759     if (!param->AdjustSectorRow(segmentID,sector,padRow))
 760      {
 761       cerr<<"AliTPC warning: invalid segment ID ! "<<segmentID<<endl;
 762       continue;
 763      }
 764     AliTPCCalROC * gainROC = gainTPC->GetCalROC(sector);  // pad gains per given sector
 765     AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(sector);  // noise per given sector
 766     digrow->SetID(segmentID);
 767
 768     Int_t nTimeBins = 0;
 769     Int_t nPads = 0;
 770
 771     Bool_t digitize = kFALSE;
 772     for (Int_t i=0;i<nInputs; i++)
 773       {
 774
 775       rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i));
 776       gime = rl->GetLoader("TPCLoader");
 777
 778       if (gime->TreeS()->GetEntryWithIndex(segmentID,segmentID) >= 0) {
 779         digarr[i]->ExpandBuffer();
 780         digarr[i]->ExpandTrackBuffer();
 781         nTimeBins = digarr[i]->GetNRows();
 782         nPads = digarr[i]->GetNCols();
 783         active[i] = kTRUE;
 784         if (!GetRegionOfInterest() || (i == 0)) digitize = kTRUE;
 785       } else {
 786         active[i] = kFALSE;
 787       }
 788       if (GetRegionOfInterest() && !digitize) break;
 789      }
 790     if (!digitize) continue;
 791
 792     digrow->Allocate(nTimeBins,nPads);
 793     digrow->AllocateTrack(3);
 794
 795     Float_t q=0;
 796     Int_t label[1000]; //stack for 300 events
 797     Int_t labptr = 0;
 798
 799     Int_t nElems = nTimeBins*nPads;
 800
 801     for (Int_t i=0;i<nInputs; i++)
 802      if (active[i]) {
 803        pdig[i] = digarr[i]->GetDigits();
 804        ptr[i]  = digarr[i]->GetTracks();
 805       }
 806
 807     Short_t *pdig1= digrow->GetDigits();
 808     Int_t   *ptr1= digrow->GetTracks() ;
 809
 810
 811
 812     for (Int_t elem=0;elem<nElems; elem++)
 813      {
 814        // padNumber=elem/nTimeBins
 815        // timeBin=elem%nTimeBins
 816        q=0;
 817        labptr=0;
 818        // looop over digits
 819         for (Int_t i=0;i<nInputs; i++) if (active[i])
 820          {
 821           //          q  += digarr[i]->GetDigitFast(rows,col);
 822             q  += *(pdig[i]);
 823
 824            for (Int_t tr=0;tr<3;tr++)
 825             {
 826              //             Int_t lab = digarr[i]->GetTrackIDFast(rows,col,tr);
 827              Int_t lab = ptr[i][tr*nElems];
 828              if ( (lab > 1) && *(pdig[i])>zerosup)
 829               {
 830                 label[labptr]=lab+masks[i];
 831                 labptr++;
 832               }
 833             }
 834            pdig[i]++;
 835            ptr[i]++;
 836          }
 837         Int_t padNumber=elem/nTimeBins;
 838
 839         q/=16.;  //conversion factor
 840         Float_t gain = gainROC->GetValue(padRow,padNumber);  // get gain for given - pad-row pad
 841         //if (gain<0.5){
 842           //printf("problem\n");
 843         //}
 844         q*= gain;
 845         // Crosstalk correction:
 846         //      Double_t qCrossTalk=crossTalkSignal(wireSegment,timeBin); // signal matrix from per sector array
 847         //
 848         //
 849         // Ion tail correction:
 850         //   padNumber=elem/nTimeBins
 851         //   timeBin=elem%nTimeBins
 852         //   elem=padNumber*nTimeBins+timeBin;
 853         //    lowerElem=elem-nIonTailBins;
 854         //    if (lowerElem<0) lowerElem=0;
 855         //    if (lowerElem in previospad) lowerElem = padNumber*nTimeBins;
 856         //
 857         // for (Int_t celem=elem-1; celem>lowerElem; celem--){
 858         //  Int_t deltaT=elem-celem
 859         //
 860         // }
 861         //
 862         Float_t noisePad = noiseROC->GetValue(padRow,padNumber);
 863         //       Float_t noise  = gRandom->Gaus(0,param->GetNoise()*param->GetNoiseNormFac());
 864         Float_t noise  = pTPC->GetNoise();
 865         q+=noise*noisePad;
 866         q=TMath::Nint(q);  // round to the nearest integer
 867         /*
 868           fill infor to check consistency of the data
 869          */
 870         if (q > zerosup)
 871          {
 872           if(q >= param->GetADCSat()) q = (Short_t)(param->GetADCSat() - 1);
 873           //digrow->SetDigitFast((Short_t)q,rows,col);
 874           *pdig1 =Short_t(q);
 875           for (Int_t tr=0;tr<3;tr++)
 876            {
 877             if (tr<labptr)
 878              ptr1[tr*nElems] = label[tr];
 879            }
 880           }
 881         pdig1++;
 882         ptr1++;
 883      }
 884     //
 885     //  glitch filter
 886     //
 887     digrow->GlitchFilter();
 888     //
 889     digrow->CompresBuffer(1,zerosup);
 890     digrow->CompresTrackBuffer(1);
 891     tree->Fill();
 892     if (fDebug>0) cerr<<sector<<"\t"<<padRow<<"\n";
 893    } //for (Int_t n=0; n<param->GetNRowsTotal(); n++)
 894
 895
 896   orl = AliRunLoader::GetRunLoader(fDigInput->GetOutputFolderName());
 897   ogime = orl->GetLoader("TPCLoader");
 898   ogime->WriteDigits("OVERWRITE");
 899
 900   //fDigInput->GetTreeDTPC()->Write(0,TObject::kOverwrite);
 901
 902   delete digrow;
 903   for (Int_t i1=0;i1<nInputs; i1++) delete digarr[i1];
 904   delete []masks;
 905   delete []pdig;
 906   delete []ptr;
 907   delete []active;
 908   delete []digarr;
 909 }