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