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