TPC/FindKrClusters.C

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