TPC/FindKrClusters.C

   1 /// \file FindKrClusters.C
   2 /// \brief Macro to create array of clusters from TPC digits
   3 ///
   4 /// Input files:
   5 ///
   6 /// * galice.root
   7 /// * digits.root: file with digits - usualy use link to galice.root
   8 ///   in splitted mode - neccesary to create link to proper file
   9 ///
  10 /// Output file:
  11 ///
  12 /// * TPC.RecPoints.root
  13 ///
  14 /// Warning - if cluster file AliTPCclusters.root already exist - macro exit and don't produce anything
  15 ///
  16 /// \author A.Matyja amatyja@cern.ch
  17
  18 #ifndef __CINT__
  19 #include <iostream.h>
  20 #include "AliRun.h"
  21 #include "AliTPCv4.h"
  22 #include "AliTPCParam.h"
  23 #include "AliTPCclusterKr.h"
  24 #include "AliTPCclustererKr.h"
  25 #include "TFile.h"
  26 #include "TStopwatch.h"
  27 #include "TTree.h"
  28 #endif
  29
  30 Int_t FindKrClusters(){
  31
  32   // Load DataBase
  33   //char *ocdbpath ="local:///afs/cern.ch/alice/tpctest/OCDB";
  34   //char *ocdbpath ="local:///home/matyja/baza/OCDB";
  35   char *ocdbpath ="local:///data/baza/OCDB";
  36   if (ocdbpath==0){
  37     ocdbpath="alien://folder=/alice/data/2007/LHC07w/OCDB/";
  38   }
  39   printf("OCDB PATH = %s\n",ocdbpath);
  40   AliCDBManager * man = AliCDBManager::Instance();
  41   man->SetDefaultStorage(ocdbpath);
  42   man->SetRun(0);
  43
  44   AliRunLoader* rl = AliRunLoader::Open("galice.root");
  45   if (rl == 0x0) {
  46     cerr<<"Can not open session"<<endl;
  47     return 1;
  48   }
  49
  50   AliTPCLoader *tpcl = (AliTPCLoader*)rl->GetLoader("TPCLoader");
  51   if (tpcl == 0x0) {
  52     cerr<<"Can not get TPC Loader"<<endl;
  53     return 1;
  54   }
  55
  56   if (tpcl->LoadDigits()) {
  57     cerr<<"Error occured while loading digits"<<endl;
  58     return 1;
  59   }
  60
  61   if (rl->LoadgAlice()) {
  62     cerr<<"Error occured while LoadgAlice"<<endl;
  63     return 1;
  64   }
  65
  66   gAlice=rl->GetAliRun();
  67   if (!gAlice) {
  68     cerr<<"Can't get gAlice !\n";
  69     return 1;
  70   }
  71
  72   TDirectory *cwd = gDirectory;
  73
  74   AliTPCv4 *tpc = (AliTPCv4*)gAlice->GetDetector("TPC");
  75   Int_t ver = tpc->IsVersion();
  76   cerr<<"TPC version "<<ver<<" has been found !\n";
  77
  78   rl->CdGAFile();
  79
  80   AliTPCParam *param=(AliTPCParamSR *)gDirectory->Get("75x40_100x60_150x60");
  81   if (!param) {cerr<<"TPC parameters have not been found !\n"; return 4;}
  82
  83   AliTPCDigitsArray *digarr=new AliTPCDigitsArray;
  84   digarr->Setup(param);
  85
  86   cerr<<"It has begun"<<endl;
  87   TStopwatch timer;
  88   timer.Start();
  89   cwd->cd();
  90
  91   TTree *output_tree;
  92
  93   AliTPCclustererKr *clusters = new AliTPCclustererKr();
  94   clusters->SetParam(param);
  95   clusters->SetOutput(output_tree);
  96
  97   clusters->SetMinAdc(3);//signal threshold (everything below is treated as 0)
  98   clusters->SetMinTimeBins(2);//number of neighbouring timebins
  99   clusters->SetMaxPadRangeCm(5.);//distance of the cluster center to the center of a pad (in cm)
 100   clusters->SetMaxRowRangeCm(5.);//distance of the cluster center to the center of a padrow (in cm)
 101   clusters->SetMaxTimeRange(7.);//distance of the cluster center to the max time bin on a pad (in tackts)
 102   //ie. fabs(centerT - time)<7
 103
 104   clusters->SetIsolCut(3);//set isolation cut threshold
 105   clusters->SetValueToSize(3.1);//cut reduce peak at 0
 106
 107   Int_t nevmax=rl->GetNumberOfEvents();//number of events in run
 108   for(Int_t nev=0;nev<nevmax /*&& nev<1*/ ;nev++){
 109     rl->GetEvent(nev);
 110
 111     TTree* input_tree= tpcl->TreeD();//tree with digits
 112     if (input_tree == 0x0){
 113       cerr << "Can not get TreeD for event " <<nev<<endl;
 114       continue;
 115     }
 116     digarr->ConnectTree(input_tree);
 117     clusters->SetInput(input_tree);
 118     clusters->SetDigArr(digarr);
 119     cout<<"Processing event "<<nev<<endl;
 120     clusters->FinderIO();
 121
 122   }
 123   delete clusters;
 124   timer.Stop(); timer.Print();
 125
 126   delete rl;//cleans everything
 127
 128   return 0;
 129 }