FMD/scripts/AliFMDAnaFlowKine.C

   1 #include <AliFMDHit.h>
   2 #include <TH1D.h>
   3 #include <TH2D.h>
   4 #include <TParticle.h>
   5 #include <TBrowser.h>
   6 #include <TArrayF.h>
   7 #include <TCanvas.h>
   8 #include <TStyle.h>
   9 #include <TLegend.h>
  10 #include <TProfile.h>
  11 #include <TProfile2D.h>
  12 #include <TROOT.h>
  13 #include <AliFMDAnaESD.h>
  14 #include "AliFMDAnaFlowRing.h"
  15 #if 1
  16 # include "background.C"
  17 #else
  18 TH2* bgFMD1i() { return 0; }
  19 TH2* bgFMD2i() { return 0; }
  20 TH2* bgFMD2o() { return 0; }
  21 TH2* bgFMD3o() { return 0; }
  22 TH2* bgFMD3o() { return 0; }
  23 #endif
  24 #if 1
  25 # include "reactionplane.C"
  26 #else
  27 Float_t reactionplane(int ev) { return 0; }
  28 #endif
  29 #include <iostream>
  30
  31
  32 //====================================================================
  33 class AliFMDAnaFlowKine : public AliFMDInput
  34 {
  35 public:
  36   //__________________________________________________________________
  37   /** Constructor */
  38   AliFMDAnaFlowKine(Int_t n=1, Bool_t segments=true,
  39                     Bool_t primary_only=true, Float_t cut0=.68)
  40     : fAxis(n, -3.7, 5.3),
  41       fV2("all", "Full FMD", 2, fAxis),
  42       fPrimaryOnly(primary_only),
  43       fSegmented(segments)
  44   {
  45     AddLoad(kHits);
  46     AddLoad(kKinematics);
  47     for (int i = 0; i < 5; i++) {
  48       UShort_t det = (i == 0 ? 1 : (i <= 2 ? 2 : 3));
  49       Char_t   rng = (i == 0 ? 'I' :  (i % 2 == 1 ? 'I' : 'O'));
  50       fRing[i] = new AliFMDAnaFlowRing(det,rng,n,EtaMin(i)-.1,EtaMax(i)+.1,
  51                                        0,cut0,Cut1(i));
  52     }
  53   }
  54   //__________________________________________________________________
  55   /** Constructor */
  56   AliFMDAnaFlowKine(Float_t deta, Bool_t segments=true,
  57                     Bool_t primary_only=true, Float_t cut0=.68)
  58     : fAxis(Int_t((5.3- -3.7)/deta), -3.7, 5.3),
  59       fV2("all", "Full FMD", 2, fAxis),
  60       fPrimaryOnly(primary_only),
  61       fSegmented(segments)
  62   {
  63     for (int i = 0; i < 5; i++) {
  64       UShort_t det = (i == 0 ? 1 : i <= 2 ? 2 : 3);
  65       Char_t   rng = (i % 2 == 0 ? 'I' : 'O');
  66       fRing[i] = new AliFMDAnaFlowRing(det,rng,fAxis,0,cut0,Cut1(i));
  67     }
  68   }
  69   //__________________________________________________________________
  70   Int_t FindRing(UShort_t det, Char_t ring) const
  71   {
  72     Int_t idx = -1;
  73     switch (det) {
  74     case 1: idx = 0; break;
  75     case 2: idx = 1 + (ring == 'O' || ring == 'o' ? 1 : 0); break;
  76     case 3: idx = 3 + (ring == 'O' || ring == 'o' ? 1 : 0); break;
  77     }
  78     return idx;
  79   }
  80   //__________________________________________________________________
  81   /** Initialize */
  82   Bool_t Init()
  83   {
  84     fPhis.Set(51200);
  85     fEtas.Set(51200);
  86     fWeights.Set(51200);
  87     fNObserv = 0;
  88     fNEvents = 0;
  89
  90     Float_t pi2 = 2 * TMath::Pi();
  91
  92     fPsi  = new TH2D("all_psi", "Event planes",
  93                      fAxis.N(), fAxis.Bins(), 40, 0, pi2);
  94     fPsi->SetXTitle("#eta");
  95     fPsi->SetYTitle("#Psi");
  96     fPsi->SetZTitle("Frequency");
  97     fPsi->SetMarkerStyle(20);
  98
  99     fDPsi = new TProfile("all_dpsi", "Distance to real #Psi_{R}",
 100                          fAxis.N(), fAxis.Bins(), -pi2, pi2);
 101     fDPsi->SetXTitle("#eta");
 102     fDPsi->SetYTitle("#Delta#Psi=|#Psi-#Psi_{R}|");
 103
 104     fDPhi = new TProfile2D("all_dphi","Distribution of phi signals",
 105                           fAxis.N(), -3.7, 5.3, 40, 0, pi2, 0, 20);
 106     fDPhi->SetXTitle("#eta");
 107     fDPhi->SetYTitle("#varphi");
 108     fDPhi->SetZTitle("signal");
 109     fDPhi->SetMarkerStyle(20);
 110
 111     for (size_t i = 0; i < 5; i++) fRing[i]->Init();
 112     return AliFMDInput::Init();
 113   }
 114   //__________________________________________________________________
 115   /** Called at the beginning of an event */
 116   Bool_t Begin(Int_t ev)
 117   {
 118     fNObserv = 0;
 119     fNEvents++;
 120     Float_t psir  = reactionplane(fNEvents-1);
 121     for (size_t i = 0; i < 5; i++) {
 122       fRing[i]->Begin();
 123       fRing[i]->fTrue.SetPsi(psir);
 124     }
 125     return AliFMDInput::Begin(ev);
 126   }
 127   //__________________________________________________________________
 128   /** Called at the beginning of an event */
 129   Bool_t ProcessHit(AliFMDHit* hit, TParticle* p)
 130   {
 131     if (fPrimaryOnly && !p->IsPrimary()) return kTRUE;
 132     if (hit->IsStop())   return kTRUE;
 133     UShort_t           det  = hit->Detector();
 134     Char_t             rng  = hit->Ring();
 135     UShort_t           sec  = hit->Sector();
 136     AliFMDAnaFlowRing* frn  = fRing[FindRing(det, rng)];
 137     Double_t           the  = p->Theta();
 138     Double_t           phi  = p->Phi();
 139     Double_t           eta  = -TMath::Log(TMath::Tan(the/2));
 140     // Phi
 141     if (fSegmented) {
 142       Double_t         dPhi =  2*TMath::Pi() / frn->NSeq();
 143       if (det == 3)    phi  =  TMath::Pi() - (sec + .5) * dPhi;
 144       else             phi  =  (sec + .5) * dPhi;
 145       if (phi < 0)     phi  += 2 * TMath::Pi();
 146     }
 147
 148     Double_t           mult = 1; // hit->Edep();
 149     frn->Fill(phi, eta, mult);
 150     Fill(phi, eta, mult);
 151     return kTRUE;
 152   }
 153
 154   //__________________________________________________________________
 155   /** Fill in one strip */
 156   void Fill(Float_t phi, Float_t eta, Float_t mult)
 157   {
 158     if (fNObserv >= fEtas.fN) {
 159       ULong_t n = Int_t(fEtas.fN * 1.5)+1;
 160       fEtas.Set(n);
 161       fPhis.Set(n);
 162       fWeights.Set(n);
 163     }
 164     fEtas[fNObserv]    = eta;
 165     fPhis[fNObserv]    = phi;
 166     fWeights[fNObserv] = mult;
 167     fNObserv++;
 168     fDPhi->Fill(eta, phi);
 169   }
 170   //__________________________________________________________________
 171   /** Called at the end of an event */
 172   Bool_t End()
 173   {
 174     for (size_t i = 0; i < 5; i++) fRing[i]->End();
 175     float       psir    = reactionplane(fNEvents-1)/2;
 176     // std::cout << "PsiR=" << psir << std::endl;
 177     fV2.Event(fPhis.fArray, fEtas.fArray, fWeights.fArray, fNObserv);
 178     for (UShort_t i = 0; i < fAxis.N(); i++) {
 179       AliFMDFlowBin* bin = fV2.GetBin(i);
 180       if (bin && bin->Counts() > 0) {
 181         fPsi->Fill(fAxis.BinCenter(i),  bin->Psi() - psir);
 182         fDPsi->Fill(fAxis.BinCenter(i), bin->Psi() - psir);
 183       }
 184     }
 185     return AliFMDInput::End();
 186   }
 187   //__________________________________________________________________
 188   /** Make a frame */
 189   TH1* MakeFrame(const char* xname, const char* yname,
 190                  Float_t ymin, Float_t ymax)
 191   {
 192     static int id = 0;
 193     id++;
 194     TH1* frame = new TH1D(Form("frame%d", id), "Frame",
 195                           fAxis.N(), fAxis.Bins()[0], fAxis.Bins()[fAxis.N()]);
 196     frame->SetMinimum(ymin);
 197     frame->SetMaximum(ymax);
 198     frame->SetXTitle(xname);
 199     frame->SetYTitle(yname);
 200     frame->SetLabelFont(132, "XYZ");
 201     frame->SetTitleFont(132, "XYZ");
 202     frame->SetLabelSize(.08, "XYZ");
 203     frame->SetTitleSize(.08, "XYZ");
 204     frame->SetTitleOffset(.5, "XYZ");
 205     frame->SetNdivisions(8, "XYZ");
 206     frame->Draw();
 207     return frame;
 208   }
 209   //__________________________________________________________________
 210   /** called at the end of a run */
 211   Bool_t Finish()
 212   {
 213     gStyle->SetPalette(1);
 214     gStyle->SetOptTitle(0);
 215     gStyle->SetOptStat(0);
 216     AliFMDAnaFlowRing* fring[5];
 217     for (size_t i = 0; i < 5; i++)
 218       fring[i] = static_cast<AliFMDAnaFlowRing*>(fRing[i]);
 219
 220     TCanvas* c = new TCanvas("c1", "Result", 600, 750);
 221     c->cd();
 222     c->SetGridy();
 223     c->SetTopMargin(0.02);
 224     c->SetRightMargin(0.02);
 225     c->SetFillColor(0);
 226     c->SetBorderSize(0);
 227     c->SetBorderMode(0);
 228     c->SetHighLightColor(0);
 229     c->Divide(1,3,0,0);
 230
 231     TVirtualPad* p = 0;
 232     TH1*         f = 0;
 233
 234     p = c->cd(1);
 235     p->SetGridx();
 236     p->SetGridy();
 237     p->SetFillColor(0);
 238     f = MakeFrame("#eta", "v_{2}", 0, .1);
 239     for (int i = 0; i < 5; i++) {
 240       fring[i]->fV2.Draw("th:same");
 241       fring[i]->fTrue.Draw("th:same");
 242     }
 243     fV2.Draw("th:same");
 244
 245     p = c->cd(2);
 246     p->SetGridx();
 247     p->SetGridy();
 248     p->SetFillColor(0);
 249     f = MakeFrame("#eta", "R_{k}", .7, 1);
 250     for (int i = 0; i < 5; i++) fring[i]->fV2.Draw("tr:same");
 251     fV2.Draw("tr:same");
 252
 253     p = c->cd(3);
 254     p->SetGridx();
 255     p->SetGridy();
 256     p->SetFillColor(0);
 257     f = MakeFrame("#eta", "Counts", 0, 2e5);
 258     f->GetYaxis()->SetNoExponent();
 259     for (int i = 0; i < 5; i++) fring[i]->fV2.Draw("tc:same");
 260     // fV2.Draw("tc:same");
 261
 262     c = new TCanvas("c2", "Event plane");
 263     c->cd();
 264     c->SetGridy();
 265     c->SetTopMargin(0.02);
 266     c->SetRightMargin(0.02);
 267     c->SetFillColor(0);
 268     c->SetBorderSize(0);
 269     c->SetBorderMode(0);
 270     c->SetHighLightColor(0);
 271     float pi2 = 2*TMath::Pi();
 272     TH2* h = new TH2D("h", "h", fAxis.N(), fAxis.Bins(), 80, -pi2, pi2);
 273     h->Draw();
 274     for (int i = 0; i < 5; i++) fring[i]->fDPsi.Draw("same col");
 275     // fDPhi->Draw("same");
 276
 277     for (int i = 0; i < 5; i++) {
 278       fring[i]->fV2.Print("s");
 279       fring[i]->fTrue.Print("s");
 280     }
 281     fV2.Print("s");
 282     // gROOT->GetListOfSpecials()->Add(this);
 283     return kTRUE;
 284   }
 285   //__________________________________________________________________
 286   Bool_t IsFolder() const { return kTRUE; }
 287   //__________________________________________________________________
 288   void Browse(TBrowser* b)
 289   {
 290     for (size_t i = 0; i < 5; i++)
 291       b->Add(fRing[i], fRing[i]->Name());
 292     b->Add(&fV2);
 293     if (fPsi)  b->Add(fPsi);
 294     if (fDPsi) b->Add(fDPsi);
 295     if (fDPhi) b->Add(fDPhi);
 296   }
 297   //__________________________________________________________________
 298   /** Get min pseudo rapidity of ring idx */
 299   Float_t EtaMin(Int_t idx) const {
 300     switch (idx) {
 301     case 0: return  3.61728; break; // FMD1I
 302     case 1: return  2.18122; break; // FMD2I
 303     case 2: return  1.79821; break; // FMD2O
 304     case 3: return -3.39913; break; // FMD3I
 305     case 4: return -2.28925; break; // FMD3O
 306     }
 307     return 0;
 308   }
 309   //__________________________________________________________________
 310   /** Get max pseudo rapidity of ring idx */
 311   Float_t EtaMax(Int_t idx) const {
 312     switch (idx) {
 313     case 0: return  5.02643; break; // FMD1I
 314     case 1: return  3.57899; break; // FMD2I
 315     case 2: return  2.39085; break; // FMD2O
 316     case 3: return -2.00644; break; // FMD3I
 317     case 4: return -1.70080; break; // FMD3O
 318     }
 319     return 0;
 320   }
 321   //__________________________________________________________________
 322   /** Get cut 1 of ring idx */
 323   Float_t Cut1(Int_t idx) const {
 324     switch (idx) {
 325     case 0: return 1.33857; break; // FMD1I
 326     case 1: return 1.30464; break; // FMD2I
 327     case 2: return 1.27263; break; // FMD2O
 328     case 3: return 1.28052; break; // FMD3I
 329     case 4: return 1.23041; break; // FMD3O
 330     }
 331     return 0;
 332   }
 333   /** Number of events */
 334   Int_t fNEvents;
 335   /** Flow rings */
 336   AliFMDAnaFlowRing* fRing[5];
 337   /** Flow axis */
 338   AliFMDFlowAxis fAxis;
 339   /** Flow object */
 340   AliFMDFlowBinned1D fV2;
 341   /** Cache of phi */
 342   TArrayD            fPhis;
 343   /** Cache of eta */
 344   TArrayD            fEtas;
 345   /** Cache of weights */
 346   TArrayD            fWeights;
 347   /** Number of observables */
 348   Int_t              fNObserv;
 349   /** Histogram of Psi */
 350   TH2D*              fPsi;       // Event planes
 351   /** Histogram of dPsi */
 352   TProfile*          fDPsi;   // Mean distance to Psi
 353   /** Histogram of dPhi */
 354   TProfile2D*        fDPhi;
 355   /** Wether to include primaries only */
 356   Bool_t             fPrimaryOnly;
 357   /** Whether to segment the phi space according to detector */
 358   Bool_t             fSegmented;
 359   ClassDef(AliFMDAnaFlowKine,1)
 360 };
 361
 362 //____________________________________________________________________
 363 //
 364 // EOF
 365 //