]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - ZDC/AliZDCv3.cxx
Correct energy calibration in case of p-p data
[u/mrichter/AliRoot.git] / ZDC / AliZDCv3.cxx
index 44a590eea2e761d69ed0b06a062227b4e9b3d477..2565e0e58c625be39631d707c4de183f92f76494 100644 (file)
@@ -38,6 +38,7 @@
 #include <TParticle.h>
 
 // --- AliRoot classes
+#include "AliLog.h"
 #include "AliConst.h"
 #include "AliMagF.h"
 #include "AliRun.h"
@@ -71,14 +72,11 @@ AliZDCv3::AliZDCv3() :
   fpLostD1C(0), 
   fpDetectedC(0),
   fnDetectedC(0),
-  fnLumiC(0),
   fpLostITA(0), 
   fpLostD1A(0), 
   fpLostTDI(0), 
   fpDetectedA(0),
-  fnDetectedA(0),
-  fnLumiA(0),
-  fnTrou(0)
+  fnDetectedA(0)
 {
   //
   // Default constructor for Zero Degree Calorimeter
@@ -106,14 +104,11 @@ AliZDCv3::AliZDCv3(const char *name, const char *title) :
   fpLostD1C(0), 
   fpDetectedC(0),
   fnDetectedC(0),
-  fnLumiC(0),
   fpLostITA(0), 
   fpLostD1A(0), 
   fpLostTDI(0), 
   fpDetectedA(0),
-  fnDetectedA(0),  
-  fnLumiA(0),
-  fnTrou(0)
+  fnDetectedA(0)  
 {
   //
   // Standard constructor for Zero Degree Calorimeter 
@@ -629,7 +624,7 @@ void AliZDCv3::CreateBeamLine()
   boxpar[0] = 8.0/2.;
   boxpar[1] = 8.0/2.;
   boxpar[2] = 15./2.;
-  gMC->Gsvolu("QLUC", "BOX ", idtmed[6], boxpar, 3);
+  gMC->Gsvolu("QLUC", "BOX ", idtmed[9], boxpar, 3);
   gMC->Gspos("QLUC", 1, "ZDCC", 0., 0.,  fPosZNC[2]+66.+boxpar[2], 0, "ONLY");
   //printf("   QLUC LUMINOMETER from z = %1.2f to z= %1.2f\n",  fPosZNC[2]+66., fPosZNC[2]+66.+2*boxpar[2]);
                 
@@ -1345,7 +1340,7 @@ void AliZDCv3::CreateBeamLine()
   boxpar[0] = 8.0/2.;
   boxpar[1] = 8.0/2.;
   boxpar[2] = 15./2.;
-  gMC->Gsvolu("QLUA", "BOX ", idtmed[6], boxpar, 3);
+  gMC->Gsvolu("QLUA", "BOX ", idtmed[9], boxpar, 3);
   gMC->Gspos("QLUA", 1, "ZDCA", 0., 0.,  fPosZNA[2]-66.-boxpar[2], 0, "ONLY");
   //printf("   QLUA LUMINOMETER from z = %1.2f to z= %1.2f\n\n",  fPosZNA[2]-66., fPosZNA[2]-66.-2*boxpar[2]);
 
@@ -1990,6 +1985,10 @@ void AliZDCv3::CreateMaterials()
   // --- Copper (energy loss taken into account)
   ubuf[0] = 1.10;
   AliMaterial(6, "COPP0", 63.54, 29., 8.96, 1.4, 0., ubuf, 1);
+
+  // --- Copper 
+  ubuf[0] = 1.10;
+  AliMaterial(9, "COPP1", 63.54, 29., 8.96, 1.4, 0., ubuf, 1);
   
   // --- Iron (energy loss taken into account)
   ubuf[0] = 1.1;
@@ -2026,8 +2025,8 @@ void AliZDCv3::CreateMaterials()
   // --- Fibers (SiO2) = 3 ; 
   // --- Fibers (SiO2) = 4 ; 
   // --- Lead = 5 ; 
-  // --- Copper (with energy loss)= 6 ;
-  // --- Copper (with energy loss)= 13 ; 
+  // --- Copper (with high thr.)= 6 ;
+  // --- Copper (with low thr.)=  9;
   // --- Iron (with energy loss) = 7 ; 
   // --- Iron (without energy loss) = 8 ; 
   // --- Vacuum (no field) = 10 
@@ -2058,10 +2057,10 @@ void AliZDCv3::CreateMaterials()
   AliMedium(6, "ZCOPP", 6, isvol, inofld, nofieldm, tmaxfd, stemax, deemax, epsil, stmin);
   AliMedium(7, "ZIRON", 7, isvol, inofld, nofieldm, tmaxfd, stemax, deemax, epsil, stmin);
   AliMedium(8, "ZIRONN",8, isvol, inofld, nofieldm, tmaxfd, stemax, deemax, epsil, stmin);
+  AliMedium(9, "ZCOPL", 6, isvol, inofld, nofieldm, tmaxfd, stemax, deemax, epsil, stmin);
   AliMedium(10,"ZVOID",10, isvol, inofld, nofieldm, tmaxfd, stemax, deemax, epsil, stmin);
-  AliMedium(12,"ZAIR", 12, isvol, inofld, nofieldm, tmaxfd, stemax, deemax, epsil, stmin);
-  //
   AliMedium(11,"ZVOIM",11, isvol, ifield, fieldm, tmaxfd, stemax, deemax, epsil, stmin);
+  AliMedium(12,"ZAIR", 12, isvolActive, inofld, nofieldm, tmaxfd, stemax, deemax, epsil, stmin);
 
 
 } 
@@ -2106,7 +2105,7 @@ void AliZDCv3::Init()
   Int_t *idtmed = fIdtmed->GetArray();  
   Int_t i;
   // Thresholds for showering in the ZDCs 
-  i = 1; //tantalum
+  i = 1; //Wa lloy
   gMC->Gstpar(idtmed[i], "CUTGAM", .001);
   gMC->Gstpar(idtmed[i], "CUTELE", .001);
   gMC->Gstpar(idtmed[i], "CUTNEU", .01);
@@ -2129,6 +2128,13 @@ void AliZDCv3::Init()
   gMC->Gstpar(idtmed[i], "CUTNEU", 1.);
   gMC->Gstpar(idtmed[i], "CUTHAD", 1.);
   
+  // Thresholds for showering in the luminometer
+  i = 9; //copper
+  gMC->Gstpar(idtmed[i], "CUTGAM", .001);
+  gMC->Gstpar(idtmed[i], "CUTELE", .001);
+  gMC->Gstpar(idtmed[i], "CUTNEU", .01);
+  gMC->Gstpar(idtmed[i], "CUTHAD", .01);
+  
   // Avoid too detailed showering along the beam line 
   i = 7; //iron with energy loss (ZIRON)
   gMC->Gstpar(idtmed[i], "CUTGAM", .1);
@@ -2137,7 +2143,7 @@ void AliZDCv3::Init()
   gMC->Gstpar(idtmed[i], "CUTHAD", 1.);
   
   // Avoid too detailed showering along the beam line 
-  i = 8; //iron with energy loss (ZIRONN)
+  i = 8; //iron without energy loss (ZIRONN)
   gMC->Gstpar(idtmed[i], "CUTGAM", .1);
   gMC->Gstpar(idtmed[i], "CUTELE", .1);
   gMC->Gstpar(idtmed[i], "CUTNEU", 1.);
@@ -2232,12 +2238,17 @@ void AliZDCv3::InitTables()
      return;
   }
   
+  int read=1;
   for(k=0; k<fNalfan; k++){
      for(j=0; j<fNben; j++){
-       fscanf(fp1,"%f",&fTablen[0][k][j]);
-       fscanf(fp2,"%f",&fTablen[1][k][j]);
-       fscanf(fp3,"%f",&fTablen[2][k][j]);
-       fscanf(fp4,"%f",&fTablen[3][k][j]);
+       read = fscanf(fp1,"%f",&fTablen[0][k][j]);
+       if(read==0) AliDebug(3, " Error in reading light table 1");
+       read = fscanf(fp2,"%f",&fTablen[1][k][j]);
+       if(read==0) AliDebug(3, " Error in reading light table 2");
+       read = fscanf(fp3,"%f",&fTablen[2][k][j]);
+       if(read==0) AliDebug(3, " Error in reading light table 3");
+       read = fscanf(fp4,"%f",&fTablen[3][k][j]);
+       if(read==0) AliDebug(3, " Error in reading light table 4");
      } 
   }
   fclose(fp1);
@@ -2269,10 +2280,14 @@ void AliZDCv3::InitTables()
   
   for(k=0; k<fNalfap; k++){
      for(j=0; j<fNbep; j++){
-       fscanf(fp5,"%f",&fTablep[0][k][j]);
-       fscanf(fp6,"%f",&fTablep[1][k][j]);
-       fscanf(fp7,"%f",&fTablep[2][k][j]);
-       fscanf(fp8,"%f",&fTablep[3][k][j]);
+       read = fscanf(fp5,"%f",&fTablep[0][k][j]);
+       if(read==0) AliDebug(3, " Error in reading light table 5");
+       read = fscanf(fp6,"%f",&fTablep[1][k][j]);
+       if(read==0) AliDebug(3, " Error in reading light table 6");
+       read = fscanf(fp7,"%f",&fTablep[2][k][j]);
+       if(read==0) AliDebug(3, " Error in reading light table 7");
+       read = fscanf(fp8,"%f",&fTablep[3][k][j]);
+       if(read==0) AliDebug(3, " Error in reading light table 8");
      } 
   }
   fclose(fp5);
@@ -2287,15 +2302,12 @@ void AliZDCv3::StepManager()
   // Routine called at every step in the Zero Degree Calorimeters
   //
   Int_t   j, vol[2]={0,0}, ibeta=0, ialfa=0, ibe=0, nphe=0;
-  Float_t hits[12], x[3], xdet[3], um[3], ud[3];
-  Float_t m=0., ekin=0., destep=0., be=0., out=0.;
-  // Parametrization for light guide uniformity
-  // NEW!!! Light guide tilted @ 51 degrees
-  Float_t guiPar[4]={0.31,-0.0006305,0.01337,0.8895};
+  Float_t hits[13], x[3], xdet[3], um[3], ud[3];
+  Float_t destep=0., be=0., out=0.;
   Double_t s[3], p[3];
   const char *knamed;
   //
-  for(j=0;j<12;j++) hits[j]=-999.;
+  for(j=0;j<13;j++) hits[j]=-999.;
   //
   // --- This part is for no shower developement in beam pipe and TDI
   // If particle interacts with beam pipe or TDI -> return
@@ -2316,7 +2328,6 @@ void AliZDCv3::StepManager()
          else fpLostD1A += 1;
          ipr=1;
        }
-       else if(!strncmp(knamed,"QAL",3)) fnTrou++;
       }
       else if(gMC->CurrentMedium() == fMedSensTDI){ 
         knamed = gMC->CurrentVolName();
@@ -2326,31 +2337,23 @@ void AliZDCv3::StepManager()
          ipr=1;
         }
        else if(!strncmp(knamed,"QTD",3)) fpLostTDI += 1;
-       else if(!strncmp(knamed,"QLU",3)){
-         if(s[2]<0) fnLumiC ++;
-         else fnLumiA++;
-         ipr=1;
-       }
       }
       //
       //gMC->TrackMomentum(p[0], p[1], p[2], p[3]);
       //printf("\t Particle: mass = %1.3f, E = %1.3f GeV, pz = %1.2f GeV -> stopped in volume %s\n", 
       //     gMC->TrackMass(), p[3], p[2], gMC->CurrentVolName());
       //
-      if(ipr!=0){
+      /*if(ipr!=0){
         printf("\n\t **********************************\n");
         printf("\t ********** Side C **********\n");
         printf("\t # of spectators in IT = %d\n",fpLostITC);
         printf("\t # of spectators in D1 = %d\n",fpLostD1C);
-        printf("\t # of spectators in luminometer = %d\n",fnLumiC);
         printf("\t ********** Side A **********\n");
         printf("\t # of spectators in IT = %d\n",fpLostITA);
         printf("\t # of spectators in D1 = %d\n",fpLostD1A);
         printf("\t # of spectators in TDI = %d\n",fpLostTDI);
-        printf("\t # of spectators in luminometer = %d\n",fnLumiA);
-        printf("\t # of spectators in trousers = %d\n",fnTrou);
         printf("\t **********************************\n");
-      }
+      }*/
       gMC->StopTrack();
     }
     return;
@@ -2395,9 +2398,6 @@ void AliZDCv3::StepManager()
         if(xdet[1]<=0.) vol[1]=2;
         else vol[1]=4;
       }
-      if((vol[1]!=1) && (vol[1]!=2) && (vol[1]!=3) && (vol[1]!=4))
-        printf("       ZDC StepManager->ERROR in ZN!!! vol[1] = %d, xdet[0] = %f,"
-       "xdet[1] = %f\n",vol[1], xdet[0], xdet[1]);
     }
     
     else if(vol[0]==2){        //Quadrant in ZPC
@@ -2414,9 +2414,6 @@ void AliZDCv3::StepManager()
           break;
         }
       }
-      if((vol[1]!=1) && (vol[1]!=2) && (vol[1]!=3) && (vol[1]!=4))
-        printf("       ZDC StepManager->ERROR in ZP!!! vol[1] = %d, xdet[0] = %f,"
-       "xdet[1] = %f\n",vol[1], xdet[0], xdet[1]);
     }
     //
     // Quadrant in ZEM: vol[1] = 1 -> particle in 1st ZEM (placed at x = 8.5 cm)
@@ -2448,9 +2445,6 @@ void AliZDCv3::StepManager()
         if(xdet[1]<=0.) vol[1]=2;
         else vol[1]=4;
       }
-      if((vol[1]!=1) && (vol[1]!=2) && (vol[1]!=3) && (vol[1]!=4))
-        printf("       ZDC StepManager->ERROR in ZNA!!! vol[1] = %d, xdet[0] = %f,"
-       "xdet[1] = %f\n",vol[1], xdet[0], xdet[1]);
     }    
     //
     else if(vol[0]==5){        //Quadrant in ZPA
@@ -2467,82 +2461,88 @@ void AliZDCv3::StepManager()
           break;
         }
       }
-      if((vol[1]!=1) && (vol[1]!=2) && (vol[1]!=3) && (vol[1]!=4))
-        printf("       ZDC StepManager->ERROR in ZPA!!! vol[1] = %d, xdet[0] = %f,"
-       "xdet[1] = %f\n",vol[1], xdet[0], xdet[1]);
     }    
+    if((vol[1]!=1) && (vol[1]!=2) && (vol[1]!=3) && (vol[1]!=4))
+      AliError(Form(" WRONG tower for det %d: tow %d with xdet=(%f, %f)\n",
+               vol[0], vol[1], xdet[0], xdet[1]));
+    // Ch. debug
+    //printf("\t *** det %d vol %d xdet(%f, %f)\n",vol[0], vol[1], xdet[0], xdet[1]);
     
     
-  // Store impact point and kinetic energy of the ENTERING particle
+    // Store impact point and kinetic energy of the ENTERING particle
     
-      if(gMC->IsTrackEntering()){
-        //Particle energy
-        gMC->TrackMomentum(p[0],p[1],p[2],p[3]);
-        hits[3] = p[3];
-        // Impact point on ZDC  
+    if(gMC->IsTrackEntering()){
+      //Particle energy
+      gMC->TrackMomentum(p[0],p[1],p[2],p[3]);
+      hits[3] = p[3];
+      // Impact point on ZDC
+      // X takes into account the LHC x-axis sign
+      // which is opposite to positive x on detcetor front face
+      // for side A detectors (ZNA and ZPA)  
+      if(vol[0]==4 || vol[0]==5){
+        hits[4] = -xdet[0];
+      }
+      else{
         hits[4] = xdet[0];
-        hits[5] = xdet[1];
-       hits[6] = 0;
-        hits[7] = 0;
-        hits[8] = 0;
-        hits[9] = 0;
-       //
-       Int_t curTrackN = gAlice->GetMCApp()->GetCurrentTrackNumber();
-        TParticle *part = (gAlice->GetMCApp())->Particle(curTrackN);
-       hits[10] = part->GetPdgCode();
-       //printf("\t PDGCode = %d\n", part->GetPdgCode());
-        //
-       hits[11] = 1.0e09*gMC->TrackTime();
-       //printf("\t TrackTime = %f\n", hits[11]);
-
-       AddHit(curTrackN, vol, hits);
-       
-       if(fNoShower==1){
-         //printf("\t VolName %s -> det %d quad %d - x = %f, y = %f, z = %f\n", 
-           //knamed, vol[0], vol[1], x[0], x[1], x[2]);
-         if(vol[0]==1){
-           fnDetectedC += 1;
-           printf("      # of particles in ZNC = %d\n\n",fnDetectedC);
-         }
-         else if(vol[0]==2){
-           fpDetectedC += 1;
-           printf("      # of particles in ZPC = %d\n\n",fpDetectedC);
-         }
-         else if(vol[0]==4){
-           fnDetectedA += 1;
-           printf("      # of particles in ZNA = %d\n\n",fnDetectedA);     
-         }
-         else if(vol[0]==5){
-           fpDetectedA += 1;
-           printf("      # of particles in ZPA = %d\n\n",fpDetectedA);      
-         }
-          //
-         //printf("\t Particle: mass = %1.3f, E = %1.3f GeV, pz = %1.2f GeV -> stopped in volume %s\n", 
-          //   gMC->TrackMass(), p[3], p[2], gMC->CurrentVolName());
-         //
-         gMC->StopTrack();
-         return;
-       }
       }
-             
-      // Charged particles -> Energy loss
-      if((destep=gMC->Edep())){
-         if(gMC->IsTrackStop()){
-           gMC->TrackMomentum(p[0],p[1],p[2],p[3]);
-          m = gMC->TrackMass();
-          ekin = p[3]-m;
-          hits[9] = ekin;
-          hits[7] = 0.;
-          hits[8] = 0.;
-          AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(), vol, hits);
-          }
-        else{
-          hits[9] = destep;
-          hits[7] = 0.;
-          hits[8] = 0.;
-          AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(), vol, hits);
-          }
+      hits[5] = xdet[1];
+      hits[6] = 0;
+      hits[7] = 0;
+      hits[8] = 0;
+      hits[9] = 0;
+      //
+      Int_t curTrackN = gAlice->GetMCApp()->GetCurrentTrackNumber();
+      TParticle *part = gAlice->GetMCApp()->Particle(curTrackN);
+      hits[10] = part->GetPdgCode();
+      //printf("\t PDGCode = %d\n", part->GetPdgCode());
+      //
+      Int_t imo = part->GetFirstMother();
+      if(imo>0){
+       TParticle * pmot = gAlice->GetMCApp()->Particle(imo);
+       hits[11] = pmot->GetPdgCode();
       }
+      else hits[11]=0;
+      //
+      hits[12] = 1.0e09*gMC->TrackTime(); // in ns!
+      //printf("\t TrackTime = %f\n", hits[12]);
+
+      AddHit(curTrackN, vol, hits);
+
+      if(fNoShower==1){
+        //printf("\t VolName %s -> det %d quad %d - x = %f, y = %f, z = %f\n", 
+          //knamed, vol[0], vol[1], x[0], x[1], x[2]);
+        if(vol[0]==1){
+          fnDetectedC += 1;
+          printf("     # of particles in ZNC = %d\n\n",fnDetectedC);
+        }
+        else if(vol[0]==2){
+          fpDetectedC += 1;
+          printf("     # of particles in ZPC = %d\n\n",fpDetectedC);
+        }
+        else if(vol[0]==4){
+          fnDetectedA += 1;
+          printf("     # of particles in ZNA = %d\n\n",fnDetectedA);     
+        }
+        else if(vol[0]==5){
+          fpDetectedA += 1;
+          printf("     # of particles in ZPA = %d\n\n",fpDetectedA);      
+        }
+       //
+        //printf("\t Particle: mass = %1.3f, E = %1.3f GeV, pz = %1.2f GeV -> stopped in volume %s\n", 
+       //   gMC->TrackMass(), p[3], p[2], gMC->CurrentVolName());
+        //
+        gMC->StopTrack();
+        return;
+      }
+    }
+          
+    // Particle energy loss
+    if(gMC->Edep() != 0){
+      hits[9] = gMC->Edep();
+      hits[7] = 0.;
+      hits[8] = 0.;
+      AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(), vol, hits);
+    }
   }
  
 
@@ -2648,6 +2648,10 @@ void AliZDCv3::StepManager()
         Float_t z = -xalic[2]+fPosZEM[2]+2*fZEMLength-xalic[1];
         //z = xalic[2]-fPosZEM[2]-fZEMLength-xalic[1]*(TMath::Tan(45.*kDegrad));
          //printf("    fPosZEM[2]+2*fZEMLength = %f", fPosZEM[2]+2*fZEMLength);
+         //
+        // Parametrization for light guide uniformity
+         // NEW!!! Light guide tilted @ 51 degrees
+         Float_t guiPar[4]={0.31,-0.0006305,0.01337,0.8895};
         Float_t guiEff = guiPar[0]*(guiPar[1]*z*z+guiPar[2]*z+guiPar[3]);
         out = out*guiEff;
         nphe = gRandom->Poisson(out);