Constant extrapolation for gain calibration

[u/mrichter/AliRoot.git] / EVGEN / AliGenParam.cxx

diff --git a/EVGEN/AliGenParam.cxx b/EVGEN/AliGenParam.cxx
index 2893051f14b9502c2890dd7176df4f5227e1d009..c3b4cf73bcbb1d5046cca4c5e77feb94aa6bcb84 100644 (file)
--- a/EVGEN/AliGenParam.cxx
+++ b/EVGEN/AliGenParam.cxx
@@ -22,6 +22,7 @@
 // andreas.morsch@cern.ch
 
 #include <TCanvas.h>
+#include <TClonesArray.h>
 #include <TDatabasePDG.h>
 #include <TF1.h>
 #include <TH1F.h>
@@ -37,6 +38,7 @@
 #include "AliGenParam.h"
 #include "AliMC.h"
 #include "AliRun.h"
+#include "AliGenEventHeader.h"
 
 ClassImp(AliGenParam)
 
@@ -67,7 +69,7 @@ ClassImp(AliGenParam)
 // Default constructor
 }
 //____________________________________________________________
-AliGenParam::AliGenParam(Int_t npart, AliGenLib * Library,  Int_t param, char* tname)
+AliGenParam::AliGenParam(Int_t npart, AliGenLib * Library,  Int_t param, const char* tname)
     :AliGenMC(npart),
      fPtParaFunc(Library->GetPt(param, tname)),
      fYParaFunc (Library->GetY (param, tname)),
@@ -129,8 +131,8 @@ AliGenParam::AliGenParam(Int_t npart, Int_t param, const char* tname, const char
 //____________________________________________________________
 
 AliGenParam::AliGenParam(Int_t npart, Int_t param,
-                         Double_t (*PtPara) (Double_t*, Double_t*),
-                         Double_t (*YPara ) (Double_t* ,Double_t*),
+                         Double_t (*PtPara) (const Double_t*, const Double_t*),
+                         Double_t (*YPara ) (const Double_t* ,const Double_t*),
                         Int_t    (*IpPara) (TRandom *))                 
     :AliGenMC(npart),
      
@@ -282,6 +284,8 @@ void AliGenParam::Generate()
   Int_t ipa=0;
   
 // Generating fNpart particles
+  fNprimaries = 0;
+  
   while (ipa<fNpart) {
       while(1) {
 //
@@ -317,7 +321,7 @@ void AliGenParam::Generate()
                    "Division by 0: Please check you rapidity range !");
          }
          
-         pl=xmt*ty/sqrt(1.-ty*ty);
+         pl=xmt*ty/sqrt((1.-ty)*(1.+ty));
          theta=TMath::ATan2(pt,pl);
 // Cut on theta
          if(theta<fThetaMin || theta>fThetaMax) continue;
@@ -344,6 +348,8 @@ void AliGenParam::Generate()
 // if fForceDecay == none Primary particle is tracked by GEANT 
 // (In the latest, make sure that GEANT actually does all the decays you want)   
 //
+         Bool_t decayed = kFALSE;
+         
 
          if (fForceDecay != kNoDecay) {
 // Using lujet to decay particle
@@ -370,7 +376,8 @@ void AliGenParam::Generate()
              }
              
              if (np >1) {
-                 TParticle* iparticle =  (TParticle *) particles->At(0);
+                 decayed = kTRUE;
+                 TParticle* iparticle =  0;
                  Int_t ipF, ipL;
                  for (i = 1; i<np ; i++) {
                      trackIt[i] = 1;
@@ -406,7 +413,7 @@ void AliGenParam::Generate()
 //
 // children
                      
-                     if (ChildSelected(TMath::Abs(kf)) || fForceDecay == kAll && trackIt[i])
+                     if ((ChildSelected(TMath::Abs(kf)) || fForceDecay == kAll) && trackIt[i])
                      {
                          if (fCutOnChild) {
                              pc[0]=iparticle->Px();
@@ -433,17 +440,22 @@ void AliGenParam::Generate()
                  ipa++;
 //
 // Parent
-                 PushTrack(0, -1, iPart, p, origin0, polar, 0, kPPrimary, nt, wgtp);
+                 
+                 
+                 PushTrack(0, -1, iPart, p, origin0, polar, 0, kPPrimary, nt, wgtp, ((decayed)? 11 : 1));
                  pParent[0] = nt;
                  KeepTrack(nt); 
+                 fNprimaries++;
+                 
 //
 // Decay Products
 //               
                  for (i = 1; i < np; i++) {
                      if (pSelected[i]) {
                          TParticle* iparticle = (TParticle *) particles->At(i);
-                         Int_t kf  = iparticle->GetPdgCode();
-                         Int_t ipa = iparticle->GetFirstMother()-1;
+                         Int_t kf   = iparticle->GetPdgCode();
+                         Int_t ksc  = iparticle->GetStatusCode();
+                         Int_t jpa  = iparticle->GetFirstMother()-1;
                          
                          och[0] = origin0[0]+iparticle->Vx()/10;
                          och[1] = origin0[1]+iparticle->Vy()/10;
@@ -452,17 +464,18 @@ void AliGenParam::Generate()
                          pc[1]  = iparticle->Py();
                          pc[2]  = iparticle->Pz();
                          
-                         if (ipa > -1) {
-                             iparent = pParent[ipa];
+                         if (jpa > -1) {
+                             iparent = pParent[jpa];
                          } else {
                              iparent = -1;
                          }
                         
-                         PushTrack(fTrackIt*trackIt[i], iparent, kf,
+                         PushTrack(fTrackIt * trackIt[i], iparent, kf,
                                           pc, och, polar,
-                                          0, kPDecay, nt, wgtch);
+                                          0, kPDecay, nt, wgtch, ksc);
                          pParent[i] = nt;
                          KeepTrack(nt); 
+                         fNprimaries++;
                      } // Selected
                  } // Particle loop 
              }  // Decays by Lujet
@@ -475,13 +488,20 @@ void AliGenParam::Generate()
          else  // nodecay option, so parent will be tracked by GEANT (pions, kaons, eta, omegas, baryons)
          {
            gAlice->GetMCApp()->
-               PushTrack(fTrackIt,-1,iPart,p,origin0,polar,0,kPPrimary,nt,wgtp);
+               PushTrack(fTrackIt,-1,iPart,p,origin0,polar,0,kPPrimary,nt,wgtp, 1);
             ipa++; 
+           fNprimaries++;
          }
          break;
     } // while
   } // event loop
+  
   SetHighWaterMark(nt);
+
+  AliGenEventHeader* header = new AliGenEventHeader("PARAM");
+  header->SetPrimaryVertex(fVertex);
+  header->SetNProduced(fNprimaries);
+  AddHeader(header);
 }
 //____________________________________________________________________________________
 Float_t AliGenParam::GetRelativeArea(Float_t ptMin, Float_t ptMax, Float_t yMin, Float_t yMax, Float_t phiMin, Float_t phiMax)