Update rawdata format for trigger (Christian)
[u/mrichter/AliRoot.git] / FLOW / AliFlowReconstruction.C
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
32 void 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
110 Double_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
146 Double_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 ////////////////////////////////////////////////////////////////////////////////