TPC/AliTPC.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, 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 //                                                                           //
  20 //  Time Projection Chamber                                                  //
  21 //  This class contains the basic functions for the Time Projection Chamber  //
  22 //  detector. Functions specific to one particular geometry are              //
  23 //  contained in the derived classes                                         //
  24 //                                                                           //
  25 //Begin_Html
  26 /*
  27 <img src="picts/AliTPCClass.gif">
  28 */
  29 //End_Html
  30 //                                                                           //
  31 //                                                                          //
  32 ///////////////////////////////////////////////////////////////////////////////
  33
  34 //
  35
  36 #include <Riostream.h>
  37 #include <stdlib.h>
  38
  39 #include <TFile.h>
  40 #include <TGeometry.h>
  41 #include <TInterpreter.h>
  42 #include <TMath.h>
  43 #include <TMatrixF.h>
  44 #include <TVector.h>
  45 #include <TNode.h>
  46 #include <TObjectTable.h>
  47 #include <TParticle.h>
  48 #include <TROOT.h>
  49 #include <TRandom.h>
  50 #include <TSystem.h>
  51 #include <TTUBS.h>
  52 #include <TTree.h>
  53 #include <TVirtualMC.h>
  54 #include <TString.h>
  55 #include <TF2.h>
  56 #include <TStopwatch.h>
  57
  58 #include "AliDigits.h"
  59 #include "AliMagF.h"
  60 #include "AliPoints.h"
  61 #include "AliRun.h"
  62 #include "AliRunLoader.h"
  63 #include "AliSimDigits.h"
  64 #include "AliTPC.h"
  65 #include "AliTPC.h"
  66 #include "AliTPCDigitsArray.h"
  67 #include "AliTPCLoader.h"
  68 #include "AliTPCPRF2D.h"
  69 #include "AliTPCParamSR.h"
  70 #include "AliTPCRF1D.h"
  71 //#include "AliTPCTrackHits.h"
  72 #include "AliTPCTrackHitsV2.h"
  73 #include "AliTrackReference.h"
  74 #include "AliMC.h"
  75 #include "AliStack.h"
  76 #include "AliTPCDigitizer.h"
  77 #include "AliTPCBuffer.h"
  78 #include "AliTPCDDLRawData.h"
  79 #include "AliLog.h"
  80 #include "AliRawReader.h"
  81 #include "AliTPCRawStream.h"
  82
  83 ClassImp(AliTPC)
  84 //_____________________________________________________________________________
  85   AliTPC::AliTPC():AliDetector(),
  86                    fDefaults(0),
  87                    fSens(0),
  88                    fNsectors(0),
  89                    fDigitsArray(0),
  90                    fTPCParam(0),
  91                    fTrackHits(0),
  92                    fHitType(0),
  93                    fDigitsSwitch(0),
  94                    fSide(0),
  95                    fPrimaryIonisation(0),
  96                    fNoiseDepth(0),
  97                    fNoiseTable(0),
  98                    fCurrentNoise(0),
  99                    fActiveSectors(0)
 100
 101
 102 {
 103   //
 104   // Default constructor
 105   //
 106   fIshunt   = 0;
 107
 108   //  fTrackHitsOld = 0;
 109 #if ROOT_VERSION_CODE >= ROOT_VERSION(4,0,1)
 110   fHitType = 4; // ROOT containers
 111 #else
 112   fHitType = 2; //default CONTAINERS - based on ROOT structure
 113 #endif
 114 }
 115
 116 //_____________________________________________________________________________
 117 AliTPC::AliTPC(const char *name, const char *title)
 118   : AliDetector(name,title),
 119                    fDefaults(0),
 120                    fSens(0),
 121                    fNsectors(0),
 122                    fDigitsArray(0),
 123                    fTPCParam(0),
 124                    fTrackHits(0),
 125                    fHitType(0),
 126                    fDigitsSwitch(0),
 127                    fSide(0),
 128                    fPrimaryIonisation(0),
 129                    fNoiseDepth(0),
 130                    fNoiseTable(0),
 131                    fCurrentNoise(0),
 132                    fActiveSectors(0)
 133
 134 {
 135   //
 136   // Standard constructor
 137   //
 138
 139   //
 140   // Initialise arrays of hits and digits
 141   fHits     = new TClonesArray("AliTPChit",  176);
 142   gAlice->GetMCApp()->AddHitList(fHits);
 143   //
 144   fTrackHits = new AliTPCTrackHitsV2;
 145   fTrackHits->SetHitPrecision(0.002);
 146   fTrackHits->SetStepPrecision(0.003);
 147   fTrackHits->SetMaxDistance(100);
 148
 149   //fTrackHitsOld = new AliTPCTrackHits;  //MI - 13.09.2000
 150   //fTrackHitsOld->SetHitPrecision(0.002);
 151   //fTrackHitsOld->SetStepPrecision(0.003);
 152   //fTrackHitsOld->SetMaxDistance(100);
 153
 154
 155 #if ROOT_VERSION_CODE >= ROOT_VERSION(4,0,1)
 156   fHitType = 4; // ROOT containers
 157 #else
 158   fHitType = 2;
 159 #endif
 160
 161
 162
 163   //
 164   fIshunt     =  0;
 165   //
 166   // Initialise color attributes
 167   //PH SetMarkerColor(kYellow);
 168
 169   //
 170   //  Set TPC parameters
 171   //
 172
 173
 174   if (!strcmp(title,"Default")) {
 175     fTPCParam = new AliTPCParamSR;
 176   } else {
 177     AliWarning("In Config.C you must set non-default parameters.");
 178     fTPCParam=0;
 179   }
 180 }
 181
 182 //_____________________________________________________________________________
 183 AliTPC::~AliTPC()
 184 {
 185   //
 186   // TPC destructor
 187   //
 188
 189   fIshunt   = 0;
 190   delete fHits;
 191   delete fDigits;
 192   delete fTPCParam;
 193   delete fTrackHits; //MI 15.09.2000
 194   //  delete fTrackHitsOld; //MI 10.12.2001
 195
 196   fDigitsArray = 0x0;
 197   delete [] fNoiseTable;
 198   delete [] fActiveSectors;
 199
 200 }
 201
 202 //_____________________________________________________________________________
 203 void AliTPC::AddHit(Int_t track, Int_t *vol, Float_t *hits)
 204 {
 205   //
 206   // Add a hit to the list
 207   //
 208   if (fHitType&1){
 209     TClonesArray &lhits = *fHits;
 210     new(lhits[fNhits++]) AliTPChit(fIshunt,track,vol,hits);
 211   }
 212   if (fHitType>1)
 213     AddHit2(track,vol,hits);
 214 }
 215
 216 //_____________________________________________________________________________
 217 void AliTPC::BuildGeometry()
 218 {
 219
 220   //
 221   // Build TPC ROOT TNode geometry for the event display
 222   //
 223   TNode *nNode, *nTop;
 224   TTUBS *tubs;
 225   Int_t i;
 226   const int kColorTPC=19;
 227   char name[5], title[25];
 228   const Double_t kDegrad=TMath::Pi()/180;
 229   const Double_t kRaddeg=180./TMath::Pi();
 230
 231
 232   Float_t innerOpenAngle = fTPCParam->GetInnerAngle();
 233   Float_t outerOpenAngle = fTPCParam->GetOuterAngle();
 234
 235   Float_t innerAngleShift = fTPCParam->GetInnerAngleShift();
 236   Float_t outerAngleShift = fTPCParam->GetOuterAngleShift();
 237
 238   Int_t nLo = fTPCParam->GetNInnerSector()/2;
 239   Int_t nHi = fTPCParam->GetNOuterSector()/2;
 240
 241   const Double_t kloAng = (Double_t)TMath::Nint(innerOpenAngle*kRaddeg);
 242   const Double_t khiAng = (Double_t)TMath::Nint(outerOpenAngle*kRaddeg);
 243   const Double_t kloAngSh = (Double_t)TMath::Nint(innerAngleShift*kRaddeg);
 244   const Double_t khiAngSh = (Double_t)TMath::Nint(outerAngleShift*kRaddeg);
 245
 246
 247   const Double_t kloCorr = 1/TMath::Cos(0.5*kloAng*kDegrad);
 248   const Double_t khiCorr = 1/TMath::Cos(0.5*khiAng*kDegrad);
 249
 250   Double_t rl,ru;
 251
 252
 253   //
 254   // Get ALICE top node
 255   //
 256
 257   nTop=gAlice->GetGeometry()->GetNode("alice");
 258
 259   //  inner sectors
 260
 261   rl = fTPCParam->GetInnerRadiusLow();
 262   ru = fTPCParam->GetInnerRadiusUp();
 263
 264
 265   for(i=0;i<nLo;i++) {
 266     sprintf(name,"LS%2.2d",i);
 267     name[4]='\0';
 268     sprintf(title,"TPC low sector %3d",i);
 269     title[24]='\0';
 270
 271     tubs = new TTUBS(name,title,"void",rl*kloCorr,ru*kloCorr,250.,
 272                      kloAng*(i-0.5)+kloAngSh,kloAng*(i+0.5)+kloAngSh);
 273     tubs->SetNumberOfDivisions(1);
 274     nTop->cd();
 275     nNode = new TNode(name,title,name,0,0,0,"");
 276     nNode->SetLineColor(kColorTPC);
 277     fNodes->Add(nNode);
 278   }
 279
 280   // Outer sectors
 281
 282   rl = fTPCParam->GetOuterRadiusLow();
 283   ru = fTPCParam->GetOuterRadiusUp();
 284
 285   for(i=0;i<nHi;i++) {
 286     sprintf(name,"US%2.2d",i);
 287     name[4]='\0';
 288     sprintf(title,"TPC upper sector %d",i);
 289     title[24]='\0';
 290     tubs = new TTUBS(name,title,"void",rl*khiCorr,ru*khiCorr,250,
 291                      khiAng*(i-0.5)+khiAngSh,khiAng*(i+0.5)+khiAngSh);
 292     tubs->SetNumberOfDivisions(1);
 293     nTop->cd();
 294     nNode = new TNode(name,title,name,0,0,0,"");
 295     nNode->SetLineColor(kColorTPC);
 296     fNodes->Add(nNode);
 297   }
 298
 299 }
 300
 301 //_____________________________________________________________________________
 302 void AliTPC::CreateMaterials()
 303 {
 304   //-----------------------------------------------
 305   // Create Materials for for TPC simulations
 306   //-----------------------------------------------
 307
 308   //-----------------------------------------------------------------
 309   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
 310   //-----------------------------------------------------------------
 311
 312    Int_t iSXFLD=gAlice->Field()->Integ();
 313   Float_t sXMGMX=gAlice->Field()->Max();
 314
 315   Float_t amat[5]; // atomic numbers
 316   Float_t zmat[5]; // z
 317   Float_t wmat[5]; // proportions
 318
 319   Float_t density;
 320
 321
 322
 323   //***************** Gases *************************
 324
 325
 326   //--------------------------------------------------------------
 327   // gases - air and CO2
 328   //--------------------------------------------------------------
 329
 330   // CO2
 331
 332   amat[0]=12.011;
 333   amat[1]=15.9994;
 334
 335   zmat[0]=6.;
 336   zmat[1]=8.;
 337
 338   wmat[0]=0.2729;
 339   wmat[1]=0.7271;
 340
 341   density=0.001977;
 342
 343
 344   AliMixture(10,"CO2",amat,zmat,density,2,wmat);
 345   //
 346   // Air
 347   //
 348   amat[0]=15.9994;
 349   amat[1]=14.007;
 350   //
 351   zmat[0]=8.;
 352   zmat[1]=7.;
 353   //
 354   wmat[0]=0.233;
 355   wmat[1]=0.767;
 356   //
 357   density=0.001205;
 358
 359   AliMixture(11,"Air",amat,zmat,density,2,wmat);
 360
 361   //----------------------------------------------------------------
 362   // drift gases
 363   //----------------------------------------------------------------
 364
 365
 366   // Drift gases 1 - nonsensitive, 2 - sensitive
 367   // Ne-CO2-N (85-10-5)
 368
 369   amat[0]= 20.18;
 370   amat[1]=12.011;
 371   amat[2]=15.9994;
 372   amat[3]=14.007;
 373
 374   zmat[0]= 10.;
 375   zmat[1]=6.;
 376   zmat[2]=8.;
 377   zmat[3]=7.;
 378
 379   wmat[0]=0.7707;
 380   wmat[1]=0.0539;
 381   wmat[2]=0.1438;
 382   wmat[3]=0.0316;
 383
 384   density=0.0010252;
 385
 386   AliMixture(12,"Ne-CO2-N-1",amat,zmat,density,4,wmat);
 387   AliMixture(13,"Ne-CO2-N-2",amat,zmat,density,4,wmat);
 388
 389   //----------------------------------------------------------------------
 390   //               solid materials
 391   //----------------------------------------------------------------------
 392
 393
 394   // Kevlar C14H22O2N2
 395
 396   amat[0] = 12.011;
 397   amat[1] = 1.;
 398   amat[2] = 15.999;
 399   amat[3] = 14.006;
 400
 401   zmat[0] = 6.;
 402   zmat[1] = 1.;
 403   zmat[2] = 8.;
 404   zmat[3] = 7.;
 405
 406   wmat[0] = 14.;
 407   wmat[1] = 22.;
 408   wmat[2] = 2.;
 409   wmat[3] = 2.;
 410
 411   density = 1.45;
 412
 413   AliMixture(14,"Kevlar",amat,zmat,density,-4,wmat);
 414
 415   // NOMEX
 416
 417   amat[0] = 12.011;
 418   amat[1] = 1.;
 419   amat[2] = 15.999;
 420   amat[3] = 14.006;
 421
 422   zmat[0] = 6.;
 423   zmat[1] = 1.;
 424   zmat[2] = 8.;
 425   zmat[3] = 7.;
 426
 427   wmat[0] = 14.;
 428   wmat[1] = 22.;
 429   wmat[2] = 2.;
 430   wmat[3] = 2.;
 431
 432   density = 0.029;
 433
 434   AliMixture(15,"NOMEX",amat,zmat,density,-4,wmat);
 435
 436   // Makrolon C16H18O3
 437
 438   amat[0] = 12.011;
 439   amat[1] = 1.;
 440   amat[2] = 15.999;
 441
 442   zmat[0] = 6.;
 443   zmat[1] = 1.;
 444   zmat[2] = 8.;
 445
 446   wmat[0] = 16.;
 447   wmat[1] = 18.;
 448   wmat[2] = 3.;
 449
 450   density = 1.2;
 451
 452   AliMixture(16,"Makrolon",amat,zmat,density,-3,wmat);
 453
 454   // Tedlar C2H3F
 455
 456   amat[0] = 12.011;
 457   amat[1] = 1.;
 458   amat[2] = 18.998;
 459
 460   zmat[0] = 6.;
 461   zmat[1] = 1.;
 462   zmat[2] = 9.;
 463
 464   wmat[0] = 2.;
 465   wmat[1] = 3.;
 466   wmat[2] = 1.;
 467
 468   density = 1.71;
 469
 470   AliMixture(17, "Tedlar",amat,zmat,density,-3,wmat);
 471
 472   // Mylar C5H4O2
 473
 474   amat[0]=12.011;
 475   amat[1]=1.;
 476   amat[2]=15.9994;
 477
 478   zmat[0]=6.;
 479   zmat[1]=1.;
 480   zmat[2]=8.;
 481
 482   wmat[0]=5.;
 483   wmat[1]=4.;
 484   wmat[2]=2.;
 485
 486   density = 1.39;
 487
 488   AliMixture(18, "Mylar",amat,zmat,density,-3,wmat);
 489   // material for "prepregs"
 490   // Epoxy - C14 H20 O3
 491   // Quartz SiO2
 492   // Carbon C
 493   // prepreg1 60% C-fiber, 40% epoxy (vol)
 494   amat[0]=12.011;
 495   amat[1]=1.;
 496   amat[2]=15.994;
 497
 498   zmat[0]=6.;
 499   zmat[1]=1.;
 500   zmat[2]=8.;
 501
 502   wmat[0]=0.923;
 503   wmat[1]=0.023;
 504   wmat[2]=0.054;
 505
 506   density=1.859;
 507
 508   AliMixture(19, "Prepreg1",amat,zmat,density,3,wmat);
 509
 510   //prepreg2 60% glass-fiber, 40% epoxy
 511
 512   amat[0]=12.01;
 513   amat[1]=1.;
 514   amat[2]=15.994;
 515   amat[3]=28.086;
 516
 517   zmat[0]=6.;
 518   zmat[1]=1.;
 519   zmat[2]=8.;
 520   zmat[3]=14.;
 521
 522   wmat[0]=0.194;
 523   wmat[1]=0.023;
 524   wmat[2]=0.443;
 525   wmat[3]=0.340;
 526
 527   density=1.82;
 528
 529   AliMixture(20, "Prepreg2",amat,zmat,density,4,wmat);
 530
 531   //prepreg3 50% glass-fiber, 50% epoxy
 532
 533   amat[0]=12.01;
 534   amat[1]=1.;
 535   amat[2]=15.994;
 536   amat[3]=28.086;
 537
 538   zmat[0]=6.;
 539   zmat[1]=1.;
 540   zmat[2]=8.;
 541   zmat[3]=14.;
 542
 543   wmat[0]=0.225;
 544   wmat[1]=0.03;
 545   wmat[2]=0.443;
 546   wmat[3]=0.3;
 547
 548   density=1.163;
 549
 550   AliMixture(21, "Prepreg3",amat,zmat,density,4,wmat);
 551
 552   // G10 60% SiO2 40% epoxy
 553
 554   amat[0]=12.01;
 555   amat[1]=1.;
 556   amat[2]=15.994;
 557   amat[3]=28.086;
 558
 559   zmat[0]=6.;
 560   zmat[1]=1.;
 561   zmat[2]=8.;
 562   zmat[3]=14.;
 563
 564   wmat[0]=0.194;
 565   wmat[1]=0.023;
 566   wmat[2]=0.443;
 567   wmat[3]=0.340;
 568
 569   density=1.7;
 570
 571   AliMixture(22, "G10",amat,zmat,density,4,wmat);
 572
 573   // Al
 574
 575   amat[0] = 26.98;
 576   zmat[0] = 13.;
 577
 578   density = 2.7;
 579
 580   AliMaterial(23,"Al",amat[0],zmat[0],density,999.,999.);
 581
 582   // Si (for electronics
 583
 584   amat[0] = 28.086;
 585   zmat[0] = 14.;
 586
 587   density = 2.33;
 588
 589   AliMaterial(24,"Si",amat[0],zmat[0],density,999.,999.);
 590
 591   // Cu
 592
 593   amat[0] = 63.546;
 594   zmat[0] = 29.;
 595
 596   density = 8.96;
 597
 598   AliMaterial(25,"Cu",amat[0],zmat[0],density,999.,999.);
 599
 600   // Epoxy - C14 H20 O3
 601
 602   amat[0]=12.011;
 603   amat[1]=1.;
 604   amat[2]=15.9994;
 605
 606   zmat[0]=6.;
 607   zmat[1]=1.;
 608   zmat[2]=8.;
 609
 610   wmat[0]=14.;
 611   wmat[1]=20.;
 612   wmat[2]=3.;
 613
 614   density=1.25;
 615
 616   AliMixture(26,"Epoxy",amat,zmat,density,-3,wmat);
 617
 618   // Plexiglas  C5H8O2
 619
 620   amat[0]=12.011;
 621   amat[1]=1.;
 622   amat[2]=15.9994;
 623
 624   zmat[0]=6.;
 625   zmat[1]=1.;
 626   zmat[2]=8.;
 627
 628   wmat[0]=5.;
 629   wmat[1]=8.;
 630   wmat[2]=2.;
 631
 632   density=1.18;
 633
 634   AliMixture(27,"Plexiglas",amat,zmat,density,-3,wmat);
 635
 636   // Carbon
 637
 638   amat[0]=12.011;
 639   zmat[0]=6.;
 640   density= 2.265;
 641
 642   AliMaterial(28,"C",amat[0],zmat[0],density,999.,999.);
 643
 644   // Fe (steel for the inner heat screen)
 645
 646   amat[0]=55.845;
 647
 648   zmat[0]=26.;
 649
 650   density=7.87;
 651
 652   AliMaterial(29,"Fe",amat[0],zmat[0],density,999.,999.);
 653
 654   //----------------------------------------------------------
 655   // tracking media for gases
 656   //----------------------------------------------------------
 657
 658   AliMedium(0, "Air", 11, 0, iSXFLD, sXMGMX, 10., 999., .1, .01, .1);
 659   AliMedium(1, "Ne-CO2-N-1", 12, 0, iSXFLD, sXMGMX, 10., 999.,.1,.001, .001);
 660   AliMedium(2, "Ne-CO2-N-2", 13, 1, iSXFLD, sXMGMX, 10., 999.,.1,.001, .001);
 661   AliMedium(3,"CO2",10,0, iSXFLD, sXMGMX, 10., 999.,.1, .001, .001);
 662   AliMedium(20, "Ne-CO2-N-3", 13, 1, iSXFLD, sXMGMX, 10., 999.,.1,.001, .001);
 663   //-----------------------------------------------------------
 664   // tracking media for solids
 665   //-----------------------------------------------------------
 666
 667   AliMedium(4,"Al",23,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
 668   AliMedium(5,"Kevlar",14,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
 669   AliMedium(6,"Nomex",15,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
 670   AliMedium(7,"Makrolon",16,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
 671   AliMedium(8,"Mylar",18,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
 672   AliMedium(9,"Tedlar",17,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
 673   //
 674   AliMedium(10,"Prepreg1",19,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
 675   AliMedium(11,"Prepreg2",20,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
 676   AliMedium(12,"Prepreg3",21,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
 677   AliMedium(13,"Epoxy",26,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
 678
 679   AliMedium(14,"Cu",25,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
 680   AliMedium(15,"Si",24,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
 681   AliMedium(16,"G10",22,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
 682   AliMedium(17,"Plexiglas",27,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
 683   AliMedium(18,"Steel",29,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
 684
 685 }
 686
 687 void AliTPC::GenerNoise(Int_t tablesize)
 688 {
 689   //
 690   //Generate table with noise
 691   //
 692   if (fTPCParam==0) {
 693     // error message
 694     fNoiseDepth=0;
 695     return;
 696   }
 697   if (fNoiseTable)  delete[] fNoiseTable;
 698   fNoiseTable = new Float_t[tablesize];
 699   fNoiseDepth = tablesize;
 700   fCurrentNoise =0; //!index of the noise in  the noise table
 701
 702   Float_t norm = fTPCParam->GetNoise()*fTPCParam->GetNoiseNormFac();
 703   for (Int_t i=0;i<tablesize;i++) fNoiseTable[i]= gRandom->Gaus(0,norm);
 704 }
 705
 706 Float_t AliTPC::GetNoise()
 707 {
 708   // get noise from table
 709   //  if ((fCurrentNoise%10)==0)
 710   //  fCurrentNoise= gRandom->Rndm()*fNoiseDepth;
 711   if (fCurrentNoise>=fNoiseDepth) fCurrentNoise=0;
 712   return fNoiseTable[fCurrentNoise++];
 713   //gRandom->Gaus(0, fTPCParam->GetNoise()*fTPCParam->GetNoiseNormFac());
 714 }
 715
 716
 717 Bool_t  AliTPC::IsSectorActive(Int_t sec) const
 718 {
 719   //
 720   // check if the sector is active
 721   if (!fActiveSectors) return kTRUE;
 722   else return fActiveSectors[sec];
 723 }
 724
 725 void    AliTPC::SetActiveSectors(Int_t * sectors, Int_t n)
 726 {
 727   // activate interesting sectors
 728   SetTreeAddress();//just for security
 729   if (!fActiveSectors) fActiveSectors = new Bool_t[fTPCParam->GetNSector()];
 730   for (Int_t i=0;i<fTPCParam->GetNSector();i++) fActiveSectors[i]=kFALSE;
 731   for (Int_t i=0;i<n;i++)
 732     if ((sectors[i]>=0) && sectors[i]<fTPCParam->GetNSector())  fActiveSectors[sectors[i]]=kTRUE;
 733
 734 }
 735
 736 void    AliTPC::SetActiveSectors(Int_t flag)
 737 {
 738   //
 739   // activate sectors which were hitted by tracks
 740   //loop over tracks
 741   SetTreeAddress();//just for security
 742   if (fHitType==0) return;  // if Clones hit - not short volume ID information
 743   if (!fActiveSectors) fActiveSectors = new Bool_t[fTPCParam->GetNSector()];
 744   if (flag) {
 745     for (Int_t i=0;i<fTPCParam->GetNSector();i++) fActiveSectors[i]=kTRUE;
 746     return;
 747   }
 748   for (Int_t i=0;i<fTPCParam->GetNSector();i++) fActiveSectors[i]=kFALSE;
 749   TBranch * branch=0;
 750   if (TreeH() == 0x0)
 751    {
 752      AliFatal("Can not find TreeH in folder");
 753      return;
 754    }
 755   if (fHitType>1) branch = TreeH()->GetBranch("TPC2");
 756   else branch = TreeH()->GetBranch("TPC");
 757   Stat_t ntracks = TreeH()->GetEntries();
 758   // loop over all hits
 759   AliDebug(1,Form("Got %d tracks",ntracks));
 760
 761   for(Int_t track=0;track<ntracks;track++) {
 762     ResetHits();
 763     //
 764     if (fTrackHits && fHitType&4) {
 765       TBranch * br1 = TreeH()->GetBranch("fVolumes");
 766       TBranch * br2 = TreeH()->GetBranch("fNVolumes");
 767       br1->GetEvent(track);
 768       br2->GetEvent(track);
 769       Int_t *volumes = fTrackHits->GetVolumes();
 770       for (Int_t j=0;j<fTrackHits->GetNVolumes(); j++) {
 771         if (volumes[j]>-1 && volumes[j]<fTPCParam->GetNSector()) {
 772           fActiveSectors[volumes[j]]=kTRUE;
 773         }
 774         else {
 775             AliError(Form("Volume %d -> sector number %d is outside (0..%d)",
 776                           j,
 777                           volumes[j],
 778                           fTPCParam->GetNSector()));
 779         }
 780       }
 781     }
 782
 783     //
 784 //     if (fTrackHitsOld && fHitType&2) {
 785 //       TBranch * br = TreeH()->GetBranch("fTrackHitsInfo");
 786 //       br->GetEvent(track);
 787 //       AliObjectArray * ar = fTrackHitsOld->fTrackHitsInfo;
 788 //       for (UInt_t j=0;j<ar->GetSize();j++){
 789 //      fActiveSectors[((AliTrackHitsInfo*)ar->At(j))->fVolumeID] =kTRUE;
 790 //       }
 791 //     }
 792   }
 793 }
 794
 795
 796
 797
 798 //_____________________________________________________________________________
 799 void AliTPC::Digits2Raw()
 800 {
 801 // convert digits of the current event to raw data
 802
 803   static const Int_t kThreshold = 0;
 804
 805   fLoader->LoadDigits();
 806   TTree* digits = fLoader->TreeD();
 807   if (!digits) {
 808     AliError("No digits tree");
 809     return;
 810   }
 811
 812   AliSimDigits digarr;
 813   AliSimDigits* digrow = &digarr;
 814   digits->GetBranch("Segment")->SetAddress(&digrow);
 815
 816   const char* fileName = "AliTPCDDL.dat";
 817   AliTPCBuffer* buffer  = new AliTPCBuffer(fileName);
 818   //Verbose level
 819   // 0: Silent
 820   // 1: cout messages
 821   // 2: txt files with digits
 822   //BE CAREFUL, verbose level 2 MUST be used only for debugging and
 823   //it is highly suggested to use this mode only for debugging digits files
 824   //reasonably small, because otherwise the size of the txt files can reach
 825   //quickly several MB wasting time and disk space.
 826   buffer->SetVerbose(0);
 827
 828   Int_t nEntries = Int_t(digits->GetEntries());
 829   Int_t previousSector = -1;
 830   Int_t subSector = 0;
 831   for (Int_t i = 0; i < nEntries; i++) {
 832     digits->GetEntry(i);
 833     Int_t sector, row;
 834     fTPCParam->AdjustSectorRow(digarr.GetID(), sector, row);
 835     if(previousSector != sector) {
 836       subSector = 0;
 837       previousSector = sector;
 838     }
 839
 840     if (sector < 36) { //inner sector [0;35]
 841       if (row != 30) {
 842         //the whole row is written into the output file
 843         buffer->WriteRowBinary(kThreshold, digrow, 0, 0, 0,
 844                                sector, subSector, row);
 845       } else {
 846         //only the pads in the range [37;48] are written into the output file
 847         buffer->WriteRowBinary(kThreshold, digrow, 37, 48, 1,
 848                                sector, subSector, row);
 849         subSector = 1;
 850         //only the pads outside the range [37;48] are written into the output file
 851         buffer->WriteRowBinary(kThreshold, digrow, 37, 48, 2,
 852                                sector, subSector, row);
 853       }//end else
 854
 855     } else { //outer sector [36;71]
 856       if (row == 54) subSector = 2;
 857       if ((row != 27) && (row != 76)) {
 858         buffer->WriteRowBinary(kThreshold, digrow, 0, 0, 0,
 859                                sector, subSector, row);
 860       } else if (row == 27) {
 861         //only the pads outside the range [43;46] are written into the output file
 862         buffer->WriteRowBinary(kThreshold, digrow, 43, 46, 2,
 863                                  sector, subSector, row);
 864         subSector = 1;
 865         //only the pads in the range [43;46] are written into the output file
 866         buffer->WriteRowBinary(kThreshold, digrow, 43, 46, 1,
 867                                  sector, subSector, row);
 868       } else if (row == 76) {
 869         //only the pads outside the range [33;88] are written into the output file
 870         buffer->WriteRowBinary(kThreshold, digrow, 33, 88, 2,
 871                                sector, subSector, row);
 872         subSector = 3;
 873         //only the pads in the range [33;88] are written into the output file
 874         buffer->WriteRowBinary(kThreshold, digrow, 33, 88, 1,
 875                                  sector, subSector, row);
 876       }
 877     }//end else
 878   }//end for
 879
 880   delete buffer;
 881   fLoader->UnloadDigits();
 882
 883   AliTPCDDLRawData rawWriter;
 884   rawWriter.SetVerbose(0);
 885
 886   rawWriter.RawData(fileName);
 887   gSystem->Unlink(fileName);
 888
 889 }
 890
 891
 892 //_____________________________________________________________________________
 893 Bool_t AliTPC::Raw2SDigits(AliRawReader* rawReader){
 894   // Converts the TPC raw data into summable digits
 895   // The method is used for merging simulated and
 896   // real data events
 897   if (fLoader->TreeS() == 0x0 ) {
 898     fLoader->MakeTree("S");
 899   }
 900
 901   if(fDefaults == 0) SetDefaults();  // check if the parameters are set
 902
 903   //setup TPCDigitsArray
 904   if(GetDigitsArray()) delete GetDigitsArray();
 905
 906   AliTPCDigitsArray *arr = new AliTPCDigitsArray;
 907   arr->SetClass("AliSimDigits");
 908   arr->Setup(fTPCParam);
 909   arr->MakeTree(fLoader->TreeS());
 910
 911   SetDigitsArray(arr);
 912
 913   // set zero suppression to "0"
 914   fTPCParam->SetZeroSup(0);
 915
 916   // Loop over sectors
 917   const Int_t kmaxTime = fTPCParam->GetMaxTBin();
 918   const Int_t kNIS = fTPCParam->GetNInnerSector();
 919   const Int_t kNOS = fTPCParam->GetNOuterSector();
 920   const Int_t kNS = kNIS + kNOS;
 921
 922   Short_t** allBins = NULL; //array which contains the data for one sector
 923
 924   for(Int_t iSector = 0; iSector < kNS; iSector++) {
 925
 926     Int_t nRows = fTPCParam->GetNRow(iSector);
 927     Int_t nDDLs = 0, indexDDL = 0;
 928     if (iSector < kNIS) {
 929       nDDLs = 2;
 930       indexDDL = iSector * 2;
 931     }
 932     else {
 933       nDDLs = 4;
 934       indexDDL = (iSector-kNIS) * 4 + kNIS * 2;
 935     }
 936
 937     // Loas the raw data for corresponding DDLs
 938     rawReader->Reset();
 939     AliTPCRawStream input(rawReader);
 940     rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
 941
 942     // Alocate and init the array with the sector data
 943     allBins = new Short_t*[nRows];
 944     for (Int_t iRow = 0; iRow < nRows; iRow++) {
 945       Int_t maxPad = fTPCParam->GetNPads(iSector,iRow);
 946       Int_t maxBin = kmaxTime*maxPad;
 947       allBins[iRow] = new Short_t[maxBin];
 948       memset(allBins[iRow],0,sizeof(Short_t)*maxBin);
 949     }
 950
 951     // Begin loop over altro data
 952     while (input.Next()) {
 953
 954       if (input.GetSector() != iSector)
 955         AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",iSector,input.GetSector()));
 956
 957       Int_t iRow = input.GetRow();
 958       if (iRow < 0 || iRow >= nRows)
 959         AliFatal(Form("Pad-row index (%d) outside the range (%d -> %d) !",
 960                       iRow, 0, nRows -1));
 961       Int_t iPad = input.GetPad();
 962
 963       Int_t maxPad = fTPCParam->GetNPads(iSector,iRow);
 964
 965       if (iPad < 0 || iPad >= maxPad)
 966         AliFatal(Form("Pad index (%d) outside the range (%d -> %d) !",
 967                       iPad, 0, maxPad -1));
 968
 969       Int_t iTimeBin = input.GetTime();
 970       if ( iTimeBin < 0 || iTimeBin >= kmaxTime)
 971         AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
 972                       iTimeBin, 0, kmaxTime -1));
 973
 974       Int_t maxBin = kmaxTime*maxPad;
 975
 976       if (((iPad*kmaxTime+iTimeBin) >= maxBin) ||
 977           ((iPad*kmaxTime+iTimeBin) < 0))
 978         AliFatal(Form("Index outside the allowed range"
 979                       " Sector=%d Row=%d Pad=%d Timebin=%d"
 980                       " (Max.index=%d)",iSector,iRow,iPad,iTimeBin,maxBin));
 981
 982       allBins[iRow][iPad*kmaxTime+iTimeBin] = input.GetSignal();
 983
 984     } // End loop over altro data
 985
 986     // Now fill the digits array
 987     if (fDigitsArray->GetTree()==0) {
 988       AliFatal("Tree not set in fDigitsArray");
 989     }
 990
 991     for (Int_t iRow = 0; iRow < nRows; iRow++) {
 992       AliDigits * dig = fDigitsArray->CreateRow(iSector,iRow);
 993
 994       Int_t maxPad = fTPCParam->GetNPads(iSector,iRow);
 995       for(Int_t iPad = 0; iPad < maxPad; iPad++) {
 996         for(Int_t iTimeBin = 0; iTimeBin < kmaxTime; iTimeBin++) {
 997           Short_t q = allBins[iRow][iPad*kmaxTime + iTimeBin];
 998           if (q <= 0) continue;
 999           q *= 16;
1000           dig->SetDigitFast((Short_t)q,iTimeBin,iPad);
1001         }
1002       }
1003       fDigitsArray->StoreRow(iSector,iRow);
1004       Int_t ndig = dig->GetDigitSize();
1005
1006       AliDebug(10,
1007                Form("*** Sector, row, compressed digits %d %d %d ***\n",
1008                     iSector,iRow,ndig));
1009
1010       fDigitsArray->ClearRow(iSector,iRow);
1011
1012     } // end of the sector digitization
1013
1014     for (Int_t iRow = 0; iRow < nRows; iRow++)
1015       delete [] allBins[iRow];
1016
1017     delete [] allBins;
1018
1019   }
1020
1021   fLoader->WriteSDigits("OVERWRITE");
1022
1023   if(GetDigitsArray()) delete GetDigitsArray();
1024   SetDigitsArray(0x0);
1025
1026   return kTRUE;
1027 }
1028
1029 //______________________________________________________________________
1030 AliDigitizer* AliTPC::CreateDigitizer(AliRunDigitizer* manager) const
1031 {
1032   return new AliTPCDigitizer(manager);
1033 }
1034 //__
1035 void AliTPC::SDigits2Digits2(Int_t /*eventnumber*/)
1036 {
1037   //create digits from summable digits
1038   GenerNoise(500000); //create teble with noise
1039
1040   //conect tree with sSDigits
1041   TTree *t = fLoader->TreeS();
1042
1043   if (t == 0x0) {
1044     fLoader->LoadSDigits("READ");
1045     t = fLoader->TreeS();
1046     if (t == 0x0) {
1047       AliError("Can not get input TreeS");
1048       return;
1049     }
1050   }
1051
1052   if (fLoader->TreeD() == 0x0) fLoader->MakeTree("D");
1053
1054   AliSimDigits digarr, *dummy=&digarr;
1055   TBranch* sdb = t->GetBranch("Segment");
1056   if (sdb == 0x0) {
1057     AliError("Can not find branch with segments in TreeS.");
1058     return;
1059   }
1060
1061   sdb->SetAddress(&dummy);
1062
1063   Stat_t nentries = t->GetEntries();
1064
1065   // set zero suppression
1066
1067   fTPCParam->SetZeroSup(2);
1068
1069   // get zero suppression
1070
1071   Int_t zerosup = fTPCParam->GetZeroSup();
1072
1073   //make tree with digits
1074
1075   AliTPCDigitsArray *arr = new AliTPCDigitsArray;
1076   arr->SetClass("AliSimDigits");
1077   arr->Setup(fTPCParam);
1078   arr->MakeTree(fLoader->TreeD());
1079
1080   AliTPCParam * par = fTPCParam;
1081
1082   //Loop over segments of the TPC
1083
1084   for (Int_t n=0; n<nentries; n++) {
1085     t->GetEvent(n);
1086     Int_t sec, row;
1087     if (!par->AdjustSectorRow(digarr.GetID(),sec,row)) {
1088       AliWarning(Form("Invalid segment ID ! %d",digarr.GetID()));
1089       continue;
1090     }
1091     if (!IsSectorActive(sec)) continue;
1092
1093     AliSimDigits * digrow =(AliSimDigits*) arr->CreateRow(sec,row);
1094     Int_t nrows = digrow->GetNRows();
1095     Int_t ncols = digrow->GetNCols();
1096
1097     digrow->ExpandBuffer();
1098     digarr.ExpandBuffer();
1099     digrow->ExpandTrackBuffer();
1100     digarr.ExpandTrackBuffer();
1101
1102
1103     Short_t * pamp0 = digarr.GetDigits();
1104     Int_t   * ptracks0 = digarr.GetTracks();
1105     Short_t * pamp1 = digrow->GetDigits();
1106     Int_t   * ptracks1 = digrow->GetTracks();
1107     Int_t  nelems =nrows*ncols;
1108     Int_t saturation = fTPCParam->GetADCSat() - 1;
1109     //use internal structure of the AliDigits - for speed reason
1110     //if you cahnge implementation
1111     //of the Alidigits - it must be rewriten -
1112     for (Int_t i= 0; i<nelems; i++){
1113       Float_t q = TMath::Nint(Float_t(*pamp0)/16.+GetNoise());
1114       if (q>zerosup){
1115         if (q>saturation) q=saturation;
1116         *pamp1=(Short_t)q;
1117
1118         ptracks1[0]=ptracks0[0];
1119         ptracks1[nelems]=ptracks0[nelems];
1120         ptracks1[2*nelems]=ptracks0[2*nelems];
1121       }
1122       pamp0++;
1123       pamp1++;
1124       ptracks0++;
1125       ptracks1++;
1126     }
1127
1128     arr->StoreRow(sec,row);
1129     arr->ClearRow(sec,row);
1130   }
1131
1132
1133   //write results
1134   fLoader->WriteDigits("OVERWRITE");
1135
1136   delete arr;
1137 }
1138 //__________________________________________________________________
1139 void AliTPC::SetDefaults(){
1140   //
1141   // setting the defaults
1142   //
1143
1144   // Set response functions
1145
1146   //
1147   AliRunLoader* rl = (AliRunLoader*)fLoader->GetEventFolder()->FindObject(AliRunLoader::GetRunLoaderName());
1148   rl->CdGAFile();
1149   AliTPCParamSR *param=(AliTPCParamSR*)gDirectory->Get("75x40_100x60");
1150   if(param){
1151     AliInfo("You are using 2 pad-length geom hits with 3 pad-lenght geom digits...");
1152     delete param;
1153     param = new AliTPCParamSR();
1154   }
1155   else {
1156     param=(AliTPCParamSR*)gDirectory->Get("75x40_100x60_150x60");
1157   }
1158   if(!param){
1159     AliFatal("No TPC parameters found");
1160   }
1161
1162
1163   AliTPCPRF2D    * prfinner   = new AliTPCPRF2D;
1164   AliTPCPRF2D    * prfouter1   = new AliTPCPRF2D;
1165   AliTPCPRF2D    * prfouter2   = new AliTPCPRF2D;
1166   AliTPCRF1D     * rf    = new AliTPCRF1D(kTRUE);
1167   rf->SetGauss(param->GetZSigma(),param->GetZWidth(),1.);
1168   rf->SetOffset(3*param->GetZSigma());
1169   rf->Update();
1170
1171   TDirectory *savedir=gDirectory;
1172   TFile *f=TFile::Open("$ALICE_ROOT/TPC/AliTPCprf2d.root");
1173   if (!f->IsOpen())
1174     AliFatal("Can't open $ALICE_ROOT/TPC/AliTPCprf2d.root !");
1175
1176   TString s;
1177   prfinner->Read("prf_07504_Gati_056068_d02");
1178   //PH Set different names
1179   s=prfinner->GetGRF()->GetName();
1180   s+="in";
1181   prfinner->GetGRF()->SetName(s.Data());
1182
1183   prfouter1->Read("prf_10006_Gati_047051_d03");
1184   s=prfouter1->GetGRF()->GetName();
1185   s+="out1";
1186   prfouter1->GetGRF()->SetName(s.Data());
1187
1188   prfouter2->Read("prf_15006_Gati_047051_d03");
1189   s=prfouter2->GetGRF()->GetName();
1190   s+="out2";
1191   prfouter2->GetGRF()->SetName(s.Data());
1192
1193   f->Close();
1194   savedir->cd();
1195
1196   param->SetInnerPRF(prfinner);
1197   param->SetOuter1PRF(prfouter1);
1198   param->SetOuter2PRF(prfouter2);
1199   param->SetTimeRF(rf);
1200
1201   // set fTPCParam
1202
1203   SetParam(param);
1204
1205
1206   fDefaults = 1;
1207
1208 }
1209 //__________________________________________________________________
1210 void AliTPC::Hits2Digits()
1211 {
1212   //
1213   // creates digits from hits
1214   //
1215   if (!fTPCParam->IsGeoRead()){
1216     //
1217     // read transformation matrices for gGeoManager
1218     //
1219     fTPCParam->ReadGeoMatrices();
1220   }
1221
1222   fLoader->LoadHits("read");
1223   fLoader->LoadDigits("recreate");
1224   AliRunLoader* runLoader = fLoader->GetRunLoader();
1225
1226   for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) {
1227     //PH    runLoader->GetEvent(iEvent);
1228     SetActiveSectors();
1229     Hits2Digits(iEvent);
1230   }
1231
1232   fLoader->UnloadHits();
1233   fLoader->UnloadDigits();
1234 }
1235 //__________________________________________________________________
1236 void AliTPC::Hits2Digits(Int_t eventnumber)
1237 {
1238  //----------------------------------------------------
1239  // Loop over all sectors for a single event
1240  //----------------------------------------------------
1241   AliRunLoader* rl = (AliRunLoader*)fLoader->GetEventFolder()->FindObject(AliRunLoader::GetRunLoaderName());
1242   rl->GetEvent(eventnumber);
1243   if (fLoader->TreeH() == 0x0) {
1244     if(fLoader->LoadHits()) {
1245       AliError("Can not load hits.");
1246     }
1247   }
1248   SetTreeAddress();
1249
1250   if (fLoader->TreeD() == 0x0 ) {
1251     fLoader->MakeTree("D");
1252     if (fLoader->TreeD() == 0x0 ) {
1253       AliError("Can not get TreeD");
1254       return;
1255     }
1256   }
1257
1258   if(fDefaults == 0) SetDefaults();  // check if the parameters are set
1259   GenerNoise(500000); //create teble with noise
1260
1261   //setup TPCDigitsArray
1262
1263   if(GetDigitsArray()) delete GetDigitsArray();
1264
1265   AliTPCDigitsArray *arr = new AliTPCDigitsArray;
1266   arr->SetClass("AliSimDigits");
1267   arr->Setup(fTPCParam);
1268
1269   arr->MakeTree(fLoader->TreeD());
1270   SetDigitsArray(arr);
1271
1272   fDigitsSwitch=0; // standard digits
1273
1274   for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++)
1275     if (IsSectorActive(isec)) {
1276       AliDebug(1,Form("Hits2Digits","Sector %d is active.",isec));
1277       Hits2DigitsSector(isec);
1278     }
1279     else {
1280       AliDebug(1,Form("Hits2Digits","Sector %d is NOT active.",isec));
1281     }
1282
1283   fLoader->WriteDigits("OVERWRITE");
1284
1285 //this line prevents the crash in the similar one
1286 //on the beginning of this method
1287 //destructor attempts to reset the tree, which is deleted by the loader
1288 //need to be redesign
1289   if(GetDigitsArray()) delete GetDigitsArray();
1290   SetDigitsArray(0x0);
1291
1292 }
1293
1294 //__________________________________________________________________
1295 void AliTPC::Hits2SDigits2(Int_t eventnumber)
1296 {
1297
1298   //-----------------------------------------------------------
1299   //   summable digits - 16 bit "ADC", no noise, no saturation
1300   //-----------------------------------------------------------
1301
1302   //----------------------------------------------------
1303   // Loop over all sectors for a single event
1304   //----------------------------------------------------
1305
1306   AliRunLoader* rl = fLoader->GetRunLoader();
1307
1308   rl->GetEvent(eventnumber);
1309   if (fLoader->TreeH() == 0x0) {
1310     if(fLoader->LoadHits()) {
1311       AliError("Can not load hits.");
1312       return;
1313     }
1314   }
1315   SetTreeAddress();
1316
1317
1318   if (fLoader->TreeS() == 0x0 ) {
1319     fLoader->MakeTree("S");
1320   }
1321
1322   if(fDefaults == 0) SetDefaults();
1323
1324   GenerNoise(500000); //create table with noise
1325   //setup TPCDigitsArray
1326
1327   if(GetDigitsArray()) delete GetDigitsArray();
1328
1329
1330   AliTPCDigitsArray *arr = new AliTPCDigitsArray;
1331   arr->SetClass("AliSimDigits");
1332   arr->Setup(fTPCParam);
1333   arr->MakeTree(fLoader->TreeS());
1334
1335   SetDigitsArray(arr);
1336
1337   fDigitsSwitch=1; // summable digits
1338
1339     // set zero suppression to "0"
1340
1341   fTPCParam->SetZeroSup(0);
1342
1343   for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++)
1344     if (IsSectorActive(isec)) {
1345       Hits2DigitsSector(isec);
1346     }
1347
1348   fLoader->WriteSDigits("OVERWRITE");
1349
1350 //this line prevents the crash in the similar one
1351 //on the beginning of this method
1352 //destructor attempts to reset the tree, which is deleted by the loader
1353 //need to be redesign
1354   if(GetDigitsArray()) delete GetDigitsArray();
1355   SetDigitsArray(0x0);
1356 }
1357 //__________________________________________________________________
1358
1359 void AliTPC::Hits2SDigits()
1360 {
1361
1362   //-----------------------------------------------------------
1363   //   summable digits - 16 bit "ADC", no noise, no saturation
1364   //-----------------------------------------------------------
1365
1366   if (!fTPCParam->IsGeoRead()){
1367     //
1368     // read transformation matrices for gGeoManager
1369     //
1370     fTPCParam->ReadGeoMatrices();
1371   }
1372
1373   fLoader->LoadHits("read");
1374   fLoader->LoadSDigits("recreate");
1375   AliRunLoader* runLoader = fLoader->GetRunLoader();
1376
1377   for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) {
1378     runLoader->GetEvent(iEvent);
1379     SetTreeAddress();
1380     SetActiveSectors();
1381     Hits2SDigits2(iEvent);
1382   }
1383
1384   fLoader->UnloadHits();
1385   fLoader->UnloadSDigits();
1386 }
1387 //_____________________________________________________________________________
1388
1389 void AliTPC::Hits2DigitsSector(Int_t isec)
1390 {
1391   //-------------------------------------------------------------------
1392   // TPC conversion from hits to digits.
1393   //-------------------------------------------------------------------
1394
1395   //-----------------------------------------------------------------
1396   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1397   //-----------------------------------------------------------------
1398
1399   //-------------------------------------------------------
1400   //  Get the access to the track hits
1401   //-------------------------------------------------------
1402
1403   // check if the parameters are set - important if one calls this method
1404   // directly, not from the Hits2Digits
1405
1406   if(fDefaults == 0) SetDefaults();
1407
1408   TTree *tH = TreeH(); // pointer to the hits tree
1409   if (tH == 0x0) {
1410     AliFatal("Can not find TreeH in folder");
1411     return;
1412   }
1413
1414   Stat_t ntracks = tH->GetEntries();
1415
1416   if( ntracks > 0){
1417
1418   //-------------------------------------------
1419   //  Only if there are any tracks...
1420   //-------------------------------------------
1421
1422     TObjArray **row;
1423
1424     Int_t nrows =fTPCParam->GetNRow(isec);
1425
1426     row= new TObjArray* [nrows+2]; // 2 extra rows for cross talk
1427
1428     MakeSector(isec,nrows,tH,ntracks,row);
1429
1430     //--------------------------------------------------------
1431     //   Digitize this sector, row by row
1432     //   row[i] is the pointer to the TObjArray of TVectors,
1433     //   each one containing electrons accepted on this
1434     //   row, assigned into tracks
1435     //--------------------------------------------------------
1436
1437     Int_t i;
1438
1439     if (fDigitsArray->GetTree()==0) {
1440       AliFatal("Tree not set in fDigitsArray");
1441     }
1442
1443     for (i=0;i<nrows;i++){
1444
1445       AliDigits * dig = fDigitsArray->CreateRow(isec,i);
1446
1447       DigitizeRow(i,isec,row);
1448
1449       fDigitsArray->StoreRow(isec,i);
1450
1451       Int_t ndig = dig->GetDigitSize();
1452
1453       AliDebug(10,
1454                Form("*** Sector, row, compressed digits %d %d %d ***\n",
1455                     isec,i,ndig));
1456
1457       fDigitsArray->ClearRow(isec,i);
1458
1459
1460     } // end of the sector digitization
1461
1462     for(i=0;i<nrows+2;i++){
1463       row[i]->Delete();
1464       delete row[i];
1465     }
1466
1467     delete [] row; // delete the array of pointers to TObjArray-s
1468
1469   } // ntracks >0
1470
1471 } // end of Hits2DigitsSector
1472
1473
1474 //_____________________________________________________________________________
1475 void AliTPC::DigitizeRow(Int_t irow,Int_t isec,TObjArray **rows)
1476 {
1477   //-----------------------------------------------------------
1478   // Single row digitization, coupling from the neighbouring
1479   // rows taken into account
1480   //-----------------------------------------------------------
1481
1482   //-----------------------------------------------------------------
1483   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1484   // Modified: Marian Ivanov GSI Darmstadt, m.ivanov@gsi.de
1485   //-----------------------------------------------------------------
1486
1487   Float_t zerosup = fTPCParam->GetZeroSup();
1488
1489   fCurrentIndex[1]= isec;
1490
1491
1492   Int_t nofPads = fTPCParam->GetNPads(isec,irow);
1493   Int_t nofTbins = fTPCParam->GetMaxTBin();
1494   Int_t indexRange[4];
1495   //
1496   //  Integrated signal for this row
1497   //  and a single track signal
1498   //
1499
1500   TMatrixF *m1 = new TMatrixF(0,nofPads,0,nofTbins); // integrated
1501   TMatrixF *m2 = new TMatrixF(0,nofPads,0,nofTbins); // single
1502   //
1503   TMatrixF &total  = *m1;
1504
1505   //  Array of pointers to the label-signal list
1506
1507   Int_t nofDigits = nofPads*nofTbins; // number of digits for this row
1508   Float_t  **pList = new Float_t* [nofDigits];
1509
1510   Int_t lp;
1511   Int_t i1;
1512   for(lp=0;lp<nofDigits;lp++)pList[lp]=0; // set all pointers to NULL
1513   //
1514   //calculate signal
1515   //
1516   Int_t row1=irow;
1517   Int_t row2=irow+2;
1518   for (Int_t row= row1;row<=row2;row++){
1519     Int_t nTracks= rows[row]->GetEntries();
1520     for (i1=0;i1<nTracks;i1++){
1521       fCurrentIndex[2]= row;
1522       fCurrentIndex[3]=irow+1;
1523       if (row==irow+1){
1524         m2->Zero();  // clear single track signal matrix
1525         Float_t trackLabel = GetSignal(rows[row],i1,m2,m1,indexRange);
1526         GetList(trackLabel,nofPads,m2,indexRange,pList);
1527       }
1528       else   GetSignal(rows[row],i1,0,m1,indexRange);
1529     }
1530   }
1531
1532   Int_t tracks[3];
1533
1534   AliDigits *dig = fDigitsArray->GetRow(isec,irow);
1535   Int_t gi=-1;
1536   Float_t fzerosup = zerosup+0.5;
1537   for(Int_t it=0;it<nofTbins;it++){
1538     for(Int_t ip=0;ip<nofPads;ip++){
1539       gi++;
1540       Float_t q=total(ip,it);
1541       if(fDigitsSwitch == 0){
1542         q+=GetNoise();
1543         if(q <=fzerosup) continue; // do not fill zeros
1544         q = TMath::Nint(q);
1545         if(q >= fTPCParam->GetADCSat()) q = fTPCParam->GetADCSat() - 1;  // saturation
1546
1547       }
1548
1549       else {
1550         if(q <= 0.) continue; // do not fill zeros
1551         if(q>2000.) q=2000.;
1552         q *= 16.;
1553         q = TMath::Nint(q);
1554       }
1555
1556       //
1557       //  "real" signal or electronic noise (list = -1)?
1558       //
1559
1560       for(Int_t j1=0;j1<3;j1++){
1561         tracks[j1] = (pList[gi]) ?(Int_t)(*(pList[gi]+j1)) : -2;
1562       }
1563
1564 //Begin_Html
1565 /*
1566   <A NAME="AliDigits"></A>
1567   using of AliDigits object
1568 */
1569 //End_Html
1570       dig->SetDigitFast((Short_t)q,it,ip);
1571       if (fDigitsArray->IsSimulated()) {
1572         ((AliSimDigits*)dig)->SetTrackIDFast(tracks[0],it,ip,0);
1573         ((AliSimDigits*)dig)->SetTrackIDFast(tracks[1],it,ip,1);
1574         ((AliSimDigits*)dig)->SetTrackIDFast(tracks[2],it,ip,2);
1575       }
1576
1577     } // end of loop over time buckets
1578   }  // end of lop over pads
1579
1580   //
1581   //  This row has been digitized, delete nonused stuff
1582   //
1583
1584   for(lp=0;lp<nofDigits;lp++){
1585     if(pList[lp]) delete [] pList[lp];
1586   }
1587
1588   delete [] pList;
1589
1590   delete m1;
1591   delete m2;
1592
1593 } // end of DigitizeRow
1594
1595 //_____________________________________________________________________________
1596
1597 Float_t AliTPC::GetSignal(TObjArray *p1, Int_t ntr,
1598              TMatrixF *m1, TMatrixF *m2,Int_t *indexRange)
1599 {
1600
1601   //---------------------------------------------------------------
1602   //  Calculates 2-D signal (pad,time) for a single track,
1603   //  returns a pointer to the signal matrix and the track label
1604   //  No digitization is performed at this level!!!
1605   //---------------------------------------------------------------
1606
1607   //-----------------------------------------------------------------
1608   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1609   // Modified: Marian Ivanov
1610   //-----------------------------------------------------------------
1611
1612   TVector *tv;
1613
1614   tv = (TVector*)p1->At(ntr); // pointer to a track
1615   TVector &v = *tv;
1616
1617   Float_t label = v(0);
1618   Int_t centralPad = (fTPCParam->GetNPads(fCurrentIndex[1],fCurrentIndex[3]-1))/2;
1619
1620   Int_t nElectrons = (tv->GetNrows()-1)/5;
1621   indexRange[0]=9999; // min pad
1622   indexRange[1]=-1; // max pad
1623   indexRange[2]=9999; //min time
1624   indexRange[3]=-1; // max time
1625
1626   TMatrixF &signal = *m1;
1627   TMatrixF &total = *m2;
1628   //
1629   //  Loop over all electrons
1630   //
1631   for(Int_t nel=0; nel<nElectrons; nel++){
1632     Int_t idx=nel*5;
1633     Float_t aval =  v(idx+4);
1634     Float_t eltoadcfac=aval*fTPCParam->GetTotalNormFac();
1635     Float_t xyz[4]={v(idx+1),v(idx+2),v(idx+3),v(idx+5)};
1636     Int_t n = ((AliTPCParamSR*)fTPCParam)->CalcResponseFast(xyz,fCurrentIndex,fCurrentIndex[3]);
1637
1638     Int_t *index = fTPCParam->GetResBin(0);
1639     Float_t *weight = & (fTPCParam->GetResWeight(0));
1640
1641     if (n>0) for (Int_t i =0; i<n; i++){
1642       Int_t pad=index[1]+centralPad;  //in digit coordinates central pad has coordinate 0
1643
1644       if (pad>=0){
1645         Int_t time=index[2];
1646         Float_t qweight = *(weight)*eltoadcfac;
1647
1648         if (m1!=0) signal(pad,time)+=qweight;
1649         total(pad,time)+=qweight;
1650         if (indexRange[0]>pad) indexRange[0]=pad;
1651         if (indexRange[1]<pad) indexRange[1]=pad;
1652         if (indexRange[2]>time) indexRange[2]=time;
1653         if (indexRange[3]<time) indexRange[3]=time;
1654
1655         index+=3;
1656         weight++;
1657
1658       }
1659     }
1660   } // end of loop over electrons
1661
1662   return label; // returns track label when finished
1663 }
1664
1665 //_____________________________________________________________________________
1666 void AliTPC::GetList(Float_t label,Int_t np,TMatrixF *m,
1667                      Int_t *indexRange, Float_t **pList)
1668 {
1669   //----------------------------------------------------------------------
1670   //  Updates the list of tracks contributing to digits for a given row
1671   //----------------------------------------------------------------------
1672
1673   //-----------------------------------------------------------------
1674   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1675   //-----------------------------------------------------------------
1676
1677   TMatrixF &signal = *m;
1678
1679   // lop over nonzero digits
1680
1681   for(Int_t it=indexRange[2];it<indexRange[3]+1;it++){
1682     for(Int_t ip=indexRange[0];ip<indexRange[1]+1;ip++){
1683
1684
1685       // accept only the contribution larger than 500 electrons (1/2 s_noise)
1686
1687       if(signal(ip,it)<0.5) continue;
1688
1689       Int_t globalIndex = it*np+ip; // globalIndex starts from 0!
1690
1691       if(!pList[globalIndex]){
1692
1693         //
1694         // Create new list (6 elements - 3 signals and 3 labels),
1695         //
1696
1697         pList[globalIndex] = new Float_t [6];
1698
1699         // set list to -1
1700
1701         *pList[globalIndex] = -1.;
1702         *(pList[globalIndex]+1) = -1.;
1703         *(pList[globalIndex]+2) = -1.;
1704         *(pList[globalIndex]+3) = -1.;
1705         *(pList[globalIndex]+4) = -1.;
1706         *(pList[globalIndex]+5) = -1.;
1707
1708         *pList[globalIndex] = label;
1709         *(pList[globalIndex]+3) = signal(ip,it);
1710       }
1711       else {
1712
1713         // check the signal magnitude
1714
1715         Float_t highest = *(pList[globalIndex]+3);
1716         Float_t middle = *(pList[globalIndex]+4);
1717         Float_t lowest = *(pList[globalIndex]+5);
1718
1719         //
1720         //  compare the new signal with already existing list
1721         //
1722
1723         if(signal(ip,it)<lowest) continue; // neglect this track
1724
1725         //
1726
1727         if (signal(ip,it)>highest){
1728           *(pList[globalIndex]+5) = middle;
1729           *(pList[globalIndex]+4) = highest;
1730           *(pList[globalIndex]+3) = signal(ip,it);
1731
1732           *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
1733           *(pList[globalIndex]+1) = *pList[globalIndex];
1734           *pList[globalIndex] = label;
1735         }
1736         else if (signal(ip,it)>middle){
1737           *(pList[globalIndex]+5) = middle;
1738           *(pList[globalIndex]+4) = signal(ip,it);
1739
1740           *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
1741           *(pList[globalIndex]+1) = label;
1742         }
1743         else{
1744           *(pList[globalIndex]+5) = signal(ip,it);
1745           *(pList[globalIndex]+2) = label;
1746         }
1747       }
1748
1749     } // end of loop over pads
1750   } // end of loop over time bins
1751
1752 }//end of GetList
1753 //___________________________________________________________________
1754 void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH,
1755                         Stat_t ntracks,TObjArray **row)
1756 {
1757
1758   //-----------------------------------------------------------------
1759   // Prepares the sector digitization, creates the vectors of
1760   // tracks for each row of this sector. The track vector
1761   // contains the track label and the position of electrons.
1762   //-----------------------------------------------------------------
1763
1764   //-----------------------------------------------------------------
1765   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1766   //-----------------------------------------------------------------
1767
1768   Float_t gasgain = fTPCParam->GetGasGain();
1769   Int_t i;
1770   Float_t xyz[5];
1771
1772   AliTPChit *tpcHit; // pointer to a sigle TPC hit
1773   //MI change
1774   TBranch * branch=0;
1775   if (fHitType>1) branch = TH->GetBranch("TPC2");
1776   else branch = TH->GetBranch("TPC");
1777
1778
1779   //----------------------------------------------
1780   // Create TObjArray-s, one for each row,
1781   // each TObjArray will store the TVectors
1782   // of electrons, one TVectors per each track.
1783   //----------------------------------------------
1784
1785   Int_t *nofElectrons = new Int_t [nrows+2]; // electron counter for each row
1786   TVector **tracks = new TVector* [nrows+2]; //pointers to the track vectors
1787
1788   for(i=0; i<nrows+2; i++){
1789     row[i] = new TObjArray;
1790     nofElectrons[i]=0;
1791     tracks[i]=0;
1792   }
1793
1794
1795
1796   //--------------------------------------------------------------------
1797   //  Loop over tracks, the "track" contains the full history
1798   //--------------------------------------------------------------------
1799
1800   Int_t previousTrack,currentTrack;
1801   previousTrack = -1; // nothing to store so far!
1802
1803   for(Int_t track=0;track<ntracks;track++){
1804     Bool_t isInSector=kTRUE;
1805     ResetHits();
1806     isInSector = TrackInVolume(isec,track);
1807     if (!isInSector) continue;
1808     //MI change
1809     branch->GetEntry(track); // get next track
1810
1811     //M.I. changes
1812
1813     tpcHit = (AliTPChit*)FirstHit(-1);
1814
1815     //--------------------------------------------------------------
1816     //  Loop over hits
1817     //--------------------------------------------------------------
1818
1819
1820     while(tpcHit){
1821
1822       Int_t sector=tpcHit->fSector; // sector number
1823       if(sector != isec){
1824         tpcHit = (AliTPChit*) NextHit();
1825         continue;
1826       }
1827
1828       // Remove hits which arrive before the TPC opening gate signal
1829       if(((fTPCParam->GetZLength(isec)-TMath::Abs(tpcHit->Z()))
1830           /fTPCParam->GetDriftV()+tpcHit->Time())<fTPCParam->GetGateDelay()) {
1831         tpcHit = (AliTPChit*) NextHit();
1832         continue;
1833       }
1834
1835       currentTrack = tpcHit->Track(); // track number
1836
1837       if(currentTrack != previousTrack){
1838
1839         // store already filled fTrack
1840
1841         for(i=0;i<nrows+2;i++){
1842           if(previousTrack != -1){
1843             if(nofElectrons[i]>0){
1844               TVector &v = *tracks[i];
1845               v(0) = previousTrack;
1846               tracks[i]->ResizeTo(5*nofElectrons[i]+1); // shrink if necessary
1847               row[i]->Add(tracks[i]);
1848             }
1849             else {
1850               delete tracks[i]; // delete empty TVector
1851               tracks[i]=0;
1852             }
1853           }
1854
1855           nofElectrons[i]=0;
1856           tracks[i] = new TVector(601); // TVectors for the next fTrack
1857
1858         } // end of loop over rows
1859
1860         previousTrack=currentTrack; // update track label
1861       }
1862
1863       Int_t qI = (Int_t) (tpcHit->fQ); // energy loss (number of electrons)
1864
1865       //---------------------------------------------------
1866       //  Calculate the electron attachment probability
1867       //---------------------------------------------------
1868
1869
1870       Float_t time = 1.e6*(fTPCParam->GetZLength(isec)-TMath::Abs(tpcHit->Z()))
1871         /fTPCParam->GetDriftV();
1872       // in microseconds!
1873       Float_t attProb = fTPCParam->GetAttCoef()*
1874         fTPCParam->GetOxyCont()*time; //  fraction!
1875
1876       //-----------------------------------------------
1877       //  Loop over electrons
1878       //-----------------------------------------------
1879       Int_t index[3];
1880       index[1]=isec;
1881       for(Int_t nel=0;nel<qI;nel++){
1882         // skip if electron lost due to the attachment
1883         if((gRandom->Rndm(0)) < attProb) continue; // electron lost!
1884         xyz[0]=tpcHit->X();
1885         xyz[1]=tpcHit->Y();
1886         xyz[2]=tpcHit->Z();
1887         //
1888         // protection for the nonphysical avalanche size (10**6 maximum)
1889         //
1890         Double_t rn=TMath::Max(gRandom->Rndm(0),1.93e-22);
1891         xyz[3]= (Float_t) (-gasgain*TMath::Log(rn));
1892         index[0]=1;
1893
1894         TransportElectron(xyz,index);
1895         Int_t rowNumber;
1896         fTPCParam->GetPadRow(xyz,index);
1897         // Electron track time (for pileup simulation)
1898         xyz[4] = tpcHit->Time()/fTPCParam->GetTSample();
1899         // row 0 - cross talk from the innermost row
1900         // row fNRow+1 cross talk from the outermost row
1901         rowNumber = index[2]+1;
1902         //transform position to local digit coordinates
1903         //relative to nearest pad row
1904         if ((rowNumber<0)||rowNumber>fTPCParam->GetNRow(isec)+1) continue;
1905         Float_t x1,y1;
1906         if (isec <fTPCParam->GetNInnerSector()) {
1907           x1 = xyz[1]*fTPCParam->GetInnerPadPitchWidth();
1908           y1 = fTPCParam->GetYInner(rowNumber);
1909         }
1910         else{
1911           x1=xyz[1]*fTPCParam->GetOuterPadPitchWidth();
1912           y1 = fTPCParam->GetYOuter(rowNumber);
1913         }
1914         // gain inefficiency at the wires edges - linear
1915         x1=TMath::Abs(x1);
1916         y1-=1.;
1917         if(x1>y1) xyz[3]*=TMath::Max(1.e-6,(y1-x1+1.));
1918
1919         nofElectrons[rowNumber]++;
1920         //----------------------------------
1921         // Expand vector if necessary
1922         //----------------------------------
1923         if(nofElectrons[rowNumber]>120){
1924           Int_t range = tracks[rowNumber]->GetNrows();
1925           if((nofElectrons[rowNumber])>(range-1)/5){
1926
1927             tracks[rowNumber]->ResizeTo(range+500); // Add 100 electrons
1928           }
1929         }
1930
1931         TVector &v = *tracks[rowNumber];
1932         Int_t idx = 5*nofElectrons[rowNumber]-4;
1933         Real_t * position = &(((TVector&)v)(idx)); //make code faster
1934         memcpy(position,xyz,5*sizeof(Float_t));
1935
1936       } // end of loop over electrons
1937
1938       tpcHit = (AliTPChit*)NextHit();
1939
1940     } // end of loop over hits
1941   } // end of loop over tracks
1942
1943     //
1944     //   store remaining track (the last one) if not empty
1945     //
1946
1947   for(i=0;i<nrows+2;i++){
1948     if(nofElectrons[i]>0){
1949       TVector &v = *tracks[i];
1950       v(0) = previousTrack;
1951       tracks[i]->ResizeTo(5*nofElectrons[i]+1); // shrink if necessary
1952       row[i]->Add(tracks[i]);
1953     }
1954     else{
1955       delete tracks[i];
1956       tracks[i]=0;
1957     }
1958   }
1959
1960   delete [] tracks;
1961   delete [] nofElectrons;
1962
1963 } // end of MakeSector
1964
1965
1966 //_____________________________________________________________________________
1967 void AliTPC::Init()
1968 {
1969   //
1970   // Initialise TPC detector after definition of geometry
1971   //
1972   AliDebug(1,"*********************************************");
1973 }
1974
1975 //_____________________________________________________________________________
1976 void AliTPC::MakeBranch(Option_t* option)
1977 {
1978   //
1979   // Create Tree branches for the TPC.
1980   //
1981   AliDebug(1,"");
1982   Int_t buffersize = 4000;
1983   char branchname[10];
1984   sprintf(branchname,"%s",GetName());
1985
1986   const char *h = strstr(option,"H");
1987
1988   if ( h && (fHitType<=1) && (fHits == 0x0)) fHits = new TClonesArray("AliTPChit", 176);//skowron 20.06.03
1989
1990   AliDetector::MakeBranch(option);
1991
1992   const char *d = strstr(option,"D");
1993
1994   if (fDigits   && fLoader->TreeD() && d) {
1995     MakeBranchInTree(gAlice->TreeD(), branchname, &fDigits, buffersize, 0);
1996   }
1997
1998   if (fHitType>1) MakeBranch2(option,0); // MI change 14.09.2000
1999 }
2000
2001 //_____________________________________________________________________________
2002 void AliTPC::ResetDigits()
2003 {
2004   //
2005   // Reset number of digits and the digits array for this detector
2006   //
2007   fNdigits   = 0;
2008   if (fDigits)   fDigits->Clear();
2009 }
2010
2011
2012
2013 //_____________________________________________________________________________
2014 void AliTPC::SetSens(Int_t sens)
2015 {
2016
2017   //-------------------------------------------------------------
2018   // Activates/deactivates the sensitive strips at the center of
2019   // the pad row -- this is for the space-point resolution calculations
2020   //-------------------------------------------------------------
2021
2022   //-----------------------------------------------------------------
2023   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2024   //-----------------------------------------------------------------
2025
2026   fSens = sens;
2027 }
2028
2029
2030 void AliTPC::SetSide(Float_t side=0.)
2031 {
2032   // choice of the TPC side
2033
2034   fSide = side;
2035
2036 }
2037 //_____________________________________________________________________________
2038
2039 void AliTPC::TransportElectron(Float_t *xyz, Int_t *index)
2040 {
2041   //
2042   // electron transport taking into account:
2043   // 1. diffusion,
2044   // 2.ExB at the wires
2045   // 3. nonisochronity
2046   //
2047   // xyz and index must be already transformed to system 1
2048   //
2049
2050   fTPCParam->Transform1to2(xyz,index);
2051
2052   //add diffusion
2053   Float_t driftl=xyz[2];
2054   if(driftl<0.01) driftl=0.01;
2055   driftl=TMath::Sqrt(driftl);
2056   Float_t sigT = driftl*(fTPCParam->GetDiffT());
2057   Float_t sigL = driftl*(fTPCParam->GetDiffL());
2058   xyz[0]=gRandom->Gaus(xyz[0],sigT);
2059   xyz[1]=gRandom->Gaus(xyz[1],sigT);
2060   xyz[2]=gRandom->Gaus(xyz[2],sigL);
2061
2062   // ExB
2063
2064   if (fTPCParam->GetMWPCReadout()==kTRUE){
2065     Float_t dx = fTPCParam->Transform2to2NearestWire(xyz,index);
2066     xyz[1]+=dx*(fTPCParam->GetOmegaTau());
2067   }
2068   //add nonisochronity (not implemented yet)
2069 }
2070
2071 ClassImp(AliTPChit)
2072   //______________________________________________________________________
2073   AliTPChit::AliTPChit()
2074             :AliHit(),
2075              fSector(0),
2076              fPadRow(0),
2077              fQ(0),
2078              fTime(0)
2079 {
2080   //
2081   // default
2082   //
2083
2084 }
2085 //_____________________________________________________________________________
2086 AliTPChit::AliTPChit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits)
2087           :AliHit(shunt,track),
2088              fSector(0),
2089              fPadRow(0),
2090              fQ(0),
2091              fTime(0)
2092 {
2093   //
2094   // Creates a TPC hit object
2095   //
2096   fSector     = vol[0];
2097   fPadRow     = vol[1];
2098   fX          = hits[0];
2099   fY          = hits[1];
2100   fZ          = hits[2];
2101   fQ          = hits[3];
2102   fTime       = hits[4];
2103 }
2104
2105 //________________________________________________________________________
2106 // Additional code because of the AliTPCTrackHitsV2
2107
2108 void AliTPC::MakeBranch2(Option_t *option,const char */*file*/)
2109 {
2110   //
2111   // Create a new branch in the current Root Tree
2112   // The branch of fHits is automatically split
2113   // MI change 14.09.2000
2114   AliDebug(1,"");
2115   if (fHitType<2) return;
2116   char branchname[10];
2117   sprintf(branchname,"%s2",GetName());
2118   //
2119   // Get the pointer to the header
2120   const char *cH = strstr(option,"H");
2121   //
2122   if (fTrackHits   && TreeH() && cH && fHitType&4) {
2123     AliDebug(1,"Making branch for Type 4 Hits");
2124     TreeH()->Branch(branchname,"AliTPCTrackHitsV2",&fTrackHits,fBufferSize,99);
2125   }
2126
2127 //   if (fTrackHitsOld   && TreeH() && cH && fHitType&2) {
2128 //     AliDebug(1,"Making branch for Type 2 Hits");
2129 //     AliObjectBranch * branch = new AliObjectBranch(branchname,"AliTPCTrackHits",&fTrackHitsOld,
2130 //                                                    TreeH(),fBufferSize,99);
2131 //     TreeH()->GetListOfBranches()->Add(branch);
2132 //   }
2133 }
2134
2135 void AliTPC::SetTreeAddress()
2136 {
2137   //Sets tree address for hits
2138   if (fHitType<=1) {
2139     if (fHits == 0x0 ) fHits = new TClonesArray("AliTPChit", 176);//skowron 20.06.03
2140     AliDetector::SetTreeAddress();
2141   }
2142   if (fHitType>1) SetTreeAddress2();
2143 }
2144
2145 void AliTPC::SetTreeAddress2()
2146 {
2147   //
2148   // Set branch address for the TrackHits Tree
2149   //
2150   AliDebug(1,"");
2151
2152   TBranch *branch;
2153   char branchname[20];
2154   sprintf(branchname,"%s2",GetName());
2155   //
2156   // Branch address for hit tree
2157   TTree *treeH = TreeH();
2158   if ((treeH)&&(fHitType&4)) {
2159     branch = treeH->GetBranch(branchname);
2160     if (branch) {
2161       branch->SetAddress(&fTrackHits);
2162       AliDebug(1,"fHitType&4 Setting");
2163     }
2164     else
2165       AliDebug(1,"fHitType&4 Failed (can not find branch)");
2166
2167   }
2168  //  if ((treeH)&&(fHitType&2)) {
2169 //     branch = treeH->GetBranch(branchname);
2170 //     if (branch) {
2171 //       branch->SetAddress(&fTrackHitsOld);
2172 //       AliDebug(1,"fHitType&2 Setting");
2173 //     }
2174 //     else
2175 //       AliDebug(1,"fHitType&2 Failed (can not find branch)");
2176 //   }
2177   //set address to TREETR
2178
2179   TTree *treeTR = TreeTR();
2180   if (treeTR && fTrackReferences) {
2181     branch = treeTR->GetBranch(GetName());
2182     if (branch) branch->SetAddress(&fTrackReferences);
2183   }
2184
2185 }
2186
2187 void AliTPC::FinishPrimary()
2188 {
2189   if (fTrackHits &&fHitType&4)      fTrackHits->FlushHitStack();
2190   //  if (fTrackHitsOld && fHitType&2)  fTrackHitsOld->FlushHitStack();
2191 }
2192
2193
2194 void AliTPC::AddHit2(Int_t track, Int_t *vol, Float_t *hits)
2195 {
2196   //
2197   // add hit to the list
2198   Int_t rtrack;
2199   if (fIshunt) {
2200     int primary = gAlice->GetMCApp()->GetPrimary(track);
2201     gAlice->GetMCApp()->Particle(primary)->SetBit(kKeepBit);
2202     rtrack=primary;
2203   } else {
2204     rtrack=track;
2205     gAlice->GetMCApp()->FlagTrack(track);
2206   }
2207   if (fTrackHits && fHitType&4)
2208     fTrackHits->AddHitKartez(vol[0],rtrack, hits[0],
2209                              hits[1],hits[2],(Int_t)hits[3],hits[4]);
2210  //  if (fTrackHitsOld &&fHitType&2 )
2211 //     fTrackHitsOld->AddHitKartez(vol[0],rtrack, hits[0],
2212 //                                 hits[1],hits[2],(Int_t)hits[3]);
2213
2214 }
2215
2216 void AliTPC::ResetHits()
2217 {
2218   if (fHitType&1) AliDetector::ResetHits();
2219   if (fHitType>1) ResetHits2();
2220 }
2221
2222 void AliTPC::ResetHits2()
2223 {
2224   //
2225   //reset hits
2226   if (fTrackHits && fHitType&4) fTrackHits->Clear();
2227   // if (fTrackHitsOld && fHitType&2) fTrackHitsOld->Clear();
2228
2229 }
2230
2231 AliHit* AliTPC::FirstHit(Int_t track)
2232 {
2233   if (fHitType>1) return FirstHit2(track);
2234   return AliDetector::FirstHit(track);
2235 }
2236 AliHit* AliTPC::NextHit()
2237 {
2238   //
2239   // gets next hit
2240   //
2241   if (fHitType>1) return NextHit2();
2242
2243   return AliDetector::NextHit();
2244 }
2245
2246 AliHit* AliTPC::FirstHit2(Int_t track)
2247 {
2248   //
2249   // Initialise the hit iterator
2250   // Return the address of the first hit for track
2251   // If track>=0 the track is read from disk
2252   // while if track<0 the first hit of the current
2253   // track is returned
2254   //
2255   if(track>=0) {
2256     gAlice->ResetHits();
2257     TreeH()->GetEvent(track);
2258   }
2259   //
2260   if (fTrackHits && fHitType&4) {
2261     fTrackHits->First();
2262     return fTrackHits->GetHit();
2263   }
2264  //  if (fTrackHitsOld && fHitType&2) {
2265 //     fTrackHitsOld->First();
2266 //     return fTrackHitsOld->GetHit();
2267 //   }
2268
2269   else return 0;
2270 }
2271
2272 AliHit* AliTPC::NextHit2()
2273 {
2274   //
2275   //Return the next hit for the current track
2276
2277
2278 //   if (fTrackHitsOld && fHitType&2) {
2279 //     fTrackHitsOld->Next();
2280 //     return fTrackHitsOld->GetHit();
2281 //   }
2282   if (fTrackHits) {
2283     fTrackHits->Next();
2284     return fTrackHits->GetHit();
2285   }
2286   else
2287     return 0;
2288 }
2289
2290 void AliTPC::LoadPoints(Int_t)
2291 {
2292   //
2293   Int_t a = 0;
2294
2295   if(fHitType==1) AliDetector::LoadPoints(a);
2296   else LoadPoints2(a);
2297 }
2298
2299
2300 void AliTPC::RemapTrackHitIDs(Int_t *map)
2301 {
2302   //
2303   // remapping
2304   //
2305   if (!fTrackHits) return;
2306
2307 //   if (fTrackHitsOld && fHitType&2){
2308 //     AliObjectArray * arr = fTrackHitsOld->fTrackHitsInfo;
2309 //     for (UInt_t i=0;i<arr->GetSize();i++){
2310 //       AliTrackHitsInfo * info = (AliTrackHitsInfo *)(arr->At(i));
2311 //       info->fTrackID = map[info->fTrackID];
2312 //     }
2313 //   }
2314 //  if (fTrackHitsOld && fHitType&4){
2315   if (fTrackHits && fHitType&4){
2316     TClonesArray * arr = fTrackHits->GetArray();;
2317     for (Int_t i=0;i<arr->GetEntriesFast();i++){
2318       AliTrackHitsParamV2 * info = (AliTrackHitsParamV2 *)(arr->At(i));
2319       info->SetTrackID(map[info->GetTrackID()]);
2320     }
2321   }
2322 }
2323
2324 Bool_t   AliTPC::TrackInVolume(Int_t id,Int_t track)
2325 {
2326   //return bool information - is track in given volume
2327   //load only part of the track information
2328   //return true if current track is in volume
2329   //
2330   //  return kTRUE;
2331  //  if (fTrackHitsOld && fHitType&2) {
2332 //     TBranch * br = TreeH()->GetBranch("fTrackHitsInfo");
2333 //     br->GetEvent(track);
2334 //     AliObjectArray * ar = fTrackHitsOld->fTrackHitsInfo;
2335 //     for (UInt_t j=0;j<ar->GetSize();j++){
2336 //       if (  ((AliTrackHitsInfo*)ar->At(j))->fVolumeID==id) return kTRUE;
2337 //     }
2338 //   }
2339
2340   if (fTrackHits && fHitType&4) {
2341     TBranch * br1 = TreeH()->GetBranch("fVolumes");
2342     TBranch * br2 = TreeH()->GetBranch("fNVolumes");
2343     br2->GetEvent(track);
2344     br1->GetEvent(track);
2345     Int_t *volumes = fTrackHits->GetVolumes();
2346     Int_t nvolumes = fTrackHits->GetNVolumes();
2347     if (!volumes && nvolumes>0) {
2348       AliWarning(Form("Problematic track\t%d\t%d",track,nvolumes));
2349       return kFALSE;
2350     }
2351     for (Int_t j=0;j<nvolumes; j++)
2352       if (volumes[j]==id) return kTRUE;
2353   }
2354
2355   if (fHitType&1) {
2356     TBranch * br = TreeH()->GetBranch("fSector");
2357     br->GetEvent(track);
2358     for (Int_t j=0;j<fHits->GetEntriesFast();j++){
2359       if (  ((AliTPChit*)fHits->At(j))->fSector==id) return kTRUE;
2360     }
2361   }
2362   return kFALSE;
2363
2364 }
2365
2366 //_____________________________________________________________________________
2367 void AliTPC::LoadPoints2(Int_t)
2368 {
2369   //
2370   // Store x, y, z of all hits in memory
2371   //
2372   //  if (fTrackHits == 0 && fTrackHitsOld==0) return;
2373   if (fTrackHits == 0 ) return;
2374   //
2375   Int_t nhits =0;
2376   if (fHitType&4) nhits = fTrackHits->GetEntriesFast();
2377   //  if (fHitType&2) nhits = fTrackHitsOld->GetEntriesFast();
2378
2379   if (nhits == 0) return;
2380   Int_t tracks = gAlice->GetMCApp()->GetNtrack();
2381   if (fPoints == 0) fPoints = new TObjArray(tracks);
2382   AliHit *ahit;
2383   //
2384   Int_t *ntrk=new Int_t[tracks];
2385   Int_t *limi=new Int_t[tracks];
2386   Float_t **coor=new Float_t*[tracks];
2387   for(Int_t i=0;i<tracks;i++) {
2388     ntrk[i]=0;
2389     coor[i]=0;
2390     limi[i]=0;
2391   }
2392   //
2393   AliPoints *points = 0;
2394   Float_t *fp=0;
2395   Int_t trk;
2396   Int_t chunk=nhits/4+1;
2397   //
2398   // Loop over all the hits and store their position
2399   //
2400   ahit = FirstHit2(-1);
2401   while (ahit){
2402     trk=ahit->GetTrack();
2403     if(ntrk[trk]==limi[trk]) {
2404       //
2405       // Initialise a new track
2406       fp=new Float_t[3*(limi[trk]+chunk)];
2407       if(coor[trk]) {
2408         memcpy(fp,coor[trk],sizeof(Float_t)*3*limi[trk]);
2409         delete [] coor[trk];
2410       }
2411       limi[trk]+=chunk;
2412       coor[trk] = fp;
2413     } else {
2414       fp = coor[trk];
2415     }
2416     fp[3*ntrk[trk]  ] = ahit->X();
2417     fp[3*ntrk[trk]+1] = ahit->Y();
2418     fp[3*ntrk[trk]+2] = ahit->Z();
2419     ntrk[trk]++;
2420     ahit = NextHit2();
2421   }
2422
2423
2424
2425   //
2426   for(trk=0; trk<tracks; ++trk) {
2427     if(ntrk[trk]) {
2428       points = new AliPoints();
2429       points->SetMarkerColor(kYellow); //PH kYellow it the default in TPC
2430       points->SetMarkerSize(1);//PH Default size=1
2431       points->SetDetector(this);
2432       points->SetParticle(trk);
2433       points->SetPolyMarker(ntrk[trk],coor[trk],1); // Default style=1
2434       fPoints->AddAt(points,trk);
2435       delete [] coor[trk];
2436       coor[trk]=0;
2437     }
2438   }
2439   delete [] coor;
2440   delete [] ntrk;
2441   delete [] limi;
2442 }
2443
2444
2445 //_____________________________________________________________________________
2446 void AliTPC::LoadPoints3(Int_t)
2447 {
2448   //
2449   // Store x, y, z of all hits in memory
2450   // - only intersection point with pad row
2451   if (fTrackHits == 0) return;
2452   //
2453   Int_t nhits = fTrackHits->GetEntriesFast();
2454   if (nhits == 0) return;
2455   Int_t tracks = gAlice->GetMCApp()->GetNtrack();
2456   if (fPoints == 0) fPoints = new TObjArray(2*tracks);
2457   fPoints->Expand(2*tracks);
2458   AliHit *ahit;
2459   //
2460   Int_t *ntrk=new Int_t[tracks];
2461   Int_t *limi=new Int_t[tracks];
2462   Float_t **coor=new Float_t*[tracks];
2463   for(Int_t i=0;i<tracks;i++) {
2464     ntrk[i]=0;
2465     coor[i]=0;
2466     limi[i]=0;
2467   }
2468   //
2469   AliPoints *points = 0;
2470   Float_t *fp=0;
2471   Int_t trk;
2472   Int_t chunk=nhits/4+1;
2473   //
2474   // Loop over all the hits and store their position
2475   //
2476   ahit = FirstHit2(-1);
2477
2478   Int_t lastrow = -1;
2479   while (ahit){
2480     trk=ahit->GetTrack();
2481     Float_t  x[3]={ahit->X(),ahit->Y(),ahit->Z()};
2482     Int_t    index[3]={1,((AliTPChit*)ahit)->fSector,0};
2483     Int_t    currentrow = fTPCParam->GetPadRow(x,index) ;
2484     if (currentrow!=lastrow){
2485       lastrow = currentrow;
2486       //later calculate intersection point
2487       if(ntrk[trk]==limi[trk]) {
2488         //
2489         // Initialise a new track
2490         fp=new Float_t[3*(limi[trk]+chunk)];
2491         if(coor[trk]) {
2492           memcpy(fp,coor[trk],sizeof(Float_t)*3*limi[trk]);
2493           delete [] coor[trk];
2494         }
2495         limi[trk]+=chunk;
2496         coor[trk] = fp;
2497       } else {
2498         fp = coor[trk];
2499       }
2500       fp[3*ntrk[trk]  ] = ahit->X();
2501       fp[3*ntrk[trk]+1] = ahit->Y();
2502       fp[3*ntrk[trk]+2] = ahit->Z();
2503       ntrk[trk]++;
2504     }
2505     ahit = NextHit2();
2506   }
2507
2508   //
2509   for(trk=0; trk<tracks; ++trk) {
2510     if(ntrk[trk]) {
2511       points = new AliPoints();
2512       points->SetMarkerColor(kMagenta); //PH kYellow + 1 is kMagenta
2513       points->SetMarkerStyle(5);
2514       points->SetMarkerSize(0.2);
2515       points->SetDetector(this);
2516       points->SetParticle(trk);
2517       points->SetPolyMarker(ntrk[trk],coor[trk],30);
2518       fPoints->AddAt(points,tracks+trk);
2519       delete [] coor[trk];
2520       coor[trk]=0;
2521     }
2522   }
2523   delete [] coor;
2524   delete [] ntrk;
2525   delete [] limi;
2526 }
2527
2528
2529
2530 AliLoader* AliTPC::MakeLoader(const char* topfoldername)
2531 {
2532   //Makes TPC loader
2533   fLoader = new AliTPCLoader(GetName(),topfoldername);
2534   return fLoader;
2535 }
2536
2537 ////////////////////////////////////////////////////////////////////////
2538 AliTPCParam* AliTPC::LoadTPCParam(TFile *file) {
2539 //
2540 // load TPC paarmeters from a given file or create new if the object
2541 // is not found there
2542 // 12/05/2003 This method should be moved to the AliTPCLoader
2543 // and one has to decide where to store the TPC parameters
2544 // M.Kowalski
2545   char paramName[50];
2546   sprintf(paramName,"75x40_100x60_150x60");
2547   AliTPCParam *paramTPC=(AliTPCParam*)file->Get(paramName);
2548   if (paramTPC) {
2549     AliDebugClass(1,Form("TPC parameters %s found.",paramName));
2550   } else {
2551     AliWarningClass("TPC parameters not found. Create new (they may be incorrect)");
2552     paramTPC = new AliTPCParamSR;
2553   }
2554   return paramTPC;
2555
2556 // the older version of parameters can be accessed with this code.
2557 // In some cases, we have old parameters saved in the file but
2558 // digits were created with new parameters, it can be distinguish
2559 // by the name of TPC TreeD. The code here is just for the case
2560 // we would need to compare with old data, uncomment it if needed.
2561 //
2562 //  char paramName[50];
2563 //  sprintf(paramName,"75x40_100x60");
2564 //  AliTPCParam *paramTPC=(AliTPCParam*)in->Get(paramName);
2565 //  if (paramTPC) {
2566 //    cout<<"TPC parameters "<<paramName<<" found."<<endl;
2567 //  } else {
2568 //    sprintf(paramName,"75x40_100x60_150x60");
2569 //    paramTPC=(AliTPCParam*)in->Get(paramName);
2570 //    if (paramTPC) {
2571 //      cout<<"TPC parameters "<<paramName<<" found."<<endl;
2572 //    } else {
2573 //      cerr<<"TPC parameters not found. Create new (they may be incorrect)."
2574 //          <<endl;
2575 //      paramTPC = new AliTPCParamSR;
2576 //    }
2577 //  }
2578 //  return paramTPC;
2579
2580 }
2581
2582