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
 256 //_________________________________________________________________
 257
 258
 259
 260 void  AliTPCFast::Hits2ExactClusters(AliRunLoader* runLoader) const{
 261
 262
 263
 264   AliLoader* loader = runLoader->GetLoader("TPCLoader");
 265   if (!loader) {
 266     AliError("No TPC loader found");
 267     return;
 268   }
 269   AliTPC* tpc = (AliTPC*) gAlice->GetDetector("TPC");
 270   if (!tpc) {
 271     AliError("Couldn't get TPC detector");
 272     return;
 273   }
 274   AliTPCParam* param = tpc->GetParam();
 275   if (!param) {
 276     AliError("No TPC parameters available");
 277     return;
 278   }
 279   //---------------------------------------------------------------
 280   //  Get the access to the tracks
 281   //---------------------------------------------------------------
 282
 283   TTree *tH = loader->TreeH();
 284   if (tH == 0x0) { runLoader->LoadHits("TPC"); loader->TreeH();}
 285
 286   tpc->SetTreeAddress();
 287
 288
 289   //Switch to the output file
 290
 291   if (loader->TreeR() == 0x0) loader->MakeTree("R");
 292
 293   AliDebug(1,Form("param->GetTitle() = %s",param->GetTitle()));
 294
 295   runLoader->CdGAFile();
 296   //param->Write(param->GetTitle());
 297
 298   AliTPCClustersArray carray;
 299   carray.Setup(param);
 300   carray.SetClusterType("AliTPCclusterMI");
 301   carray.MakeTree(loader->TreeR());
 302
 303   //------------------------------------------------------------
 304   // Loop over all sectors (72 sectors for 20 deg)
 305   //------------------------------------------------------------
 306
 307   for(Int_t isec=0;isec<param->GetNSector();isec++){
 308     Hits2ExactClustersSector(runLoader, &carray, isec);
 309   } // end of loop over sectors
 310   loader->WriteRecPoints("OVERWRITE");
 311 }
 312
 313
 314
 315
 316 //_________________________________________________________________
 317 void AliTPCFast::Hits2ExactClustersSector(AliRunLoader* runLoader,
 318                                           AliTPCClustersArray* clustersArray,
 319                                           Int_t isec) const
 320 {
 321   //--------------------------------------------------------
 322   //calculate exact cross point of track and given pad row
 323   //resulting values are expressed in "local" coordinata
 324   //the sigmaY2 and sigmaZ2 of cluster are the shape of cluster parameters
 325   //                        - thay are used later on for error parameterization
 326   //--------------------------------------------------------
 327
 328   //-----------------------------------------------------------------
 329   // Origin: Marian Ivanov  GSI Darmstadt, m.ivanov@gsi.de
 330   //-----------------------------------------------------------------
 331   //
 332   if (clustersArray==0){
 333     return;
 334   }
 335   AliLoader* loader = runLoader->GetLoader("TPCLoader");
 336   if (!loader) {
 337     AliError("No TPC loader found");
 338     return;
 339   }
 340   if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
 341   AliRun* aliRun = runLoader->GetAliRun();
 342   if (!aliRun) {
 343     AliError("Couldn't get AliRun object");
 344     return;
 345   }
 346   AliTPC* tpc = (AliTPC*) aliRun->GetDetector("TPC");
 347   if (!tpc) {
 348     AliError("Couldn't get TPC detector");
 349     return;
 350   }
 351   AliTPCParam* param = tpc->GetParam();
 352   if (!param) {
 353     AliError("No TPC parameters available");
 354     return;
 355   }
 356   //
 357   //
 358   AliTPChit *tpcHit; // pointer to a sigle TPC hit
 359   //  Int_t sector,nhits;
 360   Int_t ipart;
 361   const Int_t kcmaxhits=30000;
 362   TVector * xxxx = new TVector(kcmaxhits*4);
 363   TVector & xxx = *xxxx;
 364   Int_t maxhits = kcmaxhits;
 365   //construct array for each padrow
 366   for (Int_t i=0; i<param->GetNRow(isec);i++)
 367     clustersArray->CreateRow(isec,i);
 368
 369   //---------------------------------------------------------------
 370   //  Get the access to the tracks
 371   //---------------------------------------------------------------
 372
 373   TTree *tH = loader->TreeH();
 374   if (tH == 0x0)
 375     AliFatal("Can not find TreeH in folder");
 376
 377   tpc->SetTreeAddress();
 378
 379   Stat_t ntracks = tH->GetEntries();
 380   //MI change
 381   TBranch * branch=0;
 382   if (tpc->GetHitType()>1) branch = tH->GetBranch("TPC2");
 383   else branch = tH->GetBranch("TPC");
 384
 385   //------------------------------------------------------------
 386   // Loop over tracks
 387   //------------------------------------------------------------
 388
 389   for(Int_t track=0;track<ntracks;track++){
 390     Bool_t isInSector=kTRUE;
 391     tpc->ResetHits();
 392     isInSector = tpc->TrackInVolume(isec,track);
 393     if (!isInSector) continue;
 394     //MI change
 395     branch->GetEntry(track); // get next track
 396     //
 397     // Loop over hits
 398     //
 399     Int_t currentIndex=0;
 400     Int_t lastrow=-1;  //last writen row
 401
 402     //M.I. changes
 403
 404     tpcHit = (AliTPChit*)tpc->FirstHit(-1);
 405     while(tpcHit){
 406
 407       Int_t sector=tpcHit->fSector; // sector number
 408       if(sector != isec){
 409         tpcHit = (AliTPChit*) tpc->NextHit();
 410         continue;
 411       }
 412
 413       ipart=tpcHit->Track();
 414
 415       //find row number
 416
 417       Float_t  x[3]={tpcHit->X(),tpcHit->Y(),tpcHit->Z()};
 418       Int_t    index[3]={1,isec,0};
 419       Int_t    currentrow = param->GetPadRow(x,index);
 420
 421       if (currentrow<0) {tpcHit = (AliTPChit*)tpc->NextHit(); continue;}
 422       if (lastrow<0) lastrow=currentrow;
 423       //
 424       //
 425       if (currentrow!=lastrow){
 426         if (currentIndex>2){
 427           Float_t sumx=0;
 428           Float_t sumx2=0;
 429           Float_t sumx3=0;
 430           Float_t sumx4=0;
 431           Float_t sumy=0;
 432           Float_t sumxy=0;
 433           Float_t sumx2y=0;
 434           Float_t sumz=0;
 435           Float_t sumxz=0;
 436           Float_t sumx2z=0;
 437           Float_t sumq=0;
 438           for (Int_t index=0;index<currentIndex;index++){
 439             Float_t x,x2,x3,x4;
 440             x=x2=x3=x4=xxx(index*4);
 441             x2*=x;
 442             x3*=x2;
 443             x4*=x3;
 444             sumx+=x;
 445             sumx2+=x2;
 446             sumx3+=x3;
 447             sumx4+=x4;
 448             sumy+=xxx(index*4+1);
 449             sumxy+=xxx(index*4+1)*x;
 450             sumx2y+=xxx(index*4+1)*x2;
 451             sumz+=xxx(index*4+2);
 452             sumxz+=xxx(index*4+2)*x;
 453             sumx2z+=xxx(index*4+2)*x2;
 454             sumq+=xxx(index*4+3);
 455           }
 456           Float_t det=currentIndex*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumx*sumx4-sumx2*sumx3)+
 457             sumx2*(sumx*sumx3-sumx2*sumx2);
 458
 459           Float_t detay=sumy*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxy*sumx4-sumx2y*sumx3)+
 460             sumx2*(sumxy*sumx3-sumx2y*sumx2);
 461           Float_t detaz=sumz*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxz*sumx4-sumx2z*sumx3)+
 462             sumx2*(sumxz*sumx3-sumx2z*sumx2);
 463
 464           Float_t detby=currentIndex*(sumxy*sumx4-sumx2y*sumx3)-sumy*(sumx*sumx4-sumx2*sumx3)+
 465             sumx2*(sumx*sumx2y-sumx2*sumxy);
 466           Float_t detbz=currentIndex*(sumxz*sumx4-sumx2z*sumx3)-sumz*(sumx*sumx4-sumx2*sumx3)+
 467             sumx2*(sumx*sumx2z-sumx2*sumxz);
 468
 469           if (TMath::Abs(det)<0.00000000000000001){
 470              tpcHit = (AliTPChit*)tpc->NextHit();
 471             continue;
 472           }
 473
 474           Float_t y=detay/det;
 475           Float_t z=detaz/det;
 476           Float_t by=detby/det; //y angle
 477           Float_t bz=detbz/det; //z angle
 478           sumy/=Float_t(currentIndex);
 479           sumz/=Float_t(currentIndex);
 480
 481
 482           //
 483           //
 484           Float_t sign = (tpcHit->Z()>0)? 1.:-1.;
 485           y = (y+0.5)*param->GetPadPitchWidth(isec);
 486           z = z*param->GetZWidth();
 487           //
 488           // y expected shape
 489           Double_t sigmay2 = z*fParam->GetDiffL();
 490           sigmay2 *= sigmay2;
 491           Float_t angulary = by*param->GetPadPitchLength(isec,lastrow);
 492           angulary*=angulary;
 493           angulary/=12;
 494           sigmay2 +=angulary+0.25*param->GetPadPitchWidth(isec)*param->GetPadPitchWidth(isec);
 495           //
 496           // z expected shape
 497           Double_t sigmaz2 = z*fParam->GetDiffT();
 498           sigmaz2 *= sigmaz2;
 499           Float_t angularz = bz*param->GetPadPitchLength(isec,lastrow);
 500           angularz*=angularz;
 501           angularz/=12;
 502           sigmaz2 +=angularz+0.25;
 503           //
 504           sigmaz2 = TMath::Min(sigmaz2,1.);
 505           sigmay2 = TMath::Min(sigmay2,1.);
 506           //
 507           //
 508           z = sign*(param->GetZLength() - z);
 509           if (TMath::Abs(z)< param->GetZLength()-1){
 510             AliTPCClustersRow * row = (clustersArray->GetRow(isec,lastrow));
 511             if (row!=0) {
 512               AliTPCclusterMI* cl = new((AliTPCclusterMI*)row->Append()) AliTPCclusterMI ;
 513               //          AliTPCclusterMI cl;
 514               cl->SetZ(z);
 515               cl->SetY(y);
 516               cl->SetQ(sumq);
 517               if (TMath::Abs(sumq)<1000){
 518                 cl->SetMax(Short_t(sumq));
 519               }
 520               else{
 521                 cl->SetMax(0);
 522               }
 523               cl->SetSigmaZ2(sigmaz2);
 524               cl->SetSigmaY2(sigmay2);
 525               cl->SetLabel(ipart,0);
 526               cl->SetLabel(-1,1);
 527               cl->SetLabel(-1,2);
 528               cl->SetType(0);
 529             }
 530           }
 531         } //end of calculating cluster for given row
 532         currentIndex=0;
 533         lastrow=currentrow;
 534       }  //  end of crossing rows
 535       //
 536       if ( currentIndex>=maxhits){
 537         maxhits+=kcmaxhits;
 538         xxx.ResizeTo(4*maxhits);
 539       }
 540       xxx(currentIndex*4)=x[0];
 541       xxx(currentIndex*4+1)=x[1];
 542       xxx(currentIndex*4+2)=x[2];
 543       xxx(currentIndex*4+3)=tpcHit->fQ;
 544       currentIndex++;
 545       tpcHit = (AliTPChit*)tpc->NextHit();
 546     } // end of loop over hits
 547   }   // end of loop over tracks
 548   //write padrows to tree
 549   for (Int_t ii=0; ii<param->GetNRow(isec);ii++) {
 550     clustersArray->StoreRow(isec,ii);
 551     clustersArray->ClearRow(isec,ii);
 552   }
 553   xxxx->Delete();
 554
 555 }
 556
 557