New function Hits2ExactClusters (M.Ivanov)
[u/mrichter/AliRoot.git] / 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