- Volume name attribute replaced with volume path
[u/mrichter/AliRoot.git] / FLOW / AliFlowReconstruction.C
CommitLineData
0321b05e 1//
2// This macro is a part of "Alice PPR fast flow analysis package"
3//
4// The macro does:
5// 1) open existing database with Flow Pico DST ("flowPicoEvent.root")
6// 2) open source data - galice.root type file with events
7// 3) reconstruct psi, psiA, psiB (sub-events in eta space) and v2.
8// 4) fill Flow Pico DST with reconstructed data
9//
10// INPUT PARAMETERS:
11// inputName: file with events
12// if (inputName == 0) macro tries to fetch system variable FILE_NAME
13// and use its value as the name of the file. (useful for bath processing)
14//
15// maxN: maximum number of events to be processed, if 0 process all events
16//
17// NOTE:
18// file "flowPicoEvent.root" with "flowData" NTuple have to exist
19// it can be created by "AliFlowCreate.C"
20//
21// Sylwester Radomski, GSI
22// mail: S.Radomski@gsi
23// 23. Oct. 2002
24//
25////////////////////////////////////////////////////////////////////////////////
26
27//#include "TTree.h"
28//#include "TNtuple.h"
29//#include "TFile.h"
30//#include "AliRun.h"
31
32void AliFlowReconstruction(const char* inputName = 0, Int_t maxN = 0) {
33
34 const char *dataName = "flowPicoEvent.root";
35
36 TTree *kine;
37 TNtuple *data;
38
39 Double_t trueV2, truePsi;
40 Double_t v2, psi, psiA, psiB;
41 Int_t N;
42 Int_t mult;
43
44 if (gAlice) delete gAlice;
45 gAlice = 0;
46
47 TFile *dataFile = new TFile(dataName,"UPDATE");
48 data = (TNtuple *) dataFile->Get("flowData");
49
50 TFile *inputFile;
51
52 if (!inputName) {
53 TFile *inputFile = new TFile(gSystem->Getenv("FILE_NAME"), "READ");
54 ::Info("AliFlowReconstruction", "Using: %s",gSystem->Getenv("FILE_NAME"));
55 } else {
56 TFile *inputFile = new TFile(inputName, "READ");
57 }
58
59 gAlice = (AliRun*) inputFile->Get("gAlice");
60
61 N = gAlice->GetEventsPerRun();
62 if (maxN != 0 && maxN < N) N = maxN;
63
64 ::Info("AliFlowReconstruction", "Number of events to be processed = %d", N);
65
66 for (Int_t i=0; i<N; i++) {
67
68 ::Info("Processing event i = %d", i);
69
70 gAlice->GetEvent(i);
71
72 AliHeader *header = gAlice->GetHeader();
73 AliGenGeVSimEventHeader *egHeader = (AliGenGeVSimEventHeader*)header->GenEventHeader();
74
75 truePsi = egHeader->GetEventPlane() * 180 / TMath::Pi();
76 while(truePsi > 90) truePsi -= 180;
77
78 trueV2 = egHeader->GetEllipticFlow(211);
79
80 mult = gAlice->GetNtrack();
81 kine = gAlice->TreeK();
82
83 psi = FlowPsi(kine, 0);
84 psiA = FlowPsi(kine, 1);
85 psiB = FlowPsi(kine, -1);
86
87 v2 = FlowV2(kine, psi);
88 data->Fill(1, i, mult, truePsi, trueV2, psi, psiA, psiB, v2);
89 }
90
91 dataFile->cd();
92 data->Write();
93
94 delete data;
95
96 delete gAlice;
97 gAlice = 0;
98
99 dataFile->Close();
100 inputFile->Close();
101
102 delete dataFile;
103 delete inputFile;
104
105 gSystem->Exit(0);
106}
107
108////////////////////////////////////////////////////////////////////////////////
109
110Double_t FlowPsi(TTree *kine, Int_t etaPart) {
111 //
112 // this function calculates psi from list of particles
113 // etaPart - division to sub-events in eta
114 // 0 all particles,
115 // -1 negative eta
116 // 1 positive eta
117 //
118
119 Int_t mult = kine->GetEntries();
120 Double_t *phi, *eta;
121 Double_t psi;
122
123 Double_t Ssin = 0, Scos = 0;
124
125 kine->Draw("Particles->Phi():Particles->Eta()", "", "goff");
126 phi = kine->GetV1();
127 eta = kine->GetV2();
128
129 for (Int_t i=0; i<mult; i++) {
130
131 if ( (etaPart * eta[i]) < 0 ) continue;
132
133 Ssin += sin( 2 * phi[i] );
134 Scos += cos( 2 * phi[i] );
135 }
136
137 //psi = atan( Ssin / Scos ) / 2;
138 psi = atan2 (Ssin, Scos) / 2;
139 psi = psi * 180. / TMath::Pi();
140
141 return psi;
142}
143
144////////////////////////////////////////////////////////////////////////////////
145
146Double_t FlowV2(TTree *kine, Double_t psi) {
147
148 Int_t mult = kine->GetEntries();
149 Double_t *phi;
150 Double_t V2;
151
152 Double_t Ssin = 0, Scos = 0;
153
154 kine->Draw("Particles->Phi()", "", "goff");
155 phi = kine->GetV1();
156
157 for (Int_t i=0; i<mult; i++) {
158
159 Ssin += sin( 2 * (phi[i] - psi) );
160 Scos += cos( 2 * (phi[i] - psi) );
161 }
162
163 V2 = sqrt(Ssin*Ssin + Scos*Scos) / mult;
164
165 return V2;
166}
167
168////////////////////////////////////////////////////////////////////////////////