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