TPC/FindKrClusters.C

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