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 ////////////////////////////////////////////////////////////////////////////////