TPC/AliTPCFast.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 //  fast TPC cluster simulation                                              //
  21 //                                                                          //
  22 ///////////////////////////////////////////////////////////////////////////////
  23
  24 #include <TParticle.h>
  25 #include <TVector.h>
  26
  27 #include "AliRunLoader.h"
  28 #include "AliRun.h"
  29 #include "AliMC.h"
  30 #include "AliTPC.h"
  31 #include "AliTPCParam.h"
  32 #include "AliTPCClustersArray.h"
  33 #include "AliTPCClustersRow.h"
  34 #include "AliTPCcluster.h"
  35 #include "AliComplexCluster.h"
  36 #include "AliTPCFast.h"
  37 #include "AliLog.h"
  38
  39 ClassImp(AliTPCFast)
  40
  41
  42 //_____________________________________________________________________________
  43 void AliTPCFast::Hits2Clusters(AliRunLoader* runLoader) const
  44 {
  45   //--------------------------------------------------------
  46   // TPC simple cluster generator from hits
  47   // obtained from the TPC Fast Simulator
  48   // The point errors are taken from the parametrization
  49   //--------------------------------------------------------
  50
  51   //-----------------------------------------------------------------
  52   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
  53   //-----------------------------------------------------------------
  54   // Adopted to Marian's cluster data structure by I.Belikov, CERN,
  55   // Jouri.Belikov@cern.ch
  56   //----------------------------------------------------------------
  57
  58   /////////////////////////////////////////////////////////////////////////////
  59   //
  60   //---------------------------------------------------------------------
  61   //   ALICE TPC Cluster Parameters
  62   //--------------------------------------------------------------------
  63
  64
  65
  66   // Cluster width in rphi
  67   const Float_t kACrphi=0.18322;
  68   const Float_t kBCrphi=0.59551e-3;
  69   const Float_t kCCrphi=0.60952e-1;
  70   // Cluster width in z
  71   const Float_t kACz=0.19081;
  72   const Float_t kBCz=0.55938e-3;
  73   const Float_t kCCz=0.30428;
  74
  75
  76   AliLoader* loader = runLoader->GetLoader("TPCLoader");
  77   if (!loader) {
  78     AliError("No TPC loader found");
  79     return;
  80   }
  81   if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
  82   AliRun* aliRun = runLoader->GetAliRun();
  83   if (!aliRun) {
  84     AliError("Couldn't get AliRun object");
  85     return;
  86   }
  87   AliTPC* tpc = (AliTPC*) aliRun->GetDetector("TPC");
  88   if (!tpc) {
  89     AliError("Couldn't get TPC detector");
  90     return;
  91   }
  92   AliTPCParam* param = tpc->GetParam();
  93   if (!param) {
  94     AliError("No TPC parameters available");
  95     return;
  96   }
  97
  98   //if(fDefaults == 0) SetDefaults();
  99
 100   Float_t sigmaRphi,sigmaZ,clRphi,clZ;
 101   //
 102   TParticle *particle; // pointer to a given particle
 103   AliTPChit *tpcHit; // pointer to a sigle TPC hit
 104   Int_t sector;
 105   Int_t ipart;
 106   Float_t xyz[5];
 107   Float_t pl,pt,tanth,rpad,ratio;
 108   Float_t cph,sph;
 109
 110   //---------------------------------------------------------------
 111   //  Get the access to the tracks
 112   //---------------------------------------------------------------
 113
 114   TTree *tH = loader->TreeH();
 115   if (tH == 0x0)
 116     AliFatal("Can not find TreeH in folder");
 117
 118   tpc->SetTreeAddress();
 119
 120   Stat_t ntracks = tH->GetEntries();
 121
 122   //Switch to the output file
 123
 124   if (loader->TreeR() == 0x0) loader->MakeTree("R");
 125
 126   AliDebug(1,Form("param->GetTitle() = %s",param->GetTitle()));
 127
 128   runLoader->CdGAFile();
 129   //param->Write(param->GetTitle());
 130
 131   AliTPCClustersArray carray;
 132   carray.Setup(param);
 133   carray.SetClusterType("AliTPCcluster");
 134   carray.MakeTree(loader->TreeR());
 135
 136   Int_t nclusters=0; //cluster counter
 137
 138   //------------------------------------------------------------
 139   // Loop over all sectors (72 sectors for 20 deg
 140   // segmentation for both lower and upper sectors)
 141   // Sectors 0-35 are lower sectors, 0-17 z>0, 17-35 z<0
 142   // Sectors 36-71 are upper sectors, 36-53 z>0, 54-71 z<0
 143   //
 144   // First cluster for sector 0 starts at "0"
 145   //------------------------------------------------------------
 146
 147   for(Int_t isec=0;isec<param->GetNSector();isec++){
 148     //MI change
 149     param->AdjustCosSin(isec,cph,sph);
 150
 151     //------------------------------------------------------------
 152     // Loop over tracks
 153     //------------------------------------------------------------
 154
 155     for(Int_t track=0;track<ntracks;track++){
 156       tpc->ResetHits();
 157       tH->GetEvent(track);
 158       //
 159       //  Get number of the TPC hits
 160       //
 161        tpcHit = (AliTPChit*)tpc->FirstHit(-1);
 162
 163       // Loop over hits
 164       //
 165        while(tpcHit){
 166
 167          if (tpcHit->fQ == 0.) {
 168            tpcHit = (AliTPChit*) tpc->NextHit();
 169            continue; //information about track (I.Belikov)
 170          }
 171         sector=tpcHit->fSector; // sector number
 172
 173        if(sector != isec){
 174          tpcHit = (AliTPChit*) tpc->NextHit();
 175          continue;
 176        }
 177         ipart=tpcHit->Track();
 178         particle=aliRun->GetMCApp()->Particle(ipart);
 179         pl=particle->Pz();
 180         pt=particle->Pt();
 181         if(pt < 1.e-9) pt=1.e-9;
 182         tanth=pl/pt;
 183         tanth = TMath::Abs(tanth);
 184         rpad=TMath::Sqrt(tpcHit->X()*tpcHit->X() + tpcHit->Y()*tpcHit->Y());
 185         ratio=0.001*rpad/pt; // pt must be in MeV/c - historical reason
 186
 187         //   space-point resolutions
 188
 189         sigmaRphi=AliTPCcluster::SigmaY2(rpad,tanth,pt);
 190         sigmaZ   =AliTPCcluster::SigmaZ2(rpad,tanth   );
 191
 192         //   cluster widths
 193
 194         clRphi=kACrphi-kBCrphi*rpad*tanth+kCCrphi*ratio*ratio;
 195         clZ=kACz-kBCz*rpad*tanth+kCCz*tanth*tanth;
 196
 197         // temporary protection
 198
 199         if(sigmaRphi < 0.) sigmaRphi=0.4e-3;
 200         if(sigmaZ < 0.) sigmaZ=0.4e-3;
 201         if(clRphi < 0.) clRphi=2.5e-3;
 202         if(clZ < 0.) clZ=2.5e-5;
 203
 204         //
 205
 206         //
 207         // smearing --> rotate to the 1 (13) or to the 25 (49) sector,
 208         // then the inaccuracy in a X-Y plane is only along Y (pad row)!
 209         //
 210         Float_t xprim= tpcHit->X()*cph + tpcHit->Y()*sph;
 211         Float_t yprim=-tpcHit->X()*sph + tpcHit->Y()*cph;
 212         xyz[0]=gRandom->Gaus(yprim,TMath::Sqrt(sigmaRphi));   // y
 213           Float_t alpha=(isec < param->GetNInnerSector()) ?
 214           param->GetInnerAngle() : param->GetOuterAngle();
 215           Float_t ymax=xprim*TMath::Tan(0.5*alpha);
 216           if (TMath::Abs(xyz[0])>ymax) xyz[0]=yprim;
 217         xyz[1]=gRandom->Gaus(tpcHit->Z(),TMath::Sqrt(sigmaZ)); // z
 218           if (TMath::Abs(xyz[1])>param->GetZLength()) xyz[1]=tpcHit->Z();
 219         xyz[2]=sigmaRphi;                                     // fSigmaY2
 220         xyz[3]=sigmaZ;                                        // fSigmaZ2
 221         xyz[4]=tpcHit->fQ;                                    // q
 222
 223         AliTPCClustersRow *clrow=carray.GetRow(sector,tpcHit->fPadRow);
 224         if (!clrow) clrow=carray.CreateRow(sector,tpcHit->fPadRow);
 225
 226         Int_t tracks[3]={tpcHit->Track(), -1, -1};
 227         AliTPCcluster cluster(tracks,xyz);
 228
 229         clrow->InsertCluster(&cluster); nclusters++;
 230
 231         tpcHit = (AliTPChit*)tpc->NextHit();
 232
 233
 234       } // end of loop over hits
 235
 236     }   // end of loop over tracks
 237
 238     Int_t nrows=param->GetNRow(isec);
 239     for (Int_t irow=0; irow<nrows; irow++) {
 240         AliTPCClustersRow *clrow=carray.GetRow(isec,irow);
 241         if (!clrow) continue;
 242         carray.StoreRow(isec,irow);
 243         carray.ClearRow(isec,irow);
 244     }
 245
 246   } // end of loop over sectors
 247
 248   //  cerr<<"Number of made clusters : "<<nclusters<<"                        \n";
 249   loader->WriteRecPoints("OVERWRITE");
 250
 251
 252 } // end of function
 253
 254 //_________________________________________________________________
 255 void AliTPCFast::Hits2ExactClustersSector(AliRunLoader* runLoader,
 256                                           AliTPCClustersArray* clustersArray,
 257                                           Int_t isec) const
 258 {
 259   //--------------------------------------------------------
 260   //calculate exact cross point of track and given pad row
 261   //resulting values are expressed in "digit" coordinata
 262   //--------------------------------------------------------
 263
 264   //-----------------------------------------------------------------
 265   // Origin: Marian Ivanov  GSI Darmstadt, m.ivanov@gsi.de
 266   //-----------------------------------------------------------------
 267   //
 268   if (clustersArray==0){
 269     return;
 270   }
 271   AliLoader* loader = runLoader->GetLoader("TPCLoader");
 272   if (!loader) {
 273     AliError("No TPC loader found");
 274     return;
 275   }
 276   if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
 277   AliRun* aliRun = runLoader->GetAliRun();
 278   if (!aliRun) {
 279     AliError("Couldn't get AliRun object");
 280     return;
 281   }
 282   AliTPC* tpc = (AliTPC*) aliRun->GetDetector("TPC");
 283   if (!tpc) {
 284     AliError("Couldn't get TPC detector");
 285     return;
 286   }
 287   AliTPCParam* param = tpc->GetParam();
 288   if (!param) {
 289     AliError("No TPC parameters available");
 290     return;
 291   }
 292   //
 293   TParticle *particle; // pointer to a given particle
 294   AliTPChit *tpcHit; // pointer to a sigle TPC hit
 295   //  Int_t sector,nhits;
 296   Int_t ipart;
 297   const Int_t kcmaxhits=30000;
 298   TVector * xxxx = new TVector(kcmaxhits*4);
 299   TVector & xxx = *xxxx;
 300   Int_t maxhits = kcmaxhits;
 301   //construct array for each padrow
 302   for (Int_t i=0; i<param->GetNRow(isec);i++)
 303     clustersArray->CreateRow(isec,i);
 304
 305   //---------------------------------------------------------------
 306   //  Get the access to the tracks
 307   //---------------------------------------------------------------
 308
 309   TTree *tH = loader->TreeH();
 310   if (tH == 0x0)
 311     AliFatal("Can not find TreeH in folder");
 312
 313   tpc->SetTreeAddress();
 314
 315   Stat_t ntracks = tH->GetEntries();
 316   Int_t npart = aliRun->GetMCApp()->GetNtrack();
 317   //MI change
 318   TBranch * branch=0;
 319   if (tpc->GetHitType()>1) branch = tH->GetBranch("TPC2");
 320   else branch = tH->GetBranch("TPC");
 321
 322   //------------------------------------------------------------
 323   // Loop over tracks
 324   //------------------------------------------------------------
 325
 326   for(Int_t track=0;track<ntracks;track++){
 327     Bool_t isInSector=kTRUE;
 328     tpc->ResetHits();
 329      isInSector = tpc->TrackInVolume(isec,track);
 330     if (!isInSector) continue;
 331     //MI change
 332     branch->GetEntry(track); // get next track
 333     //
 334     //  Get number of the TPC hits and a pointer
 335     //  to the particles
 336     // Loop over hits
 337     //
 338     Int_t currentIndex=0;
 339     Int_t lastrow=-1;  //last writen row
 340
 341     //M.I. changes
 342
 343     tpcHit = (AliTPChit*)tpc->FirstHit(-1);
 344     while(tpcHit){
 345
 346       Int_t sector=tpcHit->fSector; // sector number
 347       if(sector != isec){
 348         tpcHit = (AliTPChit*) tpc->NextHit();
 349         continue;
 350       }
 351
 352       ipart=tpcHit->Track();
 353       if (ipart<npart) particle=aliRun->GetMCApp()->Particle(ipart);
 354
 355       //find row number
 356
 357       Float_t  x[3]={tpcHit->X(),tpcHit->Y(),tpcHit->Z()};
 358       Int_t    index[3]={1,isec,0};
 359       Int_t    currentrow = param->GetPadRow(x,index) ;
 360       if (currentrow<0) {tpcHit = (AliTPChit*)tpc->NextHit(); continue;}
 361       if (lastrow<0) lastrow=currentrow;
 362       if (currentrow==lastrow){
 363         if ( currentIndex>=maxhits){
 364           maxhits+=kcmaxhits;
 365           xxx.ResizeTo(4*maxhits);
 366         }
 367         xxx(currentIndex*4)=x[0];
 368         xxx(currentIndex*4+1)=x[1];
 369         xxx(currentIndex*4+2)=x[2];
 370         xxx(currentIndex*4+3)=tpcHit->fQ;
 371         currentIndex++;
 372       }
 373       else
 374         if (currentIndex>2){
 375           Float_t sumx=0;
 376           Float_t sumx2=0;
 377           Float_t sumx3=0;
 378           Float_t sumx4=0;
 379           Float_t sumy=0;
 380           Float_t sumxy=0;
 381           Float_t sumx2y=0;
 382           Float_t sumz=0;
 383           Float_t sumxz=0;
 384           Float_t sumx2z=0;
 385           Float_t sumq=0;
 386           for (Int_t index=0;index<currentIndex;index++){
 387             Float_t x,x2,x3,x4;
 388             x=x2=x3=x4=xxx(index*4);
 389             x2*=x;
 390             x3*=x2;
 391             x4*=x3;
 392             sumx+=x;
 393             sumx2+=x2;
 394             sumx3+=x3;
 395             sumx4+=x4;
 396             sumy+=xxx(index*4+1);
 397             sumxy+=xxx(index*4+1)*x;
 398             sumx2y+=xxx(index*4+1)*x2;
 399             sumz+=xxx(index*4+2);
 400             sumxz+=xxx(index*4+2)*x;
 401             sumx2z+=xxx(index*4+2)*x2;
 402             sumq+=xxx(index*4+3);
 403           }
 404           Float_t centralPad = (param->GetNPads(isec,lastrow)-1)/2;
 405           Float_t det=currentIndex*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumx*sumx4-sumx2*sumx3)+
 406             sumx2*(sumx*sumx3-sumx2*sumx2);
 407
 408           Float_t detay=sumy*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxy*sumx4-sumx2y*sumx3)+
 409             sumx2*(sumxy*sumx3-sumx2y*sumx2);
 410           Float_t detaz=sumz*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxz*sumx4-sumx2z*sumx3)+
 411             sumx2*(sumxz*sumx3-sumx2z*sumx2);
 412
 413           Float_t detby=currentIndex*(sumxy*sumx4-sumx2y*sumx3)-sumy*(sumx*sumx4-sumx2*sumx3)+
 414             sumx2*(sumx*sumx2y-sumx2*sumxy);
 415           Float_t detbz=currentIndex*(sumxz*sumx4-sumx2z*sumx3)-sumz*(sumx*sumx4-sumx2*sumx3)+
 416             sumx2*(sumx*sumx2z-sumx2*sumxz);
 417
 418           if (TMath::Abs(det)<0.00001){
 419              tpcHit = (AliTPChit*)tpc->NextHit();
 420             continue;
 421           }
 422
 423           Float_t y=detay/det+centralPad;
 424           Float_t z=detaz/det;
 425           Float_t by=detby/det; //y angle
 426           Float_t bz=detbz/det; //z angle
 427           sumy/=Float_t(currentIndex);
 428           sumz/=Float_t(currentIndex);
 429
 430           AliTPCClustersRow * row = (clustersArray->GetRow(isec,lastrow));
 431           if (row!=0) {
 432             AliComplexCluster* cl = new((AliComplexCluster*)row->Append()) AliComplexCluster ;
 433             //    AliComplexCluster cl;
 434             cl->fX=z;
 435             cl->fY=y;
 436             cl->fQ=sumq;
 437             cl->fSigmaX2=bz;
 438             cl->fSigmaY2=by;
 439             cl->fTracks[0]=ipart;
 440           }
 441           currentIndex=0;
 442           lastrow=currentrow;
 443         } //end of calculating cluster for given row
 444
 445
 446       tpcHit = (AliTPChit*)tpc->NextHit();
 447     } // end of loop over hits
 448   }   // end of loop over tracks
 449   //write padrows to tree
 450   for (Int_t ii=0; ii<param->GetNRow(isec);ii++) {
 451     clustersArray->StoreRow(isec,ii);
 452     clustersArray->ClearRow(isec,ii);
 453   }
 454   xxxx->Delete();
 455
 456 }
 457
 458
 459 void AliTPCFast::FindTrackHitsIntersection(AliRunLoader* runLoader,
 460                                            AliTPCClustersArray* clustersArray) const
 461 {
 462
 463   //
 464   //fill clones array with intersection of current point with the
 465   //middle of the row
 466   AliLoader* loader = runLoader->GetLoader("TPCLoader");
 467   if (!loader) {
 468     AliError("No TPC loader found");
 469     return;
 470   }
 471   if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
 472   AliRun* aliRun = runLoader->GetAliRun();
 473   if (!aliRun) {
 474     AliError("Couldn't get AliRun object");
 475     return;
 476   }
 477   AliTPC* tpc = (AliTPC*) aliRun->GetDetector("TPC");
 478   if (!tpc) {
 479     AliError("Couldn't get TPC detector");
 480     return;
 481   }
 482   AliTPCParam* param = tpc->GetParam();
 483   if (!param) {
 484     AliError("No TPC parameters available");
 485     return;
 486   }
 487
 488   Int_t sector;
 489   Int_t ipart;
 490
 491   const Int_t kcmaxhits=30000;
 492   TVector * xxxx = new TVector(kcmaxhits*4);
 493   TVector & xxx = *xxxx;
 494   Int_t maxhits = kcmaxhits;
 495
 496   //
 497   AliTPChit * tpcHit=0;
 498   tpcHit = (AliTPChit*)tpc->FirstHit2(-1);
 499   Int_t currentIndex=0;
 500   Int_t lastrow=-1;  //last writen row
 501
 502   while (tpcHit){
 503     if (tpcHit==0) continue;
 504     sector=tpcHit->fSector; // sector number
 505     ipart=tpcHit->Track();
 506
 507     //find row number
 508
 509     Float_t  x[3]={tpcHit->X(),tpcHit->Y(),tpcHit->Z()};
 510     Int_t    index[3]={1,sector,0};
 511     Int_t    currentrow = param->GetPadRow(x,index) ;
 512     if (currentrow<0) continue;
 513     if (lastrow<0) lastrow=currentrow;
 514     if (currentrow==lastrow){
 515       if ( currentIndex>=maxhits){
 516         maxhits+=kcmaxhits;
 517         xxx.ResizeTo(4*maxhits);
 518       }
 519       xxx(currentIndex*4)=x[0];
 520       xxx(currentIndex*4+1)=x[1];
 521       xxx(currentIndex*4+2)=x[2];
 522       xxx(currentIndex*4+3)=tpcHit->fQ;
 523       currentIndex++;
 524     }
 525     else
 526       if (currentIndex>2){
 527         Float_t sumx=0;
 528         Float_t sumx2=0;
 529         Float_t sumx3=0;
 530         Float_t sumx4=0;
 531         Float_t sumy=0;
 532         Float_t sumxy=0;
 533         Float_t sumx2y=0;
 534         Float_t sumz=0;
 535         Float_t sumxz=0;
 536         Float_t sumx2z=0;
 537         Float_t sumq=0;
 538         for (Int_t index=0;index<currentIndex;index++){
 539           Float_t x,x2,x3,x4;
 540           x=x2=x3=x4=xxx(index*4);
 541           x2*=x;
 542           x3*=x2;
 543           x4*=x3;
 544           sumx+=x;
 545           sumx2+=x2;
 546           sumx3+=x3;
 547           sumx4+=x4;
 548           sumy+=xxx(index*4+1);
 549           sumxy+=xxx(index*4+1)*x;
 550           sumx2y+=xxx(index*4+1)*x2;
 551           sumz+=xxx(index*4+2);
 552           sumxz+=xxx(index*4+2)*x;
 553           sumx2z+=xxx(index*4+2)*x2;
 554           sumq+=xxx(index*4+3);
 555         }
 556         Float_t centralPad = (param->GetNPads(sector,lastrow)-1)/2;
 557         Float_t det=currentIndex*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumx*sumx4-sumx2*sumx3)+
 558           sumx2*(sumx*sumx3-sumx2*sumx2);
 559
 560         Float_t detay=sumy*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxy*sumx4-sumx2y*sumx3)+
 561           sumx2*(sumxy*sumx3-sumx2y*sumx2);
 562         Float_t detaz=sumz*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxz*sumx4-sumx2z*sumx3)+
 563           sumx2*(sumxz*sumx3-sumx2z*sumx2);
 564
 565         Float_t detby=currentIndex*(sumxy*sumx4-sumx2y*sumx3)-sumy*(sumx*sumx4-sumx2*sumx3)+
 566           sumx2*(sumx*sumx2y-sumx2*sumxy);
 567         Float_t detbz=currentIndex*(sumxz*sumx4-sumx2z*sumx3)-sumz*(sumx*sumx4-sumx2*sumx3)+
 568           sumx2*(sumx*sumx2z-sumx2*sumxz);
 569
 570         Float_t y=detay/det+centralPad;
 571         Float_t z=detaz/det;
 572         Float_t by=detby/det; //y angle
 573         Float_t bz=detbz/det; //z angle
 574         sumy/=Float_t(currentIndex);
 575         sumz/=Float_t(currentIndex);
 576
 577         AliComplexCluster cl;
 578         cl.fX=z;
 579         cl.fY=y;
 580         cl.fQ=sumq;
 581         cl.fSigmaX2=bz;
 582         cl.fSigmaY2=by;
 583         cl.fTracks[0]=ipart;
 584
 585         AliTPCClustersRow * row = (clustersArray->GetRow(sector,lastrow));
 586         if (row!=0) row->InsertCluster(&cl);
 587         currentIndex=0;
 588         lastrow=currentrow;
 589       } //end of calculating cluster for given row
 590
 591   } // end of loop over hits
 592   xxxx->Delete();
 593
 594 }
 595