RICH/AliRICHv3.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16
  17 #include <Riostream.h>
  18
  19 #include <TBRIK.h>
  20 #include <TGeometry.h>
  21 #include <TLorentzVector.h>
  22 #include <TNode.h>
  23 #include <TParticle.h>
  24 #include <TVector3.h>
  25 #include <TVirtualMC.h>
  26 #include <TPDGCode.h> //for kNuetron
  27 #include <TCanvas.h>
  28 #include <TF1.h>
  29 #include <TH1.h>
  30 #include <TH2.h>
  31 #include <TStyle.h>
  32
  33 #include "AliConst.h"
  34 #include "AliMagF.h"
  35 #include "AliPDG.h"
  36 #include "AliRICHGeometry.h"
  37 #include "AliRICHResponseV0.h"
  38 #include "AliRICHSegmentationV1.h"
  39 #include "AliRICHv3.h"
  40 #include "AliRun.h"
  41 #include "AliRICHRecHit3D.h"
  42 #include "AliRICHRawCluster.h"
  43 #include "AliRICHDigit.h"
  44 #include "AliRICHRecHit1D.h"
  45
  46
  47 ClassImp(AliRICHv3)
  48
  49 //______________________________________________________________
  50 //    Implementation of the RICH version 3 with azimuthal rotation
  51
  52
  53 AliRICHv3::AliRICHv3(const char *sName, const char *sTitle)
  54           :AliRICH(sName,sTitle)
  55 {
  56 // The named ctor currently creates a single copy of
  57 // AliRICHGeometry AliRICHSegmentationV1 AliRICHResponseV0
  58 // and initialises the corresponding models of all 7 chambers with these stuctures.
  59 // Note: all chambers share the single copy of models. MUST be changed later (???).
  60   if(GetDebug())Info("named ctor","Start.");
  61
  62    fCkovNumber=fFreonProd=0;
  63
  64    AliRICHGeometry       *pRICHGeometry    =new AliRICHGeometry;           // ??? to be moved to AlRICHChamber::named ctor
  65    AliRICHSegmentationV1 *pRICHSegmentation=new AliRICHSegmentationV1;     // ??? to be moved to AlRICHChamber::named ctor
  66    AliRICHResponseV0     *pRICHResponse    =new AliRICHResponseV0;         // ??? to be moved to AlRICHChamber::named ctor
  67
  68    for (Int_t i=1; i<=kNCH; i++){
  69       SetGeometryModel(i,pRICHGeometry);
  70       SetSegmentationModel(i,pRICHSegmentation);
  71       SetResponseModel(i,pRICHResponse);
  72       C(i)->Init(i); // ??? to be removed
  73    }
  74   if(GetDebug())Info("named ctor","Stop.");
  75 }//AliRICHv3::ctor(const char *pcName, const char *pcTitle)
  76
  77 AliRICHv3::~AliRICHv3()
  78 {
  79 // Dtor deletes RICH models. In future (???) AliRICHChamber will be responsible for that.
  80    if(GetDebug()) cout<<ClassName()<<"::dtor()>\n";
  81
  82    if(fChambers) {
  83      AliRICHChamber *ch =C(1);
  84      if(ch) {
  85        delete ch->GetGeometryModel();
  86        delete ch->GetResponseModel();
  87        delete ch->GetSegmentationModel();
  88      }
  89    }
  90 }//AliRICHv3::dtor()
  91
  92
  93 void AliRICHv3::CreateGeometry()
  94 {
  95 // Provides geometry structure for simulation (currently GEANT volumes tree)
  96    if(GetDebug()) cout<<ClassName()<<"::CreateGeometry()>\n";
  97
  98   AliRICH *pRICH = (AliRICH *) gAlice->GetDetector("RICH");
  99   AliRICHSegmentationV0*  segmentation;
 100   AliRICHGeometry*  geometry;
 101   AliRICHChamber*       iChamber;
 102
 103   iChamber = &(pRICH->Chamber(0));
 104   segmentation=(AliRICHSegmentationV0*) iChamber->GetSegmentationModel();
 105   geometry=iChamber->GetGeometryModel();
 106
 107   Float_t distance;
 108   distance = geometry->GetFreonThickness()/2 + geometry->GetQuartzThickness() + geometry->GetGapThickness();
 109   geometry->SetRadiatorToPads(distance);
 110
 111   //Opaque quartz thickness
 112   Float_t oqua_thickness = .5;
 113   //CsI dimensions
 114
 115
 116   Float_t csi_width = segmentation->Npx()*segmentation->Dpx() + segmentation->DeadZone();
 117   Float_t csi_length = segmentation->Npy()*segmentation->Dpy() + 2*segmentation->DeadZone();
 118
 119
 120   Int_t *idtmed = fIdtmed->GetArray()-999;
 121
 122     Int_t i;
 123     Float_t zs;
 124     Int_t idrotm[1099];
 125     Float_t par[3];
 126
 127     // --- Define the RICH detector
 128     //     External aluminium box
 129     par[0] = 68.8;
 130     par[1] = 13;                 //Original Settings
 131     par[2] = 70.86;
 132     gMC->Gsvolu("RICH", "BOX ", idtmed[1009], par, 3);
 133
 134     //     Air
 135     par[0] = 66.3;
 136     par[1] = 13;                 //Original Settings
 137     par[2] = 68.35;
 138     gMC->Gsvolu("SRIC", "BOX ", idtmed[1000], par, 3);
 139
 140     //    Air 2 (cutting the lower part of the box)
 141
 142     par[0] = 1.25;
 143     par[1] = 3;                 //Original Settings
 144     par[2] = 70.86;
 145     gMC->Gsvolu("AIR2", "BOX ", idtmed[1000], par, 3);
 146
 147     //    Air 3 (cutting the lower part of the box)
 148
 149     par[0] = 66.3;
 150     par[1] = 3;                 //Original Settings
 151     par[2] = 1.2505;
 152     gMC->Gsvolu("AIR3", "BOX ", idtmed[1000], par, 3);
 153
 154     //     Honeycomb
 155     par[0] = 66.3;
 156     par[1] = .188;                 //Original Settings
 157     par[2] = 68.35;
 158     gMC->Gsvolu("HONE", "BOX ", idtmed[1001], par, 3);
 159
 160     //     Aluminium sheet
 161     par[0] = 66.3;
 162     par[1] = .025;                 //Original Settings
 163     par[2] = 68.35;
 164     /*par[0] = 66.5;
 165     par[1] = .025;
 166     par[2] = 63.1;*/
 167     gMC->Gsvolu("ALUM", "BOX ", idtmed[1009], par, 3);
 168
 169     //     Quartz
 170     par[0] = geometry->GetQuartzWidth()/2;
 171     par[1] = geometry->GetQuartzThickness()/2;
 172     par[2] = geometry->GetQuartzLength()/2;
 173     gMC->Gsvolu("QUAR", "BOX ", idtmed[1002], par, 3);
 174
 175     //     Spacers (cylinders)
 176     par[0] = 0.;
 177     par[1] = .5;
 178     par[2] = geometry->GetFreonThickness()/2;
 179     gMC->Gsvolu("SPAC", "TUBE", idtmed[1002], par, 3);
 180
 181     //     Feet (freon slabs supports)
 182
 183     par[0] = .7;
 184     par[1] = .3;
 185     par[2] = 1.9;
 186     gMC->Gsvolu("FOOT", "BOX", idtmed[1009], par, 3);
 187
 188     //     Opaque quartz
 189     par[0] = geometry->GetQuartzWidth()/2;
 190     par[1] = .2;
 191     par[2] = geometry->GetQuartzLength()/2;
 192     gMC->Gsvolu("OQUA", "BOX ", idtmed[1007], par, 3);
 193
 194     //     Frame of opaque quartz
 195     par[0] = geometry->GetOuterFreonWidth()/2;
 196     par[1] = geometry->GetFreonThickness()/2;
 197     par[2] = geometry->GetOuterFreonLength()/2;
 198     gMC->Gsvolu("OQF1", "BOX ", idtmed[1007], par, 3);
 199
 200     par[0] = geometry->GetInnerFreonWidth()/2;
 201     par[1] = geometry->GetFreonThickness()/2;
 202     par[2] = geometry->GetInnerFreonLength()/2;
 203     gMC->Gsvolu("OQF2", "BOX ", idtmed[1007], par, 3);
 204
 205
 206     //     Freon
 207     par[0] = geometry->GetOuterFreonWidth()/2 - oqua_thickness;
 208     par[1] = geometry->GetFreonThickness()/2;
 209     par[2] = geometry->GetOuterFreonLength()/2 - 2*oqua_thickness;
 210     gMC->Gsvolu("FRE1", "BOX ", idtmed[1003], par, 3);
 211
 212     par[0] = geometry->GetInnerFreonWidth()/2 - oqua_thickness;
 213     par[1] = geometry->GetFreonThickness()/2;
 214     par[2] = geometry->GetInnerFreonLength()/2 - 2*oqua_thickness;
 215     gMC->Gsvolu("FRE2", "BOX ", idtmed[1003], par, 3);
 216
 217     //     Methane
 218     par[0] = csi_width/2;
 219     par[1] = geometry->GetGapThickness()/2;
 220     par[2] = csi_length/2;
 221     gMC->Gsvolu("META", "BOX ", idtmed[1004], par, 3);
 222
 223     //     Methane gap
 224     par[0] = csi_width/2;
 225     par[1] = geometry->GetProximityGapThickness()/2;
 226     par[2] = csi_length/2;
 227     gMC->Gsvolu("GAP ", "BOX ", idtmed[1008], par, 3);
 228
 229     //     CsI photocathode
 230     par[0] = csi_width/2;
 231     par[1] = .25;
 232     par[2] = csi_length/2;
 233     gMC->Gsvolu("CSI ", "BOX ", idtmed[1005], par, 3);
 234
 235     //     Anode grid
 236     par[0] = 0.;
 237     par[1] = .001;
 238     par[2] = 20.;
 239     gMC->Gsvolu("GRID", "TUBE", idtmed[1006], par, 3);
 240
 241     // Wire supports
 242     // Bar of metal
 243
 244     par[0] = csi_width/2;
 245     par[1] = 1.05;
 246     par[2] = 1.05;
 247     gMC->Gsvolu("WSMe", "BOX ", idtmed[1009], par, 3);
 248
 249     // Ceramic pick up (base)
 250
 251     par[0] =  csi_width/2;
 252     par[1] = .25;
 253     par[2] = 1.05;
 254     gMC->Gsvolu("WSG1", "BOX ", idtmed[1010], par, 3);
 255
 256     // Ceramic pick up (head)
 257
 258     par[0] = csi_width/2;
 259     par[1] = .1;
 260     par[2] = .1;
 261     gMC->Gsvolu("WSG2", "BOX ", idtmed[1010], par, 3);
 262
 263     // Aluminium supports for methane and CsI
 264     // Short bar
 265
 266     par[0] = csi_width/2;
 267     par[1] = geometry->GetGapThickness()/2 + .25;
 268     par[2] = (68.35 - csi_length/2)/2;
 269     gMC->Gsvolu("SMSH", "BOX", idtmed[1009], par, 3);
 270
 271     // Long bar
 272
 273     par[0] = (66.3 - csi_width/2)/2;
 274     par[1] = geometry->GetGapThickness()/2 + .25;
 275     par[2] = csi_length/2 + 68.35 - csi_length/2;
 276     gMC->Gsvolu("SMLG", "BOX", idtmed[1009], par, 3);
 277
 278     // Aluminium supports for freon
 279     // Short bar
 280
 281     par[0] = geometry->GetQuartzWidth()/2;
 282     par[1] = .3;
 283     par[2] = (68.35 - geometry->GetQuartzLength()/2)/2;
 284     gMC->Gsvolu("SFSH", "BOX", idtmed[1009], par, 3);
 285
 286     // Long bar
 287
 288     par[0] = (66.3 - geometry->GetQuartzWidth()/2)/2;
 289     par[1] = .3;
 290     par[2] = geometry->GetQuartzLength()/2 + 68.35 - geometry->GetQuartzLength()/2;
 291     gMC->Gsvolu("SFLG", "BOX", idtmed[1009], par, 3);
 292
 293     // PCB backplane
 294
 295     par[0] = csi_width/2;
 296     par[1] = .25;
 297     par[2] = csi_length/4 -.5025;
 298     gMC->Gsvolu("PCB ", "BOX", idtmed[1011], par, 3);
 299
 300
 301     // Backplane supports
 302
 303     // Aluminium slab
 304
 305     par[0] = 33.15;
 306     par[1] = 2;
 307     par[2] = 21.65;
 308     gMC->Gsvolu("BACK", "BOX", idtmed[1009], par, 3);
 309
 310     // Big hole
 311
 312     par[0] = 9.05;
 313     par[1] = 2;
 314     par[2] = 4.4625;
 315     gMC->Gsvolu("BKHL", "BOX", idtmed[1000], par, 3);
 316
 317     // Small hole
 318
 319     par[0] = 5.7;
 320     par[1] = 2;
 321     par[2] = 4.4625;
 322     gMC->Gsvolu("BKHS", "BOX", idtmed[1000], par, 3);
 323
 324     // Place holes inside backplane support
 325
 326     gMC->Gspos("BKHS", 1, "BACK", .8 + 5.7,0., .6 + 4.4625, 0, "ONLY");
 327     gMC->Gspos("BKHS", 2, "BACK", -.8 - 5.7,0., .6 + 4.4625, 0, "ONLY");
 328     gMC->Gspos("BKHS", 3, "BACK", .8 + 5.7,0., -.6 - 4.4625, 0, "ONLY");
 329     gMC->Gspos("BKHS", 4, "BACK", -.8 - 5.7,0., -.6 - 4.4625, 0, "ONLY");
 330     gMC->Gspos("BKHS", 5, "BACK", .8 + 5.7,0., .6 + 8.925 + 1.2 + 4.4625, 0, "ONLY");
 331     gMC->Gspos("BKHS", 6, "BACK", -.8 - 5.7,0., .6 + 8.925 + 1.2 + 4.4625, 0, "ONLY");
 332     gMC->Gspos("BKHS", 7, "BACK", .8 + 5.7,0., -.6 - 8.925 - 1.2 - 4.4625, 0, "ONLY");
 333     gMC->Gspos("BKHS", 8, "BACK", -.8 - 5.7,0., -.6 - 8.925 - 1.2 - 4.4625, 0, "ONLY");
 334     gMC->Gspos("BKHL", 1, "BACK", .8 + 11.4 + 1.6 + 9.05, 0., .6 + 4.4625, 0, "ONLY");
 335     gMC->Gspos("BKHL", 2, "BACK", -.8 - 11.4 - 1.6 - 9.05, 0., .6 + 4.4625, 0, "ONLY");
 336     gMC->Gspos("BKHL", 3, "BACK", .8 + 11.4 + 1.6 + 9.05, 0., -.6 - 4.4625, 0, "ONLY");
 337     gMC->Gspos("BKHL", 4, "BACK", -.8 - 11.4 - 1.6 - 9.05, 0., -.6 - 4.4625, 0, "ONLY");
 338     gMC->Gspos("BKHL", 5, "BACK", .8 + 11.4+ 1.6 + 9.05, 0., .6 + 8.925 + 1.2 + 4.4625, 0, "ONLY");
 339     gMC->Gspos("BKHL", 6, "BACK", -.8 - 11.4 - 1.6 - 9.05, 0., .6 + 8.925 + 1.2 + 4.4625, 0, "ONLY");
 340     gMC->Gspos("BKHL", 7, "BACK", .8 + 11.4 + 1.6 + 9.05, 0., -.6 - 8.925 - 1.2 - 4.4625, 0, "ONLY");
 341     gMC->Gspos("BKHL", 8, "BACK", -.8 - 11.4 - 1.6 - 9.05, 0., -.6 - 8.925 - 1.2 - 4.4625, 0, "ONLY");
 342
 343
 344
 345     // --- Places the detectors defined with GSVOLU
 346     //     Place material inside RICH
 347     gMC->Gspos("SRIC", 1, "RICH", 0.,0., 0., 0, "ONLY");
 348     gMC->Gspos("AIR2", 1, "RICH", 66.3 + 1.2505, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.35, 0., 0, "ONLY");
 349     gMC->Gspos("AIR2", 2, "RICH", -66.3 - 1.2505, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.35, 0., 0, "ONLY");
 350     gMC->Gspos("AIR3", 1, "RICH", 0.,  1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.35, -68.35 - 1.25, 0, "ONLY");
 351     gMC->Gspos("AIR3", 2, "RICH", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.35,  68.35 + 1.25, 0, "ONLY");
 352
 353
 354     gMC->Gspos("ALUM", 1, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.025, 0., 0, "ONLY");
 355     gMC->Gspos("HONE", 1, "SRIC", 0., 1.276- geometry->GetGapThickness()/2  - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .188, 0., 0, "ONLY");
 356     gMC->Gspos("ALUM", 2, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .025, 0., 0, "ONLY");
 357     gMC->Gspos("FOOT", 1, "SRIC", 64.95, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, 36.9, 0, "ONLY");
 358     gMC->Gspos("FOOT", 2, "SRIC", 21.65, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3 , 36.9, 0, "ONLY");
 359     gMC->Gspos("FOOT", 3, "SRIC", -21.65, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, 36.9, 0, "ONLY");
 360     gMC->Gspos("FOOT", 4, "SRIC", -64.95, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, 36.9, 0, "ONLY");
 361     gMC->Gspos("FOOT", 5, "SRIC", 64.95, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, -36.9, 0, "ONLY");
 362     gMC->Gspos("FOOT", 6, "SRIC", 21.65, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, -36.9, 0, "ONLY");
 363     gMC->Gspos("FOOT", 7, "SRIC", -21.65, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, -36.9, 0, "ONLY");
 364     gMC->Gspos("FOOT", 8, "SRIC", -64.95, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, -36.9, 0, "ONLY");
 365     gMC->Gspos("OQUA", 1, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .2, 0., 0, "ONLY");
 366
 367     // Supports placing
 368
 369     // Methane supports
 370     gMC->Gspos("SMLG", 1, "SRIC", csi_width/2 + (66.3 - csi_width/2)/2, 1.276 + .25, 0., 0, "ONLY");
 371     gMC->Gspos("SMLG", 2, "SRIC", - csi_width/2 - (66.3 - csi_width/2)/2, 1.276 + .25, 0., 0, "ONLY");
 372     gMC->Gspos("SMSH", 1, "SRIC", 0., 1.276 + .25, csi_length/2 + (68.35 - csi_length/2)/2, 0, "ONLY");
 373     gMC->Gspos("SMSH", 2, "SRIC", 0., 1.276 + .25, - csi_length/2 - (68.35 - csi_length/2)/2, 0, "ONLY");
 374
 375     //Freon supports
 376
 377     Float_t supp_y = 1.276 - geometry->GetGapThickness()/2- geometry->GetQuartzThickness() -geometry->GetFreonThickness() - .2 + .3; //y position of freon supports
 378
 379     gMC->Gspos("SFLG", 1, "SRIC", geometry->GetQuartzWidth()/2 + (66.3 - geometry->GetQuartzWidth()/2)/2, supp_y, 0., 0, "ONLY");
 380     gMC->Gspos("SFLG", 2, "SRIC", - geometry->GetQuartzWidth()/2 - (66.3 - geometry->GetQuartzWidth()/2)/2, supp_y, 0., 0, "ONLY");
 381     gMC->Gspos("SFSH", 1, "SRIC", 0., supp_y, geometry->GetQuartzLength()/2 + (68.35 - geometry->GetQuartzLength()/2)/2, 0, "ONLY");
 382     gMC->Gspos("SFSH", 2, "SRIC", 0., supp_y, - geometry->GetQuartzLength()/2 - (68.35 - geometry->GetQuartzLength()/2)/2, 0, "ONLY");
 383
 384     AliMatrix(idrotm[1019], 0., 0., 90., 0., 90., 90.);
 385
 386
 387     Int_t nspacers = 30;
 388
 389     for (i = 0; i < nspacers/3; i++) {
 390         zs = -11.6/2 + (TMath::Abs(nspacers/6) - i) * 12.2;
 391         gMC->Gspos("SPAC", i, "FRE1", 10.5, 0., zs, idrotm[1019], "ONLY");  //Original settings
 392     }
 393
 394     for (i = nspacers/3; i < (nspacers*2)/3; i++) {
 395         zs = -11.6/2 + (nspacers/3 + TMath::Abs(nspacers/6) - i) * 12.2;
 396         gMC->Gspos("SPAC", i, "FRE1", 0, 0., zs, idrotm[1019], "ONLY");  //Original settings
 397     }
 398
 399     for (i = (nspacers*2)/3; i < nspacers; ++i) {
 400         zs = -11.6/2 + ((nspacers*2)/3 + TMath::Abs(nspacers/6) - i) * 12.2;
 401         gMC->Gspos("SPAC", i, "FRE1", -10.5, 0., zs, idrotm[1019], "ONLY"); //Original settings
 402     }
 403
 404     for (i = 0; i < nspacers/3; i++) {
 405         zs = -11.6/2 + (TMath::Abs(nspacers/6) - i) * 12.2;
 406         gMC->Gspos("SPAC", i, "FRE2", 10.5, 0., zs, idrotm[1019], "ONLY");  //Original settings
 407     }
 408
 409     for (i = nspacers/3; i < (nspacers*2)/3; i++) {
 410         zs = -11.6/2 + (nspacers/3 + TMath::Abs(nspacers/6) - i) * 12.2;
 411         gMC->Gspos("SPAC", i, "FRE2", 0, 0., zs, idrotm[1019], "ONLY");  //Original settings
 412     }
 413
 414     for (i = (nspacers*2)/3; i < nspacers; ++i) {
 415         zs = -11.6/2 + ((nspacers*2)/3 + TMath::Abs(nspacers/6) - i) * 12.2;
 416         gMC->Gspos("SPAC", i, "FRE2", -10.5, 0., zs, idrotm[1019], "ONLY"); //Original settings
 417     }
 418
 419
 420     gMC->Gspos("FRE1", 1, "OQF1", 0., 0., 0., 0, "ONLY");
 421     gMC->Gspos("FRE2", 1, "OQF2", 0., 0., 0., 0, "ONLY");
 422     gMC->Gspos("OQF1", 1, "SRIC", geometry->GetOuterFreonWidth()/2 + geometry->GetInnerFreonWidth()/2 + 2, 1.276 - geometry->GetGapThickness()/2- geometry->GetQuartzThickness() -geometry->GetFreonThickness()/2, 0., 0, "ONLY"); //Original settings (31.3)
 423     gMC->Gspos("OQF2", 2, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()/2, 0., 0, "ONLY");          //Original settings
 424     gMC->Gspos("OQF1", 3, "SRIC", - (geometry->GetOuterFreonWidth()/2 + geometry->GetInnerFreonWidth()/2) - 2, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()/2, 0., 0, "ONLY");       //Original settings (-31.3)
 425     gMC->Gspos("QUAR", 1, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness()/2, 0., 0, "ONLY");
 426     gMC->Gspos("GAP ", 1, "META", 0., geometry->GetGapThickness()/2 - geometry->GetProximityGapThickness()/2 - 0.0001, 0., 0, "ONLY");
 427     gMC->Gspos("META", 1, "SRIC", 0., 1.276, 0., 0, "ONLY");
 428     gMC->Gspos("CSI ", 1, "SRIC", 0., 1.276 + geometry->GetGapThickness()/2 + .25, 0., 0, "ONLY");
 429     printf("CSI pos: %f\n",1.276 + geometry->GetGapThickness()/2 + .25);
 430
 431     // Wire support placing
 432
 433     gMC->Gspos("WSG2", 1, "GAP ", 0., geometry->GetProximityGapThickness()/2 - .1, 0., 0, "ONLY");
 434     gMC->Gspos("WSG1", 1, "CSI ", 0., 0., 0., 0, "ONLY");
 435     gMC->Gspos("WSMe", 1, "SRIC ", 0., 1.276 + geometry->GetGapThickness()/2 + .5 + 1.05, 0., 0, "ONLY");
 436
 437     // Backplane placing
 438
 439     gMC->Gspos("BACK", 1, "SRIC ", -33.15, 1.276 + geometry->GetGapThickness()/2 + .5 + 2.1 + 2, 43.3, 0, "ONLY");
 440     gMC->Gspos("BACK", 2, "SRIC ", 33.15, 1.276 + geometry->GetGapThickness()/2 + .5 + 2.1 + 2 , 43.3, 0, "ONLY");
 441     gMC->Gspos("BACK", 3, "SRIC ", -33.15, 1.276 + geometry->GetGapThickness()/2 + .5 + 2.1 + 2, 0., 0, "ONLY");
 442     gMC->Gspos("BACK", 4, "SRIC ", 33.15, 1.276 + geometry->GetGapThickness()/2 + .5 + 2.1 + 2, 0., 0, "ONLY");
 443     gMC->Gspos("BACK", 5, "SRIC ", 33.15, 1.276 + geometry->GetGapThickness()/2 + .5 + 2.1 + 2, -43.3, 0, "ONLY");
 444     gMC->Gspos("BACK", 6, "SRIC ", -33.15, 1.276 + geometry->GetGapThickness()/2 + .5 + 2.1 + 2, -43.3, 0, "ONLY");
 445
 446     // PCB placing
 447
 448     gMC->Gspos("PCB ", 1, "SRIC ", 0.,  1.276 + geometry->GetGapThickness()/2 + .5 + 1.05, csi_width/4 + .5025 + 2.5, 0, "ONLY");
 449     gMC->Gspos("PCB ", 2, "SRIC ", 0.,  1.276 + geometry->GetGapThickness()/2 + .5 + 1.05, -csi_width/4 - .5025 - 2.5, 0, "ONLY");
 450
 451 // Place chambers into mother volume ALIC
 452
 453    Double_t dOffset        = geometry->GetOffset() - geometry->GetGapThickness()/2;  // distance from center of mother volume ALIC to methane
 454
 455    Double_t dAlpha         = geometry->GetAlphaAngle(); // angle between centers of chambers - y-z plane
 456    Double_t dAlphaRad      = dAlpha*kDegrad;
 457
 458    Double_t dBeta          = geometry->GetBetaAngle();   // angle between center of chambers - y-x plane
 459    Double_t dBetaRad       = dBeta*kDegrad;
 460
 461    Double_t dRotAngle      = geometry->GetRotationAngle();     // the whole RICH is to be rotated in x-y plane + means clockwise rotation
 462    Double_t dRotAngleRad   = dRotAngle*kDegrad;
 463
 464
 465    TRotMatrix *pRotMatrix; // tmp pointer
 466
 467    TVector3 vector(0,dOffset,0); // Position of chamber 2 without rotation
 468
 469 // Chamber 0  standalone (no other chambers in this row)
 470    pRotMatrix = new TRotMatrix("rot993","rot993", 0., 0., 0.,0.,0.,0.);
 471    const Double_t* r   = pRotMatrix->SetAngles(90., 0., 90.-dAlpha , 90.,  dAlpha, -90.);
 472    Double_t* rr  = RotateXY(r, -dRotAngleRad);
 473    AliMatrix(idrotm[1000], rr[0], rr[1], rr[2], rr[3], rr[4], rr[5]);
 474    pRotMatrix->SetAngles(rr[0], rr[1], rr[2], rr[3], rr[4], rr[5]);
 475
 476    vector.SetXYZ(0,dOffset,0);  vector.RotateX(dAlphaRad);
 477    vector.RotateZ(-dRotAngleRad);
 478
 479    gMC->Gspos("RICH",1,"ALIC",vector.X(),vector.Y(),vector.Z(),idrotm[1000], "ONLY");
 480    Chamber(0).SetChamberTransform(vector.X(),vector.Y(),vector.Z(),pRotMatrix);
 481   if(GetDebug()) Info("CreateGeometry 0","%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f",rr[0],rr[1],rr[2],rr[3],rr[4],rr[5]);
 482   if(GetDebug()) Info("CreateGeometry 0","x=%8.3f y=%8.3f z=%8.3f",vector.X(),vector.Y(),vector.Z());
 483 // Chamber 1
 484    pRotMatrix = new TRotMatrix("rot994","rot994", 0., 0., 0.,0.,0.,0.);
 485    r   = pRotMatrix->SetAngles(90., -dBeta, 90., 90.-dBeta,  0., 0.);
 486    rr  = RotateXY(r, -dRotAngleRad);
 487    AliMatrix(idrotm[1001], rr[0], rr[1], rr[2], rr[3], rr[4], rr[5]);
 488    pRotMatrix->SetAngles(rr[0], rr[1], rr[2], rr[3], rr[4], rr[5]);
 489    vector.SetXYZ(0,dOffset,0);  vector.RotateZ(-dBetaRad);
 490    vector.RotateZ(-dRotAngleRad);
 491
 492    gMC->Gspos("RICH",2,"ALIC",vector.X(),vector.Y(),vector.Z(),idrotm[1001], "ONLY");
 493    Chamber(1).SetChamberTransform(vector.X(),vector.Y(),vector.Z(),pRotMatrix);
 494   if(GetDebug()) Info("CreateGeometry 1","%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f",rr[0],rr[1],rr[2],rr[3],rr[4],rr[5]);
 495   if(GetDebug()) Info("CreateGeometry 1","x=%8.3f y=%8.3f z=%8.3f",vector.X(),vector.Y(),vector.Z());
 496 // Chamber 2   the top one with no Alpha-Beta rotation
 497    pRotMatrix = new TRotMatrix("rot995","rot995", 0., 0., 0.,0.,0.,0.);
 498    r   = pRotMatrix->SetAngles(90., 0., 90., 90.,  0., 0.);
 499    rr  = RotateXY(r, -dRotAngleRad);
 500    AliMatrix(idrotm[1002], rr[0], rr[1], rr[2], rr[3], rr[4], rr[5]);
 501    pRotMatrix->SetAngles(rr[0], rr[1], rr[2], rr[3], rr[4], rr[5]);
 502    vector.SetXYZ(0,dOffset,0);
 503    vector.RotateZ(-dRotAngleRad);
 504    gMC->Gspos("RICH",3,"ALIC",vector.X(),vector.Y(),vector.Z(),idrotm[1002], "ONLY");
 505    Chamber(2).SetChamberTransform(vector.X(),vector.Y(),vector.Z(),pRotMatrix);
 506   if(GetDebug()) Info("CreateGeometry 2","%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f",rr[0],rr[1],rr[2],rr[3],rr[4],rr[5]);
 507   if(GetDebug()) Info("CreateGeometry 2","x=%8.3f y=%8.3f z=%8.3f",vector.X(),vector.Y(),vector.Z());
 508 // Chamber 3
 509    pRotMatrix = new TRotMatrix("rot996","rot996", 0., 0., 0.,0.,0.,0.);
 510    r   = pRotMatrix->SetAngles(90., dBeta, 90., 90.+dBeta,  0., 0.);
 511    rr  = RotateXY(r, -dRotAngleRad);
 512    AliMatrix(idrotm[1003], rr[0], rr[1], rr[2], rr[3], rr[4], rr[5]);
 513    pRotMatrix->SetAngles(rr[0], rr[1], rr[2], rr[3], rr[4], rr[5]);
 514    vector.SetXYZ(0,dOffset,0);  vector.RotateZ(dBetaRad);
 515    vector.RotateZ(-dRotAngleRad);
 516
 517    gMC->Gspos("RICH",4,"ALIC",vector.X(),vector.Y(),vector.Z(),idrotm[1003], "ONLY");
 518    Chamber(3).SetChamberTransform(vector.X(),vector.Y(),vector.Z(),pRotMatrix);
 519   if(GetDebug()) Info("CreateGeometry 3","%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f",rr[0],rr[1],rr[2],rr[3],rr[4],rr[5]);
 520   if(GetDebug()) Info("CreateGeometry 3","x=%8.3f y=%8.3f z=%8.3f",vector.X(),vector.Y(),vector.Z());
 521 // Chamber 4
 522    pRotMatrix = new TRotMatrix("rot997","rot997", 0., 0., 0.,0.,0.,0.);
 523    r   = pRotMatrix->SetAngles(90., 360.-dBeta, 108.2, 90.-dBeta,  18.2, 90.-dBeta);
 524    rr  = RotateXY(r, -dRotAngleRad);
 525    AliMatrix(idrotm[1004], rr[0], rr[1], rr[2], rr[3], rr[4], rr[5]);
 526    pRotMatrix->SetAngles(rr[0], rr[1], rr[2], rr[3], rr[4], rr[5]);
 527    vector.SetXYZ(0,dOffset,0);  vector.RotateZ(-dBetaRad); vector.RotateX(-dAlphaRad);
 528    vector.RotateZ(-dRotAngleRad);
 529
 530    gMC->Gspos("RICH",5,"ALIC",vector.X(),vector.Y(),vector.Z(),idrotm[1004], "ONLY");
 531    Chamber(4).SetChamberTransform(vector.X(),vector.Y(),vector.Z(),pRotMatrix);
 532   if(GetDebug()) Info("CreateGeometry 4","%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f",rr[0],rr[1],rr[2],rr[3],rr[4],rr[5]);
 533   if(GetDebug()) Info("CreateGeometry 4","x=%8.3f y=%8.3f z=%8.3f",vector.X(),vector.Y(),vector.Z());
 534 // Chamber 5
 535    pRotMatrix = new TRotMatrix("rot998","rot998", 0., 0., 0.,0.,0.,0.);
 536    r   = pRotMatrix->SetAngles(90., 0., 90.+dAlpha, 90.,  dAlpha, 90.);
 537    rr  = RotateXY(r, -dRotAngleRad);
 538    AliMatrix(idrotm[1005], rr[0], rr[1], rr[2], rr[3], rr[4], rr[5]);
 539    pRotMatrix->SetAngles(rr[0], rr[1], rr[2], rr[3], rr[4], rr[5]);
 540    vector.SetXYZ(0,dOffset,0); vector.RotateX(-dAlphaRad);
 541    vector.RotateZ(-dRotAngleRad);
 542
 543    gMC->Gspos("RICH",6,"ALIC",vector.X(),vector.Y(),vector.Z(),idrotm[1005], "ONLY");
 544    Chamber(5).SetChamberTransform(vector.X(),vector.Y(),vector.Z(),pRotMatrix);
 545   if(GetDebug()) Info("CreateGeometry 5","%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f",rr[0],rr[1],rr[2],rr[3],rr[4],rr[5]);
 546   if(GetDebug()) Info("CreateGeometry 5","x=%8.3f y=%8.3f z=%8.3f",vector.X(),vector.Y(),vector.Z());
 547 // Chamber 6
 548    pRotMatrix = new TRotMatrix("rot999","rot999", 0., 0., 0.,0.,0.,0.);
 549    r   = pRotMatrix->SetAngles(90., dBeta, 108.2, 90.+dBeta,  18.2, 90.+dBeta);
 550    rr  = RotateXY(r, -dRotAngleRad);
 551    AliMatrix(idrotm[1006], rr[0], rr[1], rr[2], rr[3], rr[4], rr[5]);
 552    pRotMatrix->SetAngles(rr[0], rr[1], rr[2], rr[3], rr[4], rr[5]);
 553    vector.SetXYZ(0,dOffset,0);  vector.RotateZ(dBetaRad); vector.RotateX(-dAlphaRad);
 554    vector.RotateZ(-dRotAngleRad);
 555
 556    gMC->Gspos("RICH",7,"ALIC",vector.X(),vector.Y(),vector.Z(),idrotm[1006], "ONLY");
 557    Chamber(6).SetChamberTransform(vector.X(),vector.Y(),vector.Z(),pRotMatrix);
 558   if(GetDebug()) Info("CreateGeometry 6","%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f",rr[0],rr[1],rr[2],rr[3],rr[4],rr[5]);
 559   if(GetDebug()) Info("CreateGeometry 6","x=%8.3f y=%8.3f z=%8.3f",vector.X(),vector.Y(),vector.Z());
 560
 561 }//void AliRICHv3::CreateGeometry()
 562 //______________________________________________________________________________
 563 void AliRICHv3::Init()
 564 {//Makes nothing for a while
 565   if(GetDebug())Info("Init","Start.");
 566   if(GetDebug())Info("Init","Stop.");
 567 }
 568 //______________________________________________________________________________
 569 Double_t* AliRICHv3::RotateXY(const Double_t* r, Double_t a)
 570 {
 571     // Rotatation in xy-plane
 572     // by angle a
 573     // The resulting rotation matrix is given back in the G3 notation.
 574     Double_t* rr = new Double_t[6];
 575     Double_t m[9];
 576     Int_t i,j,k;
 577
 578     for (i = 0; i < 3; i++) {
 579         j = 3*i;
 580         m[j]   = r[j] * TMath::Cos(a) - r[j+1] * TMath::Sin(a);
 581         m[j+1] = r[j] * TMath::Sin(a) + r[j+1] * TMath::Cos(a);
 582         m[j+2] = r[j+2];
 583     }
 584
 585     for (i = 0; i < 3; i++) {
 586             j = 3*i;
 587             k = 2*i;
 588             rr[k]    = TMath::ACos(m[j+2])        * kRaddeg;
 589             rr[k+1]  = TMath::ATan2(m[j+1], m[j]) * kRaddeg;
 590     }
 591     return rr;
 592 }//Double_t* AliRICHv3::RotateXY(const Double_t* r, Double_t a)
 593 //______________________________________________________________________________
 594 void AliRICHv3::StepManager()
 595 {//Full Step Manager
 596
 597     Int_t          copy, id;
 598     static Int_t   idvol;
 599     static Int_t   vol[2];
 600     Int_t          ipart;
 601     static Float_t hits[22];
 602     static Float_t ckovData[19];
 603     TLorentzVector position;
 604     TLorentzVector momentum;
 605     Float_t        pos[3];
 606     Float_t        mom[4];
 607     Float_t        localPos[3];
 608     Float_t        localMom[4];
 609     Float_t        localTheta,localPhi;
 610     Float_t        theta,phi;
 611     Float_t        destep, step;
 612     Double_t        ranf[2];
 613     Float_t        coscerenkov;
 614     static Float_t eloss, xhit, yhit, tlength;
 615     const  Float_t kBig=1.e10;
 616
 617     TClonesArray &lhits = *fHits;
 618     TParticle *current = (TParticle*)(*gAlice->Particles())[gAlice->GetCurrentTrackNumber()];
 619
 620  //if (current->Energy()>1)
 621    //{
 622
 623     // Only gas gap inside chamber
 624     // Tag chambers and record hits when track enters
 625
 626
 627     id=gMC->CurrentVolID(copy);
 628     idvol = copy-1;
 629     Float_t cherenkovLoss=0;
 630     //gAlice->KeepTrack(gAlice->GetCurrentTrackNumber());
 631
 632     gMC->TrackPosition(position);
 633     pos[0]=position(0);
 634     pos[1]=position(1);
 635     pos[2]=position(2);
 636     //bzero((char *)ckovData,sizeof(ckovData)*19);
 637     ckovData[1] = pos[0];                 // X-position for hit
 638     ckovData[2] = pos[1];                 // Y-position for hit
 639     ckovData[3] = pos[2];                 // Z-position for hit
 640     ckovData[6] = 0;                      // dummy track length
 641     //ckovData[11] = gAlice->GetCurrentTrackNumber();
 642
 643     //printf("\n+++++++++++\nTrack: %d\n++++++++++++\n",gAlice->GetCurrentTrackNumber());
 644
 645     //AliRICH *RICH = (AliRICH *) gAlice->GetDetector("RICH");
 646
 647     /********************Store production parameters for Cerenkov photons************************/
 648 //is it a Cerenkov photon?
 649     if (gMC->TrackPid() == 50000050) {
 650
 651       //if (gMC->VolId("GAP ")==gMC->CurrentVolID(copy))
 652         //{
 653           Float_t ckovEnergy = current->Energy();
 654           //energy interval for tracking
 655           if  (ckovEnergy > 5.6e-09 && ckovEnergy < 7.8e-09 )
 656             //if (ckovEnergy > 0)
 657             {
 658               if (gMC->IsTrackEntering()){        //is track entering?
 659                 //printf("Track entered (1)\n");
 660                 if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
 661                   {                                                          //is it in freo?
 662                     if (gMC->IsNewTrack()){                          //is it the first step?
 663                       //printf("I'm in!\n");
 664                       Int_t mother = current->GetFirstMother();
 665
 666                       //printf("Second Mother:%d\n",current->GetSecondMother());
 667
 668                       ckovData[10] = mother;
 669                       ckovData[11] = gAlice->GetCurrentTrackNumber();
 670                       ckovData[12] = 1;             //Media where photon was produced 1->Freon, 2->Quarz
 671                       //printf("Produced in FREO\n");
 672                       fCkovNumber++;
 673                       fFreonProd=1;
 674                       //printf("Index: %d\n",fCkovNumber);
 675                     }    //first step question
 676                   }        //freo question
 677
 678                 if (gMC->IsNewTrack()){                                  //is it first step?
 679                   if (gMC->VolId("QUAR")==gMC->CurrentVolID(copy))             //is it in quarz?
 680                     {
 681                       ckovData[12] = 2;
 682                       //printf("Produced in QUAR\n");
 683                     }    //quarz question
 684                 }        //first step question
 685
 686                 //printf("Before %d\n",fFreonProd);
 687               }   //track entering question
 688
 689               if (ckovData[12] == 1)                                        //was it produced in Freon?
 690                 //if (fFreonProd == 1)
 691                 {
 692                   if (gMC->IsTrackEntering()){                                     //is track entering?
 693                     //printf("Track entered (2)\n");
 694                     //printf("Current volume (should be META): %s\n",gMC->CurrentVolName());
 695                     //printf("VolId: %d, CurrentVolID: %d\n",gMC->VolId("META"),gMC->CurrentVolID(copy));
 696                     if (gMC->VolId("META")==gMC->CurrentVolID(copy))                //is it in gap?
 697                       {
 698                         //printf("Got in META\n");
 699                         gMC->TrackMomentum(momentum);
 700                         mom[0]=momentum(0);
 701                         mom[1]=momentum(1);
 702                         mom[2]=momentum(2);
 703                         mom[3]=momentum(3);
 704
 705                         gMC->Gmtod(mom,localMom,2);
 706                         Float_t cophi = TMath::Cos(TMath::ATan2(localMom[0], localMom[1]));
 707                         Float_t t = (1. - .025 / cophi) * (1. - .05 /  cophi);
 708                         /**************** Photons lost in second grid have to be calculated by hand************/
 709                         gMC->GetRandom()->RndmArray(1,ranf);
 710                         if (ranf[0] > t) {
 711                           gMC->StopTrack();
 712                           ckovData[13] = 5;
 713                           AddCerenkov(gAlice->GetCurrentTrackNumber(),vol,ckovData);
 714                           //printf("Added One (1)!\n");
 715                           //printf("Lost one in grid\n");
 716                         }
 717                         /**********************************************************************************/
 718                       }    //gap
 719
 720                     //printf("Current volume (should be CSI) (1): %s\n",gMC->CurrentVolName());
 721                     //printf("VolId: %d, CurrentVolID: %d\n",gMC->VolId("CSI "),gMC->CurrentVolID(copy));
 722                     if (gMC->VolId("CSI ")==gMC->CurrentVolID(copy))             //is it in csi?
 723                       {
 724                         //printf("Got in CSI\n");
 725                         gMC->TrackMomentum(momentum);
 726                         mom[0]=momentum(0);
 727                         mom[1]=momentum(1);
 728                         mom[2]=momentum(2);
 729                         mom[3]=momentum(3);
 730
 731                         gMC->Gmtod(mom,localMom,2);
 732                         /********* Photons lost by Fresnel reflection have to be calculated by hand********/
 733                         /***********************Cerenkov phtons (always polarised)*************************/
 734                         Double_t localTc = localMom[0]*localMom[0]+localMom[2]*localMom[2];
 735                         Double_t localRt = TMath::Sqrt(localTc);
 736                         localTheta   = Float_t(TMath::ATan2(localRt,Double_t(localMom[1])));
 737                         Double_t cotheta = TMath::Abs(cos(localTheta));
 738                         Float_t t = Fresnel(ckovEnergy*1e9,cotheta,1);
 739                             gMC->GetRandom()->RndmArray(1,ranf);
 740                             if (ranf[0] < t) {
 741                               gMC->StopTrack();
 742                               ckovData[13] = 6;
 743                               AddCerenkov(gAlice->GetCurrentTrackNumber(),vol,ckovData);
 744
 745                               //printf("Added One (2)!\n");
 746                               //printf("Lost by Fresnel\n");
 747                             }
 748                             /**********************************************************************************/
 749                       }
 750                   } //track entering?
 751
 752
 753                   /********************Evaluation of losses************************/
 754                   /******************still in the old fashion**********************/
 755
 756                   TArrayI procs;
 757                   Int_t i1 = gMC->StepProcesses(procs);            //number of physics mechanisms acting on the particle
 758                   for (Int_t i = 0; i < i1; ++i) {
 759                     //        Reflection loss
 760                     if (procs[i] == kPLightReflection) {        //was it reflected
 761                       ckovData[13]=10;
 762                       if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
 763                         ckovData[13]=1;
 764                       if (gMC->CurrentVolID(copy) == gMC->VolId("QUAR"))
 765                         ckovData[13]=2;
 766                       //gMC->StopTrack();
 767                       //AddCerenkov(gAlice->GetCurrentTrackNumber(),vol,ckovData);
 768                     } //reflection question
 769
 770                     //        Absorption loss
 771                     else if (procs[i] == kPLightAbsorption) {              //was it absorbed?
 772                       //printf("Got in absorption\n");
 773                       ckovData[13]=20;
 774                       if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
 775                         ckovData[13]=11;
 776                       if (gMC->CurrentVolID(copy) == gMC->VolId("QUAR"))
 777                         ckovData[13]=12;
 778                       if (gMC->CurrentVolID(copy) == gMC->VolId("META"))
 779                         ckovData[13]=13;
 780                       if (gMC->CurrentVolID(copy) == gMC->VolId("GAP "))
 781                         ckovData[13]=13;
 782
 783                       if (gMC->CurrentVolID(copy) == gMC->VolId("SRIC"))
 784                         ckovData[13]=15;
 785
 786                       //        CsI inefficiency
 787                       if (gMC->CurrentVolID(copy) == gMC->VolId("CSI ")) {
 788                         ckovData[13]=16;
 789                       }
 790                       gMC->StopTrack();
 791                       AddCerenkov(gAlice->GetCurrentTrackNumber(),vol,ckovData);
 792                       //printf("Added One (3)!\n");
 793                       //printf("Added cerenkov %d\n",fCkovNumber);
 794                     } //absorption question
 795
 796
 797                     //        Photon goes out of tracking scope
 798                     else if (procs[i] == kPStop) {                 //is it below energy treshold?
 799                       ckovData[13]=21;
 800                       gMC->StopTrack();
 801                       AddCerenkov(gAlice->GetCurrentTrackNumber(),vol,ckovData);
 802                       //printf("Added One (4)!\n");
 803                     }   // energy treshold question
 804                   }  //number of mechanisms cycle
 805                   /**********************End of evaluation************************/
 806                 } //freon production question
 807             } //energy interval question
 808         //}//inside the proximity gap question
 809     } //cerenkov photon question
 810
 811     /**************************************End of Production Parameters Storing*********************/
 812
 813
 814     /*******************************Treat photons that hit the CsI (Ckovs and Feedbacks)************/
 815
 816     if (gMC->TrackPid() == 50000050 || gMC->TrackPid() == 50000051) {
 817       //printf("Cerenkov\n");
 818
 819       //if (gMC->TrackPid() == 50000051)
 820         //printf("Tracking a feedback\n");
 821
 822       if (gMC->VolId("CSI ")==gMC->CurrentVolID(copy))
 823         {
 824           //printf("Current volume (should be CSI) (2): %s\n",gMC->CurrentVolName());
 825           //printf("VolId: %d, CurrentVolID: %d\n",gMC->VolId("CSI "),gMC->CurrentVolID(copy));
 826           //printf("Got in CSI\n");
 827           //printf("Tracking a %d\n",gMC->TrackPid());
 828           if (gMC->Edep() > 0.){
 829                 gMC->TrackPosition(position);
 830                 gMC->TrackMomentum(momentum);
 831                 pos[0]=position(0);
 832                 pos[1]=position(1);
 833                 pos[2]=position(2);
 834                 mom[0]=momentum(0);
 835                 mom[1]=momentum(1);
 836                 mom[2]=momentum(2);
 837                 mom[3]=momentum(3);
 838                 Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
 839                 Double_t rt = TMath::Sqrt(tc);
 840                 theta   = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
 841                 phi     = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
 842
 843                 gMC->CurrentVolOffID(2,copy);
 844                 vol[0]=copy;
 845                 idvol=vol[0]-1;
 846
 847
 848                 gMC->Gmtod(pos,localPos,1);
 849
 850                 //Chamber(idvol).GlobaltoLocal(pos,localPos);
 851
 852                 gMC->Gmtod(mom,localMom,2);
 853
 854                 //Chamber(idvol).GlobaltoLocal(mom,localMom);
 855
 856                 gMC->CurrentVolOffID(2,copy);
 857                 vol[0]=copy;
 858                 idvol=vol[0]-1;
 859
 860                 //Int_t sector=((AliRICHChamber*) (*fChambers)[idvol])
 861                         //->Sector(localPos[0], localPos[2]);
 862                 //printf("Sector:%d\n",sector);
 863
 864                 /*if (gMC->TrackPid() == 50000051){
 865                   fFeedbacks++;
 866                   printf("Feedbacks:%d\n",fFeedbacks);
 867                 }*/
 868
 869         //PH            ((AliRICHChamber*) (*fChambers)[idvol])
 870                 ((AliRICHChamber*)fChambers->At(idvol))
 871                     ->SigGenInit(localPos[0], localPos[2], localPos[1]);
 872                 if(idvol<kNCH) {
 873                     ckovData[0] = gMC->TrackPid();        // particle type
 874                     ckovData[1] = pos[0];                 // X-position for hit
 875                     ckovData[2] = pos[1];                 // Y-position for hit
 876                     ckovData[3] = pos[2];                 // Z-position for hit
 877                     ckovData[4] = theta;                      // theta angle of incidence
 878                     ckovData[5] = phi;                      // phi angle of incidence
 879                     ckovData[8] = (Float_t) fNsdigits;      // first sdigit
 880                     ckovData[9] = -1;                       // last pad hit
 881                     ckovData[13] = 4;                       // photon was detected
 882                     ckovData[14] = mom[0];
 883                     ckovData[15] = mom[1];
 884                     ckovData[16] = mom[2];
 885
 886                     destep = gMC->Edep();
 887                     gMC->SetMaxStep(kBig);
 888                     cherenkovLoss  += destep;
 889                     ckovData[7]=cherenkovLoss;
 890
 891                     ckovData[17] = Hits2SDigits(localPos[0],localPos[2],cherenkovLoss,idvol,kPhoton);//for photons in CsI
 892
 893                     if (fNsdigits > (Int_t)ckovData[8]) {
 894                         ckovData[8]= ckovData[8]+1;
 895                         ckovData[9]= (Float_t) fNsdigits;
 896                     }
 897
 898
 899                     //TClonesArray *Hits = RICH->Hits();
 900                     AliRICHhit *mipHit =  (AliRICHhit*) (fHits->UncheckedAt(0));
 901                     if (mipHit)
 902                       {
 903                         mom[0] = current->Px();
 904                         mom[1] = current->Py();
 905                         mom[2] = current->Pz();
 906                         Float_t mipPx = mipHit->MomX();
 907                         Float_t mipPy = mipHit->MomY();
 908                         Float_t mipPz = mipHit->MomZ();
 909
 910                         Float_t r = mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2];
 911                         Float_t rt = TMath::Sqrt(r);
 912                         Float_t mipR = mipPx*mipPx + mipPy*mipPy + mipPz*mipPz;
 913                         Float_t mipRt = TMath::Sqrt(mipR);
 914                         if ((rt*mipRt) > 0)
 915                           {
 916                             coscerenkov = (mom[0]*mipPx + mom[1]*mipPy + mom[2]*mipPz)/(rt*mipRt);
 917                           }
 918                         else
 919                           {
 920                             coscerenkov = 0;
 921                           }
 922                         Float_t cherenkov = TMath::ACos(coscerenkov);
 923                         ckovData[18]=cherenkov;
 924                       }
 925                     //if (sector != -1)
 926                     //{
 927                     AddHit(gAlice->GetCurrentTrackNumber(),vol,ckovData);
 928                     AddCerenkov(gAlice->GetCurrentTrackNumber(),vol,ckovData);
 929                     //printf("Added One (5)!\n");
 930                     //}
 931                 }
 932             }
 933         }
 934     }
 935
 936     /***********************************************End of photon hits*********************************************/
 937
 938
 939     /**********************************************Charged particles treatment*************************************/
 940
 941     else if (gMC->TrackCharge()){
 942 //If MIP
 943         /*if (gMC->IsTrackEntering())
 944           {
 945             hits[13]=20;//is track entering?
 946           }*/
 947         if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
 948           {
 949             gMC->TrackMomentum(momentum);
 950             mom[0]=momentum(0);
 951             mom[1]=momentum(1);
 952             mom[2]=momentum(2);
 953             mom[3]=momentum(3);
 954             hits [19] = mom[0];
 955             hits [20] = mom[1];
 956             hits [21] = mom[2];
 957             fFreonProd=1;
 958           }
 959
 960         if (gMC->VolId("GAP ")== gMC->CurrentVolID(copy)) {//is in GAP?
 961 // Get current particle id (ipart), track position (pos)  and momentum (mom)
 962
 963             gMC->CurrentVolOffID(3,copy);
 964             vol[0]=copy;
 965             idvol=vol[0]-1;
 966
 967             //Int_t sector=((AliRICHChamber*) (*fChambers)[idvol])
 968                         //->Sector(localPos[0], localPos[2]);
 969             //printf("Sector:%d\n",sector);
 970
 971             gMC->TrackPosition(position);
 972             gMC->TrackMomentum(momentum);
 973             pos[0]=position(0);
 974             pos[1]=position(1);
 975             pos[2]=position(2);
 976             mom[0]=momentum(0);
 977             mom[1]=momentum(1);
 978             mom[2]=momentum(2);
 979             mom[3]=momentum(3);
 980
 981             gMC->Gmtod(pos,localPos,1);
 982
 983             //Chamber(idvol).GlobaltoLocal(pos,localPos);
 984
 985             gMC->Gmtod(mom,localMom,2);
 986
 987             //Chamber(idvol).GlobaltoLocal(mom,localMom);
 988
 989             ipart  = gMC->TrackPid();
 990             //
 991             // momentum loss and steplength in last step
 992             destep = gMC->Edep();
 993             step   = gMC->TrackStep();
 994
 995             //
 996             // record hits when track enters ...
 997             if( gMC->IsTrackEntering()) {
 998 //              gMC->SetMaxStep(fMaxStepGas);
 999                 Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
1000                 Double_t rt = TMath::Sqrt(tc);
1001                 theta   = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
1002                 phi     = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
1003
1004
1005                 Double_t localTc = localMom[0]*localMom[0]+localMom[2]*localMom[2];
1006                 Double_t localRt = TMath::Sqrt(localTc);
1007                 localTheta   = Float_t(TMath::ATan2(localRt,Double_t(localMom[1])))*kRaddeg;
1008                 localPhi     = Float_t(TMath::ATan2(Double_t(localMom[2]),Double_t(localMom[0])))*kRaddeg;
1009
1010                 hits[0] = Float_t(ipart);         // particle type
1011                 hits[1] = localPos[0];                 // X-position for hit
1012                 hits[2] = localPos[1];                 // Y-position for hit
1013                 hits[3] = localPos[2];                 // Z-position for hit
1014                 hits[4] = localTheta;                  // theta angle of incidence
1015                 hits[5] = localPhi;                    // phi angle of incidence
1016                 hits[8] = (Float_t) fNsdigits;    // first sdigit
1017                 hits[9] = -1;                     // last pad hit
1018                 hits[13] = fFreonProd;           // did id hit the freon?
1019                 hits[14] = mom[0];
1020                 hits[15] = mom[1];
1021                 hits[16] = mom[2];
1022                 hits[18] = 0;               // dummy cerenkov angle
1023
1024                 tlength = 0;
1025                 eloss   = 0;
1026                 fFreonProd = 0;
1027
1028                 Chamber(idvol).LocaltoGlobal(localPos,hits+1);
1029
1030
1031                 //To make chamber coordinates x-y had to pass localPos[0], localPos[2]
1032                 xhit    = localPos[0];
1033                 yhit    = localPos[2];
1034                 // Only if not trigger chamber
1035                 if(idvol<kNCH) {
1036                     //
1037                     //  Initialize hit position (cursor) in the segmentation model
1038           //PH              ((AliRICHChamber*) (*fChambers)[idvol])
1039                     ((AliRICHChamber*)fChambers->At(idvol))
1040                         ->SigGenInit(localPos[0], localPos[2], localPos[1]);
1041                 }
1042             }
1043
1044             //
1045             // Calculate the charge induced on a pad (disintegration) in case
1046             //
1047             // Mip left chamber ...
1048             if( gMC->IsTrackExiting() || gMC->IsTrackStop() || gMC->IsTrackDisappeared()){
1049                 gMC->SetMaxStep(kBig);
1050                 eloss   += destep;
1051                 tlength += step;
1052
1053
1054                 // Only if not trigger chamber
1055                 if(idvol<kNCH) {
1056                   if (eloss > 0)
1057                     {
1058                       if(gMC->TrackPid() == kNeutron)
1059                         printf("\n\n\n\n\n Neutron Making Pad Hit!!! \n\n\n\n");
1060                       hits[17] = Hits2SDigits(xhit,yhit,eloss,idvol,kMip); //for MIP
1061                     }
1062                 }
1063
1064                 hits[6]=tlength;
1065                 hits[7]=eloss;
1066                 if (fNsdigits > (Int_t)hits[8]) {
1067                     hits[8]= hits[8]+1;
1068                     hits[9]= (Float_t) fNsdigits;
1069                 }
1070
1071                 //if(sector !=-1)
1072                 new(lhits[fNhits++]) AliRICHhit(fIshunt,gAlice->GetCurrentTrackNumber(),vol,hits);
1073                 eloss = 0;
1074                 //
1075                 // Check additional signal generation conditions
1076                 // defined by the segmentation
1077                 // model (boundary crossing conditions)
1078             }else if(((AliRICHChamber*)fChambers->At(idvol))->SigGenCond(localPos[0], localPos[2], localPos[1])){
1079                 ((AliRICHChamber*)fChambers->At(idvol))->SigGenInit(localPos[0], localPos[2], localPos[1]);
1080                 if (eloss > 0)
1081                   {
1082                     if(gMC->TrackPid() == kNeutron)
1083                       printf("\n\n\n\n\n Neutron Making Pad Hit!!! \n\n\n\n");
1084                     hits[17] = Hits2SDigits(xhit,yhit,eloss,idvol,kMip);//for n
1085                   }
1086                 xhit     = localPos[0];
1087                 yhit     = localPos[2];
1088                 eloss    = destep;
1089                 tlength += step ;
1090                 //
1091                 // nothing special  happened, add up energy loss
1092             } else {
1093                 eloss   += destep;
1094                 tlength += step ;
1095             }
1096         }//is in GAP?
1097       }//is MIP?
1098     /*************************************************End of MIP treatment**************************************/
1099 }//void AliRICHv3::StepManager()
1100 //__________________________________________________________________________________________________
1101 Int_t AliRICHv3::Hits2SDigits(Float_t xhit,Float_t yhit,Float_t eloss, Int_t idvol, ResponseType res)
1102 {//calls the charge disintegration method of the current chamber and adds all generated sdigits to the list of digits
1103
1104    Float_t newclust[4][500];
1105    Int_t clhits[5];
1106    Int_t iNdigits;
1107    clhits[0]=fNhits+1;
1108
1109   ((AliRICHChamber*)fChambers->At(idvol))->DisIntegration(eloss, xhit, yhit, iNdigits,newclust, res);
1110
1111   for (Int_t i=0; i<iNdigits; i++) {
1112     if (Int_t(newclust[0][i]) > 0) {
1113             clhits[1] = Int_t(newclust[0][i]);//  Cluster Charge
1114             clhits[2] = Int_t(newclust[1][i]);//  Pad: ix
1115             clhits[3] = Int_t(newclust[2][i]);//  Pad: iy
1116             clhits[4] = Int_t(newclust[3][i]);//  Pad: chamber sector
1117             AddSpecialOld(clhits);
1118         }
1119     }
1120   return iNdigits;
1121 }//Int_t AliRICHv3::Hits2SDigits(Float_t xhit,Float_t yhit,Float_t eloss, Int_t idvol, ResponseType res)
1122 //__________________________________________________________________________________________________
1123 void AliRICHv3::DiagnosticsFE(Int_t evNumber1,Int_t evNumber2)
1124 {
1125
1126   Int_t NpadX = 162;                 // number of pads on X
1127   Int_t NpadY = 162;                 // number of pads on Y
1128
1129   Int_t Pad[162][162];
1130   for (Int_t i=0;i<NpadX;i++) {
1131     for (Int_t j=0;j<NpadY;j++) {
1132       Pad[i][j]=0;
1133     }
1134   }
1135
1136   //  Create some histograms
1137
1138   TH1F *pionspectra1 = new TH1F("pionspectra1","Pion Spectra",200,-4,2);
1139   TH1F *pionspectra2 = new TH1F("pionspectra2","Pion Spectra",200,-4,2);
1140   TH1F *pionspectra3 = new TH1F("pionspectra3","Pion Spectra",200,-4,2);
1141   TH1F *protonspectra1 = new TH1F("protonspectra1","Proton Spectra",200,-4,2);
1142   TH1F *protonspectra2 = new TH1F("protonspectra2","Proton Spectra",200,-4,2);
1143   TH1F *protonspectra3 = new TH1F("protonspectra3","Proton Spectra",200,-4,2);
1144   TH1F *kaonspectra1 = new TH1F("kaonspectra1","Kaon Spectra",100,-4,2);
1145   TH1F *kaonspectra2 = new TH1F("kaonspectra2","Kaon Spectra",100,-4,2);
1146   TH1F *kaonspectra3 = new TH1F("kaonspectra3","Kaon Spectra",100,-4,2);
1147   TH1F *electronspectra1 = new TH1F("electronspectra1","Electron Spectra",100,-4,2);
1148   TH1F *electronspectra2 = new TH1F("electronspectra2","Electron Spectra",100,-4,2);
1149   TH1F *electronspectra3 = new TH1F("electronspectra3","Electron Spectra",100,-4,2);
1150   TH1F *muonspectra1 = new TH1F("muonspectra1","Muon Spectra",100,-4,2);
1151   TH1F *muonspectra2 = new TH1F("muonspectra2","Muon Spectra",100,-4,2);
1152   TH1F *muonspectra3 = new TH1F("muonspectra3","Muon Spectra",100,-4,2);
1153   TH1F *neutronspectra1 = new TH1F("neutronspectra1","Neutron Spectra",100,-4,2);
1154   TH1F *neutronspectra2 = new TH1F("neutronspectra2","Neutron Spectra",100,-4,2);
1155   TH1F *neutronspectra3 = new TH1F("neutronspectra2","Neutron Spectra",100,-4,2);
1156   TH1F *chargedspectra1 = new TH1F("chargedspectra1","Charged particles above 1 GeV Spectra",100,-1,3);
1157   TH1F *chargedspectra2 = new TH1F("chargedspectra2","Charged particles above 1 GeV Spectra",100,-1,3);
1158   TH1F *chargedspectra3 = new TH1F("chargedspectra2","Charged particles above 1 GeV Spectra",100,-1,3);
1159   TH1F *pionptspectrafinal = new TH1F("pionptspectrafinal","Primary Pions Transverse Momenta at HMPID",20,0,5);
1160   TH1F *pionptspectravertex = new TH1F("pionptspectravertex","Primary Pions Transverse Momenta at vertex",20,0,5);
1161   TH1F *kaonptspectrafinal = new TH1F("kaonptspectrafinal","Primary Kaons Transverse Momenta at HMPID",20,0,5);
1162   TH1F *kaonptspectravertex = new TH1F("kaonptspectravertex","Primary Kaons Transverse Momenta at vertex",20,0,5);
1163   //TH1F *hitsPhi = new TH1F("hitsPhi","Distribution of phi angle of incidence",100,-180,180);
1164   TH1F *hitsTheta = new TH1F("hitsTheta","Distribution of Theta angle of incidence, all tracks",100,0,50);
1165   TH1F *hitsTheta500MeV = new TH1F("hitsTheta500MeV","Distribution of Theta angle of incidence, 0.5-1 GeV primary tracks",100,0,50);
1166   TH1F *hitsTheta1GeV = new TH1F("hitsTheta1GeV","Distribution of Theta angle of incidence, 1-2 GeV primary tracks",100,0,50);
1167   TH1F *hitsTheta2GeV = new TH1F("hitsTheta2GeV","Distribution of Theta angle of incidence, 2-3 GeV primary tracks",100,0,50);
1168   TH1F *hitsTheta3GeV = new TH1F("hitsTheta3GeV","Distribution of Theta angle of incidence, >3 GeV primary tracks",100,0,50);
1169   TH2F *production = new TH2F("production","Mother production vertices",100,-300,300,100,0,600);
1170
1171
1172
1173
1174 //   Start loop over events
1175
1176   Int_t pion=0, kaon=0, proton=0, electron=0, positron=0, neutron=0, highneutrons=0, muon=0;
1177   Int_t chargedpions=0,primarypions=0,highprimarypions=0,chargedkaons=0,primarykaons=0,highprimarykaons=0;
1178   Int_t photons=0, primaryphotons=0, highprimaryphotons=0;
1179   TRandom* random=0;
1180
1181    for (int nev=0; nev<= evNumber2; nev++) {
1182        Int_t nparticles = gAlice->GetEvent(nev);
1183
1184
1185        if (nev < evNumber1) continue;
1186        if (nparticles <= 0) return;
1187
1188 // Get pointers to RICH detector and Hits containers
1189
1190        AliRICH *pRICH = (AliRICH *) gAlice->GetDetector("RICH");
1191
1192        TTree *treeH = TreeH();
1193        Int_t ntracks =(Int_t) treeH->GetEntries();
1194
1195 // Start loop on tracks in the hits containers
1196
1197        for (Int_t track=0; track<ntracks;track++) {
1198            printf ("Processing Track: %d\n",track);
1199            gAlice->ResetHits();
1200            treeH->GetEvent(track);
1201
1202            for(AliRICHhit* mHit=(AliRICHhit*)pRICH->FirstHit(-1);
1203                mHit;
1204                mHit=(AliRICHhit*)pRICH->NextHit())
1205              {
1206                //Int_t nch  = mHit->fChamber;              // chamber number
1207                //Float_t x  = mHit->X();                    // x-pos of hit
1208                //Float_t y  = mHit->Z();                    // y-pos
1209                //Float_t z  = mHit->Y();
1210                //Float_t phi = mHit->Phi();                 //Phi angle of incidence
1211                Float_t theta = mHit->Theta();             //Theta angle of incidence
1212                Float_t px = mHit->MomX();
1213                Float_t py = mHit->MomY();
1214                Int_t index = mHit->Track();
1215                Int_t particle = (Int_t)(mHit->Particle());
1216                Float_t R;
1217                Float_t PTfinal;
1218                Float_t PTvertex;
1219
1220               TParticle *current = gAlice->Particle(index);
1221
1222               //Float_t energy=current->Energy();
1223
1224               R=TMath::Sqrt(current->Vx()*current->Vx() + current->Vy()*current->Vy());
1225               PTfinal=TMath::Sqrt(px*px + py*py);
1226               PTvertex=TMath::Sqrt(current->Px()*current->Px() + current->Py()*current->Py());
1227
1228
1229
1230               if (TMath::Abs(particle) < 10000000)
1231                 {
1232                   hitsTheta->Fill(theta,(float) 1);
1233                   if (R<5)
1234                     {
1235                       if (PTvertex>.5 && PTvertex<=1)
1236                         {
1237                           hitsTheta500MeV->Fill(theta,(float) 1);
1238                         }
1239                       if (PTvertex>1 && PTvertex<=2)
1240                         {
1241                           hitsTheta1GeV->Fill(theta,(float) 1);
1242                         }
1243                       if (PTvertex>2 && PTvertex<=3)
1244                         {
1245                           hitsTheta2GeV->Fill(theta,(float) 1);
1246                         }
1247                       if (PTvertex>3)
1248                         {
1249                           hitsTheta3GeV->Fill(theta,(float) 1);
1250                         }
1251                     }
1252
1253                 }
1254
1255               //if (nch == 3)
1256                 //{
1257
1258               if (TMath::Abs(particle) < 50000051)
1259                 {
1260                   //if (TMath::Abs(particle) == 50000050 || TMath::Abs(particle) == 2112)
1261                   if (TMath::Abs(particle) == 2112 || TMath::Abs(particle) == 50000050)
1262                     {
1263                       //gMC->Rndm(&random, 1);
1264                       if (random->Rndm() < .1)
1265                         production->Fill(current->Vz(),R,(float) 1);
1266                       if (TMath::Abs(particle) == 50000050)
1267                         //if (TMath::Abs(particle) > 50000000)
1268                         {
1269                           photons +=1;
1270                           if (R<5)
1271                             {
1272                               primaryphotons +=1;
1273                               if (current->Energy()>0.001)
1274                                 highprimaryphotons +=1;
1275                             }
1276                         }
1277                       if (TMath::Abs(particle) == 2112)
1278                         {
1279                           neutron +=1;
1280                           if (current->Energy()>0.0001)
1281                             highneutrons +=1;
1282                         }
1283                     }
1284                   if (TMath::Abs(particle) < 50000000)
1285                     {
1286                       production->Fill(current->Vz(),R,(float) 1);
1287                     }
1288                   //mip->Fill(x,y,(float) 1);
1289                 }
1290
1291               if (TMath::Abs(particle)==211 || TMath::Abs(particle)==111)
1292                 {
1293                   if (R<5)
1294                     {
1295                       pionptspectravertex->Fill(PTvertex,(float) 1);
1296                       pionptspectrafinal->Fill(PTfinal,(float) 1);
1297                     }
1298                 }
1299
1300               if (TMath::Abs(particle)==321 || TMath::Abs(particle)==130 || TMath::Abs(particle)==310
1301                   || TMath::Abs(particle)==311)
1302                 {
1303                   if (R<5)
1304                     {
1305                       kaonptspectravertex->Fill(PTvertex,(float) 1);
1306                       kaonptspectrafinal->Fill(PTfinal,(float) 1);
1307                     }
1308                 }
1309
1310
1311               if (TMath::Abs(particle)==211 || TMath::Abs(particle)==111)
1312                 {
1313                   pionspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
1314                   if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
1315                     pionspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
1316                   if (R>250 && R<450)
1317                     {
1318                       pionspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
1319                     }
1320                   pion +=1;
1321                   if (TMath::Abs(particle)==211)
1322                     {
1323                       chargedpions +=1;
1324                       if (R<5)
1325                         {
1326                           primarypions +=1;
1327                           if (current->Energy()>1)
1328                             highprimarypions +=1;
1329                         }
1330                     }
1331                 }
1332               if (TMath::Abs(particle)==2212)
1333                 {
1334                   protonspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
1335                   //ptspectra->Fill(Pt,(float) 1);
1336                   if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
1337                     protonspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
1338                   if (R>250 && R<450)
1339                     protonspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
1340                   proton +=1;
1341                 }
1342               if (TMath::Abs(particle)==321 || TMath::Abs(particle)==130 || TMath::Abs(particle)==310
1343                   || TMath::Abs(particle)==311)
1344                 {
1345                   kaonspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
1346                   //ptspectra->Fill(Pt,(float) 1);
1347                   if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
1348                     kaonspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
1349                   if (R>250 && R<450)
1350                     kaonspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
1351                   kaon +=1;
1352                   if (TMath::Abs(particle)==321)
1353                     {
1354                       chargedkaons +=1;
1355                       if (R<5)
1356                         {
1357                           primarykaons +=1;
1358                           if (current->Energy()>1)
1359                             highprimarykaons +=1;
1360                         }
1361                     }
1362                 }
1363               if (TMath::Abs(particle)==11)
1364                 {
1365                   electronspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
1366                   //ptspectra->Fill(Pt,(float) 1);
1367                   if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
1368                     electronspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
1369                   if (R>250 && R<450)
1370                     electronspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
1371                   if (particle == 11)
1372                     electron +=1;
1373                   if (particle == -11)
1374                     positron +=1;
1375                 }
1376               if (TMath::Abs(particle)==13)
1377                 {
1378                   muonspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
1379                   //ptspectra->Fill(Pt,(float) 1);
1380                   if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
1381                     muonspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
1382                   if (R>250 && R<450)
1383                     muonspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
1384                   muon +=1;
1385                 }
1386               if (TMath::Abs(particle)==2112)
1387                 {
1388                   neutronspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
1389                   //ptspectra->Fill(Pt,(float) 1);
1390                   if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
1391                     neutronspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
1392                   if (R>250 && R<450)
1393                     {
1394                       neutronspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
1395                     }
1396                   neutron +=1;
1397                 }
1398               if(TMath::Abs(particle)==211 || TMath::Abs(particle)==2212 || TMath::Abs(particle)==321)
1399                 {
1400                   if (current->Energy()-current->GetCalcMass()>1)
1401                     {
1402                       chargedspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
1403                       if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
1404                         chargedspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
1405                       if (R>250 && R<450)
1406                         chargedspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
1407                     }
1408                 }
1409               // Fill the histograms
1410               //Nh1+=nhits;
1411               //h->Fill(x,y,(float) 1);
1412               //}
1413               //}
1414            }
1415
1416        }
1417
1418    }
1419    //   }
1420
1421    TStyle *mystyle=new TStyle("Plain","mystyle");
1422    mystyle->SetPalette(1,0);
1423    mystyle->cd();
1424
1425    //Create canvases, set the view range, show histograms
1426
1427     TCanvas *c2 = new TCanvas("c2","Angles of incidence",150,150,100,150);
1428     c2->Divide(2,2);
1429     //c2->SetFillColor(42);
1430
1431     c2->cd(1);
1432     hitsTheta500MeV->SetFillColor(5);
1433     hitsTheta500MeV->Draw();
1434     c2->cd(2);
1435     hitsTheta1GeV->SetFillColor(5);
1436     hitsTheta1GeV->Draw();
1437     c2->cd(3);
1438     hitsTheta2GeV->SetFillColor(5);
1439     hitsTheta2GeV->Draw();
1440     c2->cd(4);
1441     hitsTheta3GeV->SetFillColor(5);
1442     hitsTheta3GeV->Draw();
1443
1444
1445
1446     TCanvas *c15 = new TCanvas("c15","Mothers Production Vertices",50,50,600,600);
1447     c15->cd();
1448     production->SetFillColor(42);
1449     production->SetXTitle("z (m)");
1450     production->SetYTitle("R (m)");
1451     production->Draw();
1452
1453     TCanvas *c10 = new TCanvas("c10","Pt Spectra",50,50,600,700);
1454     c10->Divide(2,2);
1455     c10->cd(1);
1456     pionptspectravertex->SetFillColor(5);
1457     pionptspectravertex->SetXTitle("Pt (GeV)");
1458     pionptspectravertex->Draw();
1459     c10->cd(2);
1460     pionptspectrafinal->SetFillColor(5);
1461     pionptspectrafinal->SetXTitle("Pt (GeV)");
1462     pionptspectrafinal->Draw();
1463     c10->cd(3);
1464     kaonptspectravertex->SetFillColor(5);
1465     kaonptspectravertex->SetXTitle("Pt (GeV)");
1466     kaonptspectravertex->Draw();
1467     c10->cd(4);
1468     kaonptspectrafinal->SetFillColor(5);
1469     kaonptspectrafinal->SetXTitle("Pt (GeV)");
1470     kaonptspectrafinal->Draw();
1471
1472
1473    TCanvas *c16 = new TCanvas("c16","Particles Spectra II",150,150,600,350);
1474    c16->Divide(2,1);
1475
1476    c16->cd(1);
1477    //TCanvas *c13 = new TCanvas("c13","Electron Spectra",400,10,600,700);
1478    electronspectra1->SetFillColor(5);
1479    electronspectra1->SetXTitle("log(GeV)");
1480    electronspectra2->SetFillColor(46);
1481    electronspectra2->SetXTitle("log(GeV)");
1482    electronspectra3->SetFillColor(10);
1483    electronspectra3->SetXTitle("log(GeV)");
1484    //c13->SetLogx();
1485    electronspectra1->Draw();
1486    electronspectra2->Draw("same");
1487    electronspectra3->Draw("same");
1488
1489    c16->cd(2);
1490    //TCanvas *c14 = new TCanvas("c14","Muon Spectra",400,10,600,700);
1491    muonspectra1->SetFillColor(5);
1492    muonspectra1->SetXTitle("log(GeV)");
1493    muonspectra2->SetFillColor(46);
1494    muonspectra2->SetXTitle("log(GeV)");
1495    muonspectra3->SetFillColor(10);
1496    muonspectra3->SetXTitle("log(GeV)");
1497    //c14->SetLogx();
1498    muonspectra1->Draw();
1499    muonspectra2->Draw("same");
1500    muonspectra3->Draw("same");
1501
1502    //c16->cd(3);
1503    //TCanvas *c16 = new TCanvas("c16","Neutron Spectra",400,10,600,700);
1504    //neutronspectra1->SetFillColor(42);
1505    //neutronspectra1->SetXTitle("log(GeV)");
1506    //neutronspectra2->SetFillColor(46);
1507    //neutronspectra2->SetXTitle("log(GeV)");
1508    //neutronspectra3->SetFillColor(10);
1509    //neutronspectra3->SetXTitle("log(GeV)");
1510    //c16->SetLogx();
1511    //neutronspectra1->Draw();
1512    //neutronspectra2->Draw("same");
1513    //neutronspectra3->Draw("same");
1514
1515    TCanvas *c9 = new TCanvas("c9","Particles Spectra",150,150,600,700);
1516    //TCanvas *c9 = new TCanvas("c9","Pion Spectra",400,10,600,700);
1517    c9->Divide(2,2);
1518
1519    c9->cd(1);
1520    pionspectra1->SetFillColor(5);
1521    pionspectra1->SetXTitle("log(GeV)");
1522    pionspectra2->SetFillColor(46);
1523    pionspectra2->SetXTitle("log(GeV)");
1524    pionspectra3->SetFillColor(10);
1525    pionspectra3->SetXTitle("log(GeV)");
1526    //c9->SetLogx();
1527    pionspectra1->Draw();
1528    pionspectra2->Draw("same");
1529    pionspectra3->Draw("same");
1530
1531    c9->cd(2);
1532    //TCanvas *c10 = new TCanvas("c10","Proton Spectra",400,10,600,700);
1533    protonspectra1->SetFillColor(5);
1534    protonspectra1->SetXTitle("log(GeV)");
1535    protonspectra2->SetFillColor(46);
1536    protonspectra2->SetXTitle("log(GeV)");
1537    protonspectra3->SetFillColor(10);
1538    protonspectra3->SetXTitle("log(GeV)");
1539    //c10->SetLogx();
1540    protonspectra1->Draw();
1541    protonspectra2->Draw("same");
1542    protonspectra3->Draw("same");
1543
1544    c9->cd(3);
1545    //TCanvas *c11 = new TCanvas("c11","Kaon Spectra",400,10,600,700);
1546    kaonspectra1->SetFillColor(5);
1547    kaonspectra1->SetXTitle("log(GeV)");
1548    kaonspectra2->SetFillColor(46);
1549    kaonspectra2->SetXTitle("log(GeV)");
1550    kaonspectra3->SetFillColor(10);
1551    kaonspectra3->SetXTitle("log(GeV)");
1552    //c11->SetLogx();
1553    kaonspectra1->Draw();
1554    kaonspectra2->Draw("same");
1555    kaonspectra3->Draw("same");
1556
1557    c9->cd(4);
1558    //TCanvas *c12 = new TCanvas("c12","Charged Particles Spectra",400,10,600,700);
1559    chargedspectra1->SetFillColor(5);
1560    chargedspectra1->SetXTitle("log(GeV)");
1561    chargedspectra2->SetFillColor(46);
1562    chargedspectra2->SetXTitle("log(GeV)");
1563    chargedspectra3->SetFillColor(10);
1564    chargedspectra3->SetXTitle("log(GeV)");
1565    //c12->SetLogx();
1566    chargedspectra1->Draw();
1567    chargedspectra2->Draw("same");
1568    chargedspectra3->Draw("same");
1569
1570
1571
1572    printf("*****************************************\n");
1573    printf("* Particle                   *  Counts  *\n");
1574    printf("*****************************************\n");
1575
1576    printf("* Pions:                     *   %4d   *\n",pion);
1577    printf("* Charged Pions:             *   %4d   *\n",chargedpions);
1578    printf("* Primary Pions:             *   %4d   *\n",primarypions);
1579    printf("* Primary Pions (p>1GeV/c):  *   %4d   *\n",highprimarypions);
1580    printf("* Kaons:                     *   %4d   *\n",kaon);
1581    printf("* Charged Kaons:             *   %4d   *\n",chargedkaons);
1582    printf("* Primary Kaons:             *   %4d   *\n",primarykaons);
1583    printf("* Primary Kaons (p>1GeV/c):  *   %4d   *\n",highprimarykaons);
1584    printf("* Muons:                     *   %4d   *\n",muon);
1585    printf("* Electrons:                 *   %4d   *\n",electron);
1586    printf("* Positrons:                 *   %4d   *\n",positron);
1587    printf("* Protons:                   *   %4d   *\n",proton);
1588    printf("* All Charged:               *   %4d   *\n",(chargedpions+chargedkaons+muon+electron+positron+proton));
1589    printf("*****************************************\n");
1590    //printf("* Photons:                   *   %3.1f   *\n",photons);
1591    //printf("* Primary Photons:           *   %3.1f   *\n",primaryphotons);
1592    //printf("* Primary Photons (p>1MeV/c):*   %3.1f   *\n",highprimaryphotons);
1593    //printf("*****************************************\n");
1594    //printf("* Neutrons:                  *   %3.1f   *\n",neutron);
1595    //printf("* Neutrons (p>100keV/c):     *   %3.1f   *\n",highneutrons);
1596    //printf("*****************************************\n");
1597
1598    if (gAlice->TreeD())
1599      {
1600        gAlice->TreeD()->GetEvent(0);
1601
1602        Float_t occ[7];
1603        Float_t sum=0;
1604        Float_t mean=0;
1605        printf("\n*****************************************\n");
1606        printf("* Chamber   * Digits      * Occupancy   *\n");
1607        printf("*****************************************\n");
1608
1609        for (Int_t ich=0;ich<7;ich++)
1610          {
1611            TClonesArray *Digits = DigitsAddress(ich);    //  Raw clusters branch
1612            Int_t ndigits = Digits->GetEntriesFast();
1613            occ[ich] = Float_t(ndigits)/(160*144);
1614            sum += Float_t(ndigits)/(160*144);
1615            printf("*   %d      *    %d      *   %3.1f%%     *\n",ich,ndigits,occ[ich]*100);
1616          }
1617        mean = sum/7;
1618        printf("*****************************************\n");
1619        printf("* Mean occupancy          *   %3.1f%%     *\n",mean*100);
1620        printf("*****************************************\n");
1621      }
1622
1623   printf("\nEnd of analysis\n");
1624
1625 }//void AliRICHv3::DiagnosticsFE(Int_t evNumber1,Int_t evNumber2)
1626 //__________________________________________________________________________________________________
1627 void AliRICHv3::DiagnosticsSE(Int_t diaglevel,Int_t evNumber1,Int_t evNumber2)
1628 {
1629
1630 AliRICH *pRICH  = (AliRICH*)gAlice->GetDetector("RICH");
1631    AliRICHSegmentationV0*  segmentation;
1632    AliRICHChamber*       chamber;
1633
1634    chamber = &(pRICH->Chamber(0));
1635    segmentation=(AliRICHSegmentationV0*) chamber->GetSegmentationModel();
1636
1637    Int_t NpadX = segmentation->Npx();                 // number of pads on X
1638    Int_t NpadY = segmentation->Npy();                 // number of pads on Y
1639
1640    Int_t xmin= -NpadX/2;
1641    Int_t xmax=  NpadX/2;
1642    Int_t ymin= -NpadY/2;
1643    Int_t ymax=  NpadY/2;
1644
1645    Float_t PTfinal = 0;
1646    Int_t pionCount = 0;
1647    Int_t kaonCount = 0;
1648    Int_t protonCount = 0;
1649
1650    TH2F *feedback = 0;
1651    TH2F *mip = 0;
1652    TH2F *cerenkov = 0;
1653    TH2F *h = 0;
1654    TH1F *hitsX = 0;
1655    TH1F *hitsY = 0;
1656
1657    TH2F *hc0 = new TH2F("hc0","Zoom on center of central chamber",150,-25,25,150,-45,5);
1658
1659    if (diaglevel == 1)
1660      {
1661        printf("Single Ring Hits\n");
1662        feedback = new TH2F("feedback","Feedback hit distribution",150,-20,20,150,-35,5);
1663        mip = new TH2F("mip","Mip hit distribution",150,-20,20,150,-35,5);
1664        cerenkov = new TH2F("cerenkov","Cerenkov hit distribution",150,-20,20,150,-35,5);
1665        h = new TH2F("h","Detector hit distribution",150,-20,20,150,-35,5);
1666        hitsX = new TH1F("hitsX","Distribution of hits along x-axis",150,-50,50);
1667        hitsY = new TH1F("hitsY","Distribution of hits along z-axis",150,-50,50);
1668      }
1669    else
1670      {
1671        printf("Full Event Hits\n");
1672
1673        feedback = new TH2F("feedback","Feedback hit distribution",150,-300,300,150,-300,300);
1674        mip = new TH2F("mip","Mip hit distribution",150,-300,300,150,-300,300);
1675        cerenkov = new TH2F("cerenkov","Cerenkov hit distribution",150,-300,300,150,-300,300);
1676        h = new TH2F("h","Detector hit distribution",150,-300,300,150,-300,300);
1677        hitsX = new TH1F("digitsX","Distribution of hits along x-axis",200,-300,300);
1678        hitsY = new TH1F("digitsY","Distribution of hits along z-axis",200,-300,300);
1679      }
1680
1681
1682
1683    TH2F *hc1 = new TH2F("hc1","Chamber 1 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
1684    TH2F *hc2 = new TH2F("hc2","Chamber 2 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
1685    TH2F *hc3 = new TH2F("hc3","Chamber 3 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
1686    TH2F *hc4 = new TH2F("hc4","Chamber 4 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
1687    TH2F *hc5 = new TH2F("hc5","Chamber 5 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
1688    TH2F *hc6 = new TH2F("hc6","Chamber 6 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
1689    TH2F *hc7 = new TH2F("hc7","Chamber 7 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
1690
1691    TH1F *Clcharge = new TH1F("Clcharge","Cluster Charge Distribution",500,0.,500.);
1692    TH1F *ckovangle = new TH1F("ckovangle","Cerenkov angle per photon",100,.35,.8);
1693    TH1F *hckphi = new TH1F("hckphi","Cerenkov phi angle per photon",620,-3.1,3.1);
1694    TH1F *mother = new TH1F("mother","Cerenkovs per Mip",75,0.,75.);
1695    TH1F *radius = new TH1F("radius","Mean distance to Mip",100,0.,20.);
1696    TH1F *phspectra1 = new TH1F("phspectra1","Detected Photon Spectra",200,5.,10.);
1697    TH1F *phspectra2 = new TH1F("phspectra2","Produced Photon Spectra",200,5.,10.);
1698    TH1F *totalphotonstrack = new TH1F("totalphotonstrack","Produced Photons per Mip",100,200,700.);
1699    TH1F *totalphotonsevent = new TH1F("totalphotonsevent","Produced Photons per Mip",100,200,700.);
1700    //TH1F *feedbacks = new TH1F("feedbacks","Produced Feedbacks per Mip",50,0.5,50.);
1701    TH1F *padnumber = new TH1F("padnumber","Number of pads per cluster",50,-0.5,50.);
1702    TH1F *padsev = new TH1F("padsev","Number of pads hit per MIP",50,0.5,100.);
1703    TH1F *clusev = new TH1F("clusev","Number of clusters per MIP",50,0.5,50.);
1704    TH1F *photev = new TH1F("photev","Number of detected photons per MIP",50,0.5,50.);
1705    TH1F *feedev = new TH1F("feedev","Number of feedbacks per MIP",50,0.5,50.);
1706    TH1F *padsmip = new TH1F("padsmip","Number of pads per event inside MIP region",50,0.5,50.);
1707    TH1F *padscl = new TH1F("padscl","Number of pads per event from cluster count",50,0.5,100.);
1708    TH1F *pionspectra = new TH1F("pionspectra","Pion Spectra",200,.5,10.);
1709    TH1F *protonspectra = new TH1F("protonspectra","Proton Spectra",200,.5,10.);
1710    TH1F *kaonspectra = new TH1F("kaonspectra","Kaon Spectra",100,.5,10.);
1711    TH1F *chargedspectra = new TH1F("chargedspectra","Charged particles above 1 GeV Spectra",100,.5,10.);
1712    TH1F *hitsPhi = new TH1F("hitsPhi","Distribution of phi angle of incidence",50,0,360);
1713    TH1F *hitsTheta = new TH1F("hitsTheta","Distribution of theta angle of incidence",50,0,15);
1714    TH1F *Omega1D = new TH1F("omega","Reconstructed Cerenkov angle per track",50,.5,1);
1715    TH1F *Theta = new TH1F("theta","Reconstructed theta incidence angle per track",100,0,15);
1716    TH1F *Phi = new TH1F("phi","Reconstructed phi incidence per track",100,0,360);
1717    TH1F *Omega3D = new TH1F("omega","Reconstructed Cerenkov angle per track",100,.35,.8);
1718    TH1F *PhotonCer = new TH1F("photoncer","Reconstructed Cerenkov angle per photon",100,.35,.8);
1719    TH2F *PadsUsed = new TH2F("padsused","Pads Used for Reconstruction",100,-30,30,100,-30,30);
1720    TH1F *MeanRadius = new TH1F("radius","Mean Radius for reconstructed track",100,0.,20.);
1721    TH2F *identification = new TH2F("identification","Particle Identification",100,1,5,100,0,.8);
1722    TH1F *OriginalOmega = new TH1F("Original Omega","Cerenkov angle per track",100,.35,.8);
1723    TH1F *OriginalPhi = new TH1F("Original Phi","Distribution of phi angle of incidence per track",100,0,360);
1724    TH1F *OriginalTheta = new TH1F("Original Theta","Distribution of theta angle per track",100,0,15);
1725    TH1F *OmegaError = new TH1F("Omega Error","Difference between original an reconstructed cerenkov angle",100,0,.2);
1726    TH1F *PhiError = new TH1F("Phi Error","Difference between original an reconstructed phi angle",100,0,360);
1727    TH1F *ThetaError = new TH1F("Theta Error","Difference between original an reconstructed phi angle",100,0,15);
1728
1729
1730 //   Start loop over events
1731
1732    Int_t Nh=0;
1733    Int_t pads=0;
1734    Int_t Nh1=0;
1735    Int_t mothers[80000];
1736    Int_t mothers2[80000];
1737    Float_t mom[3];
1738    Int_t nraw=0;
1739    Int_t phot=0;
1740    Int_t feed=0;
1741    Int_t padmip=0;
1742    Float_t x=0,y=0;
1743
1744    Float_t chiSquareOmega = 0;
1745    Float_t chiSquareTheta = 0;
1746    Float_t chiSquarePhi = 0;
1747
1748    Float_t recEffEvent = 0;
1749    Float_t recEffTotal = 0;
1750
1751    Float_t trackglob[3];
1752    Float_t trackloc[3];
1753
1754
1755    for (Int_t i=0;i<100;i++) mothers[i]=0;
1756
1757    for (int nev=0; nev<= evNumber2; nev++) {
1758        Int_t nparticles = gAlice->GetEvent(nev);
1759
1760
1761        //cout<<"nev  "<<nev<<endl;
1762        printf ("\n**********************************\nProcessing Event: %d\n",nev);
1763        //cout<<"nparticles  "<<nparticles<<endl;
1764        printf ("Particles       : %d\n\n",nparticles);
1765        if (nev < evNumber1) continue;
1766        if (nparticles <= 0) return;
1767
1768 // Get pointers to RICH detector and Hits containers
1769
1770
1771        TTree *TH = TreeH();
1772        Stat_t ntracks = TH->GetEntries();
1773
1774        // Start loop on tracks in the hits containers
1775        //Int_t Nc=0;
1776        for (Int_t track=0; track<ntracks;track++) {
1777
1778          printf ("\nProcessing Track: %d\n",track);
1779          gAlice->ResetHits();
1780          TH->GetEvent(track);
1781          Int_t nhits = pRICH->Hits()->GetEntriesFast();
1782          if (nhits) Nh+=nhits;
1783          printf("Hits            : %d\n",nhits);
1784          for(AliRICHhit* mHit=(AliRICHhit*)pRICH->FirstHit(-1);
1785              mHit;
1786              mHit=(AliRICHhit*)pRICH->NextHit())
1787            {
1788              Int_t nch  = mHit->Chamber();              // chamber number
1789              trackglob[0] = mHit->X();                 // x-pos of hit
1790              trackglob[1] = mHit->Y();
1791              trackglob[2] = mHit->Z();                 // y-pos of hit
1792              //x  = mHit->X();                           // x-pos of hit
1793              //y  = mHit->Z();                           // y-pos
1794              Float_t phi = mHit->Phi();                 //Phi angle of incidence
1795              Float_t theta = mHit->Theta();             //Theta angle of incidence
1796              Int_t index = mHit->Track();
1797              Int_t particle = (Int_t)(mHit->Particle());
1798              //Int_t freon = (Int_t)(mHit->fLoss);
1799              Float_t px = mHit->MomX();
1800              Float_t py = mHit->MomY();
1801
1802              if (TMath::Abs(particle) < 10000000)
1803                {
1804                  PTfinal=TMath::Sqrt(px*px + py*py);
1805                }
1806
1807              chamber = &(pRICH->Chamber(nch-1));
1808
1809
1810              chamber->GlobaltoLocal(trackglob,trackloc);
1811
1812              chamber->LocaltoGlobal(trackloc,trackglob);
1813
1814
1815              x=trackloc[0];
1816              y=trackloc[2];
1817
1818              hitsX->Fill(x,(float) 1);
1819              hitsY->Fill(y,(float) 1);
1820
1821
1822               TParticle *current = (TParticle*)gAlice->Particle(index);
1823
1824               hitsTheta->Fill(theta,(float) 1);
1825
1826               if (current->GetPdgCode() < 10000000)
1827                 {
1828                   mip->Fill(x,y,(float) 1);
1829                   hitsPhi->Fill(TMath::Abs(phi),(float) 1);
1830                 }
1831
1832               if (TMath::Abs(particle)==211 || TMath::Abs(particle)==111)
1833                 {
1834                   pionspectra->Fill(current->Energy() - current->GetCalcMass(),(float) 1);
1835                 }
1836               if (TMath::Abs(particle)==2212)
1837                 {
1838                   protonspectra->Fill(current->Energy() - current->GetCalcMass(),(float) 1);
1839                 }
1840               if (TMath::Abs(particle)==321 || TMath::Abs(particle)==130 || TMath::Abs(particle)==310
1841                   || TMath::Abs(particle)==311)
1842                 {
1843                   kaonspectra->Fill(current->Energy() - current->GetCalcMass(),(float) 1);
1844                 }
1845               if(TMath::Abs(particle)==211 || TMath::Abs(particle)==2212 || TMath::Abs(particle)==321)
1846                 {
1847                   if (current->Energy() - current->GetCalcMass()>1)
1848                     chargedspectra->Fill(current->Energy() - current->GetCalcMass(),(float) 1);
1849                 }
1850               //printf("Hits:%d\n",hit);
1851               //printf ("Chamber number:%d x:%f y:%f\n",nch,x,y);
1852               // Fill the histograms
1853               Nh1+=nhits;
1854               h->Fill(x,y,(float) 1);
1855                   //}
1856               //}
1857            }
1858
1859            Int_t ncerenkovs = pRICH->Cerenkovs()->GetEntriesFast();
1860            //if (current->GetPdgCode() < 50000051 && current->GetPdgCode() > 50000040)
1861            //totalphotonsevent->Fill(ncerenkovs,(float) 1);
1862
1863            if (ncerenkovs) {
1864              printf("Cerenkovs       : %d\n",ncerenkovs);
1865              totalphotonsevent->Fill(ncerenkovs,(float) 1);
1866              for (Int_t hit=0;hit<ncerenkovs;hit++) {
1867                AliRICHCerenkov* cHit = (AliRICHCerenkov*) pRICH->Cerenkovs()->UncheckedAt(hit);
1868                Int_t nchamber = cHit->fChamber;     // chamber number
1869                Int_t index =    cHit->Track();
1870                //Int_t pindex =   (Int_t)(cHit->fIndex);
1871                trackglob[0] = cHit->X();                 // x-pos of hit
1872                trackglob[1] = cHit->Y();
1873                trackglob[2] = cHit->Z();                 // y-pos of hit
1874                //Float_t cx  =      cHit->X();                // x-position
1875                //Float_t cy  =      cHit->Z();                // y-position
1876                Int_t cmother =  cHit->fCMother;      // Index of mother particle
1877                Int_t closs =    (Int_t)(cHit->fLoss);           // How did the particle get lost?
1878                Float_t cherenkov = cHit->fCerenkovAngle;   //production cerenkov angle
1879
1880                chamber = &(pRICH->Chamber(nchamber-1));
1881
1882                //printf("Nch:%d\n",nch);
1883
1884                chamber->GlobaltoLocal(trackglob,trackloc);
1885
1886                chamber->LocaltoGlobal(trackloc,trackglob);
1887
1888
1889                Float_t cx=trackloc[0];
1890                Float_t cy=trackloc[2];
1891
1892                //printf ("Cerenkov hit number %d/%d, X:%f, Y:%f\n",hit,ncerenkovs,cx,cy);
1893
1894
1895                //printf("Particle:%9d\n",index);
1896
1897                TParticle *current = (TParticle*)gAlice->Particle(index);
1898                Float_t energyckov = current->Energy();
1899
1900                if (current->GetPdgCode() == 50000051)
1901                  {
1902                    if (closs==4)
1903                      {
1904                        feedback->Fill(cx,cy,(float) 1);
1905                        feed++;
1906                      }
1907                  }
1908                if (current->GetPdgCode() == 50000050)
1909                  {
1910
1911                    if (closs !=4)
1912                      {
1913                        phspectra2->Fill(energyckov*1e9,(float) 1);
1914                      }
1915
1916                    if (closs==4)
1917                      {
1918                        cerenkov->Fill(cx,cy,(float) 1);
1919
1920                        //printf ("Cerenkov hit number %d/%d, X:%d, Y:%d\n",hit,ncerenkovs,cx,cy);
1921
1922                        //TParticle *MIP = (TParticle*)gAlice->Particle(cmother);
1923                        AliRICHhit* mipHit = (AliRICHhit*) pRICH->Hits()->UncheckedAt(0);
1924                        mom[0] = current->Px();
1925                        mom[1] = current->Py();
1926                        mom[2] = current->Pz();
1927                        //mom[0] = cHit->fMomX;
1928                        // mom[1] = cHit->fMomZ;
1929                        //mom[2] = cHit->fMomY;
1930                        //Float_t energymip = MIP->Energy();
1931                        //Float_t Mip_px = mipHit->fMomFreoX;
1932                        //Float_t Mip_py = mipHit->fMomFreoY;
1933                        //Float_t Mip_pz = mipHit->fMomFreoZ;
1934                        //Float_t Mip_px = MIP->Px();
1935                        //Float_t Mip_py = MIP->Py();
1936                        //Float_t Mip_pz = MIP->Pz();
1937
1938
1939
1940                        //Float_t r = mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2];
1941                        //Float_t rt = TMath::Sqrt(r);
1942                        //Float_t Mip_r = Mip_px*Mip_px + Mip_py*Mip_py + Mip_pz*Mip_pz;
1943                        //Float_t Mip_rt = TMath::Sqrt(Mip_r);
1944                        //Float_t coscerenkov = (mom[0]*Mip_px + mom[1]*Mip_py + mom[2]*Mip_pz)/(rt*Mip_rt+0.0000001);
1945                        //Float_t cherenkov = TMath::ACos(coscerenkov);
1946                        ckovangle->Fill(cherenkov,(float) 1);                           //Cerenkov angle calculus
1947                        //printf("Cherenkov: %f\n",cherenkov);
1948                        Float_t ckphi=TMath::ATan2(mom[0], mom[2]);
1949                        hckphi->Fill(ckphi,(float) 1);
1950
1951
1952                        //Float_t mix = MIP->Vx();
1953                        //Float_t miy = MIP->Vy();
1954                        Float_t mx = mipHit->X();
1955                        Float_t my = mipHit->Z();
1956                        //printf("FX %e, FY %e, VX %e, VY %e\n",cx,cy,mx,my);
1957                        Float_t dx = trackglob[0] - mx;
1958                        Float_t dy = trackglob[2] - my;
1959                        //printf("Dx:%f, Dy:%f\n",dx,dy);
1960                        Float_t final_radius = TMath::Sqrt(dx*dx+dy*dy);
1961                        //printf("Final radius:%f\n",final_radius);
1962                        radius->Fill(final_radius,(float) 1);
1963
1964                        phspectra1->Fill(energyckov*1e9,(float) 1);
1965                        phot++;
1966                      }
1967                    for (Int_t nmothers=0;nmothers<=ntracks;nmothers++){
1968                      if (cmother == nmothers){
1969                        if (closs == 4)
1970                          mothers2[cmother]++;
1971                        mothers[cmother]++;
1972                      }
1973                    }
1974                  }
1975              }
1976            }
1977
1978
1979            if(gAlice->TreeR())
1980              {
1981                Int_t nent=(Int_t)gAlice->TreeR()->GetEntries();
1982                gAlice->TreeR()->GetEvent(nent-1);
1983                TClonesArray *Rawclusters = pRICH->RawClustAddress(2);    //  Raw clusters branch
1984                //printf ("Rawclusters:%p",Rawclusters);
1985                Int_t nrawclusters = Rawclusters->GetEntriesFast();
1986
1987                if (nrawclusters) {
1988                  printf("Raw Clusters    : %d\n",nrawclusters);
1989                  for (Int_t hit=0;hit<nrawclusters;hit++) {
1990                    AliRICHRawCluster* rcHit = (AliRICHRawCluster*) pRICH->RawClustAddress(2)->UncheckedAt(hit);
1991                    //Int_t nchamber = rcHit->fChamber;     // chamber number
1992                    //Int_t nhit = cHit->fHitNumber;        // hit number
1993                    Int_t qtot = rcHit->fQ;                 // charge
1994                    Float_t fx  =  rcHit->fX;                 // x-position
1995                    Float_t fy  =  rcHit->fY;                 // y-position
1996                    //Int_t type = rcHit->fCtype;             // cluster type ?
1997                    Int_t mult = rcHit->fMultiplicity;      // How many pads form the cluster
1998                    pads += mult;
1999                    if (qtot > 0) {
2000                      //printf ("fx: %d, fy: %d\n",fx,fy);
2001                      if (fx>(x-4) && fx<(x+4)  && fy>(y-4) && fy<(y+4)) {
2002                        //printf("There %d \n",mult);
2003                        padmip+=mult;
2004                      } else {
2005                        padnumber->Fill(mult,(float) 1);
2006                        nraw++;
2007                        if (mult<4) Clcharge->Fill(qtot,(float) 1);
2008                      }
2009
2010                    }
2011                  }
2012                }
2013
2014
2015                TClonesArray *RecHits1D = pRICH->RecHitsAddress1D(2);
2016                Int_t nrechits1D = RecHits1D->GetEntriesFast();
2017                //printf (" nrechits:%d\n",nrechits);
2018
2019                if(nrechits1D)
2020                  {
2021                    for (Int_t hit=0;hit<nrechits1D;hit++) {
2022                      AliRICHRecHit1D* recHit1D = (AliRICHRecHit1D*) pRICH->RecHitsAddress1D(2)->UncheckedAt(hit);
2023                      Float_t r_omega = recHit1D->fOmega;                  // Cerenkov angle
2024                      Float_t *cer_pho = recHit1D->fCerPerPhoton;        // Cerenkov angle per photon
2025                      Int_t *padsx = recHit1D->fPadsUsedX;           // Pads Used fo reconstruction (x)
2026                      Int_t *padsy = recHit1D->fPadsUsedY;           // Pads Used fo reconstruction (y)
2027                      Int_t goodPhotons = recHit1D->fGoodPhotons;    // Number of pads used for reconstruction
2028
2029                      Omega1D->Fill(r_omega,(float) 1);
2030
2031                      for (Int_t i=0; i<goodPhotons; i++)
2032                        {
2033                          PhotonCer->Fill(cer_pho[i],(float) 1);
2034                          PadsUsed->Fill(padsx[i],padsy[i],1);
2035                          //printf("Angle:%f, pad: %d %d\n",cer_pho[i],padsx[i],padsy[i]);
2036                        }
2037
2038                      //printf("Omega: %f, Theta: %f, Phi: %f\n",r_omega,r_theta,r_phi);
2039                    }
2040                  }
2041
2042
2043                TClonesArray *RecHits3D = pRICH->RecHitsAddress3D(2);
2044                Int_t nrechits3D = RecHits3D->GetEntriesFast();
2045                //printf (" nrechits:%d\n",nrechits);
2046
2047                if(nrechits3D)
2048                  {
2049                    recEffEvent = 0;
2050
2051                    //for (Int_t hit=0;hit<nrechits3D;hit++) {
2052                    AliRICHRecHit3D* recHit3D = (AliRICHRecHit3D*) pRICH->RecHitsAddress3D(2)->UncheckedAt(track);
2053                    Float_t r_omega    = recHit3D->fOmega;                  // Cerenkov angle
2054                    Float_t r_theta    = recHit3D->fTheta;                  // Theta angle of incidence
2055                    Float_t r_phi      = recHit3D->fPhi;                    // Phi angle if incidence
2056                    Float_t meanradius = recHit3D->fMeanRadius;              // Mean radius for reconstructed point
2057                    Float_t originalOmega = recHit3D->fOriginalOmega;       // Real Cerenkov angle
2058                    Float_t originalTheta = recHit3D->fOriginalTheta;       // Real incidence angle
2059                    Float_t originalPhi = recHit3D->fOriginalPhi;           // Real azimuthal angle
2060
2061
2062                    //correction to track cerenkov angle
2063                    originalOmega = (Float_t) ckovangle->GetMean();
2064
2065                    if(diaglevel == 4)
2066                      {
2067                        printf("\nMean cerenkov angle: %f\n", originalOmega);
2068                        printf("Reconstructed cerenkov angle: %f\n",r_omega);
2069                      }
2070
2071                    Float_t omegaError = TMath::Abs(originalOmega - r_omega);
2072                    Float_t thetaError = TMath::Abs(originalTheta - r_theta);
2073                    Float_t phiError   = TMath::Abs(originalPhi - r_phi);
2074
2075
2076                    if(TMath::Abs(omegaError) < 0.015)
2077                      recEffEvent += 1;
2078
2079                    Omega3D->Fill(r_omega,(float) 1);
2080                    Theta->Fill(r_theta*180/TMath::Pi(),(float) 1);
2081                    Phi->Fill(r_phi*180/TMath::Pi()-180,(float) 1);
2082                    MeanRadius->Fill(meanradius,(float) 1);
2083                    identification->Fill(PTfinal, r_omega,1);
2084                    OriginalOmega->Fill(originalOmega, (float) 1);
2085                    OriginalTheta->Fill(originalTheta, (float) 1);
2086                    OriginalPhi->Fill(TMath::Abs(originalPhi), (float) 1);
2087                    OmegaError->Fill(omegaError, (float) 1);
2088                    ThetaError->Fill(thetaError, (float) 1);
2089                    PhiError->Fill(phiError, (float) 1);
2090
2091                    recEffEvent = recEffEvent;
2092                    recEffTotal += recEffEvent;
2093
2094                    Float_t pioncer = acos(sqrt((.139*.139+PTfinal*PTfinal)/(PTfinal*PTfinal*1.285*1.285)));
2095                    Float_t kaoncer = acos(sqrt((.439*.439+PTfinal*PTfinal)/(PTfinal*PTfinal*1.285*1.285)));
2096                    Float_t protoncer = acos(sqrt((.938*.938+PTfinal*PTfinal)/(PTfinal*PTfinal*1.285*1.285)));
2097
2098                    Float_t piondist = TMath::Abs(r_omega - pioncer);
2099                    Float_t kaondist = TMath::Abs(r_omega - kaoncer);
2100                    Float_t protondist = TMath::Abs(r_omega - protoncer);
2101
2102                    if(diaglevel == 4)
2103                      {
2104                        if(pioncer<r_omega)
2105                          {
2106                            printf("Identified as a PION!\n");
2107                            pionCount += 1;
2108                          }
2109                        if(kaoncer<r_omega && pioncer>r_omega)
2110                          {
2111                            if(kaondist>piondist)
2112                              {
2113                                printf("Identified as a PION!\n");
2114                                pionCount += 1;
2115                              }
2116                            else
2117                              {
2118                                printf("Identified as a KAON!\n");
2119                                kaonCount += 1;
2120                              }
2121                          }                       }
2122                        if(protoncer<r_omega && kaoncer>r_omega)
2123                          {
2124                            if(kaondist>protondist)
2125                              {
2126                                printf("Identified as a PROTON!\n");
2127                                protonCount += 1;
2128                              }
2129                            else
2130                              {
2131                                printf("Identified as a KAON!\n");
2132                                pionCount += 1;
2133                              }
2134                          }
2135                        if(protoncer>r_omega)
2136                          {
2137                            printf("Identified as a PROTON!\n");
2138                            protonCount += 1;
2139                          }
2140
2141                        printf("\nReconstruction efficiency: %5.2f%%\n", recEffEvent*100);
2142                  }
2143              }
2144        }
2145
2146
2147        for (Int_t nmothers=0;nmothers<ntracks;nmothers++){
2148          totalphotonstrack->Fill(mothers[nmothers],(float) 1);
2149          mother->Fill(mothers2[nmothers],(float) 1);
2150        }
2151
2152        clusev->Fill(nraw,(float) 1);
2153        photev->Fill(phot,(float) 1);
2154        feedev->Fill(feed,(float) 1);
2155        padsmip->Fill(padmip,(float) 1);
2156        padscl->Fill(pads,(float) 1);
2157        phot = 0;
2158        feed = 0;
2159        pads = 0;
2160        nraw=0;
2161        padmip=0;
2162
2163
2164
2165        gAlice->ResetDigits();
2166        gAlice->TreeD()->GetEvent(0);
2167
2168        if (diaglevel < 4)
2169          {
2170
2171
2172            TClonesArray *Digits  = pRICH->DigitsAddress(2);
2173            Int_t ndigits = Digits->GetEntriesFast();
2174            printf("Digits          : %d\n",ndigits);
2175            padsev->Fill(ndigits,(float) 1);
2176            for (Int_t hit=0;hit<ndigits;hit++) {
2177              AliRICHDigit* dHit = (AliRICHDigit*) Digits->UncheckedAt(hit);
2178              Int_t qtot = dHit->Signal();                // charge
2179              Int_t ipx  = dHit->PadX();               // pad number on X
2180              Int_t ipy  = dHit->PadY();               // pad number on Y
2181              //printf("%d, %d\n",ipx,ipy);
2182              if( ipx<=100 && ipy <=100) hc0->Fill(ipx,ipy,(float) qtot);
2183            }
2184          }
2185
2186        if (diaglevel == 5)
2187          {
2188            for (Int_t ich=0;ich<7;ich++)
2189              {
2190                TClonesArray *Digits = pRICH->DigitsAddress(ich);    //  Raw clusters branch
2191                Int_t ndigits = Digits->GetEntriesFast();
2192                //printf("Digits:%d\n",ndigits);
2193                padsev->Fill(ndigits,(float) 1);
2194                if (ndigits) {
2195                  for (Int_t hit=0;hit<ndigits;hit++) {
2196                    AliRICHDigit* dHit = (AliRICHDigit*) Digits->UncheckedAt(hit);
2197                    Int_t qtot = dHit->Signal();                // charge
2198                    Int_t ipx  = dHit->PadX();               // pad number on X
2199                    Int_t ipy  = dHit->PadY();               // pad number on Y
2200                    if( ipx<=100 && ipy <=100 && ich==2) hc0->Fill(ipx,ipy,(float) qtot);
2201                    if( ipx<=162 && ipy <=162 && ich==0) hc1->Fill(ipx,ipy,(float) qtot);
2202                    if( ipx<=162 && ipy <=162 && ich==1) hc2->Fill(ipx,ipy,(float) qtot);
2203                    if( ipx<=162 && ipy <=162 && ich==2) hc3->Fill(ipx,ipy,(float) qtot);
2204                    if( ipx<=162 && ipy <=162 && ich==3) hc4->Fill(ipx,ipy,(float) qtot);
2205                    if( ipx<=162 && ipy <=162 && ich==4) hc5->Fill(ipx,ipy,(float) qtot);
2206                    if( ipx<=162 && ipy <=162 && ich==5) hc6->Fill(ipx,ipy,(float) qtot);
2207                    if( ipx<=162 && ipy <=162 && ich==6) hc7->Fill(ipx,ipy,(float) qtot);
2208                  }
2209                }
2210              }
2211          }
2212    }
2213
2214    if(diaglevel == 4)
2215      {
2216
2217        Stat_t omegaE;
2218        Stat_t thetaE;
2219        Stat_t phiE;
2220
2221        Stat_t omegaO;
2222        Stat_t thetaO;
2223        Stat_t phiO;
2224
2225        for(Int_t i=0;i<99;i++)
2226          {
2227            omegaE = OriginalOmega->GetBinContent(i);
2228            if(omegaE != 0)
2229              {
2230                omegaO = Omega3D->GetBinContent(i);
2231                chiSquareOmega += (TMath::Power(omegaE,2) - TMath::Power(omegaO,2))/omegaO;
2232              }
2233
2234            thetaE = OriginalTheta->GetBinContent(i);
2235            if(thetaE != 0)
2236              {
2237                thetaO = Theta->GetBinContent(i);
2238                chiSquareTheta += (TMath::Power(thetaE,2) - TMath::Power(thetaO,2))/thetaO;
2239              }
2240
2241            phiE = OriginalPhi->GetBinContent(i);
2242            if(phiE != 0)
2243              {
2244                phiO = Phi->GetBinContent(i);
2245                chiSquarePhi += (TMath::Power(phiE,2) - TMath::Power(phiO,2))/phiO;
2246              }
2247          }
2248
2249
2250
2251        printf("\nChi square test values:   Omega - %f\n", chiSquareOmega);
2252        printf("                          Theta - %f\n", chiSquareTheta);
2253        printf("                          Phi   - %f\n", chiSquarePhi);
2254
2255        printf("\nKolmogorov test values:   Omega - %5.4f\n", Omega3D->KolmogorovTest(OriginalOmega));
2256        printf("                          Theta - %5.4f\n", Theta->KolmogorovTest(OriginalTheta));
2257        printf("                          Phi   - %5.4f\n", Phi->KolmogorovTest(OriginalPhi));
2258
2259        recEffTotal = recEffTotal/evNumber2;
2260        printf("\nTotal reconstruction efficiency: %5.2f%%\n", recEffTotal*100);
2261        printf("\n Pions: %d\n Kaons: %d\n Protons:%d\n",pionCount, kaonCount, protonCount);
2262
2263      }
2264
2265
2266    //Create canvases, set the view range, show histograms
2267
2268    TCanvas *c1 = 0;
2269    TCanvas *c2 = 0;
2270    TCanvas *c3 = 0;
2271    TCanvas *c4 = 0;
2272    TCanvas *c5 = 0;
2273    TCanvas *c6 = 0;
2274    TCanvas *c7 = 0;
2275    TCanvas *c8 = 0;
2276    TCanvas *c9 = 0;
2277    TCanvas *c10 = 0;
2278    TCanvas *c11 = 0;
2279    TCanvas *c12 = 0;
2280    TCanvas *c13 = 0;
2281
2282
2283    TStyle *mystyle=new TStyle("Plain","mystyle");
2284    mystyle->SetPalette(1,0);
2285    mystyle->SetFuncColor(2);
2286    mystyle->SetDrawBorder(0);
2287    mystyle->SetTitleBorderSize(0);
2288    mystyle->SetOptFit(1111);
2289    mystyle->cd();
2290
2291
2292    TClonesArray *RecHits3D = pRICH->RecHitsAddress3D(2);
2293    Int_t nrechits3D = RecHits3D->GetEntriesFast();
2294    TClonesArray *RecHits1D = pRICH->RecHitsAddress1D(2);
2295    Int_t nrechits1D = RecHits1D->GetEntriesFast();
2296
2297   switch(diaglevel)
2298      {
2299      case 1:
2300
2301        c1 = new TCanvas("c1","Alice RICH digits",50,50,300,350);
2302        hc0->SetXTitle("ix (npads)");
2303        hc0->Draw("colz");
2304
2305        c2 = new TCanvas("c2","Hits per type",100,100,600,700);
2306        c2->Divide(2,2);
2307        //c4->SetFillColor(42);
2308
2309        c2->cd(1);
2310        feedback->SetXTitle("x (cm)");
2311        feedback->SetYTitle("y (cm)");
2312        feedback->Draw("colz");
2313
2314        c2->cd(2);
2315        //mip->SetFillColor(5);
2316        mip->SetXTitle("x (cm)");
2317        mip->SetYTitle("y (cm)");
2318        mip->Draw("colz");
2319
2320        c2->cd(3);
2321        //cerenkov->SetFillColor(5);
2322        cerenkov->SetXTitle("x (cm)");
2323        cerenkov->SetYTitle("y (cm)");
2324        cerenkov->Draw("colz");
2325
2326        c2->cd(4);
2327        //h->SetFillColor(5);
2328        h->SetXTitle("x (cm)");
2329        h->SetYTitle("y (cm)");
2330        h->Draw("colz");
2331
2332        c3 = new TCanvas("c3","Hits distribution",150,150,600,350);
2333        c3->Divide(2,1);
2334        //c10->SetFillColor(42);
2335
2336        c3->cd(1);
2337        hitsX->SetFillColor(5);
2338        hitsX->SetXTitle("(cm)");
2339        hitsX->Draw();
2340
2341        c3->cd(2);
2342        hitsY->SetFillColor(5);
2343        hitsY->SetXTitle("(cm)");
2344        hitsY->Draw();
2345
2346
2347        break;
2348      case 2:
2349
2350        c4 = new TCanvas("c4","Photon Spectra",50,50,600,350);
2351        c4->Divide(2,1);
2352
2353        c4->cd(1);
2354        phspectra2->SetFillColor(5);
2355        phspectra2->SetXTitle("energy (eV)");
2356        phspectra2->Draw();
2357        c4->cd(2);
2358        phspectra1->SetFillColor(5);
2359        phspectra1->SetXTitle("energy (eV)");
2360        phspectra1->Draw();
2361
2362        c5 = new TCanvas("c5","Particles Spectra",100,100,600,700);
2363        c5->Divide(2,2);
2364
2365        c5->cd(1);
2366        pionspectra->SetFillColor(5);
2367        pionspectra->SetXTitle("(GeV)");
2368        pionspectra->Draw();
2369
2370        c5->cd(2);
2371        protonspectra->SetFillColor(5);
2372        protonspectra->SetXTitle("(GeV)");
2373        protonspectra->Draw();
2374
2375        c5->cd(3);
2376        kaonspectra->SetFillColor(5);
2377        kaonspectra->SetXTitle("(GeV)");
2378        kaonspectra->Draw();
2379
2380        c5->cd(4);
2381        chargedspectra->SetFillColor(5);
2382        chargedspectra->SetXTitle("(GeV)");
2383        chargedspectra->Draw();
2384
2385        break;
2386
2387      case 3:
2388
2389
2390        if(gAlice->TreeR())
2391          {
2392            c6=new TCanvas("c6","Clusters Statistics",50,50,600,700);
2393            c6->Divide(2,2);
2394
2395            c6->cd(1);
2396            Clcharge->SetFillColor(5);
2397            Clcharge->SetXTitle("ADC counts");
2398            if (evNumber2>10)
2399              {
2400                Clcharge->Fit("expo");
2401              }
2402            Clcharge->Draw();
2403
2404            c6->cd(2);
2405            padnumber->SetFillColor(5);
2406            padnumber->SetXTitle("(counts)");
2407            padnumber->Draw();
2408
2409            c6->cd(3);
2410            clusev->SetFillColor(5);
2411            clusev->SetXTitle("(counts)");
2412            if (evNumber2>10)
2413              {
2414                clusev->Fit("gaus");
2415                //gaus->SetLineColor(2);
2416                //gaus->SetLineWidth(3);
2417              }
2418            clusev->Draw();
2419
2420            c6->cd(4);
2421            padsmip->SetFillColor(5);
2422            padsmip->SetXTitle("(counts)");
2423            padsmip->Draw();
2424          }
2425
2426        if(evNumber2<1)
2427          {
2428            c11 = new TCanvas("c11","Cherenkov per Mip",400,10,600,700);
2429            mother->SetFillColor(5);
2430            mother->SetXTitle("counts");
2431            mother->Draw();
2432          }
2433
2434        c7 = new TCanvas("c7","Production Statistics",100,100,600,700);
2435        c7->Divide(2,2);
2436        //c7->SetFillColor(42);
2437
2438        c7->cd(1);
2439        totalphotonsevent->SetFillColor(5);
2440        totalphotonsevent->SetXTitle("Photons (counts)");
2441        if (evNumber2>10)
2442            {
2443              totalphotonsevent->Fit("gaus");
2444              //gaus->SetLineColor(2);
2445              //gaus->SetLineWidth(3);
2446            }
2447        totalphotonsevent->Draw();
2448
2449        c7->cd(2);
2450        photev->SetFillColor(5);
2451        photev->SetXTitle("(counts)");
2452        if (evNumber2>10)
2453          {
2454            photev->Fit("gaus");
2455            //gaus->SetLineColor(2);
2456            //gaus->SetLineWidth(3);
2457          }
2458        photev->Draw();
2459
2460        c7->cd(3);
2461        feedev->SetFillColor(5);
2462        feedev->SetXTitle("(counts)");
2463        if (evNumber2>10)
2464          {
2465            feedev->Fit("gaus");
2466          }
2467        feedev->Draw();
2468
2469        c7->cd(4);
2470        padsev->SetFillColor(5);
2471        padsev->SetXTitle("(counts)");
2472        if (evNumber2>10)
2473          {
2474            padsev->Fit("gaus");
2475          }
2476        padsev->Draw();
2477
2478        break;
2479
2480      case 4:
2481
2482
2483        if(nrechits3D)
2484          {
2485            c8 = new TCanvas("c8","3D reconstruction of Phi angle",50,50,300,1050);
2486            c8->Divide(1,3);
2487            //c2->SetFillColor(42);
2488
2489
2490            // data per hit
2491            c8->cd(1);
2492            hitsPhi->SetFillColor(5);
2493            if (evNumber2>10)
2494              hitsPhi->Fit("gaus");
2495            hitsPhi->Draw();
2496
2497             //data per track
2498            c8->cd(2);
2499            OriginalPhi->SetFillColor(5);
2500            if (evNumber2>10)
2501              OriginalPhi->Fit("gaus");
2502            OriginalPhi->Draw();
2503
2504            //recontructed data
2505            c8->cd(3);
2506            Phi->SetFillColor(5);
2507            if (evNumber2>10)
2508              Phi->Fit("gaus");
2509            Phi->Draw();
2510
2511            c9 = new TCanvas("c9","3D reconstruction of theta angle",75,75,300,1050);
2512            c9->Divide(1,3);
2513
2514            // data per hit
2515            c9->cd(1);
2516            hitsTheta->SetFillColor(5);
2517            if (evNumber2>10)
2518              hitsTheta->Fit("gaus");
2519            hitsTheta->Draw();
2520
2521            //data per track
2522            c9->cd(2);
2523            OriginalTheta->SetFillColor(5);
2524            if (evNumber2>10)
2525              OriginalTheta->Fit("gaus");
2526            OriginalTheta->Draw();
2527
2528            //recontructed data
2529            c9->cd(3);
2530            Theta->SetFillColor(5);
2531            if (evNumber2>10)
2532              Theta->Fit("gaus");
2533            Theta->Draw();
2534
2535            c10 = new TCanvas("c10","3D reconstruction of cherenkov angle",100,100,300,1050);
2536            c10->Divide(1,3);
2537
2538            // data per hit
2539            c10->cd(1);
2540            ckovangle->SetFillColor(5);
2541            ckovangle->SetXTitle("angle (radians)");
2542            if (evNumber2>10)
2543              ckovangle->Fit("gaus");
2544            ckovangle->Draw();
2545
2546            //data per track
2547            c10->cd(2);
2548            OriginalOmega->SetFillColor(5);
2549            OriginalOmega->SetXTitle("angle (radians)");
2550            if (evNumber2>10)
2551              OriginalOmega->Fit("gaus");
2552            OriginalOmega->Draw();
2553
2554            //recontructed data
2555            c10->cd(3);
2556            Omega3D->SetFillColor(5);
2557            Omega3D->SetXTitle("angle (radians)");
2558            if (evNumber2>10)
2559              Omega3D->Fit("gaus");
2560            Omega3D->Draw();
2561
2562
2563            c11 = new TCanvas("c11","3D reconstruction of mean radius",125,125,300,700);
2564            c11->Divide(1,2);
2565
2566            // data per hit
2567            c11->cd(1);
2568            radius->SetFillColor(5);
2569            radius->SetXTitle("radius (cm)");
2570            radius->Draw();
2571
2572            //recontructed data
2573            c11->cd(2);
2574            MeanRadius->SetFillColor(5);
2575            MeanRadius->SetXTitle("radius (cm)");
2576            MeanRadius->Draw();
2577
2578
2579            c12 = new TCanvas("c12","Cerenkov angle vs. Momentum",150,150,550,350);
2580
2581            c12->cd(1);
2582            identification->SetFillColor(5);
2583            identification->SetXTitle("Momentum (GeV/c)");
2584            identification->SetYTitle("Cherenkov angle (radians)");
2585
2586            TF1 *pionplot = new TF1("pion","acos(sqrt((.139*.139+x*x)/(x*x*1.285*1.285)))",1,5);
2587            TF1 *kaonplot = new TF1("kaon","acos(sqrt((.439*.439+x*x)/(x*x*1.285*1.285)))",1,5);
2588            TF1 *protonplot = new TF1("proton","acos(sqrt((.938*.938+x*x)/(x*x*1.285*1.285)))",1,5);
2589
2590            identification->Draw();
2591
2592            pionplot->SetLineColor(5);
2593            pionplot->Draw("same");
2594
2595            kaonplot->SetLineColor(4);
2596            kaonplot->Draw("same");
2597
2598            protonplot->SetLineColor(3);
2599            protonplot->Draw("same");
2600
2601            c13 = new TCanvas("c13","Reconstruction Errors",200,200,900,350);
2602            c13->Divide(3,1);
2603
2604            c13->cd(1);
2605            PhiError->SetFillColor(5);
2606            if (evNumber2>10)
2607              PhiError->Fit("gaus");
2608            PhiError->Draw();
2609            c13->cd(2);
2610            ThetaError->SetFillColor(5);
2611            if (evNumber2>10)
2612              ThetaError->Fit("gaus");
2613            ThetaError->Draw();
2614            c13->cd(3);
2615            OmegaError->SetFillColor(5);
2616            OmegaError->SetXTitle("angle (radians)");
2617            if (evNumber2>10)
2618              OmegaError->Fit("gaus");
2619            OmegaError->Draw();
2620
2621          }
2622
2623        if(nrechits1D)
2624          {
2625            c9 = new TCanvas("c9","1D Reconstruction",100,100,1100,700);
2626            c9->Divide(3,2);
2627            //c5->SetFillColor(42);
2628
2629            c9->cd(1);
2630            ckovangle->SetFillColor(5);
2631            ckovangle->SetXTitle("angle (radians)");
2632            ckovangle->Draw();
2633
2634            c9->cd(2);
2635            radius->SetFillColor(5);
2636            radius->SetXTitle("radius (cm)");
2637            radius->Draw();
2638
2639            c9->cd(3);
2640            hc0->SetXTitle("pads");
2641            hc0->Draw("box");
2642
2643            c9->cd(5);
2644            Omega1D->SetFillColor(5);
2645            Omega1D->SetXTitle("angle (radians)");
2646            Omega1D->Draw();
2647
2648            c9->cd(4);
2649            PhotonCer->SetFillColor(5);
2650            PhotonCer->SetXTitle("angle (radians)");
2651            PhotonCer->Draw();
2652
2653            c9->cd(6);
2654            PadsUsed->SetXTitle("pads");
2655            PadsUsed->Draw("box");
2656          }
2657
2658        break;
2659
2660      case 5:
2661
2662        printf("Drawing histograms.../n");
2663
2664        c10 = new TCanvas("c10","Alice RICH digits",50,50,1200,700);
2665        c1->Divide(4,2);
2666
2667        c10->cd(1);
2668        hc1->SetXTitle("ix (npads)");
2669        hc1->Draw("box");
2670        c10->cd(2);
2671        hc2->SetXTitle("ix (npads)");
2672        hc2->Draw("box");
2673        c10->cd(3);
2674        hc3->SetXTitle("ix (npads)");
2675        hc3->Draw("box");
2676        c10->cd(4);
2677        hc4->SetXTitle("ix (npads)");
2678        hc4->Draw("box");
2679        c10->cd(5);
2680        hc5->SetXTitle("ix (npads)");
2681        hc5->Draw("box");
2682        c10->cd(6);
2683        hc6->SetXTitle("ix (npads)");
2684        hc6->Draw("box");
2685        c10->cd(7);
2686        hc7->SetXTitle("ix (npads)");
2687        hc7->Draw("box");
2688        c10->cd(8);
2689        hc0->SetXTitle("ix (npads)");
2690        hc0->Draw("box");
2691        c11 = new TCanvas("c11","Hits per type",100,100,600,700);
2692        c11->Divide(2,2);
2693
2694        c11->cd(1);
2695        feedback->SetXTitle("x (cm)");
2696        feedback->SetYTitle("y (cm)");
2697        feedback->Draw();
2698
2699        c11->cd(2);
2700        mip->SetXTitle("x (cm)");
2701        mip->SetYTitle("y (cm)");
2702        mip->Draw();
2703
2704        c11->cd(3);
2705        cerenkov->SetXTitle("x (cm)");
2706        cerenkov->SetYTitle("y (cm)");
2707        cerenkov->Draw();
2708
2709        c11->cd(4);
2710        h->SetXTitle("x (cm)");
2711        h->SetYTitle("y (cm)");
2712        h->Draw();
2713
2714        c12 = new TCanvas("c12","Hits distribution",150,150,600,350);
2715        c12->Divide(2,1);
2716
2717        c12->cd(1);
2718        hitsX->SetFillColor(5);
2719        hitsX->SetXTitle("(cm)");
2720        hitsX->Draw();
2721
2722        c12->cd(2);
2723        hitsY->SetFillColor(5);
2724        hitsY->SetXTitle("(cm)");
2725        hitsY->Draw();
2726
2727        break;
2728
2729      }
2730
2731
2732    printf("\nEnd of analysis\n");
2733    printf("**********************************\n");
2734 }//void AliRICHv3::DiagnosticsSE(Int_t diaglevel,Int_t evNumber1,Int_t evNumber2)
2735