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   //
  33   //Load DataBase
  34   //
  35   //char *ocdbpath ="local:///afs/cern.ch/alice/tpctest/OCDB";
  36   //char *ocdbpath ="local:///home/matyja/baza/OCDB";
  37   char *ocdbpath ="local:///data/baza/OCDB";
  38   if (ocdbpath==0){
  39     ocdbpath="alien://folder=/alice/data/2007/LHC07w/OCDB/";
  40   }
  41   printf("OCDB PATH = %s\n",ocdbpath);
  42   AliCDBManager * man = AliCDBManager::Instance();
  43   man->SetDefaultStorage(ocdbpath);
  44   man->SetRun(0);
  45
  46   AliRunLoader* rl = AliRunLoader::Open("galice.root");
  47   if (rl == 0x0) {
  48     cerr<<"Can not open session"<<endl;
  49     return 1;
  50   }
  51
  52   AliTPCLoader *tpcl = (AliTPCLoader*)rl->GetLoader("TPCLoader");
  53   if (tpcl == 0x0) {
  54     cerr<<"Can not get TPC Loader"<<endl;
  55     return 1;
  56   }
  57
  58   if (tpcl->LoadDigits()) {
  59     cerr<<"Error occured while loading digits"<<endl;
  60     return 1;
  61   }
  62
  63   if (rl->LoadgAlice()) {
  64     cerr<<"Error occured while LoadgAlice"<<endl;
  65     return 1;
  66   }
  67
  68   gAlice=rl->GetAliRun();
  69   if (!gAlice) {
  70     cerr<<"Can't get gAlice !\n";
  71     return 1;
  72   }
  73
  74   TDirectory *cwd = gDirectory;
  75
  76   AliTPCv4 *tpc = (AliTPCv4*)gAlice->GetDetector("TPC");
  77   Int_t ver = tpc->IsVersion();
  78   cerr<<"TPC version "<<ver<<" has been found !\n";
  79
  80   rl->CdGAFile();
  81
  82   AliTPCParam *param=(AliTPCParamSR *)gDirectory->Get("75x40_100x60_150x60");
  83   if (!param) {cerr<<"TPC parameters have not been found !\n"; return 4;}
  84
  85   AliTPCDigitsArray *digarr=new AliTPCDigitsArray;
  86   digarr->Setup(param);
  87
  88   cerr<<"It has begun"<<endl;
  89   TStopwatch timer;
  90   timer.Start();
  91   cwd->cd();
  92
  93   TTree *output_tree;
  94
  95   AliTPCclustererKr *clusters = new AliTPCclustererKr();
  96   clusters->SetParam(param);
  97   clusters->SetOutput(output_tree);
  98
  99   clusters->SetMinAdc(3);//signal threshold (everything below is treated as 0)
 100   clusters->SetMinTimeBins(2);//number of neighbouring timebins
 101   clusters->SetMaxPadRangeCm(5.);//distance of the cluster center to the center of a pad (in cm)
 102   clusters->SetMaxRowRangeCm(5.);//distance of the cluster center to the center of a padrow (in cm)
 103   clusters->SetMaxTimeRange(7.);//distance of the cluster center to the max time bin on a pad (in tackts)
 104   //ie. fabs(centerT - time)<7
 105
 106   clusters->SetIsolCut(3);//set isolation cut threshold
 107   clusters->SetValueToSize(3.1);//cut reduce peak at 0
 108
 109   Int_t nevmax=rl->GetNumberOfEvents();//number of events in run
 110   for(Int_t nev=0;nev<nevmax /*&& nev<1*/ ;nev++){
 111     rl->GetEvent(nev);
 112
 113     TTree* input_tree= tpcl->TreeD();//tree with digits
 114     if (input_tree == 0x0){
 115       cerr << "Can not get TreeD for event " <<nev<<endl;
 116       continue;
 117     }
 118     digarr->ConnectTree(input_tree);
 119     clusters->SetInput(input_tree);
 120     clusters->SetDigArr(digarr);
 121     cout<<"Processing event "<<nev<<endl;
 122     clusters->FinderIO();
 123
 124   }
 125   delete clusters;
 126   timer.Stop(); timer.Print();
 127
 128   delete rl;//cleans everything
 129
 130   return 0;
 131 }