add h-jet jet mass analysis class
[u/mrichter/AliRoot.git] / TPC / FindKrClusters.C
CommitLineData
a65a7e70 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
34Int_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}