JETAN/jetan2004/anaAliJets.C

   1 // $Id$
   2
   3 #if !defined(__CINT__) || defined(__MAKECINT__)
   4 #include <stdio.h>
   5 #include <iostream.h>
   6 #include <time.h>
   7 #include <TROOT.h>
   8 #include <TCanvas.h>
   9 #include <TH1F.h>
  10 #include <TH2F.h>
  11 #include <TProfile.h>
  12 #include <TProfile2D.h>
  13 #include <TGraph.h>
  14 #include <TNtuple.h>
  15 #include <TRandom.h>
  16 #include <TSystem.h>
  17 #include <TStopwatch.h>
  18 #include <TFile.h>
  19 #include <TChain.h>
  20 #include <TSystem.h>
  21 #include <TStyle.h>
  22 #include <TTimeStamp.h>
  23 #include "AliTkConeJetEvent.h"
  24 #include "AliTkConeJet.h"
  25 #include "AliTkConeJetFinderV2.h"
  26 #include <AliJetParticle.h>
  27 #include <AliJetParticlesReader.h>
  28 #include <AliJetParticlesReaderKine.h>
  29 #include <AliJetParticlesReaderESD.h>
  30 #include <AliJetParticlesReaderHLT.h>
  31 #include <AliJetEventParticles.h>
  32 #endif
  33
  34 Float_t relphi(Float_t phi1, Float_t phi2);
  35 Float_t addphi(Float_t phi1, Float_t phi2);
  36 Float_t diffphi(Float_t phi1, Float_t phi2);
  37 Int_t eventindex(Float_t phi1, Float_t phi2);
  38 void convert(Float_t pjet[4], Float_t &pt, Float_t &theta, Float_t &eta, Float_t &phi);
  39
  40 void anaAliJets(Char_t *filename,Char_t *savefilename,
  41                 Int_t mEnergy=100,Int_t mBackEn=-1,Int_t nMaxEvents=-1,
  42                 Char_t *evfilename=0,Char_t *sigevfilename=0,
  43                 Char_t *signalfilename=0,Char_t *monteconefilename=0)
  44   /*
  45     filename       = cone finder event
  46     evfilename     = background event or jetevent only (if 0 take from jetevent)
  47     sigevfilename  = signal event if background event is used (otherwice == 0)
  48     signalfilename = original signal event (eg. for real tracking to get trigger jets)
  49     monteconefilename = reconstructed jets from signal event (esd=0,esd=10)
  50    */
  51
  52 {
  53   const Float_t minEtRecoCut=mEnergy*0.95;
  54   const Float_t maxEtRecoCut=mEnergy*1.05;
  55   Float_t minJetEt;
  56   if(mBackEn) minJetEt=mBackEn;
  57   else minJetEt=mEnergy/4.;
  58   const Float_t minPartPt=0.5;
  59   //const Char_t *figprefix=0;
  60   //const Char_t *figdirname=".";
  61
  62   const Int_t nclasses=22;
  63   const Float_t cletmin[nclasses] = {0,10,20,30,40,50,60,70,80,90,100,
  64                                      110,120,130,140,150,160,170,180,190,200,0};
  65   const Float_t cletmax[nclasses] = {10,20,30,40,50,60,70,80,90,100,
  66                                      110,120,130,140,150,160,170,180,190,200,350,350};
  67
  68   //differential shape
  69   const Float_t deltaR=0.1/2;
  70
  71   Float_t corrfac=0.;
  72 #ifdef APPLYCORRECTION
  73   if(!allparts) corrfac=2./3;
  74   Float_t conefluc=1.;
  75   if(Radius<=0.3) conefluc=0.8;
  76   else if(Radius<=0.5) conefluc=0.9;
  77   else if(Radius<=0.7) conefluc=0.9;
  78   corrfac*=conefluc;
  79 #endif
  80
  81   Char_t dummy[1024];
  82   Char_t name[1024];
  83   Char_t title[1024];
  84
  85   TH1F *hJetEt = new TH1F("hJetEt","E_{T} (jet)",350,0,350);
  86   hJetEt->GetXaxis()->SetTitle("E_{T} (jet) [GeV]");
  87   hJetEt->GetYaxis()->SetTitle("dN/dE_{T} [GeV^{-1}]");
  88
  89   TH1F *hBackJetEt = new TH1F("hBackJetEt","E_{T} (jet)",350,0,350);
  90   hBackJetEt->GetXaxis()->SetTitle("E_{T} (jet) [GeV]");
  91   hBackJetEt->GetYaxis()->SetTitle("dN/dE_{T} [GeV^{-1}]");
  92
  93   TH1F *hJetEtall = new TH1F("hAllJetEt","E_{T} (jet)",350,0,350);
  94   hJetEtall->GetXaxis()->SetTitle("E_{T} (jet) [GeV]");
  95   hJetEtall->GetYaxis()->SetTitle("dN/dE_{T} [GeV^{-1}]");
  96
  97   TProfile *hJetEttrue = new TProfile("hJetEttrue","E_{T} (jet)",350,0,350,0,1.5);
  98   hJetEttrue->GetXaxis()->SetTitle("E_{T} (jet) [GeV]");
  99   hJetEttrue->GetYaxis()->SetTitle("E^{true}_{T}/E_{T}");
 100
 101   TProfile *hBackJetEttrue = new TProfile("hBackJetEttrue","E_{T} (jet)",350,0,350,0,1.5);
 102   hBackJetEttrue->GetXaxis()->SetTitle("E_{T} (jet) [GeV]");
 103   hBackJetEttrue->GetYaxis()->SetTitle("E^{true}_{T}/E_{T}");
 104
 105   TProfile *hJetEtalltrue = new TProfile("hAllJetEttrue","E_{T} (jet)",350,0,350,0,1.5);
 106   hJetEtalltrue->GetXaxis()->SetTitle("E_{T} (jet) [GeV]");
 107   hJetEtalltrue->GetYaxis()->SetTitle("E^{true}_{T}/E_{T}");
 108
 109   TH1F *hJetEtUQTrigger = new TH1F("hJetEtUQTrigger","E_{T} (trigger jet)",350,0,350);
 110   hJetEtUQTrigger->GetXaxis()->SetTitle("E_{T} (jet) [GeV]");
 111   hJetEtUQTrigger->GetYaxis()->SetTitle("dN/dE_{T} [GeV^{-1}]");
 112
 113   TH1F *hJetEtTrigger = new TH1F("hJetEtTrigger","E_{T} (trigger jet)",350,0,350);
 114   hJetEtTrigger->GetXaxis()->SetTitle("E_{T} (jet) [GeV]");
 115   hJetEtTrigger->GetYaxis()->SetTitle("dN/dE_{T} [GeV^{-1}]");
 116
 117   TH2F *hJetEtvsEll = new TH2F("hJetEtvsEll","E_{T} (jet) versus Length",350,0,350,100,0,10);
 118   hJetEtvsEll->GetXaxis()->SetTitle("E_{T} (jet) [GeV]");
 119   hJetEtvsEll->GetYaxis()->SetTitle("L [fm]");
 120
 121   TH2F *hJetEtallvsEll = new TH2F("hAllJetEtallvsEll","E_{T} (jet) versus Length",350,0,350,100,0,10);
 122   hJetEtallvsEll->GetXaxis()->SetTitle("E_{T} (jet) [GeV]");
 123   hJetEtallvsEll->GetYaxis()->SetTitle("L [fm]");
 124
 125   TH2F *hJetEtvsTrigger = new TH2F("hJetEtvsTrigger","",350,0,350,350,0,350);
 126   hJetEtvsTrigger->GetXaxis()->SetTitle("E_{T} (jet) [GeV]");
 127   hJetEtvsTrigger->GetYaxis()->SetTitle("E_{T} (jettrigger) [GeV]");
 128
 129   TH2F *hJetEtvsUQTrigger = new TH2F("hJetEtvsUQTrigger","",350,0,350,350,0,350);
 130   hJetEtvsUQTrigger->GetXaxis()->SetTitle("E_{T} (jet) [GeV]");
 131   hJetEtvsUQTrigger->GetYaxis()->SetTitle("E_{T} (jettrigger) [GeV]");
 132
 133   TH2F *hJetEtvsEt = new TH2F("hJetEtvsEt","",350,0,350,350,0,350);
 134   hJetEtvsEt->GetXaxis()->SetTitle("E_{T} (jet) [GeV]");
 135   hJetEtvsEt->GetYaxis()->SetTitle("E_{T} (jet) [GeV]");
 136
 137   TH1F *hJetZ     =  new TH1F("hjetZ","Z distribution",100,0,1);
 138   hJetZ->GetXaxis()->SetTitle("Z");
 139   hJetZ->GetYaxis()->SetTitle("dN/dZ");
 140
 141   TH1F *hJet1    =  new TH1F("hjet1","E_{t} distribution",350,0,350);
 142   hJet1->GetXaxis()->SetTitle("E_{T} (parton) [GeV]");
 143   hJet1->GetYaxis()->SetTitle("dN/dE_{T} [GeV^{-1}]");
 144
 145   TH1F *hJet2    =  new TH1F("hjet2","E_{t} distribution",350,0,350);
 146   hJet2->GetXaxis()->SetTitle("E_{T} (parton) [GeV]");
 147   hJet2->GetYaxis()->SetTitle("dN/dE_{T} [GeV^{-1}]");
 148
 149   TH1F **hJettype=  new TH1F*[3];
 150   for(Int_t i=0;i<3;i++){
 151     Char_t t[100];
 152     sprintf(t,"hJettype%d",i);
 153     Char_t tit[100];
 154     if(i==0)
 155       sprintf(tit,"Unknown");
 156     else if(i==1)
 157       sprintf(tit,"Quark");
 158     else
 159       sprintf(tit,"Gluon");
 160     hJettype[i]=new TH1F(t,tit,350,0,350);
 161     hJettype[i]->GetXaxis()->SetTitle("E_{T} (jet) [GeV]");
 162     hJettype[i]->GetYaxis()->SetTitle("dN/dE_{T} [GeV^{-1}]");
 163   }
 164
 165   TH2F *hAxesDiff = new TH2F("hAxesDiff","",120,0,TMath::Pi(),40,0,2);
 166   hAxesDiff->GetXaxis()->SetTitle("#Delta #phi");
 167   hAxesDiff->GetYaxis()->SetTitle("#Delta #eta");
 168
 169   //---
 170   TH1F *hJetEtres = new TH1F("hJetEtres","E_{T} (jet)",350,0,350);
 171   hJetEtres->GetXaxis()->SetTitle("E_{T} (jet) [GeV]");
 172   hJetEtres->GetYaxis()->SetTitle("dN/dE_{T} [GeV^{-1}]");
 173
 174   TH1F *hJetEtratio = new TH1F("hJetEtratio","E_{T} (jet)",200,0,2);
 175   hJetEtratio->GetXaxis()->SetTitle("E_{T} (jet)/E_{T} (trigger)");
 176   hJetEtratio->GetYaxis()->SetTitle("Number of jets");
 177
 178   TProfile *hJetEtrestrue = new TProfile("hJetEttestrue","E_{T} (jet)",35,0,350,-2,2);
 179   hJetEtrestrue->GetXaxis()->SetTitle("E_{T} (jet) [GeV]");
 180   hJetEtrestrue->GetYaxis()->SetTitle("E^{true}_{T}/E_{T}");
 181
 182   TH1F *hJetEtresTrigger = new TH1F("hJetEtresTrigger","E_{T} (trigger jet)",350,0,350);
 183   hJetEtresTrigger->GetXaxis()->SetTitle("E_{T} (jet) [GeV]");
 184   hJetEtresTrigger->GetYaxis()->SetTitle("dN/dE_{T} [GeV^{-1}]");
 185
 186   TProfile2D *hAxesDiffres = new TProfile2D("hAxesDiffres","",240,0,TMath::Pi(),100,0,2);
 187   hAxesDiffres->GetXaxis()->SetTitle("#Delta #phi");
 188   hAxesDiffres->GetYaxis()->SetTitle("#Delta #eta");
 189
 190   TH1F *hPhires = new TH1F("hPhires","",250,-1,1);
 191   hPhires->GetXaxis()->SetTitle("#Delta #phi [rad]");
 192   hPhires->GetYaxis()->SetTitle("Number of jets");
 193
 194   TH1F *hEtares = new TH1F("hEtares","",250,-1,1);
 195   hEtares->GetXaxis()->SetTitle("#Delta #eta [rad]");
 196   hEtares->GetYaxis()->SetTitle("Number of jets");
 197
 198   TH1F *hmJetEtres = new TH1F("hmJetEtres","E_{T} (jet)",350,0,350);
 199   hmJetEtres->GetXaxis()->SetTitle("E_{T} (jet) [GeV]");
 200   hmJetEtres->GetYaxis()->SetTitle("dN/dE_{T} [GeV^{-1}]");
 201
 202   TH1F *hmJetEtratio = new TH1F("hmJetEtratio","E_{T} (jet)",200,0,2);
 203   hmJetEtratio->GetXaxis()->SetTitle("E_{T} (jet)/E_{T} (trigger)");
 204   hmJetEtratio->GetYaxis()->SetTitle("Number of jets");
 205
 206   TProfile *hmJetEtrestrue = new TProfile("hmJetEttestrue","E_{T} (jet)",35,0,350,-2,2);
 207   hmJetEtrestrue->GetXaxis()->SetTitle("E_{T} (jet) [GeV]");
 208   hmJetEtrestrue->GetYaxis()->SetTitle("E^{true}_{T}/E_{T}");
 209
 210   TProfile2D *hmAxesDiffres = new TProfile2D("hmAxesDiffres","",240,0,TMath::Pi(),100,0,2);
 211   hmAxesDiffres->GetXaxis()->SetTitle("#Delta #phi");
 212   hmAxesDiffres->GetYaxis()->SetTitle("#Delta #eta");
 213
 214   TH1F *hmPhires = new TH1F("hmPhires","",250,-1,1);
 215   hmPhires->GetXaxis()->SetTitle("#Delta #phi [rad]");
 216   hmPhires->GetYaxis()->SetTitle("Number of jets");
 217
 218   TH1F *hmEtares = new TH1F("hmEtares","",250,-1,1);
 219   hmEtares->GetXaxis()->SetTitle("#Delta #eta [rad]");
 220   hmEtares->GetYaxis()->SetTitle("Number of jets");
 221
 222   TH1F *hPhiMonteres = new TH1F("hPhiMonteres","",250,-1,1);
 223   hPhiMonteres->GetXaxis()->SetTitle("#Delta #phi [rad]");
 224   hPhiMonteres->GetYaxis()->SetTitle("Number of jets");
 225
 226   TH1F *hEtaMonteres = new TH1F("hEtaMonteres","",250,-1,1);
 227   hEtaMonteres->GetXaxis()->SetTitle("#Delta #eta [rad]");
 228   hEtaMonteres->GetYaxis()->SetTitle("Number of jets");
 229
 230   TH1F *hEtMonteres = new TH1F("hEtMonteres","",250,-125,125);
 231   hEtMonteres->GetXaxis()->SetTitle("#Delta E_{T} [GeV]");
 232   hEtMonteres->GetYaxis()->SetTitle("Number of jets");
 233
 234   TH1F *hEtMonteratio = new TH1F("hEtMonteratio","E_{T} (jet)",200,0,2);
 235   hEtMonteratio->GetXaxis()->SetTitle("E^{rec}_{T} / E^{monte}_{T}");
 236   hEtMonteratio->GetYaxis()->SetTitle("Number of jets");
 237
 238   //---
 239
 240   TH1F **hJetEttrigger=new TH1F*[9];
 241   TH1F **hJetEttrigger2=new TH1F*[9];
 242   for(Int_t i=0;i<9;i++){
 243     sprintf(dummy,"hJetEttrigger%d",i);
 244     hJetEttrigger[i] = new TH1F(dummy,"E_{T} (jet)",350,0,350);
 245     hJetEttrigger[i]->GetXaxis()->SetTitle("E_{T} (jet) [GeV]");
 246     hJetEttrigger[i]->GetYaxis()->SetTitle("dN/dE_{T} [GeV^{-1}]");
 247     sprintf(dummy,"hJetEttrigger2%d",i);
 248     hJetEttrigger2[i] = new TH1F(dummy,"E_{T} (jet)",350,0,350);
 249     hJetEttrigger2[i]->GetXaxis()->SetTitle("E_{T} (jet) [GeV]");
 250     hJetEttrigger2[i]->GetYaxis()->SetTitle("dN/dE_{T} [GeV^{-1}]");
 251   }
 252   sprintf(dummy,"hJetEttriggernorm");
 253   TH1F *hJetEttriggernorm = new TH1F(dummy,"E_{T} (jet)",350,0,350);
 254   hJetEttriggernorm->GetXaxis()->SetTitle("E_{T} (jet) [GeV]");
 255   hJetEttriggernorm->GetYaxis()->SetTitle("dN/dE_{T} [GeV^{-1}]");
 256
 257   //---
 258
 259   TH1F **hJetLeadingPt = new TH1F*[nclasses];
 260   for(Int_t i=0;i<nclasses;i++){
 261     sprintf(dummy,"hJetLeadingPt%d",i);
 262     sprintf(title,"Transverse Momentum Fraction of Leading Particle (%.0f GeV<E_{T}<%.0f GeV)",cletmin[i],cletmax[i]);
 263     hJetLeadingPt[i]=new TH1F(dummy,title,100,0,1);
 264     hJetLeadingPt[i]->GetXaxis()->SetTitle("z=p_{T}^{lead.part.}/P_{T}^{jet}");
 265     hJetLeadingPt[i]->GetYaxis()->SetTitle("dN/dz");
 266   }
 267
 268   TH1F **hJetFragLeadingPt = new TH1F*[nclasses];
 269   for(Int_t i=0;i<nclasses;i++){
 270     sprintf(dummy,"hJetFragLeadingPt%d",i);
 271     sprintf(title,"Fragmentation Function (%.0f GeV<E_{T}<%.0f GeV)",cletmin[i],cletmax[i]);
 272     hJetFragLeadingPt[i]=new TH1F(dummy,title,100,0,1);
 273     hJetFragLeadingPt[i]->GetXaxis()->SetTitle("z=p_{T}/p_{T}^{lead.part.}");
 274     hJetFragLeadingPt[i]->GetYaxis()->SetTitle("dN/dz");
 275   }
 276
 277   TH1F **hJetLeadingPtDist = new TH1F*[nclasses];
 278   for(Int_t i=0;i<nclasses;i++){
 279     sprintf(dummy,"hJetLeadingPtDist%d",i);
 280     sprintf(title,"Transverse Momentum Fraction of Leading Particle (%.0f GeV<E_{T}<%.0f GeV)",cletmin[i],cletmax[i]);
 281     hJetLeadingPtDist[i]=new TH1F(dummy,title,600,0,150);
 282     hJetLeadingPtDist[i]->GetXaxis()->SetTitle("p_{T}^{lead.part.} [GeV]");
 283     hJetLeadingPtDist[i]->GetYaxis()->SetTitle("dN/dp_{T}^{lead.part.} [GeV^{-1}]");
 284   }
 285
 286   TH1F **hJetFragL = new TH1F*[nclasses];
 287   for(Int_t i=0;i<nclasses;i++){
 288     sprintf(dummy,"hJetFragL%d",i);
 289     sprintf(title,"Longitudinal Fragmentation Function (%.0f GeV<E_{T}<%.0f GeV)",cletmin[i],cletmax[i]);
 290     hJetFragL[i]=new TH1F(dummy,title,100,0,1);
 291     hJetFragL[i]->GetXaxis()->SetTitle("z=p_{L}/P^{jet}");
 292     hJetFragL[i]->GetYaxis()->SetTitle("dN/dz");
 293   }
 294
 295   TH1F **hJetFragPL = new TH1F*[nclasses];
 296   for(Int_t i=0;i<nclasses;i++){
 297     sprintf(dummy,"hJetFragPL%d",i);
 298     sprintf(title,"Longitudinal Momentum (%.0f GeV<E_{T}<%.0f GeV)",cletmin[i],cletmax[i]);
 299     hJetFragPL[i]=new TH1F(dummy,title,600,0,150);
 300     hJetFragPL[i]->GetXaxis()->SetTitle("p_{L}");
 301     hJetFragPL[i]->GetYaxis()->SetTitle("dN/dp_{L} [GeV^{-1}]");
 302   }
 303
 304   TH1F **hJetFragT = new TH1F*[nclasses];
 305   for(Int_t i=0;i<nclasses;i++){
 306     sprintf(dummy,"hJetFragT%d",i);
 307     sprintf(title,"Transversal Fragmentation Function (%.0f GeV<E_{T}<%.0f GeV)",cletmin[i],cletmax[i]);
 308     hJetFragT[i]=new TH1F(dummy,title,250,0,25);
 309     hJetFragT[i]->GetXaxis()->SetTitle("j_{T} [GeV]");
 310     hJetFragT[i]->GetYaxis()->SetTitle("dN/dj_{T}");
 311   }
 312
 313   TH1F **hJetFragPt = new TH1F*[nclasses];
 314   for(Int_t i=0;i<nclasses;i++){
 315     sprintf(dummy,"hJetFragPt%d",i);
 316     sprintf(title,"Particle Transversal Momentum Distribution (%.0f GeV<E_{T}<%.0f GeV)",cletmin[i],cletmax[i]);
 317     hJetFragPt[i]=new TH1F(dummy,title,600,0,150);
 318     hJetFragPt[i]->GetXaxis()->SetTitle("p_{T} [GeV]");
 319     hJetFragPt[i]->GetYaxis()->SetTitle("dN/dp_{T}");
 320   }
 321
 322   TH1F **hJetN = new TH1F*[nclasses];
 323   for(Int_t i=0;i<nclasses;i++){
 324     sprintf(dummy,"hJetN%d",i);
 325     sprintf(title,"Particle Multiplicity (%.0f GeV<E_{T}<%.0f GeV)",cletmin[i],cletmax[i]);
 326     hJetN[i]=new TH1F(dummy,title,150+1,-0.5,150+0.5);
 327     hJetN[i]->GetXaxis()->SetTitle("n = Number of Particles within Jet");
 328     hJetN[i]->GetYaxis()->SetTitle("dN/dn");
 329   }
 330
 331   TH1F **hJetMeanPt = new TH1F*[nclasses];
 332   for(Int_t i=0;i<nclasses;i++){
 333     sprintf(dummy,"hJetMeanPt%d",i);
 334     sprintf(title,"Mean Transverse Jet Energy per Particle (%.0f GeV<E_{T}<%.0f GeV)",cletmin[i],cletmax[i]);
 335     hJetMeanPt[i]=new TH1F(dummy,title,125,0,25);
 336     hJetMeanPt[i]->GetXaxis()->SetTitle("#bar{p} = P_{T}^{jet}/n [GeV]");
 337     hJetMeanPt[i]->GetYaxis()->SetTitle("dN/d#bar{p}");
 338   }
 339
 340   //intshape
 341   Float_t retall[nclasses][10];
 342   Float_t retlow[nclasses][10];
 343   Float_t retlow1[nclasses][10];
 344   Float_t retlow2[nclasses][10];
 345   Float_t retlow3[nclasses][10];
 346   Float_t retlow4[nclasses][10];
 347   Float_t rxet[10];
 348   for(Int_t i=0;i<10;i++) {
 349     for(Int_t j=0;j<nclasses;j++) {
 350       retall[j][i]=0.;
 351       retlow[j][i]=0.;
 352       retlow1[j][i]=0.;
 353       retlow2[j][i]=0.;
 354       retlow3[j][i]=0.;
 355       retlow4[j][i]=0.;
 356     }
 357     rxet[i]=(i+1.)/10;
 358   }
 359
 360   //diffshape
 361   Float_t dretall[nclasses][11];
 362   Float_t dretlow[nclasses][11];
 363   Float_t dretlow1[nclasses][11];
 364   Float_t dretlow2[nclasses][11];
 365   Float_t dretlow3[nclasses][11];
 366   Float_t dretlow4[nclasses][11];
 367   Float_t drxet[11];
 368   for(Int_t i=0;i<11;i++) {
 369     for(Int_t j=0;j<nclasses;j++) {
 370       dretall[j][i]=0.;
 371       dretlow[j][i]=0.;
 372       dretlow1[j][i]=0.;
 373       dretlow2[j][i]=0.;
 374       dretlow3[j][i]=0.;
 375       dretlow4[j][i]=0.;
 376     }
 377     drxet[i]=1.*i/10.;
 378   }
 379
 380   TH1F **hPhiCorr = new TH1F*[nclasses];
 381   for(Int_t i=0;i<nclasses;i++){
 382     sprintf(dummy,"hPartPhiCorr%d",i);
 383     sprintf(title,"Azimuthal Correlation (%.0f GeV<E_{T}<%.0f GeV)",cletmin[i],cletmax[i]);
 384     hPhiCorr[i]= new TH1F(dummy,title,120,0,TMath::TwoPi());
 385     hPhiCorr[i]->GetXaxis()->SetTitle("#phi");
 386     hPhiCorr[i]->GetYaxis()->SetTitle("dN/d#phi");
 387   }
 388
 389   //
 390   //global event properties (ue)
 391   // [0]=toward, [1]=away, [2]=transverse,
 392   //
 393   AliTkEtaPhiVector **centers=new AliTkEtaPhiVector*[4];
 394
 395   TH1F ***hgJetFragLeadingPt = new TH1F**[3];
 396   TH1F ***hgJetFragL = new TH1F**[3];
 397   TH1F ***hgJetFragPL = new TH1F**[3];
 398   TH1F ***hgJetFragT = new TH1F**[3];
 399   TH1F ***hgJetFragPt = new TH1F**[3];
 400   TH1F ***hgDiJetFragLeadingPt = new TH1F**[3];
 401   TH1F ***hgDiJetFragL = new TH1F**[3];
 402   TH1F ***hgDiJetFragPL = new TH1F**[3];
 403   TH1F ***hgDiJetFragT = new TH1F**[3];
 404   TH1F ***hgDiJetFragPt = new TH1F**[3];
 405   for(Int_t k=0;k<3;k++){
 406     hgJetFragLeadingPt[k] = new TH1F*[nclasses];
 407     hgJetFragL[k] = new TH1F*[nclasses];
 408     hgJetFragPL[k] = new TH1F*[nclasses];
 409     hgJetFragT[k] = new TH1F*[nclasses];
 410     hgJetFragPt[k] = new TH1F*[nclasses];
 411     hgDiJetFragLeadingPt[k] = new TH1F*[nclasses];
 412     hgDiJetFragL[k] = new TH1F*[nclasses];
 413     hgDiJetFragPL[k] = new TH1F*[nclasses];
 414     hgDiJetFragT[k] = new TH1F*[nclasses];
 415     hgDiJetFragPt[k] = new TH1F*[nclasses];
 416     if(k==0){
 417       sprintf(name,"toward");
 418     } else if (k==1) {
 419       sprintf(name,"away");
 420     } else {
 421       sprintf(name,"transverse");
 422     }
 423     for(Int_t i=0;i<nclasses;i++){
 424       sprintf(dummy,"h%s-JetFragLeadingPt%d",name,i);
 425       sprintf(title,"Fragmentation Function (%.0f GeV<E_{T}<%.0f GeV)",cletmin[i],cletmax[i]);
 426       hgJetFragLeadingPt[k][i]=new TH1F(dummy,title,100,0,1);
 427       hgJetFragLeadingPt[k][i]->GetXaxis()->SetTitle("z=p_{T}/p_{T}^{lead.part.}");
 428       hgJetFragLeadingPt[k][i]->GetYaxis()->SetTitle("dN/dz");
 429
 430       sprintf(dummy,"h%s-JetFragL%d",name,i);
 431       sprintf(title,"Longitudinal Fragmentation Function (%.0f GeV<E_{T}<%.0f GeV)",cletmin[i],cletmax[i]);
 432       hgJetFragL[k][i]=new TH1F(dummy,title,100,0,1);
 433       hgJetFragL[k][i]->GetXaxis()->SetTitle("z=p_{L}/P^{jet}");
 434       hgJetFragL[k][i]->GetYaxis()->SetTitle("dN/dz");
 435
 436       sprintf(dummy,"h%s-JetFragPL%d",name,i);
 437       sprintf(title,"Longitudinal Momentum (%.0f GeV<E_{T}<%.0f GeV)",cletmin[i],cletmax[i]);
 438       hgJetFragPL[k][i]=new TH1F(dummy,title,600,0,150);
 439       hgJetFragPL[k][i]->GetXaxis()->SetTitle("p_{L}");
 440       hgJetFragPL[k][i]->GetYaxis()->SetTitle("dN/dp_{L} [GeV^{-1}]");
 441
 442       sprintf(dummy,"h%s-JetFragT%d",name,i);
 443       sprintf(title,"Transversal Fragmentation Function (%.0f GeV<E_{T}<%.0f GeV)",cletmin[i],cletmax[i]);
 444       hgJetFragT[k][i]=new TH1F(dummy,title,250,0,25);
 445       hgJetFragT[k][i]->GetXaxis()->SetTitle("j_{T} [GeV]");
 446       hgJetFragT[k][i]->GetYaxis()->SetTitle("dN/dj_{T}");
 447
 448       sprintf(dummy,"h%s-JetFragPt%d",name,i);
 449       sprintf(title,"Particle Transversal Momentum Distribution (%.0f GeV<E_{T}<%.0f GeV)",cletmin[i],cletmax[i]);
 450       hgJetFragPt[k][i]=new TH1F(dummy,title,600,0,150);
 451       hgJetFragPt[k][i]->GetXaxis()->SetTitle("p_{T} [GeV]");
 452       hgJetFragPt[k][i]->GetYaxis()->SetTitle("dN/dp_{T}");
 453
 454       sprintf(dummy,"h%s-DiJetFragLeadingPt%d",name,i);
 455       sprintf(title,"Fragmentation Function (%.0f GeV<E_{T}<%.0f GeV)",cletmin[i],cletmax[i]);
 456       hgDiJetFragLeadingPt[k][i]=new TH1F(dummy,title,100,0,1);
 457       hgDiJetFragLeadingPt[k][i]->GetXaxis()->SetTitle("z=p_{T}/p_{T}^{lead.part.}");
 458       hgDiJetFragLeadingPt[k][i]->GetYaxis()->SetTitle("dN/dz");
 459
 460       sprintf(dummy,"h%s-DiJetFragL%d",name,i);
 461       sprintf(title,"Longitudinal Fragmentation Function (%.0f GeV<E_{T}<%.0f GeV)",cletmin[i],cletmax[i]);
 462       hgDiJetFragL[k][i]=new TH1F(dummy,title,100,0,1);
 463       hgDiJetFragL[k][i]->GetXaxis()->SetTitle("z=p_{L}/P^{jet}");
 464       hgDiJetFragL[k][i]->GetYaxis()->SetTitle("dN/dz");
 465
 466       sprintf(dummy,"h%s-DiJetFragPL%d",name,i);
 467       sprintf(title,"Longitudinal Momentum (%.0f GeV<E_{T}<%.0f GeV)",cletmin[i],cletmax[i]);
 468       hgDiJetFragPL[k][i]=new TH1F(dummy,title,600,0,150);
 469       hgDiJetFragPL[k][i]->GetXaxis()->SetTitle("p_{L}");
 470       hgDiJetFragPL[k][i]->GetYaxis()->SetTitle("dN/dp_{L} [GeV^{-1}]");
 471
 472       sprintf(dummy,"h%s-DiJetFragT%d",name,i);
 473       sprintf(title,"Transversal Fragmentation Function (%.0f GeV<E_{T}<%.0f GeV)",cletmin[i],cletmax[i]);
 474       hgDiJetFragT[k][i]=new TH1F(dummy,title,250,0,25);
 475       hgDiJetFragT[k][i]->GetXaxis()->SetTitle("j_{T} [GeV]");
 476       hgDiJetFragT[k][i]->GetYaxis()->SetTitle("dN/dj_{T}");
 477
 478       sprintf(dummy,"h%s-DiJetFragPt%d",name,i);
 479       sprintf(title,"Particle Transversal Momentum Distribution (%.0f GeV<E_{T}<%.0f GeV)",cletmin[i],cletmax[i]);
 480       hgDiJetFragPt[k][i]=new TH1F(dummy,title,600,0,150);
 481       hgDiJetFragPt[k][i]->GetXaxis()->SetTitle("p_{T} [GeV]");
 482       hgDiJetFragPt[k][i]->GetYaxis()->SetTitle("dN/dp_{T}");
 483     }
 484   }
 485
 486   //intshape
 487   Float_t gretall[3][nclasses][10];
 488   Float_t gretlow[3][nclasses][10];
 489   Float_t gretlow1[3][nclasses][10];
 490   Float_t gretlow2[3][nclasses][10];
 491   Float_t gretlow3[3][nclasses][10];
 492   Float_t gretlow4[3][nclasses][10];
 493   Float_t grxet[10];
 494   for(Int_t i=0;i<10;i++) {
 495     for(Int_t j=0;j<nclasses;j++) {
 496       for(Int_t k=0;k<3;k++){
 497         gretall[k][j][i]=0.;
 498         gretlow[k][j][i]=0.;
 499         gretlow1[k][j][i]=0.;
 500         gretlow2[k][j][i]=0.;
 501         gretlow3[k][j][i]=0.;
 502         gretlow4[k][j][i]=0.;
 503       }
 504     }
 505     grxet[i]=(i+1.)/10;
 506   }
 507
 508   //diffshape
 509   Float_t gdretall[3][nclasses][11];
 510   Float_t gdretlow[3][nclasses][11];
 511   Float_t gdretlow1[3][nclasses][11];
 512   Float_t gdretlow2[3][nclasses][11];
 513   Float_t gdretlow3[3][nclasses][11];
 514   Float_t gdretlow4[3][nclasses][11];
 515   Float_t gdrxet[11];
 516   for(Int_t i=0;i<11;i++) {
 517     for(Int_t j=0;j<nclasses;j++) {
 518       for(Int_t k=0;k<3;k++){
 519         gdretall[k][j][i]=0.;
 520         gdretlow[k][j][i]=0.;
 521         gdretlow1[k][j][i]=0.;
 522         gdretlow2[k][j][i]=0.;
 523         gdretlow3[k][j][i]=0.;
 524         gdretlow4[k][j][i]=0.;
 525       }
 526     }
 527     gdrxet[i]=1.*i/10.;
 528   }
 529
 530   //intshape
 531   Float_t gdiretall[3][nclasses][10];
 532   Float_t gdiretlow[3][nclasses][10];
 533   Float_t gdiretlow1[3][nclasses][10];
 534   Float_t gdiretlow2[3][nclasses][10];
 535   Float_t gdiretlow3[3][nclasses][10];
 536   Float_t gdiretlow4[3][nclasses][10];
 537   for(Int_t i=0;i<10;i++) {
 538     for(Int_t j=0;j<nclasses;j++) {
 539       for(Int_t k=0;k<3;k++){
 540         gdiretall[k][j][i]=0.;
 541         gdiretlow[k][j][i]=0.;
 542         gdiretlow1[k][j][i]=0.;
 543         gdiretlow2[k][j][i]=0.;
 544         gdiretlow3[k][j][i]=0.;
 545         gdiretlow4[k][j][i]=0.;
 546       }
 547     }
 548   }
 549
 550   //diffshape
 551   Float_t gdidretall[3][nclasses][11];
 552   Float_t gdidretlow[3][nclasses][11];
 553   Float_t gdidretlow1[3][nclasses][11];
 554   Float_t gdidretlow2[3][nclasses][11];
 555   Float_t gdidretlow3[3][nclasses][11];
 556   Float_t gdidretlow4[3][nclasses][11];
 557   for(Int_t i=0;i<11;i++) {
 558     for(Int_t j=0;j<nclasses;j++) {
 559       for(Int_t k=0;k<3;k++){
 560         gdidretall[k][j][i]=0.;
 561         gdidretlow[k][j][i]=0.;
 562         gdidretlow1[k][j][i]=0.;
 563         gdidretlow2[k][j][i]=0.;
 564         gdidretlow3[k][j][i]=0.;
 565         gdidretlow4[k][j][i]=0.;
 566       }
 567     }
 568   }
 569
 570   TH1F *hPartPtDist = new TH1F("hPartPtDist","Transverse Momentum Distribution",600,0,150);
 571   hPartPtDist->GetXaxis()->SetTitle("p_{T} [GeV]");
 572   hPartPtDist->GetYaxis()->SetTitle("dN/dp_{T} [GeV^{-1}]");
 573
 574   TH1F *hPartPhiDist = new TH1F("hPartPhiDist","Azimuthal Distribution",120,0,TMath::TwoPi());
 575   hPartPhiDist->GetXaxis()->SetTitle("#phi");
 576   hPartPhiDist->GetYaxis()->SetTitle("dN/d#phi");
 577
 578   TH1F *hPartEtaDist = new TH1F("hPartEtaDist","Pseudo-Rapidity Distribution",100,-1,1);
 579   hPartEtaDist->GetXaxis()->SetTitle("#eta");
 580   hPartEtaDist->GetYaxis()->SetTitle("dN/d#eta");
 581
 582   TH1F *hPartPhiCorr = new TH1F("hPartPhiCorr","Azimuthal Correlation",120,0,TMath::TwoPi());
 583   hPartPhiCorr->GetXaxis()->SetTitle("#phi");
 584   hPartPhiCorr->GetYaxis()->SetTitle("dN/d#phi");
 585
 586   TH1F *hPartDiPhiCorr = new TH1F("hPartDiPhiCorr","Azimuthal Correlation",120,0,TMath::TwoPi());
 587   hPartDiPhiCorr->GetXaxis()->SetTitle("#phi");
 588   hPartDiPhiCorr->GetYaxis()->SetTitle("dN/d#phi");
 589
 590   TH1F *hPartACorr = new TH1F("hPartACorr","Alpha Correlation",100,-1.1,1.1);
 591   hPartACorr->GetXaxis()->SetTitle("#alpha");
 592   hPartACorr->GetYaxis()->SetTitle("dN/d#alpha");
 593
 594   TH1F *hPartDiACorr = new TH1F("hPartDiACorr","Alpha Correlation",100,-1.1,1.1);
 595   hPartDiACorr->GetXaxis()->SetTitle("#alpha");
 596   hPartDiACorr->GetYaxis()->SetTitle("dN/d#alpha");
 597
 598   Int_t nTotalEvents = 0;
 599   Int_t nclGoodEvents[nclasses];
 600   for(Int_t i=0;i<nclasses;i++) nclGoodEvents[i]=0;
 601   Int_t nclLeadEvents[nclasses];
 602   for(Int_t i=0;i<nclasses;i++) nclLeadEvents[i]=0;
 603   Int_t nclDiEvents[nclasses];
 604   for(Int_t i=0;i<nclasses;i++) nclDiEvents[i]=0;
 605
 606   //connect to jets
 607   TChain *theTree = new TChain("jets");
 608   theTree->Add(filename);
 609   AliTkConeJetEvent *event = new AliTkConeJetEvent();
 610   theTree->SetBranchAddress("ConeFinder",&event);
 611
 612   Int_t treeentries=(Int_t)theTree->GetEntries();
 613   if((nMaxEvents<0) || (nMaxEvents>treeentries))
 614     nMaxEvents=treeentries;
 615
 616   cout << "Found " << nMaxEvents << " in " << filename << endl;
 617
 618   TChain *theEvTree=0;
 619   AliJetEventParticles *jetevent=0;
 620   TChain *theEvTree_sig=0;
 621   AliJetEventParticles *jetevent_sig=0;
 622   Int_t backtreeentries=0;
 623   Int_t nPerBackground=0;
 624   if(nPerBackground==0) nPerBackground=1;
 625
 626   if(evfilename){
 627     theEvTree = new TChain("AJEPtree");
 628     theEvTree->Add(evfilename);
 629     jetevent=new AliJetEventParticles();
 630     theEvTree->SetBranchAddress("particles",&jetevent);
 631     Int_t evtreeentries=(Int_t)theEvTree->GetEntries();
 632
 633     if(sigevfilename){
 634       theEvTree_sig = new TChain("AJEPtree");
 635       theEvTree_sig->Add(sigevfilename);
 636       jetevent_sig=new AliJetEventParticles();
 637       theEvTree_sig->SetBranchAddress("particles",&jetevent_sig);
 638       evtreeentries=(Int_t)theEvTree_sig->GetEntries();
 639       backtreeentries=(Int_t)theEvTree->GetEntries();
 640       nPerBackground=nMaxEvents/backtreeentries;
 641       if(nPerBackground==0) nPerBackground=1;
 642     }
 643
 644     if(evtreeentries!=treeentries){
 645       cerr << "WARNING: ";
 646       cerr << "Total jet event number not equals event number: "
 647            << evtreeentries << " " << treeentries <<endl;
 648       exit(1);
 649     }
 650   }
 651
 652   TChain *theSignalTree=0;
 653   AliJetEventParticles *jetsigevent=0;
 654   if(signalfilename){
 655     theSignalTree = new TChain("AJEPtree");
 656     theSignalTree->Add(signalfilename);
 657     jetsigevent=new AliJetEventParticles();
 658     theSignalTree->SetBranchAddress("particles",&jetsigevent);
 659
 660     Int_t sigtreeentries=(Int_t)theSignalTree->GetEntries();
 661     if(sigtreeentries!=treeentries){
 662       cerr << "WARNING: ";
 663       cerr << "Total jet signal event number not equals event number: "
 664            << sigtreeentries << " " << treeentries <<endl;
 665       //exit(1);
 666     }
 667   }
 668
 669   //connect to monte jets (if wanted)
 670   TChain *theTreeMonte=0;
 671   AliTkConeJetEvent *monteevent = 0;
 672   if(monteconefilename){
 673     theTreeMonte = new TChain("jets");
 674     theTreeMonte->Add(monteconefilename);
 675     monteevent= new AliTkConeJetEvent();
 676     theTreeMonte->SetBranchAddress("ConeFinder",&monteevent);
 677     Int_t montetreeentries=(Int_t)theTreeMonte->GetEntries();
 678     if(montetreeentries!=treeentries){
 679       cerr << "WARNING: ";
 680       cerr << "Total monte jet event number not equals jet event number: "
 681            << montetreeentries << " " << treeentries <<endl;
 682       //exit(1);
 683       if(theSignalTree){
 684         Int_t sigtreeentries=(Int_t)theSignalTree->GetEntries();
 685         if(montetreeentries!=sigtreeentries){
 686           cerr << "WARNING: ";
 687           cerr << "Total signal cone jet event number not equals signal event number: "
 688                << sigtreeentries << " " << montetreeentries <<endl;
 689           exit(1);
 690         }
 691       }
 692     }
 693   }
 694
 695   //=========================================================================
 696   // start the event loop
 697   //=========================================================================
 698   Int_t nEvent = -1;
 699   Int_t nEventSig = -1;
 700   Int_t nEventHijing = -1;
 701   Int_t nEventHijingCounter = nPerBackground;
 702   while(nEvent<nMaxEvents-1){
 703     nEvent++;
 704     nEventSig++;
 705
 706     if ((nEvent % 100) == 0) {
 707       cout << "Analysing event " << nEvent << endl;
 708     }
 709
 710     //connect the cone event/jets
 711     event->Clear();
 712     theTree->GetEvent(nEvent);
 713     event->sortJets();
 714     //event->Print("");
 715     Float_t ptcutused=event->getPtCut();
 716     if(ptcutused<minPartPt) ptcutused=minPartPt;
 717
 718     TClonesArray *jets=event->getJets();
 719
 720     if(theEvTree_sig){ // need to mix
 721       if(nEventHijingCounter==nPerBackground){
 722         jetevent->Reset();
 723         nEventHijing++;
 724         theEvTree->GetEvent(nEventHijing);
 725         //jetevent->Print();
 726         if(nEventHijing==backtreeentries) nEventHijing=0;
 727         nEventHijingCounter=0;
 728       }
 729       jetevent_sig->Reset();
 730       theEvTree_sig->GetEvent(nEvent);
 731       //jetevent_sig->Print();
 732       jetevent->AddSignal(*jetevent_sig);
 733       TString dummy="Counter: ";
 734       dummy+=nEvent;
 735       dummy+=" ";
 736       dummy+="(Pythia ";dummy+=nEvent;
 737       dummy+=" ";
 738       dummy+=", Hijing ";dummy+=nEventHijing;
 739       dummy+=")";
 740       jetevent->SetHeader(dummy);
 741       nEventHijingCounter++;
 742     } else if(theEvTree){
 743       theEvTree->GetEvent(nEvent);
 744     }
 745     else{
 746       jetevent=event->getJetParticles();
 747     }
 748
 749     if(theSignalTree){ //get MC event containing signal
 750       jetsigevent->Reset();
 751       theSignalTree->GetEvent(nEventSig);
 752 #if 1
 753       if(jetsigevent->GetEventNr()!=jetevent->GetEventNr()){
 754         Int_t js=jetsigevent->GetEventNr();
 755         Int_t es=jetevent->GetEventNr();
 756         cerr << "Need to skip event: " << nEvent << ": " <<  es  << " != " << js << endl;
 757         if(es>js) nEventSig++;
 758         else nEvent++;
 759         continue;
 760       }
 761 #endif
 762     }
 763     else {
 764       jetsigevent=jetevent;
 765     }
 766
 767     //connect the monte cone jets
 768     TClonesArray *montejets=0;
 769     if(theTreeMonte){
 770       monteevent->Clear();
 771       theTreeMonte->GetEvent(nEventSig);
 772       monteevent->sortJets();
 773       //monteevent->Print("");
 774       montejets=monteevent->getJets();
 775     }
 776
 777     //
 778     //jetevent->Print();
 779     //
 780
 781     Int_t njets=jets->GetEntries();
 782     if(!jets || !njets){
 783       cerr << "No Cone jet found in event " << nEvent << ", continuing..." <<endl;
 784       continue;
 785     }
 786
 787     Int_t nhard=0;    //get hard partons
 788     TClonesArray *chard_jets=new TClonesArray("AliTkConeJet",0);
 789     Float_t phard[4];
 790     Float_t type;
 791     //TString header=jetsigevent->getHeader();
 792     //Int_t ntrials=jetsigevent->Trials();
 793     jetsigevent->Hard(0,phard,type);
 794     if(type!=0){
 795       Float_t ptj,thj,etaj,phj;
 796       convert(phard,ptj,thj,etaj,phj);
 797       new ((*chard_jets)[nhard]) AliTkConeJet(ptj,etaj,phj,(Int_t)type);
 798       nhard++;
 799     }
 800     jetsigevent->Hard(1,phard,type);
 801     if(type!=0){
 802       Float_t ptj,thj,etaj,phj;
 803       convert(phard,ptj,thj,etaj,phj);
 804       new ((*chard_jets)[nhard]) AliTkConeJet(ptj,etaj,phj,(Int_t)type);
 805       nhard++;
 806     }
 807     chard_jets->Sort();
 808     //chard_jets->Print();
 809
 810     Int_t ntr=jetsigevent->NTriggerJets();
 811     Float_t tr_jets[ntr][3];    // trigger jets
 812     TClonesArray *ctr_jets=new TClonesArray("AliTkConeJet",ntr);
 813     Float_t lead_tr_pt=-1.,lead_tr_phi=0,lead_tr_eta=0.;
 814     for(Int_t j=0;j<ntr;j++){
 815       Float_t pjet[4];
 816       Float_t ptj,thj,etaj,phj;
 817       jetsigevent->TriggerJet(j,pjet);
 818       convert(pjet,ptj,thj,etaj,phj);
 819       tr_jets[j][0]=ptj;
 820       tr_jets[j][1]=phj;
 821       tr_jets[j][2]=etaj;
 822       new ((*ctr_jets)[j]) AliTkConeJet(ptj,etaj,phj);
 823       if(lead_tr_pt<ptj) {
 824          lead_tr_pt=ptj;
 825          lead_tr_phi=phj;
 826          lead_tr_eta=etaj;
 827       }
 828       for(Int_t i=0;i<nhard;i++){
 829         Float_t diff,etdiff,etadiff,phidiff;
 830         AliTkConeJet *jet=(AliTkConeJet*)ctr_jets->At(j);
 831         diff=AliTkConeJet::Diff(jet,(AliTkConeJet*)chard_jets->At(i),etdiff,phidiff,etadiff);
 832         //cout << diff << " " << etdiff << " " << phidiff << " " << etadiff << endl;
 833         if(diff<0.25) {
 834           jet->setType((Int_t)type);
 835           break;
 836         }
 837       }
 838     }
 839     ctr_jets->Sort();
 840     //ctr_jets->Print();
 841     Int_t ntruq=jetsigevent->NUQTriggerJets();
 842     Float_t uq_jets[ntruq][3];  // unquenched jets
 843     TClonesArray *cuq_jets=new TClonesArray("AliTkConeJet",ntruq);
 844     Float_t lead_uq_pt=-1.,lead_uq_phi=0,lead_uq_eta=0.;
 845     for(Int_t j=0;j<ntruq;j++){
 846       Float_t pjet[4];
 847       Float_t ptj,thj,etaj,phj;
 848       jetsigevent->UQJet(j,pjet);
 849       convert(pjet,ptj,thj,etaj,phj);
 850       uq_jets[j][0]=ptj;
 851       uq_jets[j][1]=phj;
 852       uq_jets[j][2]=etaj;
 853       new ((*cuq_jets)[j]) AliTkConeJet(ptj,etaj,phj);
 854       if(lead_uq_pt<ptj) {
 855         lead_uq_pt=ptj;
 856         lead_uq_phi=phj;
 857         lead_uq_eta=etaj;
 858       }
 859     }
 860     cuq_jets->Sort();
 861     //cuq_jets->Print();
 862
 863     Double_t x0=jetsigevent->GetXJet();
 864     Double_t y0=jetsigevent->GetYJet();
 865     Double_t prodlength=TMath::Sqrt(x0*x0+y0*y0);
 866     Double_t zquench[4];
 867     jetsigevent->GetZQuench(zquench);
 868
 869     //fill trigger histos (without cuts)
 870     hJetEtUQTrigger->Fill(lead_uq_pt);
 871     hJetEtTrigger->Fill(lead_tr_pt);
 872     for(Int_t i=0;i<2;i++){
 873       hJetZ->Fill(zquench[i]);
 874     }
 875     //reconstruction efficiency plots
 876     if(ntr) for(Int_t j=0;j<1/*ntr*/;j++){ //loop over MC jets
 877       Float_t diff,etdiff,etadiff,phidiff;
 878       Float_t mindiff=100;
 879       Int_t index=-1;
 880       AliTkConeJet *jt=(AliTkConeJet*)ctr_jets->At(j);
 881       if(jt->getEt()<minEtRecoCut||jt->getEt()>maxEtRecoCut) continue;
 882       //if(TMath:Abs(jt->getEta())>0.1) continue,
 883       //jets in event
 884       Int_t njets=jets->GetEntries();
 885       AliTkConeJet *jet=0;
 886       for(Int_t i=0;i<njets;i++){
 887         jet=(AliTkConeJet*)jets->At(i);
 888         if(!jet) continue;
 889         jet->calculateValues();
 890         diff=AliTkConeJet::Diff(jt,jet,etdiff,phidiff,etadiff);
 891         //_cut_ here if wanted (prob. for mixed events)
 892         //if(phidiff>0.25||etadiff>0.25) continue;
 893         //if(TMath::Abs(etdiff)/jt->getEt()>0.1) continue;
 894         //cout << diff << " " << etdiff << " " << phidiff << " " << etadiff << endl;
 895         if(mindiff>diff){
 896           mindiff=diff;
 897           index=i;
 898         }
 899       }
 900       if(index>-1) {
 901         jet=(AliTkConeJet*)jets->At(index);
 902         hJetEtres->Fill(jet->getEt());
 903         hJetEtresTrigger->Fill(jt->getEt());
 904         hJetEtratio->Fill(jet->getEt()/jt->getEt());
 905         Float_t phidiff=(jt->getPhi()-jet->getPhi());
 906         if(phidiff>TMath::Pi()) phidiff=TMath::TwoPi()-phidiff;
 907         Float_t etadiff=(jt->getEta()-jet->getEta());
 908         Float_t etfract=jet->getEtMarkedFrac();
 909         hAxesDiffres->Fill(phidiff,etadiff,TMath::Abs(jet->getEt()-jt->getEt()));
 910         hJetEtrestrue->Fill(jet->getEt(),etfract);
 911         hPhires->Fill(phidiff);
 912         hEtares->Fill(etadiff);
 913         //cout << "test" << nEvent << " " << phidiff << " " << etadiff << endl;
 914       }
 915     }
 916     if(theTreeMonte){
 917       //compare leading jets from Monte and Reconstruction
 918       AliTkConeJet *jt=(AliTkConeJet*)jets->At(0);
 919       AliTkConeJet *mjt=(AliTkConeJet*)montejets->At(0);
 920       if(jt && mjt){
 921         Float_t phidiff=(jt->getPhi()-mjt->getPhi());
 922         if(phidiff>TMath::Pi()) phidiff=TMath::TwoPi()-phidiff;
 923         Float_t etadiff=(jt->getEta()-mjt->getEta());
 924         Float_t etdiff=jt->getEt()-mjt->getEt();
 925         hPhiMonteres->Fill(phidiff);
 926         hEtaMonteres->Fill(etadiff);
 927         hEtMonteres->Fill(etdiff);
 928         hEtMonteratio->Fill(jt->getEt()/mjt->getEt());
 929       }
 930       for(Int_t j=0;j<1/*montejets->Entries()*/;j++){
 931         Float_t diff,etdiff,etadiff,phidiff;
 932         Float_t mindiff=100;
 933         Int_t index=-1;
 934         AliTkConeJet *jt=(AliTkConeJet*)montejets->At(j);
 935         if(!jt) continue;
 936         //if(jt->getEt()<minEtRecoCut||jt->getEt()>maxEtRecoCut) continue;
 937         //jets in event
 938         Int_t njets=jets->GetEntries();
 939         AliTkConeJet *jet=0;
 940         for(Int_t i=0;i<njets;i++){
 941           jet=(AliTkConeJet*)jets->At(i);
 942           if(!jet) continue;
 943           jet->calculateValues();
 944           diff=AliTkConeJet::Diff(jt,jet,etdiff,phidiff,etadiff);
 945           //_cut_ here if wanted (prob. for mixed events)
 946           //if(phidiff>0.25||etadiff>0.25) continue;
 947           //if(TMath::Abs(etdiff)/jt->getEt()>0.1) continue;
 948           //cout << diff << " " << etdiff << " " << phidiff << " " << etadiff << endl;
 949           if(mindiff>diff){
 950             mindiff=diff;
 951             index=i;
 952           }
 953         }
 954         if(index>-1) {
 955           jet=(AliTkConeJet*)jets->At(index);
 956           hmJetEtres->Fill(jet->getEt());
 957           hmJetEtratio->Fill(jet->getEt()/jt->getEt());
 958           Float_t phidiff=(jt->getPhi()-jet->getPhi());
 959           if(phidiff>TMath::Pi()) phidiff=TMath::TwoPi()-phidiff;
 960           Float_t etadiff=(jt->getEta()-jet->getEta());
 961           Float_t etfract=jet->getEtMarkedFrac();
 962           hmAxesDiffres->Fill(phidiff,etadiff,TMath::Abs(jet->getEt()-jt->getEt()));
 963           hmJetEtrestrue->Fill(jet->getEt(),etfract);
 964           hmPhires->Fill(phidiff);
 965           hmEtares->Fill(etadiff);
 966         }
 967       }
 968     }
 969
 970     //could _cut_ on event
 971 #if 0
 972     //cout << lead_tr_pt << " " << lead_tr_eta << " " << lead_tr_phi << endl;
 973     if((lead_tr_eta>0.5) || (lead_tr_eta<-0.5)) continue;
 974     if(zquench[0]<0.1||zquench[1]<0.1) continue;
 975 #endif
 976
 977     //particles in event
 978     const TClonesArray *parts=jetevent->GetParticles();
 979
 980     //jets in event
 981     AliTkConeJet *lead_jet=0;
 982     AliTkConeJet *back_jet=0;
 983     for(Int_t i=0;i<njets;i++){
 984       AliTkConeJet* jet=0;
 985       Int_t whichjet=0;
 986       jet=(AliTkConeJet*)jets->At(i);
 987       if(!jet) continue;
 988       jet->calculateValues();
 989
 990       hJetEtall->Fill(jet->getEt()); //without any cuts
 991       Float_t etfract=jet->getEtMarkedFrac();
 992       hJetEtalltrue->Fill(jet->getEt(),etfract); //without any cuts
 993
 994       //here could _cut_ on jet
 995       //-------
 996       Float_t et=jet->getEt();
 997       Float_t corret=et;
 998       if(corrfac>0) corret/=corrfac;
 999       if(corret<minJetEt) continue;
1000
1001       Int_t clindex=Int_t(corret/10);
1002       if(clindex>nclasses-2) clindex=nclasses-2;
1003       nclGoodEvents[clindex]++;
1004       nclGoodEvents[nclasses-1]++;
1005       //-------
1006
1007       //set jet type
1008       for(Int_t i=0;i<nhard;i++){
1009         Float_t diff,etdiff,etadiff,phidiff;
1010         AliTkConeJet *jh=(AliTkConeJet*)chard_jets->At(i);
1011         diff=AliTkConeJet::Diff(jet,jh,etdiff,phidiff,etadiff);
1012         //cout << diff << " " << etdiff << " " << phidiff << " " << etadiff << endl;
1013         if(diff<0.25) {
1014           jet->setType(jh->getType());
1015           break;
1016         }
1017       }
1018
1019       //check leading jet
1020       if(!lead_jet){
1021         lead_jet=jet;
1022         whichjet=1;
1023         hJetEt->Fill(lead_jet->getEt());
1024         hJetEttrue->Fill(lead_jet->getEt(),etfract);
1025         Float_t mindiff=100;
1026         Float_t minetdiff=0.,minphidiff=0.,minetadiff=0.;
1027         Int_t index;
1028         for(Int_t j=0;j<ntr;j++){
1029           Float_t diff,etdiff,etadiff,phidiff;
1030           AliTkConeJet *jt=(AliTkConeJet*)ctr_jets->At(j);
1031           diff=AliTkConeJet::Diff(lead_jet,jt,etdiff,phidiff,etadiff);
1032           if(mindiff>diff){
1033             mindiff=diff;
1034             index=j;
1035             minphidiff=phidiff;
1036             minetadiff=etadiff;
1037             minetdiff=etdiff;
1038           }
1039         }
1040         if(TMath::Abs(minetdiff)/lead_jet->getEt()<0.15)
1041           hAxesDiff->Fill(minphidiff,minetadiff);
1042
1043         //trigger
1044         Float_t triget=lead_jet->getEt();
1045         Float_t leadet=((AliTkConeJet*)ctr_jets->At(0))->getEt();
1046         for(Int_t i=0;i<9;i++){
1047           Float_t minet=i*10+10;
1048           if(triget>=minet) hJetEttrigger[i]->Fill(leadet,1);
1049           else hJetEttrigger2[i]->Fill(leadet,1);
1050         }
1051         hJetEttriggernorm->Fill(leadet,1);
1052       }// check the back-to-back jet
1053       else if ((jet->getEt()/lead_jet->getEt()>0.75) &&
1054                 (relphi(jet->getPhi(),lead_jet->getPhi())>5/6*TMath::Pi())){
1055         if(!back_jet) {
1056           whichjet=2;
1057           back_jet=jet;
1058           hBackJetEt->Fill(back_jet->getEt());
1059           hBackJetEttrue->Fill(back_jet->getEt(),etfract);
1060           hJetEtvsEt->Fill(lead_jet->getEt(),back_jet->getEt());
1061
1062         } else{
1063           cerr << "Already found one back jet, disregarding the new one." << endl;
1064         }
1065       }
1066
1067       //fill properties
1068       hJetEtvsTrigger->Fill(corret,lead_tr_pt);
1069       hJetEtvsUQTrigger->Fill(corret,lead_uq_pt);
1070       hJetEtallvsEll->Fill(corret,prodlength);
1071       hJetEtvsEll->Fill(corret,prodlength);
1072       hJettype[jet->getType()]->Fill(jet->getEt());
1073
1074       Int_t njetparts=jet->getNParts();
1075       hJetN[clindex]->Fill(njetparts);
1076       hJetN[nclasses-1]->Fill(njetparts);
1077       Float_t leadPartPt=jet->getPtLead();
1078       Float_t ptRatio=0;
1079       if(corret>0) ptRatio=leadPartPt/corret;
1080       Float_t meanpt=0.;
1081       if(njetparts>0) meanpt=corret/njetparts;
1082       hJetMeanPt[clindex]->Fill(meanpt);
1083       hJetMeanPt[nclasses-1]->Fill(meanpt);
1084       hJetLeadingPt[clindex]->Fill(ptRatio);
1085       hJetLeadingPt[nclasses-1]->Fill(ptRatio);
1086       hJetLeadingPtDist[clindex]->Fill(leadPartPt);
1087       hJetLeadingPtDist[nclasses-1]->Fill(leadPartPt);
1088
1089       Float_t jetAxisX=jet->getXAxis();
1090       Float_t jetAxisY=jet->getYAxis();
1091       Float_t jetAxisZ=jet->getZAxis();
1092       Float_t jetAxisLength=jet->getPLength();
1093       if(jetAxisLength) {
1094         TClonesArray *particles = jet->getParts();
1095         TIterator *partit = particles->MakeIterator();
1096         TParticle *particle = NULL;
1097         Int_t firstval=0;
1098         while ((particle = (TParticle *) partit->Next()) != NULL) {
1099           Float_t al=particle->Px()*jetAxisX+particle->Py()*jetAxisY+particle->Pz()*jetAxisZ;
1100           if(al<0){
1101             //cout << "Should not happen!" << al << endl;
1102             continue;
1103           }
1104           Float_t at=TMath::Sqrt(TMath::Abs(particle->P()*particle->P()-al*al));
1105           hJetFragL[clindex]->Fill(al/jetAxisLength);
1106           hJetFragT[clindex]->Fill(at);
1107           hJetFragPL[clindex]->Fill(al);
1108           hJetFragPt[clindex]->Fill(particle->Pt());
1109           hJetFragL[nclasses-1]->Fill(al/jetAxisLength);
1110           hJetFragT[nclasses-1]->Fill(at);
1111           hJetFragPL[nclasses-1]->Fill(al);
1112           hJetFragPt[nclasses-1]->Fill(particle->Pt());
1113           if(leadPartPt>0){
1114             Float_t z=particle->Pt()/leadPartPt;
1115             if(firstval||z<1){
1116               hJetFragLeadingPt[clindex]->Fill(z);
1117               hJetFragLeadingPt[nclasses-1]->Fill(z);
1118             } else if(z==1) firstval=1;
1119           }
1120         }
1121         delete partit;
1122       }
1123
1124       //calculate int/diff shapes
1125       Float_t jeteta=jet->getCEta();
1126       Float_t jetphi=jet->getCPhi();
1127       AliTkEtaPhiVector center(jeteta,jetphi);
1128       TIterator *partit = parts->MakeIterator();
1129       AliJetParticle *particle = NULL;
1130       while ((particle = (AliJetParticle *) partit->Next()) != NULL) {
1131         Float_t pt=particle->Pt();
1132         if(pt<ptcutused) continue;
1133         AliTkEtaPhiVector centerp;
1134         centerp.setVector(particle->Eta(),particle->Phi());
1135         Float_t diff=center.diff(centerp);
1136         for(Int_t loop=Int_t(diff*10)+1;loop<=10;loop++){
1137           if(pt<1.) retlow1[clindex][loop-1]+=pt;
1138           if(pt<2.) retlow[clindex][loop-1]+=pt;
1139           if(pt<3.) retlow2[clindex][loop-1]+=pt;
1140           if(pt<4.) retlow3[clindex][loop-1]+=pt;
1141           if(pt<5.) retlow4[clindex][loop-1]+=pt;
1142           retall[clindex][loop-1]+=pt;
1143           if(pt<1.) retlow1[nclasses-1][loop-1]+=pt;
1144           if(pt<2.) retlow[nclasses-1][loop-1]+=pt;
1145           if(pt<3.) retlow2[nclasses-1][loop-1]+=pt;
1146           if(pt<4.) retlow3[nclasses-1][loop-1]+=pt;
1147           if(pt<5.) retlow4[nclasses-1][loop-1]+=pt;
1148           retall[nclasses-1][loop-1]+=pt;
1149         }
1150         if(diff<=1){
1151           Int_t index=0;
1152           if(diff<=deltaR) index=0;
1153           else index=Int_t((diff-deltaR)*10)+1;
1154           if(pt<1.) dretlow1[clindex][index]+=pt;
1155           if(pt<2.) dretlow[clindex][index]+=pt;
1156           if(pt<3.) dretlow2[clindex][index]+=pt;
1157           if(pt<4.) dretlow3[clindex][index]+=pt;
1158           if(pt<5.) dretlow4[clindex][index]+=pt;
1159           dretall[clindex][index]+=pt;
1160           if(pt<1.) dretlow1[nclasses-1][index]+=pt;
1161           if(pt<2.) dretlow[nclasses-1][index]+=pt;
1162           if(pt<3.) dretlow2[nclasses-1][index]+=pt;
1163           if(pt<4.) dretlow3[nclasses-1][index]+=pt;
1164           if(pt<5.) dretlow4[nclasses-1][index]+=pt;
1165           dretall[nclasses-1][index]+=pt;
1166         }
1167         if(pt<ptcutused) continue;
1168         Float_t phi=particle->Phi();
1169
1170         Float_t dphi=diffphi(jet->getPhi(),phi);
1171         hPhiCorr[clindex]->Fill(dphi);
1172         hPhiCorr[nclasses-1]->Fill(dphi);
1173       }
1174       delete partit;
1175
1176     } //loop over cone jets
1177
1178     //global event studies
1179     TIterator *partit = parts->MakeIterator();
1180     AliJetParticle *particle = NULL;
1181     Float_t leta=0.;
1182     Float_t lphi=0.;
1183     Int_t clindex=-1;
1184     Float_t jetAxisX=-1;
1185     Float_t jetAxisY=-1;
1186     Float_t jetAxisZ=-1;
1187     Float_t jetAxisLength=-1;
1188     Float_t leadPartPt=-1;
1189     if(lead_jet){
1190       leta=lead_jet->getCEta();
1191       lphi=lead_jet->getCPhi();
1192       jetAxisX=lead_jet->getXAxis();
1193       jetAxisY=lead_jet->getYAxis();
1194       jetAxisZ=lead_jet->getZAxis();
1195       jetAxisLength=lead_jet->getPLength();
1196       leadPartPt=lead_jet->getPtLead();
1197       centers[0]=new AliTkEtaPhiVector(leta,lphi);
1198       centers[1]=new AliTkEtaPhiVector(leta,addphi(lphi,TMath::Pi()));
1199       centers[2]=new AliTkEtaPhiVector(leta,addphi(lphi,TMath::Pi()/2));
1200       centers[3]=new AliTkEtaPhiVector(leta,addphi(lphi,3*TMath::Pi()/2));
1201       Float_t et=lead_jet->getEt();
1202       Float_t corret=et;
1203       if(corrfac>0) corret/=corrfac;
1204       clindex=Int_t(corret/10);
1205       if(clindex>nclasses-2) clindex=nclasses-2;
1206       nclLeadEvents[clindex]++;
1207       nclLeadEvents[nclasses-1]++;
1208       if(back_jet){
1209         nclDiEvents[clindex]++;
1210         nclDiEvents[nclasses-1]++;
1211       }
1212     }
1213
1214     //loop over particles in event
1215     Int_t firstval=0;
1216     Int_t difirstval=0;
1217     while ((particle = (AliJetParticle *) partit->Next()) != NULL) {
1218       Float_t pt=particle->Pt();
1219       if(pt<ptcutused) continue;
1220       Float_t eta=particle->Eta();
1221       Float_t phi=particle->Phi();
1222       hPartPtDist->Fill(pt);
1223       hPartPhiDist->Fill(phi);
1224       hPartEtaDist->Fill(eta);
1225       if(lead_jet){
1226         Float_t dphi=diffphi(lphi,phi);
1227         Float_t al=particle->Px()*jetAxisX+particle->Py()*jetAxisY+particle->Pz()*jetAxisZ;
1228         al/=particle->P();
1229         hPartPhiCorr->Fill(dphi);
1230         hPartACorr->Fill(al);
1231         if(back_jet){
1232           hPartDiPhiCorr->Fill(dphi);
1233           hPartDiACorr->Fill(al);
1234         }
1235
1236         //ue studies
1237         for(Int_t ceindex=0;ceindex<4;ceindex++){
1238           Int_t heindex=ceindex;
1239           if(heindex==3) heindex=2; //store both sides of trans plane in one histo
1240
1241           AliTkEtaPhiVector centerp;
1242           centerp.setVector(eta,phi);
1243           Float_t diff=centers[ceindex]->diff(centerp);
1244           for(Int_t loop=Int_t(diff*10)+1;loop<=10;loop++){
1245             if(pt<1.) gretlow1[heindex][clindex][loop-1]+=pt;
1246             if(pt<2.) gretlow[heindex][clindex][loop-1]+=pt;
1247             if(pt<3.) gretlow2[heindex][clindex][loop-1]+=pt;
1248             if(pt<4.) gretlow3[heindex][clindex][loop-1]+=pt;
1249             if(pt<5.) gretlow4[heindex][clindex][loop-1]+=pt;
1250             gretall[heindex][clindex][loop-1]+=pt;
1251             if(pt<1.) gretlow1[heindex][nclasses-1][loop-1]+=pt;
1252             if(pt<2.) gretlow[heindex][nclasses-1][loop-1]+=pt;
1253             if(pt<3.) gretlow2[heindex][nclasses-1][loop-1]+=pt;
1254             if(pt<4.) gretlow3[heindex][nclasses-1][loop-1]+=pt;
1255             if(pt<5.) gretlow4[heindex][nclasses-1][loop-1]+=pt;
1256             gretall[heindex][nclasses-1][loop-1]+=pt;
1257           }
1258           if(diff<=1){
1259             Int_t index=0;
1260             if(diff<=deltaR) index=0;
1261             else index=Int_t((diff-deltaR)*10)+1;
1262             if(pt<1.) gdretlow1[heindex][clindex][index]+=pt;
1263             if(pt<2.) gdretlow[heindex][clindex][index]+=pt;
1264             if(pt<3.) gdretlow2[heindex][clindex][index]+=pt;
1265             if(pt<4.) gdretlow3[heindex][clindex][index]+=pt;
1266             if(pt<5.) gdretlow4[heindex][clindex][index]+=pt;
1267             gdretall[heindex][clindex][index]+=pt;
1268             if(pt<1.) gdretlow1[heindex][nclasses-1][index]+=pt;
1269             if(pt<2.) gdretlow[heindex][nclasses-1][index]+=pt;
1270             if(pt<3.) gdretlow2[heindex][nclasses-1][index]+=pt;
1271             if(pt<4.) gdretlow3[heindex][nclasses-1][index]+=pt;
1272             if(pt<5.) gdretlow4[heindex][nclasses-1][index]+=pt;
1273             gdretall[heindex][nclasses-1][index]+=pt;
1274           }
1275           if(back_jet){
1276             for(Int_t loop=Int_t(diff*10)+1;loop<=10;loop++){
1277               if(pt<1.) gdiretlow1[heindex][clindex][loop-1]+=pt;
1278               if(pt<2.) gdiretlow[heindex][clindex][loop-1]+=pt;
1279               if(pt<3.) gdiretlow2[heindex][clindex][loop-1]+=pt;
1280               if(pt<4.) gdiretlow3[heindex][clindex][loop-1]+=pt;
1281               if(pt<5.) gdiretlow4[heindex][clindex][loop-1]+=pt;
1282               gdiretall[heindex][clindex][loop-1]+=pt;
1283               if(pt<1.) gdiretlow1[heindex][nclasses-1][loop-1]+=pt;
1284               if(pt<2.) gdiretlow[heindex][nclasses-1][loop-1]+=pt;
1285               if(pt<3.) gdiretlow2[heindex][nclasses-1][loop-1]+=pt;
1286               if(pt<4.) gdiretlow3[heindex][nclasses-1][loop-1]+=pt;
1287               if(pt<5.) gdiretlow4[heindex][nclasses-1][loop-1]+=pt;
1288               gdiretall[heindex][nclasses-1][loop-1]+=pt;
1289             }
1290             if(diff<=1){
1291               Int_t index=0;
1292               if(diff<=deltaR) index=0;
1293               else index=Int_t((diff-deltaR)*10)+1;
1294               if(pt<1.) gdidretlow1[heindex][clindex][index]+=pt;
1295               if(pt<2.) gdidretlow[heindex][clindex][index]+=pt;
1296               if(pt<3.) gdidretlow2[heindex][clindex][index]+=pt;
1297               if(pt<4.) gdidretlow3[heindex][clindex][index]+=pt;
1298               if(pt<5.) gdidretlow4[heindex][clindex][index]+=pt;
1299               gdidretall[heindex][clindex][index]+=pt;
1300               if(pt<1.) gdretlow1[heindex][nclasses-1][index]+=pt;
1301               if(pt<2.) gdretlow[heindex][nclasses-1][index]+=pt;
1302               if(pt<3.) gdidretlow2[heindex][nclasses-1][index]+=pt;
1303               if(pt<4.) gdidretlow3[heindex][nclasses-1][index]+=pt;
1304               if(pt<5.) gdidretlow4[heindex][nclasses-1][index]+=pt;
1305               gdidretall[heindex][nclasses-1][index]+=pt;
1306             }
1307           }
1308         }//ceindex
1309
1310         Int_t ceindex=eventindex(lphi,phi);
1311         Int_t heindex=ceindex;
1312         if(heindex==3) heindex=2; //store both sides of trans plane in one histo
1313
1314         if(jetAxisLength) {
1315           Float_t jetX=jetAxisX;
1316           Float_t jetY=jetAxisY;
1317           if(ceindex==1) {
1318             jetX=-jetAxisX;
1319             jetY=-jetAxisY;
1320           } else if(ceindex==2) {
1321             jetX=-jetAxisY;
1322             jetY=jetAxisX;
1323           }
1324           else if(ceindex==3) {
1325             jetX=jetAxisY;
1326             jetY=-jetAxisX;
1327           }
1328           Float_t al=particle->Px()*jetX+particle->Py()*jetY+particle->Pz()*jetAxisZ;
1329           Float_t at=TMath::Sqrt(TMath::Abs(particle->P()*particle->P()-al*al));
1330           hgJetFragL[heindex][clindex]->Fill(al/jetAxisLength);
1331           hgJetFragT[heindex][clindex]->Fill(at);
1332           hgJetFragPL[heindex][clindex]->Fill(al);
1333           hgJetFragPt[heindex][clindex]->Fill(particle->Pt());
1334           hgJetFragL[heindex][nclasses-1]->Fill(al/jetAxisLength);
1335           hgJetFragT[heindex][nclasses-1]->Fill(at);
1336           hgJetFragPL[heindex][nclasses-1]->Fill(al);
1337           hgJetFragPt[heindex][nclasses-1]->Fill(particle->Pt());
1338           if(leadPartPt>0){
1339             Float_t z=particle->Pt()/leadPartPt;
1340             if(firstval||z<1){
1341               hgJetFragLeadingPt[heindex][clindex]->Fill(z);
1342               hgJetFragLeadingPt[heindex][nclasses-1]->Fill(z);
1343             } else if(z==1) firstval=1;
1344           }
1345           if(back_jet){
1346             hgDiJetFragL[heindex][clindex]->Fill(al/jetAxisLength);
1347             hgDiJetFragT[heindex][clindex]->Fill(at);
1348             hgDiJetFragPL[heindex][clindex]->Fill(al);
1349             hgDiJetFragPt[heindex][clindex]->Fill(particle->Pt());
1350             hgDiJetFragL[heindex][nclasses-1]->Fill(al/jetAxisLength);
1351             hgDiJetFragT[heindex][nclasses-1]->Fill(at);
1352             hgDiJetFragPL[heindex][nclasses-1]->Fill(al);
1353             hgDiJetFragPt[heindex][nclasses-1]->Fill(particle->Pt());
1354             if(leadPartPt>0){
1355               Float_t z=particle->Pt()/leadPartPt;
1356               if(difirstval||z<1){
1357                 hgDiJetFragLeadingPt[heindex][clindex]->Fill(z);
1358                 hgDiJetFragLeadingPt[heindex][nclasses-1]->Fill(z);
1359               } else if(z==1) difirstval=1;
1360             }
1361           }
1362         }
1363       } //lead_jet
1364
1365     }
1366     delete partit;
1367     if(lead_jet)
1368       for(Int_t i=0;i<4;i++) delete centers[i];
1369
1370     nTotalEvents++;
1371   } //end of nev loop
1372
1373   cout << "Finished analysing events " << nTotalEvents << endl;
1374
1375   delete centers;
1376   delete event;
1377   delete theTree;
1378   if(theEvTree){
1379     delete jetevent;
1380     delete theEvTree;
1381   }
1382   if(theSignalTree){
1383     delete jetsigevent;
1384     delete theSignalTree;
1385   }
1386   if(theTreeMonte){
1387     delete monteevent;
1388     delete theTreeMonte;
1389   }
1390
1391   //========================================================================
1392   // draw histograms
1393   //========================================================================
1394
1395   TTimeStamp tstamp;
1396   Char_t timestamp[255];
1397   sprintf(timestamp,"%d",tstamp.GetDate(0,0));
1398
1399   //store all objects in root file (you never know, when you need it)
1400   Char_t rootfilename[1024];
1401   if(savefilename!=NULL)
1402     sprintf(rootfilename,"%s",savefilename);
1403   else
1404     sprintf(rootfilename,"%s/anaAliJets-%s.root",gSystem->Getenv("JF_PLOTSDIR"),timestamp);
1405   TFile *rootfile=new TFile(rootfilename,"RECREATE");
1406
1407   // let's start with the drawing...
1408   Int_t nents=(Int_t)hJetEt->GetEntries();
1409   if(nents){
1410     cout << "hJetEt: " << nents << " entries in histogram" << endl;
1411     hJetEt->SetLineWidth(3);
1412     hJetEt->Write();
1413   }
1414
1415   nents=(Int_t)hJetEttrue->GetEntries();
1416   if(nents){
1417     cout << "hJetEttrue: " << nents << " entries in histogram" << endl;
1418     hJetEttrue->SetLineWidth(3);
1419     hJetEttrue->Write();
1420   }
1421
1422   nents=(Int_t)hBackJetEt->GetEntries();
1423   if(nents){
1424     cout << "hBackJetEt: " << nents << " entries in histogram" << endl;
1425     hBackJetEt->SetLineWidth(3);
1426     hBackJetEt->Write();
1427   }
1428
1429   nents=(Int_t)hBackJetEttrue->GetEntries();
1430   if(nents){
1431     cout << "hBackJetEttrue: " << nents << " entries in histogram" << endl;
1432     hBackJetEttrue->SetLineWidth(3);
1433     hBackJetEttrue->Write();
1434   }
1435
1436   nents=(Int_t)hJetEtall->GetEntries();
1437   if(nents){
1438     cout << "hJetEtall: " << nents << " entries in histogram" << endl;
1439     hJetEtall->SetLineWidth(3);
1440     hJetEtall->Write();
1441   }
1442
1443   nents=(Int_t)hJetEtalltrue->GetEntries();
1444   if(nents){
1445     cout << "hJetEtalltrue: " << nents << " entries in histogram" << endl;
1446     hJetEtalltrue->SetLineWidth(3);
1447     hJetEtalltrue->Write();
1448   }
1449
1450   nents=(Int_t)hJetEtUQTrigger->GetEntries();
1451   if(nents){
1452     cout << "hJetEtUQTrigger: " << nents << " entries in histogram" << endl;
1453     hJetEtUQTrigger->SetLineWidth(3);
1454     hJetEtUQTrigger->Write();
1455   }
1456   nents=(Int_t)hJetEtTrigger->GetEntries();
1457   if(nents){
1458     cout << "hJetEtTrigger: " << nents << " entries in histogram" << endl;
1459     hJetEtTrigger->SetLineWidth(3);
1460     hJetEtTrigger->Write();
1461   }
1462
1463   rootfile->cd();
1464   rootfile->mkdir("LeadingPt");
1465   rootfile->cd("LeadingPt");
1466   for(Int_t i=0;i<nclasses;i++){
1467     nents=(Int_t)hJetLeadingPt[i]->GetEntries();
1468     if(nents){
1469       //cout << "hJetLeadingPt " << i << ": "  << nents << " entries in histogram" << endl;
1470       hJetLeadingPt[i]->SetLineWidth(3);
1471       hJetLeadingPt[i]->Write();
1472     }
1473   }
1474
1475   rootfile->cd();
1476   rootfile->mkdir("FragLeadingPt");
1477   rootfile->cd("FragLeadingPt");
1478   for(Int_t i=0;i<nclasses;i++){
1479     nents=(Int_t)hJetFragLeadingPt[i]->GetEntries();
1480     if(nents){
1481       //cout << "hJetFragLeadingPt " << i << ": "  << nents << " entries in histogram" << endl;
1482       hJetFragLeadingPt[i]->SetLineWidth(3);
1483       hJetFragLeadingPt[i]->Write();
1484     }
1485   }
1486
1487   rootfile->cd();
1488   rootfile->mkdir("LeadingPtDist");
1489   rootfile->cd("LeadingPtDist");
1490   for(Int_t i=0;i<nclasses;i++){
1491     nents=(Int_t)hJetLeadingPtDist[i]->GetEntries();
1492     if(nents){
1493       //cout << "hJetLeadingPtDist " << i << ": " << nents << " entries in histogram" << endl;
1494       hJetLeadingPtDist[i]->SetLineWidth(3);
1495       hJetLeadingPtDist[i]->Write();
1496     }
1497   }
1498
1499   rootfile->cd();
1500   rootfile->mkdir("FragLong");
1501   rootfile->cd("FragLong");
1502   for(Int_t i=0;i<nclasses;i++){
1503     nents=(Int_t)hJetFragL[i]->GetEntries();
1504     if(nents){
1505       //cout << "hJetFragL " << i << ": " << nents << " entries in histogram " << endl;
1506       hJetFragL[i]->SetLineWidth(3);
1507       hJetFragL[i]->Write();
1508     }
1509   }
1510
1511   rootfile->cd();
1512   rootfile->mkdir("FragPL");
1513   rootfile->cd("FragPL");
1514   for(Int_t i=0;i<nclasses;i++){
1515     nents=(Int_t)hJetFragPL[i]->GetEntries();
1516     if(nents){
1517       //cout << "hJetFragPL " << i << ": " << nents << " entries in histogram " << endl;
1518       hJetFragPL[i]->SetLineWidth(3);
1519       hJetFragPL[i]->Write();
1520     }
1521   }
1522
1523   rootfile->cd();
1524   rootfile->mkdir("FragTrans");
1525   rootfile->cd("FragTrans");
1526   for(Int_t i=0;i<nclasses;i++){
1527     nents=(Int_t)hJetFragT[i]->GetEntries();
1528     if(nents){
1529       //cout << "hJetFragT " << i << ": " << nents << " entries in histogram" << endl;
1530       hJetFragT[i]->SetLineWidth(3);
1531       hJetFragT[i]->Write();
1532     }
1533   }
1534
1535   rootfile->cd();
1536   rootfile->mkdir("Multiplicity");
1537   rootfile->cd("Multiplicity");
1538   for(Int_t i=0;i<nclasses;i++){
1539     nents=(Int_t)hJetN[i]->GetEntries();
1540     //cout <<"hJetN " << i << ": " << nents << " entries in histogram " << endl;
1541     if(nents){
1542       hJetN[i]->SetLineWidth(3);
1543       hJetN[i]->Write();
1544     }
1545   }
1546
1547   rootfile->cd();
1548   rootfile->mkdir("PtDistribution");
1549   rootfile->cd("PtDistribution");
1550   for(Int_t i=0;i<nclasses;i++){
1551     nents=(Int_t)hJetFragPt[i]->GetEntries();
1552     if(nents){
1553       //cout << "hJetFragPt " << i << ": " << nents << " entries in histogram " << endl;
1554       hJetFragPt[i]->SetLineWidth(3);
1555       hJetFragPt[i]->Write();
1556     }
1557   }
1558
1559   rootfile->cd();
1560   rootfile->mkdir("MeanPt");
1561   rootfile->cd("MeanPt");
1562   for(Int_t i=0;i<nclasses;i++){
1563     nents=(Int_t)hJetMeanPt[i]->GetEntries();
1564     if(nents){
1565       //cout << "hJetMeanPt " << i << ": " << nents << " entries in histogram " << endl;
1566       hJetMeanPt[i]->SetLineWidth(3);
1567       hJetMeanPt[i]->Write();
1568     }
1569   }
1570
1571   rootfile->cd();
1572   rootfile->mkdir("PhiCorr");
1573   rootfile->cd("PhiCorr");
1574   for(Int_t i=0;i<nclasses;i++){
1575     nents=(Int_t)hPhiCorr[i]->GetEntries();
1576     if(nents){
1577       //cout << "hPhiCorr " << i << ": " << nents << " entries in histogram " << endl;
1578       hPhiCorr[i]->SetLineWidth(3);
1579       hPhiCorr[i]->Write();
1580     }
1581   }
1582
1583   rootfile->cd();//intshape
1584   rootfile->mkdir("ConeFluc");
1585   rootfile->cd("ConeFluc");
1586   for(Int_t j=0;j<nclasses;j++){
1587     TGraph *graphall=new TGraph(10);
1588     TGraph *graphlow=new TGraph(10);
1589     TGraph *graphlow1=new TGraph(10);
1590     TGraph *graphlow2=new TGraph(10);
1591     TGraph *graphlow3=new TGraph(10);
1592     TGraph *graphlow4=new TGraph(10);
1593     for(Int_t i=0;i<10;i++) {
1594       graphall->SetPoint(i,rxet[i],retall[j][i]);
1595       graphlow->SetPoint(i,rxet[i],retlow[j][i]);
1596       graphlow1->SetPoint(i,rxet[i],retlow1[j][i]);
1597       graphlow2->SetPoint(i,rxet[i],retlow2[j][i]);
1598       graphlow3->SetPoint(i,rxet[i],retlow3[j][i]);
1599       graphlow4->SetPoint(i,rxet[i],retlow4[j][i]);
1600     }
1601     sprintf(dummy,"gretall%d",j);
1602     graphall->Write(dummy);
1603     sprintf(dummy,"gretlow%d",j);
1604     graphlow->Write(dummy);
1605     sprintf(dummy,"gret2low%d",j);
1606     graphlow->Write(dummy);
1607     sprintf(dummy,"gret1low%d",j);
1608     graphlow1->Write(dummy);
1609     sprintf(dummy,"gret3low%d",j);
1610     graphlow2->Write(dummy);
1611     sprintf(dummy,"gret4low%d",j);
1612     graphlow3->Write(dummy);
1613     sprintf(dummy,"gret5low%d",j);
1614     graphlow4->Write(dummy);
1615     delete graphall;
1616     delete graphlow;
1617     delete graphlow1;
1618     delete graphlow2;
1619     delete graphlow3;
1620   }
1621
1622   rootfile->cd();
1623   rootfile->mkdir("Shape");
1624   rootfile->cd("Shape");
1625   for(Int_t j=0;j<nclasses;j++){
1626     TGraph *graphall=new TGraph(11);
1627     TGraph *graphlow=new TGraph(11);
1628     TGraph *graphlow1=new TGraph(11);
1629     TGraph *graphlow2=new TGraph(11);
1630     TGraph *graphlow3=new TGraph(11);
1631     TGraph *graphlow4=new TGraph(11);
1632     for(Int_t i=0;i<11;i++) {
1633       graphall->SetPoint(i,drxet[i],dretall[j][i]/deltaR);
1634       graphlow->SetPoint(i,drxet[i],dretlow[j][i]/deltaR);
1635       graphlow1->SetPoint(i,drxet[i],dretlow1[j][i]/deltaR);
1636       graphlow2->SetPoint(i,drxet[i],dretlow2[j][i]/deltaR);
1637       graphlow3->SetPoint(i,drxet[i],dretlow3[j][i]/deltaR);
1638       graphlow4->SetPoint(i,drxet[i],dretlow4[j][i]/deltaR);
1639     }
1640     sprintf(dummy,"gretall%d",j);
1641     graphall->Write(dummy);
1642     sprintf(dummy,"gretlow%d",j);
1643     graphlow->Write(dummy);
1644     sprintf(dummy,"gret2low%d",j);
1645     graphlow->Write(dummy);
1646     sprintf(dummy,"gret1low%d",j);
1647     graphlow1->Write(dummy);
1648     sprintf(dummy,"gret3low%d",j);
1649     graphlow2->Write(dummy);
1650     sprintf(dummy,"gret4low%d",j);
1651     graphlow3->Write(dummy);
1652     sprintf(dummy,"gret5low%d",j);
1653     graphlow4->Write(dummy);
1654     delete graphall;
1655     delete graphlow;
1656     delete graphlow1;
1657     delete graphlow2;
1658     delete graphlow3;
1659   }
1660
1661   rootfile->cd();
1662   rootfile->mkdir("Extended");
1663   rootfile->cd("Extended");
1664
1665   nents=(Int_t)hJetEtvsTrigger->GetEntries();
1666   if(nents){
1667     cout << "hJetEtvsTrigger" << nents << " entries in histogram" << endl;
1668     hJetEtvsTrigger->SetLineWidth(3);
1669     hJetEtvsTrigger->Write();
1670   }
1671
1672   nents=(Int_t)hJetEtvsEt->GetEntries();
1673   if(nents){
1674     cout << "hJetEtvsEt" << nents << " entries in histogram" << endl;
1675     hJetEtvsEt->SetLineWidth(3);
1676     hJetEtvsEt->Write();
1677   }
1678
1679   nents=(Int_t)hAxesDiff->GetEntries();
1680   if(nents){
1681     cout << "hAxesDiff" << nents << " entries in histogram" << endl;
1682     hAxesDiff->SetLineWidth(3);
1683     hAxesDiff->Write();
1684   }
1685
1686   nents=(Int_t)hJet1->GetEntries();
1687   if(nents){
1688     cout << "hJet1" << nents << " entries in histogram" << endl;
1689     hJet1->SetLineWidth(3);
1690     hJet1->Write();
1691   }
1692   nents=(Int_t)hJet2->GetEntries();
1693   if(nents){
1694     cout << "hJet2" << nents << " entries in histogram" << endl;
1695     hJet2->SetLineWidth(3);
1696     hJet2->Write();
1697   }
1698
1699   for(Int_t i=0;i<3;i++){
1700     nents=(Int_t)hJettype[i]->GetEntries();
1701     if(nents){
1702       hJettype[i]->SetLineWidth(3);
1703       hJettype[i]->Write();
1704     }
1705   }
1706
1707   nents=(Int_t)hJetEtvsEll->GetEntries();
1708   if(nents){
1709     cout << "hJetEtvsEll" << nents << " entries in histogram " << endl;
1710     hJetEtvsEll->SetLineWidth(3);
1711     hJetEtvsEll->Write();
1712   }
1713
1714   nents=(Int_t)hJetEtallvsEll->GetEntries();
1715   if(nents){
1716     cout << "hJetEtallvsEll" << nents << " entries in histogram " << endl;
1717     hJetEtallvsEll->SetLineWidth(3);
1718     hJetEtallvsEll->Write();
1719   }
1720
1721   nents=(Int_t)hJetZ->GetEntries();
1722   if(nents){
1723     cout << "hJetZ: " << nents << " entries in histogram " << endl;
1724     hJetZ->SetLineWidth(3);
1725     hJetZ->Write();
1726   }
1727
1728   rootfile->cd();
1729   rootfile->mkdir("Global");
1730   rootfile->cd("Global");
1731
1732   nents=(Int_t)hPartPtDist->GetEntries();
1733   if(nents){
1734     cout << "hPartPtDist: " << nents << " entries in histogram" << endl;
1735     hPartPtDist->SetLineWidth(3);
1736     hPartPtDist->Write();
1737   }
1738
1739   nents=(Int_t)hPartEtaDist->GetEntries();
1740   if(nents){
1741     cout << "hPartEtaDist: " << nents << " entries in histogram" << endl;
1742     hPartEtaDist->SetLineWidth(3);
1743     hPartEtaDist->Write();
1744   }
1745
1746   nents=(Int_t)hPartPhiDist->GetEntries();
1747   if(nents){
1748     cout << "hPartPhiDist: " << nents << " entries in histogram" << endl;
1749     hPartPhiDist->SetLineWidth(3);
1750     hPartPhiDist->Write();
1751   }
1752
1753   nents=(Int_t)hPartPhiCorr->GetEntries();
1754   if(nents){
1755     cout << "hPartPhiCorr: " << nents << " entries in histogram" << endl;
1756     hPartPhiCorr->SetLineWidth(3);
1757     hPartPhiCorr->Write();
1758   }
1759
1760   nents=(Int_t)hPartDiPhiCorr->GetEntries();
1761   if(nents){
1762     cout << "hPartDiPhiCorr: " << nents << " entries in histogram" << endl;
1763     hPartDiPhiCorr->SetLineWidth(3);
1764     hPartDiPhiCorr->Write();
1765   }
1766
1767   nents=(Int_t)hPartACorr->GetEntries();
1768   if(nents){
1769     cout << "hPartACorr: " << nents << " entries in histogram" << endl;
1770     hPartACorr->SetLineWidth(3);
1771     hPartACorr->Write();
1772   }
1773
1774   nents=(Int_t)hPartDiACorr->GetEntries();
1775   if(nents){
1776     cout << "hPartDiACorr: " << nents << " entries in histogram" << endl;
1777     hPartDiACorr->SetLineWidth(3);
1778     hPartDiACorr->Write();
1779   }
1780
1781   rootfile->cd();
1782   rootfile->mkdir("gConeFluc");
1783   rootfile->cd("gConeFluc");
1784   for(Int_t k=0;k<3;k++){
1785     if(k==0){
1786       sprintf(name,"toward");
1787     } else if (k==1) {
1788       sprintf(name,"away");
1789     } else {
1790       sprintf(name,"transverse");
1791     }
1792
1793     for(Int_t j=0;j<nclasses;j++){
1794       TGraph *graphall=new TGraph(10);
1795       TGraph *graphlow=new TGraph(10);
1796       TGraph *graphlow1=new TGraph(10);
1797       TGraph *graphlow2=new TGraph(10);
1798       TGraph *graphlow3=new TGraph(10);
1799       TGraph *graphlow4=new TGraph(10);
1800       for(Int_t i=0;i<10;i++) {
1801         graphall->SetPoint(i,grxet[i],gretall[k][j][i]);
1802         graphlow->SetPoint(i,grxet[i],gretlow[k][j][i]);
1803         graphlow1->SetPoint(i,grxet[i],gretlow1[k][j][i]);
1804         graphlow2->SetPoint(i,grxet[i],gretlow2[k][j][i]);
1805         graphlow3->SetPoint(i,grxet[i],gretlow3[k][j][i]);
1806         graphlow4->SetPoint(i,grxet[i],gretlow4[k][j][i]);
1807       }
1808       sprintf(dummy,"%s-gretall%d",name,j);
1809       graphall->Write(dummy);
1810       sprintf(dummy,"%s-gretlow%d",name,j);
1811       graphlow->Write(dummy);
1812       sprintf(dummy,"%s-gret2low%d",name,j);
1813       graphlow->Write(dummy);
1814       sprintf(dummy,"%s-gret1low%d",name,j);
1815       graphlow1->Write(dummy);
1816       sprintf(dummy,"%s-gret3low%d",name,j);
1817       graphlow2->Write(dummy);
1818       sprintf(dummy,"%s-gret4low%d",name,j);
1819       graphlow3->Write(dummy);
1820       sprintf(dummy,"%s-gret5low%d",name,j);
1821       graphlow4->Write(dummy);
1822       delete graphall;
1823       delete graphlow;
1824       delete graphlow1;
1825       delete graphlow2;
1826       delete graphlow3;
1827     }
1828   }
1829
1830   rootfile->cd();
1831   rootfile->mkdir("gShape");
1832   rootfile->cd("gShape");
1833   for(Int_t k=0;k<3;k++){
1834     if(k==0){
1835       sprintf(name,"toward");
1836     } else if (k==1) {
1837       sprintf(name,"away");
1838     } else {
1839       sprintf(name,"transverse");
1840     }
1841
1842     for(Int_t j=0;j<nclasses;j++){
1843       TGraph *graphall=new TGraph(10);
1844       TGraph *graphlow=new TGraph(10);
1845       TGraph *graphlow1=new TGraph(10);
1846       TGraph *graphlow2=new TGraph(10);
1847       TGraph *graphlow3=new TGraph(10);
1848       TGraph *graphlow4=new TGraph(10);
1849       for(Int_t i=0;i<10;i++) {
1850         graphall->SetPoint(i,gdrxet[i],gdretall[k][j][i]/deltaR);
1851         graphlow->SetPoint(i,gdrxet[i],gdretlow[k][j][i]/deltaR);
1852         graphlow1->SetPoint(i,gdrxet[i],gdretlow1[k][j][i]/deltaR);
1853         graphlow2->SetPoint(i,gdrxet[i],gdretlow2[k][j][i]/deltaR);
1854         graphlow3->SetPoint(i,gdrxet[i],gdretlow3[k][j][i]/deltaR);
1855         graphlow4->SetPoint(i,gdrxet[i],gdretlow4[k][j][i]/deltaR);
1856       }
1857       sprintf(dummy,"%s-gretall%d",name,j);
1858       graphall->Write(dummy);
1859       sprintf(dummy,"%s-gretlow%d",name,j);
1860       graphlow->Write(dummy);
1861       sprintf(dummy,"%s-gret2low%d",name,j);
1862       graphlow->Write(dummy);
1863       sprintf(dummy,"%s-gret1low%d",name,j);
1864       graphlow1->Write(dummy);
1865       sprintf(dummy,"%s-gret3low%d",name,j);
1866       graphlow2->Write(dummy);
1867       sprintf(dummy,"%s-gret4low%d",name,j);
1868       graphlow3->Write(dummy);
1869       sprintf(dummy,"%s-gret5low%d",name,j);
1870       graphlow4->Write(dummy);
1871       delete graphall;
1872       delete graphlow;
1873       delete graphlow1;
1874       delete graphlow2;
1875       delete graphlow3;
1876     }
1877   }
1878
1879   rootfile->cd();
1880   rootfile->mkdir("gDiConeFluc");
1881   rootfile->cd("gDiConeFluc");
1882   for(Int_t k=0;k<3;k++){
1883     if(k==0){
1884       sprintf(name,"toward");
1885     } else if (k==1) {
1886       sprintf(name,"away");
1887     } else {
1888       sprintf(name,"transverse");
1889     }
1890
1891     for(Int_t j=0;j<nclasses;j++){
1892       TGraph *graphall=new TGraph(10);
1893       TGraph *graphlow=new TGraph(10);
1894       TGraph *graphlow1=new TGraph(10);
1895       TGraph *graphlow2=new TGraph(10);
1896       TGraph *graphlow3=new TGraph(10);
1897       TGraph *graphlow4=new TGraph(10);
1898       for(Int_t i=0;i<10;i++) {
1899         graphall->SetPoint(i,grxet[i],gdiretall[k][j][i]);
1900         graphlow->SetPoint(i,grxet[i],gdiretlow[k][j][i]);
1901         graphlow1->SetPoint(i,grxet[i],gdiretlow1[k][j][i]);
1902         graphlow2->SetPoint(i,grxet[i],gdiretlow2[k][j][i]);
1903         graphlow3->SetPoint(i,grxet[i],gdiretlow3[k][j][i]);
1904         graphlow4->SetPoint(i,grxet[i],gdiretlow4[k][j][i]);
1905       }
1906       sprintf(dummy,"%s-gretall%d",name,j);
1907       graphall->Write(dummy);
1908       sprintf(dummy,"%s-gretlow%d",name,j);
1909       graphlow->Write(dummy);
1910       sprintf(dummy,"%s-gret2low%d",name,j);
1911       graphlow->Write(dummy);
1912       sprintf(dummy,"%s-gret1low%d",name,j);
1913       graphlow1->Write(dummy);
1914       sprintf(dummy,"%s-gret3low%d",name,j);
1915       graphlow2->Write(dummy);
1916       sprintf(dummy,"%s-gret4low%d",name,j);
1917       graphlow3->Write(dummy);
1918       sprintf(dummy,"%s-gret5low%d",name,j);
1919       graphlow4->Write(dummy);
1920       delete graphall;
1921       delete graphlow;
1922       delete graphlow1;
1923       delete graphlow2;
1924       delete graphlow3;
1925     }
1926   }
1927
1928   rootfile->cd();
1929   rootfile->mkdir("gDiShape");
1930   rootfile->cd("gDiShape");
1931   for(Int_t k=0;k<3;k++){
1932     if(k==0){
1933       sprintf(name,"toward");
1934     } else if (k==1) {
1935       sprintf(name,"away");
1936     } else {
1937       sprintf(name,"transverse");
1938     }
1939
1940     for(Int_t j=0;j<nclasses;j++){
1941       TGraph *graphall=new TGraph(10);
1942       TGraph *graphlow=new TGraph(10);
1943       TGraph *graphlow1=new TGraph(10);
1944       TGraph *graphlow2=new TGraph(10);
1945       TGraph *graphlow3=new TGraph(10);
1946       TGraph *graphlow4=new TGraph(10);
1947       for(Int_t i=0;i<10;i++) {
1948         graphall->SetPoint(i,gdrxet[i],gdidretall[k][j][i]/deltaR);
1949         graphlow->SetPoint(i,gdrxet[i],gdidretlow[k][j][i]/deltaR);
1950         graphlow1->SetPoint(i,gdrxet[i],gdidretlow1[k][j][i]/deltaR);
1951         graphlow2->SetPoint(i,gdrxet[i],gdidretlow2[k][j][i]/deltaR);
1952         graphlow3->SetPoint(i,gdrxet[i],gdidretlow3[k][j][i]/deltaR);
1953         graphlow4->SetPoint(i,gdrxet[i],gdidretlow4[k][j][i]/deltaR);
1954       }
1955       sprintf(dummy,"%s-gretall%d",name,j);
1956       graphall->Write(dummy);
1957       sprintf(dummy,"%s-gretlow%d",name,j);
1958       graphlow->Write(dummy);
1959       sprintf(dummy,"%s-gret2low%d",name,j);
1960       graphlow->Write(dummy);
1961       sprintf(dummy,"%s-gret1low%d",name,j);
1962       graphlow1->Write(dummy);
1963       sprintf(dummy,"%s-gret3low%d",name,j);
1964       graphlow2->Write(dummy);
1965       sprintf(dummy,"%s-gret4low%d",name,j);
1966       graphlow3->Write(dummy);
1967       sprintf(dummy,"%s-gret5low%d",name,j);
1968       graphlow4->Write(dummy);
1969       delete graphall;
1970       delete graphlow;
1971       delete graphlow1;
1972       delete graphlow2;
1973       delete graphlow3;
1974     }
1975   }
1976
1977   rootfile->cd();
1978   rootfile->mkdir("gFragLeadingPt");
1979   rootfile->cd("gFragLeadingPt");
1980   for(Int_t k=0;k<3;k++){
1981     for(Int_t i=0;i<nclasses;i++){
1982       nents=(Int_t)hgJetFragLeadingPt[k][i]->GetEntries();
1983       if(nents){
1984         //cout << "hJetFragLeadingPt " << i << ": "  << nents << " entries in histogram" << endl;
1985         hgJetFragLeadingPt[k][i]->SetLineWidth(3);
1986         hgJetFragLeadingPt[k][i]->Write();
1987       }
1988     }
1989   }
1990
1991   rootfile->cd();
1992   rootfile->mkdir("gFragLong");
1993   rootfile->cd("gFragLong");
1994   for(Int_t k=0;k<3;k++){
1995     for(Int_t i=0;i<nclasses;i++){
1996       nents=(Int_t)hgJetFragL[k][i]->GetEntries();
1997       if(nents){
1998         //cout << "hgJetFragL " << i << ": " << nents << " entries in histogram " << endl;
1999         hgJetFragL[k][i]->SetLineWidth(3);
2000         hgJetFragL[k][i]->Write();
2001       }
2002     }
2003   }
2004
2005   rootfile->cd();
2006   rootfile->mkdir("gFragPL");
2007   rootfile->cd("gFragPL");
2008   for(Int_t k=0;k<3;k++){
2009     for(Int_t i=0;i<nclasses;i++){
2010       nents=(Int_t)hgJetFragPL[k][i]->GetEntries();
2011       if(nents){
2012         //cout << "hJetFragPL " << i << ": " << nents << " entries in histogram " << endl;
2013         hgJetFragPL[k][i]->SetLineWidth(3);
2014         hgJetFragPL[k][i]->Write();
2015       }
2016     }
2017   }
2018
2019   rootfile->cd();
2020   rootfile->mkdir("gFragTrans");
2021   rootfile->cd("gFragTrans");
2022   for(Int_t k=0;k<3;k++){
2023     for(Int_t i=0;i<nclasses;i++){
2024       nents=(Int_t)hgJetFragT[k][i]->GetEntries();
2025       if(nents){
2026         //cout << "hJetFragT " << i << ": " << nents << " entries in histogram" << endl;
2027         hgJetFragT[k][i]->SetLineWidth(3);
2028         hgJetFragT[k][i]->Write();
2029       }
2030     }
2031   }
2032
2033   rootfile->cd();
2034   rootfile->mkdir("gDiFragLeadingPt");
2035   rootfile->cd("gDiFragLeadingPt");
2036   for(Int_t k=0;k<3;k++){
2037     for(Int_t i=0;i<nclasses;i++){
2038       nents=(Int_t)hgDiJetFragLeadingPt[k][i]->GetEntries();
2039       if(nents){
2040         //cout << "hJetFragLeadingPt " << i << ": "  << nents << " entries in histogram" << endl;
2041         hgDiJetFragLeadingPt[k][i]->SetLineWidth(3);
2042         hgDiJetFragLeadingPt[k][i]->Write();
2043       }
2044     }
2045   }
2046
2047   rootfile->cd();
2048   rootfile->mkdir("gDiFragLong");
2049   rootfile->cd("gDiFragLong");
2050   for(Int_t k=0;k<3;k++){
2051     for(Int_t i=0;i<nclasses;i++){
2052       nents=(Int_t)hgDiJetFragL[k][i]->GetEntries();
2053       if(nents){
2054         //cout << "hgDiJetFragL " << i << ": " << nents << " entries in histogram " << endl;
2055         hgDiJetFragL[k][i]->SetLineWidth(3);
2056         hgDiJetFragL[k][i]->Write();
2057       }
2058     }
2059   }
2060
2061   rootfile->cd();
2062   rootfile->mkdir("gDiFragPL");
2063   rootfile->cd("gDiFragPL");
2064   for(Int_t k=0;k<3;k++){
2065     for(Int_t i=0;i<nclasses;i++){
2066       nents=(Int_t)hgDiJetFragPL[k][i]->GetEntries();
2067       if(nents){
2068         //cout << "hgDJetFragPL " << i << ": " << nents << " entries in histogram " << endl;
2069         hgDiJetFragPL[k][i]->SetLineWidth(3);
2070         hgDiJetFragPL[k][i]->Write();
2071       }
2072     }
2073   }
2074
2075   rootfile->cd();
2076   rootfile->mkdir("gDiFragTrans");
2077   rootfile->cd("gDiFragTrans");
2078   for(Int_t k=0;k<3;k++){
2079     for(Int_t i=0;i<nclasses;i++){
2080       nents=(Int_t)hgDiJetFragT[k][i]->GetEntries();
2081       if(nents){
2082         //cout << "hJetFragT " << i << ": " << nents << " entries in histogram" << endl;
2083         hgDiJetFragT[k][i]->SetLineWidth(3);
2084         hgDiJetFragT[k][i]->Write();
2085       }
2086     }
2087   }
2088
2089   //reconstruction
2090   rootfile->cd();
2091   rootfile->mkdir("reconstruction");
2092   rootfile->cd("reconstruction");
2093
2094   nents=(Int_t)hJetEtres->GetEntries();
2095   if(nents){
2096     cout << "hJetEtres " << nents << " entries in histogram" << endl;
2097     hJetEtres->SetLineWidth(3);
2098     hJetEtres->Write();
2099   }
2100
2101   nents=(Int_t)hJetEtratio->GetEntries();
2102   if(nents){
2103     cout << "hJetEtratio " << nents << " entries in histogram" << endl;
2104     hJetEtratio->SetLineWidth(3);
2105     hJetEtratio->Write();
2106   }
2107
2108   nents=(Int_t)hJetEtrestrue->GetEntries();
2109   if(nents){
2110     cout << "hJetEtrestrue " << nents << " entries in histogram" << endl;
2111     hJetEtrestrue->SetLineWidth(3);
2112     hJetEtrestrue->Write();
2113   }
2114
2115   nents=(Int_t)hJetEtresTrigger->GetEntries();
2116   if(nents){
2117     cout << "hJetEtresTrigger " << nents << " entries in histogram" << endl;
2118     hJetEtresTrigger->SetLineWidth(3);
2119     hJetEtresTrigger->Write();
2120   }
2121
2122   nents=(Int_t)hAxesDiffres->GetEntries();
2123   if(nents){
2124     cout << "hAxesDiffres " << nents << " entries in histogram" << endl;
2125     hAxesDiffres->SetLineWidth(3);
2126     hAxesDiffres->Write();
2127   }
2128
2129   nents=(Int_t)hPhires->GetEntries();
2130   if(nents){
2131     cout << "hPhires " << nents << " entries in histogram" << endl;
2132     hPhires->SetLineWidth(3);
2133     hPhires->Write();
2134   }
2135
2136   nents=(Int_t)hEtares->GetEntries();
2137   if(nents){
2138     cout << "hEtares " << nents << " entries in histogram" << endl;
2139     hEtares->SetLineWidth(3);
2140     hEtares->Write();
2141   }
2142
2143   nents=(Int_t)hEtMonteres->GetEntries();
2144   if(nents){
2145     cout << "hEtMonteres " << nents << " entries in histogram" << endl;
2146     hEtMonteres->SetLineWidth(3);
2147     hEtMonteres->Write();
2148   }
2149
2150   nents=(Int_t)hEtaMonteres->GetEntries();
2151   if(nents){
2152     cout << "hEtaMonteres " << nents << " entries in histogram" << endl;
2153     hEtaMonteres->SetLineWidth(3);
2154     hEtaMonteres->Write();
2155   }
2156
2157   nents=(Int_t)hPhiMonteres->GetEntries();
2158   if(nents){
2159     cout << "hPhiMonteres " << nents << " entries in histogram" << endl;
2160     hPhiMonteres->SetLineWidth(3);
2161     hPhiMonteres->Write();
2162   }
2163
2164   nents=(Int_t)hEtMonteratio->GetEntries();
2165   if(nents){
2166     cout << "hEtMonteratio " << nents << " entries in histogram" << endl;
2167     hEtMonteratio->SetLineWidth(3);
2168     hEtMonteratio->Write();
2169   }
2170
2171   nents=(Int_t)hmJetEtres->GetEntries();
2172   if(nents){
2173     cout << "hmJetEtres " << nents << " entries in histogram" << endl;
2174     hmJetEtres->SetLineWidth(3);
2175     hmJetEtres->Write();
2176   }
2177
2178   nents=(Int_t)hmJetEtratio->GetEntries();
2179   if(nents){
2180     cout << "hmJetEtratio " << nents << " entries in histogram" << endl;
2181     hmJetEtratio->SetLineWidth(3);
2182     hmJetEtratio->Write();
2183   }
2184
2185   nents=(Int_t)hmJetEtrestrue->GetEntries();
2186   if(nents){
2187     cout << "hmJetEtrestrue " << nents << " entries in histogram" << endl;
2188     hmJetEtrestrue->SetLineWidth(3);
2189     hmJetEtrestrue->Write();
2190   }
2191
2192   nents=(Int_t)hmAxesDiffres->GetEntries();
2193   if(nents){
2194     cout << "hmAxesDiffres " << nents << " entries in histogram" << endl;
2195     hmAxesDiffres->SetLineWidth(3);
2196     hmAxesDiffres->Write();
2197   }
2198
2199   nents=(Int_t)hmPhires->GetEntries();
2200   if(nents){
2201     cout << "hmPhires " << nents << " entries in histogram" << endl;
2202     hmPhires->SetLineWidth(3);
2203     hmPhires->Write();
2204   }
2205
2206   nents=(Int_t)hmEtares->GetEntries();
2207   if(nents){
2208     cout << "hmEtares " << nents << " entries in histogram" << endl;
2209     hmEtares->SetLineWidth(3);
2210     hmEtares->Write();
2211   }
2212
2213   rootfile->cd();
2214   rootfile->mkdir("trigger");
2215   rootfile->cd("trigger");
2216   for(Int_t i=0;i<9;i++){
2217     hJetEttrigger[i]->Write();
2218     hJetEttrigger2[i]->Write();
2219   }
2220   hJetEttriggernorm->Write();
2221
2222   //store event info
2223   rootfile->cd();
2224   TGraph *graph=new TGraph(nclasses+1);
2225   cout << "Good jets counters" << endl;
2226   for(Int_t i=0;i<nclasses;i++){
2227     cout << i << " " << nclGoodEvents[i] << endl;
2228     graph->SetPoint(i,i,nclGoodEvents[i]);
2229   }
2230   graph->SetPoint(nclasses,nclasses,nTotalEvents);
2231   cout << nclasses << " " << nTotalEvents << endl;
2232   graph->Write("ggoodevents");
2233   delete graph;
2234   graph=new TGraph(nclasses+1);
2235   for(Int_t i=0;i<nclasses;i++){
2236     //cout << i << " " << nclLeadEvents[i] << endl;
2237     graph->SetPoint(i,i,nclLeadEvents[i]);
2238   }
2239   graph->SetPoint(nclasses,nclasses,nTotalEvents);
2240   graph->Write("ggoodlead");
2241   delete graph;
2242   graph=new TGraph(nclasses+1);
2243   for(Int_t i=0;i<nclasses;i++){
2244     //cout << i << " " << nclDiEvents[i] << endl;
2245     graph->SetPoint(i,i,nclDiEvents[i]);
2246   }
2247   graph->SetPoint(nclasses,nclasses,nTotalEvents);
2248   graph->Write("ggooddilead");
2249   delete graph;
2250
2251   //close the result file
2252   rootfile->Close();
2253   //
2254
2255   delete hJetEt;
2256   delete hBackJetEt;
2257   delete hJetEtall;
2258   delete hJetEttrue;
2259   delete hBackJetEttrue;
2260   delete hJetEtalltrue;
2261   delete hJetEtTrigger;
2262   delete hJetEtUQTrigger;
2263   delete hJetEtvsTrigger;
2264   delete hJetEtvsUQTrigger;
2265   delete hJetEtvsEt;
2266   delete hAxesDiff;
2267   delete hPartPtDist;
2268   delete hPartEtaDist;
2269   delete hPartPhiDist;
2270   delete hPartPhiCorr;
2271   delete hPartACorr;
2272   delete hPartDiACorr;
2273   delete hPartDiPhiCorr;
2274   delete hJet1;
2275   delete hJet2;
2276   for(Int_t i=0;i<3;i++) delete hJettype[i];
2277   delete[] hJettype;
2278   delete hJetEtvsEll;
2279   delete hJetEtallvsEll;
2280   delete hJetZ;
2281
2282   delete hJetEtres;
2283   delete hJetEtratio;
2284   delete hJetEtrestrue;
2285   delete hJetEtresTrigger;
2286   delete hAxesDiffres;
2287   delete hPhires;
2288   delete hEtares;
2289   delete hmJetEtres;
2290   delete hmJetEtratio;
2291   delete hmJetEtrestrue;
2292   delete hmAxesDiffres;
2293   delete hmPhires;
2294   delete hmEtares;
2295   delete hPhiMonteres;
2296   delete hEtaMonteres;
2297   delete hEtMonteres;
2298   delete hEtMonteratio;
2299
2300   for(Int_t i=0;i<nclasses;i++) delete hJetLeadingPt[i];
2301   delete[] hJetLeadingPt;
2302   for(Int_t i=0;i<nclasses;i++) delete hJetFragLeadingPt[i];
2303   delete[] hJetFragLeadingPt;
2304   for(Int_t i=0;i<nclasses;i++) delete hJetLeadingPtDist[i];
2305   delete[] hJetLeadingPtDist;
2306   for(Int_t i=0;i<nclasses;i++) delete hJetFragL[i];
2307   delete[] hJetFragL;
2308   for(Int_t i=0;i<nclasses;i++) delete hJetFragPL[i];
2309   delete[] hJetFragPL;
2310   for(Int_t i=0;i<nclasses;i++) delete hJetFragT[i];
2311   delete[] hJetFragT;
2312   for(Int_t i=0;i<nclasses;i++) delete hJetFragPt[i];
2313   delete[] hJetFragPt;
2314   for(Int_t i=0;i<nclasses;i++) delete hJetN[i];
2315   delete[] hJetN;
2316   for(Int_t i=0;i<nclasses;i++) delete hJetMeanPt[i];
2317   delete[] hJetMeanPt;
2318   for(Int_t i=0;i<nclasses;i++) delete hPhiCorr[i];
2319   delete[] hPhiCorr;
2320
2321   for(Int_t k=0;k<3;k++){
2322     for(Int_t i=0;i<nclasses;i++){
2323       delete hgJetFragLeadingPt[k][i];
2324       delete hgJetFragL[k][i];
2325       delete hgJetFragPL[k][i];
2326       delete hgJetFragT[k][i];
2327       delete hgJetFragPt[k][i];
2328     }
2329     delete[] hgJetFragLeadingPt[k];
2330     delete[] hgJetFragL[k];
2331     delete[] hgJetFragPL[k];
2332     delete[] hgJetFragT[k];
2333     delete[] hgJetFragPt[k];
2334   }
2335   delete[] hgJetFragLeadingPt;
2336   delete[] hgJetFragL;
2337   delete[] hgJetFragPL;
2338   delete[] hgJetFragT;
2339   delete[] hgJetFragPt;
2340
2341   for(Int_t k=0;k<3;k++){
2342     for(Int_t i=0;i<nclasses;i++){
2343       delete hgDiJetFragLeadingPt[k][i];
2344       delete hgDiJetFragL[k][i];
2345       delete hgDiJetFragPL[k][i];
2346       delete hgDiJetFragT[k][i];
2347       delete hgDiJetFragPt[k][i];
2348     }
2349     delete[] hgDiJetFragLeadingPt[k];
2350     delete[] hgDiJetFragL[k];
2351     delete[] hgDiJetFragPL[k];
2352     delete[] hgDiJetFragT[k];
2353     delete[] hgDiJetFragPt[k];
2354   }
2355   delete[] hgDiJetFragLeadingPt;
2356   delete[] hgDiJetFragL;
2357   delete[] hgDiJetFragPL;
2358   delete[] hgDiJetFragT;
2359   delete[] hgDiJetFragPt;
2360
2361   for(Int_t i=0;i<9;i++){
2362     delete hJetEttrigger[i];
2363     delete hJetEttrigger2[i];
2364   }
2365   delete[] hJetEttrigger;
2366   delete[] hJetEttrigger2;
2367   delete hJetEttriggernorm;
2368   delete rootfile;
2369 }
2370
2371 //------------------------------------------------------------------
2372
2373 Float_t relphi(Float_t phi1, Float_t phi2)
2374 { //rel to phi1
2375   Float_t ret=TMath::Abs(phi1-phi2);
2376   if(ret>TMath::Pi()) ret=TMath::TwoPi()-ret;
2377   return ret;
2378 }
2379
2380 Float_t addphi(Float_t phi1, Float_t phi2)
2381 {
2382   Float_t addphi=phi1+phi2;
2383   if(addphi>TMath::TwoPi()) addphi-=TMath::TwoPi();
2384   else if(addphi<0) addphi+=TMath::TwoPi();
2385   return addphi;
2386 }
2387
2388 Float_t diffphi(Float_t phi1, Float_t phi2)
2389 { //and move correlation to pi/2
2390   Float_t diffphi=TMath::Pi()/2+phi1-phi2;
2391   if(diffphi>TMath::TwoPi()) diffphi-=TMath::TwoPi();
2392   else if(diffphi<0) diffphi+=TMath::TwoPi();
2393   return diffphi;
2394 }
2395
2396 Int_t eventindex(Float_t phi1, Float_t phi2)
2397 { //rel to phi1
2398   Int_t ret=0; //toward (300 - 60)
2399   Float_t dphi=addphi(phi1,-phi2);
2400   const Float_t slice=TMath::Pi()/3.;
2401   if (dphi>1*slice && dphi < 2*slice)
2402     ret=2; //transverse (60-120)
2403   else if (dphi>4*slice && dphi < 5*slice)
2404     ret=3; //transverse (240-300)
2405   else if(dphi>=2*slice && dphi<=4*slice)
2406     ret=1; //away side (120-200)
2407
2408   return ret;
2409 }
2410
2411 void convert(Float_t pjet[4], Float_t &pt, Float_t &theta, Float_t &eta, Float_t &phi)
2412 {
2413   pt=TMath::Sqrt(pjet[0]*pjet[0]+pjet[1]*pjet[1]);
2414   theta=TMath::ATan2(pt,pjet[2]);
2415   eta=-TMath::Log(TMath::Tan(theta/2));
2416   phi=TMath::Pi()+TMath::ATan2(-pjet[1],-pjet[0]);
2417 }
2418