VZERO/V0_efficiencies.C

   1 void V0_efficiencies()
   2
   3 {
   4 // Caution: time windows and cuts on signals are HARD-coded
   5 //          they need to be adjusted to the current configuration
   6
   7         gROOT->Reset();
   8         rl = AliRunLoader::Open("galice.root");
   9         rl->LoadgAlice();
  10         gAlice = rl->GetAliRun();
  11         Int_t nevent = rl->GetNumberOfEvents();
  12         cout<<" --- Number of events in file = "<< nevent <<" ---"<<endl;
  13
  14 //___________________________________________________
  15 // Book HISTOGRAMS
  16
  17         TH2F *hMapV0L   = new TH2F("hMapV0L","V0L",161,-0.8,0.8,161,-0.8,0.8);
  18         TH2F *hMapV0R   = new TH2F("hMapV0R","V0R",161,-0.8,0.8,161,-0.8,0.8);
  19         TH1F *hpdgV0    = new TH1F("hpdgV0","pdgV0",1001,-500,500);
  20         TH1F *hvertexZ  = new TH1F("hvertexZ","vertex Z",972,-40.,40.);
  21         TH1F *hTimC     = new TH1F("hTimC","Time of Flight in V0C",500,0,49);
  22         TH1F *hTimA     = new TH1F("hTimA","Time of Flight in V0A",500,0,49);
  23         TH1F *hCell     = new TH1F("hCell","Cell Number",100,0,99);
  24         TH2F *hCellADC  = new TH2F("hCellADC","ADC vs Cell",100,0,99,1000,0,3999);
  25
  26         TH1F *hMul0     = new TH1F("hMul0 ","Multiplicity in V0",80,0,79);
  27         TH1F *hMulC0    = new TH1F("hMulC0","Multiplicity in V0C",50,0,49);
  28         TH1F *hMulA0    = new TH1F("hMulA0","Multiplicity in V0A",50,0,49);
  29         TH1F *hMulAnd0  = new TH1F("hMulAnd0","Trigger and",50,0,49);
  30
  31 //___________________________________________________
  32
  33         AliLoader *ld = rl->GetLoader("VZEROLoader");
  34         ld->GetHitsDataLoader()->Load("READ");
  35
  36         rl->LoadHeader();
  37
  38         // to access the particle  Stack
  39         rl->LoadKinematics("READ");
  40
  41         AliVZERO  *v0 = (AliVZERO*)gAlice->GetDetector("VZERO");
  42
  43         bool    flagV0L  = false;
  44         bool    flagV0R  = false;
  45         Float_t timeV0L  = 1e12;
  46         Float_t timeV0R  = 1e12;
  47
  48         Int_t   fNdigits = 0;
  49         Int_t   fMulA    = 0;
  50         Int_t   fMulC    = 0;
  51         Int_t   fDigits  = 0;
  52         Float_t fPhotoCathodeEfficiency =   0.18;
  53
  54
  55         Int_t   nVOL     = 0;
  56         Int_t   nVOR     = 0;
  57         Int_t   nVOLetR  = 0;
  58         Int_t   nVOLouR  = 0;
  59
  60         TDatabasePDG * pdgdb = TDatabasePDG::Instance();
  61
  62
  63         Float_t fPMVoltage =  768.0;
  64         Float_t fPMGain1   = TMath::Power((fPMVoltage / 112.5) ,7.04277);
  65         Float_t fPMGain[64];
  66         Float_t cPM[64];
  67
  68         for(Int_t ii=0; ii<64; ii++){
  69         fPMGain[ii] = gRandom->Gaus(fPMGain1, fPMGain1/5); // 20% uncertainty on the PM gain
  70         cPM[ii] = fPhotoCathodeEfficiency * fPMGain[ii];}
  71
  72         Float_t kC     =  1e-11;
  73         Float_t kthau  =  2.1*1e-9;
  74         Float_t ktheta =  50.0 * kC;
  75         Float_t kQe    =  1.6e-19;
  76         Float_t coef   =  200.0;       //  p-p
  77 //      Float_t coef   =    1.0;       // Pb-Pb
  78
  79         Int_t map[80];
  80         for(Int_t i=0; i<80; i++) map[i] = 0;
  81         Int_t hit[80];
  82         for(Int_t i=0; i<80; i++) hit[i] = 0;
  83
  84         for(Int_t event = 0; event < nevent; event++){
  85         Float_t timeV0L = 1e12;
  86         Float_t timeV0R = 1e12;
  87
  88             rl->GetEvent(event);
  89             for(Int_t i=0; i<80; i++) map[i] = 0;
  90             for(Int_t i=0; i<80; i++) hit[i] = 0;
  91             ld->LoadHits();
  92             v0->SetTreeAddress();
  93             TTree *treeH = ld->TreeH();
  94             Int_t ntracks = treeH->GetEntries();
  95
  96             for (Int_t itrack=0; itrack<ntracks;itrack++) {
  97
  98                 v0->ResetHits();
  99                 treeH->GetEvent(itrack);
 100                 if(v0){
 101                   for (AliVZEROhit *vhit=(AliVZEROhit*)v0->FirstHit(itrack);
 102                           vhit; vhit=(AliVZEROhit*)v0->NextHit())
 103                   {
 104 //                      vhit->Dump();
 105                         hpdgV0->Fill(vhit->TrackPiD());
 106                         hvertexZ->Fill(vhit->Vz()/100.);
 107
 108                         Float_t dt_scintillator = gRandom->Gaus(0,0.7);
 109                         Float_t time = dt_scintillator + 1e9*vhit->Tof();
 110
 111                         if(vhit->Z() > 0){
 112                         if (time < 9 || time > 13) {continue;}          //time window V0A : 11 +/- 2 ns
 113
 114                         flagV0L = true;
 115                         hTimA->Fill(time);
 116                         if(time < timeV0L) timeV0L = time;
 117                         hMapV0L->Fill(vhit->X()/100.,vhit->Y()/100.);}
 118
 119                         if(vhit->Z() < 0){
 120                         if (time < 1 ||time > 5) {continue;}            //time window V0C : 3 +/- 2 ns
 121
 122                         flagV0R = true;
 123                         hTimC->Fill(time);
 124                         if(time < timeV0R) timeV0R = time;
 125                         hMapV0R->Fill(vhit->X()/100.,vhit->Y()/100.);}
 126
 127                         Int_t nPhot = vhit->Nphot();
 128                         Int_t cell  = vhit->Cell();
 129                         map[cell] += nPhot;
 130                         hit[cell] ++;
 131
 132                   }
 133                 }
 134              }
 135
 136
 137         Int_t map2[64];      // cell to digits
 138         Int_t hit2[64];      // cell to digits
 139         Int_t j;
 140         Float_t time1[80];
 141         Float_t time2[64];
 142
 143         for (j=0; j<16; j++){
 144         map2[j] = map [j];
 145         hit2[j] = hit [j];
 146         //time2[j]=time1[j];
 147         }
 148
 149         for (j=0; j<16; j++){
 150         map2[j+16] = map [2*j+16]+map [2*j+17];
 151         hit2[j+16] = hit [2*j+16]+hit [2*j+17];
 152         //time2[j+16]= TMath::Min(time1 [16 + 2*j], time1 [16 + 2*j + 1]);
 153         }
 154
 155         for (j=32; j<64; j++){
 156         map2[j] =map[j+16];
 157         hit2[j] =hit[j+16];
 158         //time2[j] =time[j+16];
 159         }
 160
 161            fNdigits = 0;
 162            fMulC    = 0;
 163            fMulA    = 0;
 164            for (Int_t i=0; i<64; i++) {
 165                  Float_t q1 = Float_t ( map2[i] )* cPM[i] * kQe;
 166                  Float_t noise = gRandom->Gaus(10.5,3.22);
 167                  Float_t pmResponse  =  q1/kC*TMath::Power(ktheta/kthau,1/(1-ktheta/kthau))
 168                               + noise*1e-3;
 169                  map2[i] = Int_t( pmResponse * coef);
 170                  int test = 0;
 171                  if(map2[i] > 10) {map2[i] = Int_t( gRandom->Gaus(map2[i], 80));}               // charge smearing of MIP/4 -> sigma = MIP/4 = 120 (MIP = 480)
 172                  if(map2[i] > 240) {                                                            // cut at MIP/2 = 240
 173                        hCell->Fill(float(i));
 174                        hCellADC->Fill(float(i),map2[i]);
 175                        fNdigits++;
 176                        if(i<32) fMulC++;
 177                        if(i>31) fMulA++;
 178                        }
 179
 180              }
 181              hMul0->Fill(fNdigits);
 182              hMulC0->Fill(fMulC);
 183              hMulA0->Fill(fMulA);
 184              hMulAnd0->Fill(TMath::Min(fMulA, fMulC));
 185
 186              if(fMulA > 0){
 187                 nVOL++;}
 188
 189              if(fMulC > 0){
 190                 nVOR++;}
 191
 192              if(fMulA > 0 && fMulC > 0){
 193                 nVOLetR++;}
 194
 195              if(fMulA > 0 || fMulC > 0){
 196                 nVOLouR++;}
 197
 198              if(event%100==0) cout <<" event    = " << event <<endl;
 199 //           cout <<" multi   = " << fNdigits <<endl;
 200         }
 201
 202         cout <<" nVOA     = " << nVOL <<endl;
 203         cout <<" nVOC     = " << nVOR <<endl;
 204         cout <<" nVOAandC = " << nVOLetR <<endl;
 205         cout <<" nVOAorC  = " << nVOLouR <<endl;
 206
 207 //__________________________________________________
 208 //      Fill root file
 209
 210         TFile *histoFile = new TFile("Efficiencies.root","RECREATE");
 211
 212         hMapV0L->Write();
 213         hMapV0R->Write();
 214         hpdgV0->Write();
 215         hvertexZ->Write();
 216         hTimC->Write();
 217         hTimA->Write();
 218         hCell->Write();
 219         hCellADC->Write();
 220
 221         hMul0->Write();
 222         hMulC0->Write();
 223         hMulA0->Write();
 224         hMulAnd0->Write();
 225         histoFile->Close();
 226
 227 }