Put back the two bellows in front of the absorber.
[u/mrichter/AliRoot.git] / TPC / AliTPCDigitsH.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 /*
17 $Log$
18 Revision 1.1.4.2  2000/04/10 11:37:42  kowal2
19
20 Digits handling in a new data structure
21
22 */
23
24 //-----------------------------------------------------------------------------
25 //
26 //
27 //  Author:   Marian Ivanov
28 //
29 //  Implementation of class TTPCDigitsH
30 //
31 //-----------------------------------------------------------------------------
32 ////////////////////////////////////////////////
33 //  Manager class for AliTPCDigitsH             //
34 ////////////////////////////////////////////////
35  
36 #include <iostream.h>
37 #include "TMath.h"
38 // other include files follow here
39 #include "TFile.h"
40 //*KEEP,TFile.
41 #include "TROOT.h"
42 #include "TGraph.h"
43 #include "AliRun.h"
44 #include "AliDisplay.h"
45 #include "TPad.h"
46 #include "TCanvas.h"
47 #include "TCanvasImp.h"
48 #include "TPaveText.h"
49 //*KEEP,TH1.
50 #include "TH1.h"
51 //*KEEP,TH2
52 #include "AliH2F.h"
53 //*KEEP,TF2.
54 #include "TF2.h"
55 //*KEEP,TClonesArray,T=C++.
56 #include "TClonesArray.h"
57 #include "TTree.h"
58 //GALICE includes
59 #include "GParticle.h"
60 #include "AliTPC.h"
61 #include "AliTPCParam.h"
62 #include "AliTPCD.h"
63 #include "AliTPCDigitsH.h"
64
65
66 R__EXTERN TSystem *  gSystem;
67 R__EXTERN AliRun *   gAlice;
68
69
70
71 ClassImp(AliTPCDigitsH)
72
73 AliTPCDigitsH::AliTPCDigitsH() 
74 {
75 //Begin_Html
76 /*
77 <img src="gif/TPCDigitsH.gif">
78 */
79 //End_Html
80  
81   fsec = 1;
82   frow = 1;
83   fpad = 1;
84   fTimeN = 500;
85   fTimeStart = 0;
86   fTimeStop = 500;  
87   fOccuN    =25;
88   fbDelHisto = kFALSE;
89   fEventN = 0;  //event nuber connected to digit tree
90   fThreshold  = 5;
91   
92   fParticles = 0;
93   fDigits= 0; 
94
95   fDParam  =0;
96   
97   fbIOState  = kFALSE;
98   fbDigState = kFALSE;
99 }
100
101 AliTPCDigitsH::~AliTPCDigitsH() 
102 {
103
104 }
105
106 AliTPCDigitsH::AliTPCDigitsH(const AliTPCDigitsH &) 
107 {
108 }
109
110 AliTPCDigitsH & AliTPCDigitsH::operator = (const AliTPCDigitsH &) 
111 {
112    return *this;
113 }
114 void  AliTPCDigitsH::SetDParam(AliTPCD * dig)
115 {
116   if (dig!=0) 
117     {
118       fDParam =dig;
119       fTPCParam = &(fDParam->GetParam());
120     }
121   else
122     {
123       fTPCParam = 0;
124       fDParam =0;
125     }
126
127
128 AliTPCParam *&  AliTPCDigitsH::GetParam()
129 {
130   return fTPCParam;
131 }
132
133
134 void AliTPCDigitsH::CloseFiles()
135 {
136  if (fin) 
137     {
138       fin->Close();
139       //try 
140       // {
141       delete fin; 
142       //  //    }
143       //  //catch(...)
144       //  //    {
145       //  cout<<"FIN Delete error. Contact autor of root \n";
146       //        };      
147       fin = 0;
148     };
149  if (fout) 
150     {
151       fout->Close();  
152       //try 
153       //{
154       //  delete fout;
155       //}
156       //      catch(...)
157       //{
158       //  cout<<"Out Delete error. Contact autor of root \n";
159       //};      
160       fout =0;
161     };
162  fbIOState  = kFALSE;
163  fbDigState = kFALSE;
164 }
165
166 Bool_t AliTPCDigitsH::SetIO(const char *  inFile, const char* outFile )
167 {
168   ///   Set input and output file 
169   ///   control permisions and so on
170   ///   if all is OK then set flag fbIOState = kTRUE
171  fbIOState = kFALSE; 
172  TString s1=inFile;
173  TString s2=outFile; 
174  //  important ---- it close previious open file if this file exist
175  if (fin) 
176    {
177      if (fin == fout) fout = NULL;
178      fin->Close();
179      //     try
180      //{
181      delete fin;
182      //}
183      //catch(...)
184      //{
185      //  cout<<"Fin  Delete error. Contact autor of root \n";
186      //};
187      fin = 0;
188    };   
189  //  important ---- it close previous open file if this file exist
190  if (fout) 
191    {
192      fout->Close(); 
193      //try
194      //{
195      delete fout;
196      //}
197      //catch(...)
198      //{
199      //  cout<<"Fout  Delete error. Contact autor of root \n";
200      //};
201      fout = 0;
202    };
203
204  //close the files if exist in root enviroment
205  TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject((char*)inFile);
206  if (file) file->Close();
207  file = (TFile*)gROOT->GetListOfFiles()->FindObject((char*)outFile);
208  if (file) file->Close();
209  //if input file is the output file 
210   if (s1 == s2 )
211     {
212       if (gSystem->AccessPathName((char*)inFile, kWritePermission ) == kFALSE)
213          fin  = new TFile((char*)inFile,"UPDATE","trees with digits");                
214       if (!(fin)) 
215        {
216          cout<<"Input file couldn't be open for writing \n";
217          cout<<"Loook if the file exiest or if you have permision to write \n";        
218          return  kFALSE;
219        }
220       else
221         { 
222           fout = fin;
223           fbIOState = kTRUE;
224           return kTRUE;
225         } 
226     }
227   //open input file 
228   if (gSystem->AccessPathName((char*)inFile, kReadPermission ) == kFALSE)      
229     fin = new TFile((char*)inFile,"UPDATE","trees with digits");
230   else
231     {
232       cout<<"Input file couldn't be open\n";
233       cout<<"Maybe the file is not open \n";        
234       return kFALSE ;
235     }          
236   if (!(fin)) 
237     {
238       cout<<"Input file couldn't be open\n";
239       cout<<"Probably not root file \n";     
240       return kFALSE ;
241     }
242
243   //open output file  
244   if (gSystem->AccessPathName((char*)outFile, kWritePermission) == kFALSE)    
245     fout = new TFile((char *)outFile,"UPDATE");   
246   else  
247     fout = new TFile((char *)outFile,"NEW");     
248    
249   if (!(fout)) 
250     {
251       cout<<"Output  file couldn't be open\n";
252       return kFALSE ;
253     }
254   //if input and output file is OK set state variable to true
255   fbIOState = kTRUE;  
256   return kTRUE;
257 }
258
259 Bool_t AliTPCDigitsH::SetEventN(Int_t EventN=0)
260 {
261   if (!(fin))
262     {
263       cout<<"Warning: Input file not open !!! \n";   
264       return kFALSE;
265     }
266   fin->cd();
267   if (EventN>-1) fEventN = EventN;
268   fParticles = 0;
269   if (gAlice) 
270     {
271       delete gAlice;
272       gAlice =0;
273     }
274   gAlice = (AliRun*) fin->Get("gAlice");   
275   if (gAlice == 0 )
276     cout<<" Warning : AliRun objects not found in input file \n.";
277   else
278     {
279       gAlice->GetEvent(fEventN);
280       fParticles = gAlice->Particles();
281     }
282   if (fParticles == 0) 
283     return kFALSE;
284   else
285     return kTRUE;
286 }
287
288
289
290 Bool_t AliTPCDigitsH::SetTree(Int_t eventn  )
291 {
292    fbDigState = kFALSE;
293    ftree = 0;
294    if (  fbIOState == kFALSE)
295      {
296        cout<<"IO files not adjusted \n FIX IT!!!";
297        return kFALSE;
298      }
299    fin->cd();
300    if (fDParam->SetTree(eventn)==0){
301      fbDigState = kFALSE;
302      cout<<"Input tree doesn't exist !!! \n";
303      return kFALSE;
304    }
305    fDigits = fDParam->GetArray();
306    ftree = fDParam->GetTree();
307    fbDigState=kTRUE;
308    return kTRUE; 
309 }
310
311 void AliTPCDigitsH::SetSecRowTime(Int_t sec , Int_t row , Int_t TimeN, Float_t TimeStart, Float_t TimeStop )
312 {
313   fsec = sec;
314   frow = row;
315   fTimeN =TimeN;
316   fTimeStart = TimeStart;
317   fTimeStop = TimeStop;
318 }
319
320 void   AliTPCDigitsH::SetParticles(Int_t sec = -1, Int_t row = -1 ,
321                        Int_t size1 = 30000,Int_t size2=300,
322                        Bool_t all=kTRUE )
323 {
324   if (sec>0) fsec = sec;
325   if (row>-1) frow =row;
326   char s[80];
327   char sh[80];
328   //  create particles histograms
329   sprintf(s,"Sector %d   Row %d\n",fsec,frow);  
330   sprintf(sh,"Particles%d_%d",fsec,frow);   
331   fHParticles = new TH1F(sh,s,size1,4,size1); 
332
333   sprintf(s,"Sector %d   Row %d\n",fsec,frow);  
334   sprintf(sh,"All particles%d_%d",fsec,frow);   
335   fHAllP = new TH1F(sh,s,200,1,25);  
336
337   sprintf(s,"Sector %d   Row %d\n",fsec,frow);  
338   sprintf(sh,"Secondary Particles%d_%d",fsec,frow);   
339   fHSecondaryP = new TH1F(sh,s,200,1,25); 
340
341   if (!(fin)) 
342   {
343     cout<<"Input  file not open, open file before \n";
344   }
345   else
346 {
347   fin->cd();
348   Int_t sectors_by_rows=(Int_t)ftree->GetEntries();
349   GParticle * particle;
350   //  loop over all sectors and rows
351   for (Int_t n=0; n<sectors_by_rows; n++) 
352     {
353     if (!ftree->GetEvent(n)) continue;
354     AliTPCdigit *dig=(AliTPCdigit*)fDigits->UncheckedAt(0);
355     if (fsec  < dig->fSector) break;
356     if (fsec != dig->fSector) continue;
357     if (frow != dig->fPadRow) continue;
358     
359     Int_t ndigits=fDigits->GetEntriesFast();
360     //loop over all digits  in sector pad
361     for (Int_t ndig=0; ndig<ndigits; ndig++) 
362         {
363       Float_t x,y;
364         dig=(AliTPCdigit*)fDigits->UncheckedAt(ndig);
365         fHParticles->Fill(dig->fTracks[0]);
366         //     get pointer to particle information and fill All and secondary histo
367         particle = (GParticle*) fParticles->UncheckedAt(dig->fTracks[0]);
368       fHAllP->Fill(particle->GetKF());
369       //                Int_t id = particle->GetKF();
370       x = (Float_t)particle->GetVx();
371       y = (Float_t)particle->GetVy();
372       if ( (x*x+y*y ) > 1 )
373           fHSecondaryP->Fill(particle->GetKF());
374         if (all==kTRUE)
375           {
376             fHParticles->Fill(dig->fTracks[1]);
377             particle = (GParticle*) fParticles->UncheckedAt(dig->fTracks[1]);
378             fHAllP->Fill(particle->GetKF()); 
379             x = (Float_t)particle->GetVx();
380             y = (Float_t)particle->GetVy();
381             if ( (x*x+y*y ) > 1 )           
382               fHSecondaryP->Fill(particle->GetKF());
383             fHParticles->Fill(dig->fTracks[2]);
384             particle = (GParticle*) fParticles->UncheckedAt(dig->fTracks[2]);
385             fHAllP->Fill(particle->GetKF());
386             x = (Float_t)particle->GetVx();
387             y = (Float_t)particle->GetVy();
388             if ( (x*x+y*y ) > 1 )           
389               fHSecondaryP->Fill(particle->GetKF());
390           }
391         }
392       }
393     //make histogram with multiplicity
394       char s[80];
395       char sh[80];
396       if (all==kTRUE)
397         sprintf(s,"Number of AliDigits over threshold  per one track in sector %d   Row %d\n (all three most important track recorded)",fsec,frow);
398       else
399         sprintf(s,"Number of AliDigits over threshold per sector %d   Row %d\n (only most important track)",fsec,frow);  
400       sprintf(sh,"His_%d_%d",fsec,frow);   
401       fHPartMultiplicity = new TH1F(sh,s,size2,1,size2);
402       for (Int_t i=1;i<size1;i++)
403         {
404           Int_t mul=Int_t(fHParticles->GetBinContent(i));
405           if (mul>0) fHPartMultiplicity->Fill(mul);
406         }
407       if (fout)
408         {
409           fout->cd();          
410           fHParticles->Write();     
411           fHPartMultiplicity->Write();     
412           //fHSecondaryP->Write(); 
413           //fHAllP->Write();
414         }
415   }
416 }
417
418
419 void AliTPCDigitsH::Anal()
420 {
421   if (fbIOState == kFALSE)
422     {
423       cout<<"Input output problem. \n Do you initialize IO files ? \n";
424       return;
425     }
426   if (fbDigState == kFALSE)
427     {
428       cout<<"Input file doesn't enhalt selected tree \n";
429       return;
430     }  
431   fout->cd();
432  
433 //if we dont want let histogram in memory then we delete old histogram 
434   if ( (fH2Digit) && (fbDelHisto == kTRUE) )  
435     {
436       //      try 
437       //{
438           // delete fH2Digit;
439       //}
440       //      catch(...)
441       //{
442       //  cout<<"Delete error. Contact autor of root \n";
443       //};
444       fH2Digit = 0;
445     }
446   char s[80];
447   char sh[80];
448   sprintf(s,"Sector %d   Row %d\n",fsec,frow);  
449   sprintf(sh,"h%d_%d",fsec,frow);   
450
451   if  ( (fout) && (fbDelHisto == kFALSE) )
452     {
453       fH2Digit = (AliH2F *) fout->Get(sh);
454       if (fH2Digit) return;     
455     } 
456
457   Int_t n_of_pads =fTPCParam->GetNPads(fsec,frow);
458      
459   fH2Digit = new AliH2F(sh, s, fTimeN, fTimeStart, fTimeStop, n_of_pads, 0, n_of_pads-1);
460  
461   if (!(fout)) 
462   {
463     cout<<"Input  file not open, open file before \n";
464   }
465   else
466   {
467     //fin->cd();
468     Int_t sectors_by_rows=(Int_t)ftree->GetEntries();
469     //loop over all sectors and rows
470     for (Int_t n=0; n<sectors_by_rows; n++) {
471       if (!ftree->GetEvent(n)) continue;
472       AliTPCdigit *dig=(AliTPCdigit*)fDigits->UncheckedAt(0);
473       if (fsec  < dig->fSector) break;
474       if (fsec != dig->fSector) continue;
475       if (frow != dig->fPadRow) continue;
476       
477       Int_t ndigits=fDigits->GetEntriesFast();
478       //loop over all digits  in sector pad
479       for (Int_t ndig=0; ndig<ndigits; ndig++) {
480         dig=(AliTPCdigit*)fDigits->UncheckedAt(ndig);
481         fH2Digit->Fill(dig->fTime,dig->fPad,dig->fSignal);
482       }
483     }
484     if (fout) fout->cd();          
485     fH2Digit->Write();          
486   }
487
488 }
489 void  AliTPCDigitsH::Draw(Option_t * opt1 ="cont1"  , Option_t * opt2 = "error",
490                         Option_t * opt3 = "L" )
491 {
492   TString o1 = opt1;
493   o1.ToLower();
494
495   TString o2 = opt2;
496   o2.ToLower();
497
498   TString o3 = opt3;
499   o3.ToLower();
500
501   fcanvas  = new TCanvas("dh","Digits Histograms",700,900);
502   
503   if (fTitle) 
504     {
505       //       try 
506       //{
507           //delete fTitle;
508       //}
509       //      catch(...)
510       //{
511       //  cout<<"Delete error. Contact autor of root \n";
512       //};
513       fTitle = 0;
514     } 
515   fTitle = new TPaveText(0.2,0.96,0.8,0.995);
516   fTitle->AddText("Occupancy calculation for TPC");
517   fTitle->Draw();
518   
519   fpad1 = new TPad("pad1","",0.05,0.7,0.95,0.95,21);
520   fpad1->Draw();
521   fpad2 = new TPad("pad2","",0.05,0.4,0.95,0.65,21);
522   fpad2->Draw();
523   fpad3 = new TPad("pad3","",0.05,0.05,0.95,0.35,21);
524   fpad3->Draw();
525
526   fpad1->cd();
527   //  pad1->TPaveText::title.SetSize(0.1);
528   if (fH2Digit) 
529     { 
530       fH2DigitBW->Draw(o1);
531       fH2DigitBW->SetXTitle("time bin");
532       fH2DigitBW->SetYTitle("pad number");
533     }
534   fpad2->cd();
535   fH1Occu->Fit("pol0");
536   if (fH1Occu) 
537     {
538       fH1Occu->Draw(o2);
539       fH1Occu->SetXTitle("time bin");
540       fH1Occu->SetYTitle("occupancy");
541     }
542   fpad3->cd();  
543   if (fH1Digit)
544     {
545       fH1Digit->Draw(o3);
546       fH1Digit->SetXTitle("time bin");
547       fH1Digit->SetYTitle("ADC amplitude ");
548     }
549
550 };
551
552 void  AliTPCDigitsH::DeleteHisto(const Text_t *namecycle)
553 {
554   if (fout) fout->Delete(namecycle);
555 }
556      
557 void  AliTPCDigitsH::SetHisto(Int_t pad = 1 )
558 {
559
560   Int_t n_of_pads = fTPCParam->GetNPads(fsec,frow);
561   if (pad > (n_of_pads-1)) 
562     {
563       cout<<"Pad number is greater then actula number of pads in thi row \n";
564       cout<<"Noch einmal \n";
565       return;
566     }
567       
568   fpad = pad;
569   Anal();
570
571   //  if ( (fH1Digit) && (fbDelHisto == kTRUE)) 
572   //    {
573   //      try 
574   //        {
575   //      // delete fH1Digit;
576   //    }
577   //      catch(...)
578   //    {
579   //      cout<<"Delete error. Contact autor of root \n";
580   //    };
581   //      fH1Digit = 0;
582   //    }
583
584   char s[80];
585   char sh[80];
586   sprintf(s,"example sector %d   Row %d  Pad %d",fsec,frow,fpad);  
587   sprintf(sh,"h%d_%d_%d",fsec,frow,fpad);   
588  fH2DigitBW = new AliH2F("bw", "", fTimeN, fTimeStart, fTimeStop, n_of_pads, 0, n_of_pads-1);
589   fH1Digit = new TH1F(sh,s,fTimeN,fTimeStart,fTimeStop);
590
591   for (Int_t i = 0;i<fTimeN;i++)
592     {
593       Int_t index = fH2Digit->GetBin(i,pad);
594       Float_t weight = fH2Digit->GetBinContent(index);
595       fH1Digit->Fill(i,weight);
596     };  
597
598  
599   sprintf(s,"Occupancy in sector %d   Row %d  threshold = %d",fsec,frow,fThreshold);  
600   sprintf(sh,"hoccu%d_%d_%d",fsec,frow,fpad);   
601   fH1Occu = new TH1F(sh,s,fOccuN,fTimeStart,fTimeStop);
602   
603   for (Int_t i = 0;i<fOccuN;i++)
604     {
605       Int_t over =0;  
606       Int_t all  =0;    
607       for (int itime = i*(fTimeN/fOccuN); itime<(i+1)*(fTimeN/fOccuN);itime++)
608         {
609           for (Int_t ipad = 0; ipad < n_of_pads; ipad++)
610             {
611               Int_t index = fH2Digit->GetBin(itime,ipad);
612               if ( (ipad>3) && ((ipad+3)<n_of_pads)){
613                 all++;
614                 if (fH2Digit->GetBinContent(index) >fThreshold) over++ ;
615               }
616               if (fH2Digit->GetBinContent(index) >fThreshold) 
617                 fH2DigitBW->Fill(itime,ipad,1);
618               else
619                 fH2DigitBW->Fill(itime,ipad,0);
620             }
621         }
622      Float_t occu = ((Float_t)over) /((Float_t) (all)); 
623      //     Float_t time =   ((fTimeStop-fTimeStart)/fOccuN)*i+fTimeStart;
624    
625      fH1Occu->SetBinContent(i,occu); 
626      // Int_t index = fH1Occu->GetBin(i);
627      Float_t error = sqrt( ((Float_t) ((over)/25+1)) )/((Float_t)(all)/25.);
628      fH1Occu->SetBinError(i,error);
629
630     };  
631  
632
633 }
634   
635
636
637 void AliTPCDigitsH::Streamer(TBuffer & R__b)
638 {
639   if (R__b.IsReading()) {
640     //      Version_t R__v = R__b.ReadVersion();
641    } else {
642       R__b.WriteVersion(AliTPCDigitsH::IsA());    
643    } 
644 }