MUON/AliMUONproto.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 $Log$
  17 Revision 1.5  2002/10/23 07:24:57  alibrary
  18 Introducing Riostream.h
  19
  20 Revision 1.4  2001/01/26 20:00:53  hristov
  21 Major upgrade of AliRoot code
  22
  23 Revision 1.3  2000/10/11 09:19:12  egangler
  24 Corrected some bugs - it should compile now
  25
  26 Revision 1.2  2000/06/15 07:58:49  morsch
  27 Code from MUON-dev joined
  28
  29 Revision 1.1.2.1  2000/04/18 09:11:15  morsch
  30 Implementation of MUON Chamber Prototype Class
  31 Both read digits from raw data or use the Monte-Carlo.
  32 Rachid GUERNANE, IPN Lyon guernane@ipnl.in2p3.fr
  33
  34 */
  35
  36 /*
  37 Implementation of MUON Chamber Prototype Class
  38 Both read digits from raw data or use the Monte-Carlo.
  39 1-Feb-2000 Rachid GUERNANE, IPN Lyon guernane@ipnl.in2p3.fr
  40 */
  41
  42
  43 ////////////////////////////////////////////////
  44 //  Manager and hits classes for set:PROTO    //
  45 ////////////////////////////////////////////////
  46
  47 #include <TTUBE.h>
  48 #include <TBRIK.h>
  49 #include <TRotMatrix.h>
  50 #include <TNode.h>
  51 #include <TTree.h>
  52 #include <TRandom.h>
  53 #include <TObject.h>
  54 #include <TVector.h>
  55 #include <TObjArray.h>
  56 #include <TMinuit.h>
  57 #include <TParticle.h>
  58 #include <TROOT.h>
  59 #include <TFile.h>
  60 #include <TNtuple.h>
  61 #include <TCanvas.h>
  62 #include <TPad.h>
  63 #include <TDirectory.h>
  64 #include <TObjectTable.h>
  65 #include <AliPDG.h>
  66 #include <TArrayI.h>
  67 #include <AliDetector.h>
  68
  69 #include "AliMUONChamber.h"
  70 #include "AliMUONproto.h"
  71 #include "AliMUONHit.h"
  72 #include "TTUBE.h"
  73 #include "AliMUONClusterFinder.h"
  74 #include "AliRun.h"
  75 #include "Riostream.h"
  76 #include "AliCallf77.h"
  77 #include "AliConst.h"
  78 //#include "chainalice2.h"
  79 #include "AliMUONSegmentationV0.h"
  80 //#include "AliMUONSegResV11.h"
  81
  82 ClassImp(AliMUONproto)
  83
  84 //___________________________________________
  85 AliMUONproto::AliMUONproto()
  86         : AliMUON()
  87 {
  88     cout << "\n Calling AliMUONproto constructor..." << endl;
  89
  90         //
  91         //
  92         //
  93 }
  94
  95 //___________________________________________
  96 AliMUONproto::AliMUONproto(const char *name, const char *title)
  97         : AliMUON(name,title)
  98 //  AliMUON defines everything, but the chamber for NCH=1
  99 {
 100
 101 //  z-Positions of Chambers
 102     const Float_t zch[1] = {975.};
 103 //
 104 //  inner diameter
 105     const Float_t dmi[1] = {0.};
 106 //
 107 //  outer diameter
 108     const Float_t dma[1] = {29.2};
 109 //
 110 //
 111 //    Default Parameters for ALICE2 prototype
 112
 113 //
 114     (*fChambers)[0] = new AliMUONChamber();
 115     AliMUONChamber* chamber = (AliMUONChamber*) (*fChambers)[0];
 116     chamber->SetGid(0);
 117     chamber->SetZ(zch[0]);
 118 //
 119     chamber->InitGeo(zch[0]);
 120     chamber->SetRInner(dmi[0]/2);
 121     chamber->SetROuter(dma[0]/2);
 122
 123     for (int i = 0; i <= 99; i++) {
 124         fThreshold[i] = 0.;
 125     }
 126
 127 }
 128
 129 #ifdef WE_FORGET_FOR_THE_MOMENT
 130 //___________________________________________
 131 void AliMUONproto::GetRawDigits(Int_t evnb, Int_t *lptr, Int_t ilen) {
 132
 133     Int_t ip = 0;
 134     Int_t equip = 0;
 135     Int_t nochnl;
 136     Int_t loop;
 137     Int_t val;
 138     Int_t itype;
 139     Int_t id;
 140     Int_t serial;
 141     Int_t equiplen;
 142
 143     Int_t digits[5];
 144
 145     AliMUON* MUON  = (AliMUON*)gAlice->GetModule("MUON");
 146     AliMUONSegmentationV0* seg = (AliMUONSegmentationV0*) Chamber(0).SegmentationModel(1);
 147
 148     Int_t tracks[10];
 149     Int_t charges[10];
 150
 151     for (Int_t i = 0; i < 10; i++) {
 152         tracks[i] = 0;
 153         charges[i] = 0;
 154     }
 155
 156     Int_t ich = 0;
 157
 158   nwtype:
 159
 160     itype = lptr[ip++];
 161     id = lptr[ip++];
 162     equiplen = lptr[ip++];
 163
 164     if (equiplen < 0 ) return;
 165
 166     if (itype != (int)(0XCACCA008)) {
 167         ip += equiplen;
 168         if (ip < ilen-2)  goto nwtype;
 169     }
 170     else {
 171         serial = id >> 16;
 172         equip = id & 0x1;
 173         if ((serial == 190) && (equip == 1)) {
 174             for (loop = 0; loop < equiplen; loop++) {
 175                 nochnl = (lptr[ip] & 0x7ff000 ) >> 12;
 176                 val = lptr[ip] & 0x3ff;
 177                 // fill digits from raw data according to cathode connexions
 178                 if (group[nochnl][2][1]!=0)
 179                     digits[0] = group[nochnl][2][1] - 12;
 180                 else if  (group[nochnl][1][1]!=0)
 181                     digits[0] = group[nochnl][1][1] - 12;
 182                 else
 183                     digits[0] = group[nochnl][0][1] - 12;
 184                 digits[1] = group[nochnl][0][0];
 185                 if (digits[0] != seg->Ix(digits[0], digits[1]))
 186                     printf("Pb pour ix=%d,iy=%d\n", digits[0], digits[1]);
 187                 digits[2] = val;
 188                 digits[3] = 0;
 189                 digits[4] = 0;
 190                 if (digits[2] >= fThreshold[nochnl])
 191                     MUON->AddDigits(ich, tracks, charges, digits);
 192
 193                 ip++;
 194             }
 195         }
 196         else
 197             ip += equiplen;
 198
 199         if (ip < ilen-2) goto nwtype;
 200     }
 201
 202     gAlice->TreeD()->Fill();
 203     MUON->ResetDigits();
 204
 205     gAlice->TreeD()->Fill();
 206     MUON->ResetDigits();
 207
 208     char hname[30];
 209     sprintf(hname, "TreeD%d", evnb);
 210     gAlice->TreeD()->Write(hname,TObject::kOverwrite);
 211     // reset tree
 212     gAlice->TreeD()->Reset();
 213
 214 }
 215 #endif
 216
 217 //___________________________________________
 218 void AliMUONproto::SetPadSize(Int_t id, Int_t isec, Float_t p1, Float_t p2)
 219 {
 220     Int_t i=2*(id-1);
 221     ((AliMUONChamber*)(*fChambers)[i])->SetPadSize(isec,p1,p2);
 222 }
 223
 224 //___________________________________________
 225 void AliMUONproto::SetChargeSlope(Int_t id, Float_t p1)
 226 {
 227     Int_t i=2*(id-1);
 228     ((AliMUONChamber*) (*fChambers)[i])->SetChargeSlope(p1);
 229 }
 230
 231 //___________________________________________
 232 void AliMUONproto::SetChargeSpread(Int_t id, Float_t p1, Float_t p2)
 233 {
 234     Int_t i=2*(id-1);
 235     ((AliMUONChamber*) (*fChambers)[i])->SetChargeSpread(p1,p2);
 236 }
 237
 238 //___________________________________________
 239 void AliMUONproto::SetMaxAdc(Int_t id, Float_t p1)
 240 {
 241     Int_t i=2*(id-1);
 242     ((AliMUONChamber*) (*fChambers)[i])->SetMaxAdc(p1);
 243 }
 244
 245 //___________________________________________
 246 void AliMUONproto::CreateGeometry()
 247 {
 248     Int_t *idtmed = fIdtmed->GetArray()-1099;
 249 //
 250     Float_t zpos;
 251     Float_t bpar[3];
 252     Int_t idGas=idtmed[1100];
 253
 254     AliMUONChamber *iChamber;
 255 //********************************************************************
 256 //                          Prototype ALICE2                        **
 257 //********************************************************************
 258     iChamber=(AliMUONChamber*) (*fChambers)[0];
 259
 260     zpos=iChamber->Z();
 261 //
 262     const Float_t X_POS = 11*.975/2; //half x
 263     const Float_t Y_POS = 18*.55/2;  //half y
 264     const Float_t Z_POS = 0.325;
 265
 266     bpar[0] = X_POS;
 267     bpar[1] = Y_POS;
 268     bpar[2] = Z_POS;
 269
 270     gMC->Gsvolu("C01G", "BOX", idGas, bpar, 3);
 271     gMC->Gspos("C01G", 1, "ALIC", -bpar[0], bpar[1], zpos, 0, "ONLY");
 272 }
 273
 274 //___________________________________________
 275 void AliMUONproto::CreateMaterials()
 276 {
 277   // *** DEFINITION OF AVAILABLE MUON MATERIALS ***
 278   //
 279   //     Ar-CO2 gas
 280     Float_t ag1[3]   = { 39.95,12.01,16. };
 281     Float_t zg1[3]   = { 18.,6.,8. };
 282     Float_t wg1[3]   = { .8,.0667,.13333 };
 283     Float_t dg1      = .001821;
 284
 285     Float_t epsil  = .001; // Tracking precision,
 286     Float_t tmaxfd = -20.; // Maximum angle due to field deflection
 287     Float_t stmin  = -.8;
 288     Int_t ISXFLD   = gAlice->Field()->Integ();
 289     Float_t SXMGMX = gAlice->Field()->Max();
 290
 291     //
 292     // --- Define the various materials for GEANT ---
 293     AliMixture(22, "ArCO2 80%$", ag1, zg1, dg1, 3, wg1);
 294
 295     //
 296     //    Ar/CO2
 297     AliMedium(1, "ARG_CO2   ", 22, 1, ISXFLD, SXMGMX, tmaxfd, fMaxStepGas,
 298               fMaxDestepAlu, epsil, stmin);
 299     //    Air
 300         //AliMedium(1, "AIR_CH_US         ", 15, 1, ISXFLD, SXMGMX, tmaxfd, -1., -.3, epsil, stmin);
 301 }
 302
 303 //___________________________________________
 304 void AliMUONproto::Init()
 305 {
 306    printf("\n\n\n Start Init for Prototype ALICE2 - CPC chamber type\n\n\n");
 307
 308    //
 309    // Initialize Tracking Chambers
 310    //
 311
 312    ( (AliMUONChamber*) (*fChambers)[0])->Init();
 313
 314    //
 315    // Set the chamber (sensitive region) GEANT identifier
 316    ((AliMUONChamber*)(*fChambers)[0])->SetGid(gMC->VolId("C01G"));
 317
 318    printf("\n\n\n Finished Init for Prototype ALICE2 - CPC chamber type\n\n\n");
 319
 320 }
 321
 322 //___________________________________________
 323 void AliMUONproto::StepManager()
 324 {
 325   Int_t          copy, id;
 326   static Int_t   idvol;
 327   static Int_t   vol[2];
 328   Int_t          ipart;
 329   TLorentzVector pos;
 330   TLorentzVector mom;
 331   Float_t        theta,phi;
 332   Float_t        destep, step;
 333
 334   static Float_t eloss, eloss2, xhit, yhit, tlength;
 335   const  Float_t big=1.e10;
 336
 337   //  modifs perso
 338   static Float_t hits[14];
 339
 340   TClonesArray &lhits = *fHits;
 341
 342   //
 343   // Set maximum step size for gas
 344   // numed=gMC->GetMedium();
 345   //
 346   // Only charged tracks
 347   if( !(gMC->TrackCharge()) ) return;
 348   //
 349   // Only gas gap inside chamber
 350   // Tag chambers and record hits when track enters
 351   idvol=-1;
 352   id=gMC->CurrentVolID(copy);
 353
 354     for (Int_t i=1; i<=kNCH; i++) {
 355       if(id==((AliMUONChamber*)(*fChambers)[i-1])->GetGid()){
 356           vol[0]=i;
 357           idvol=i-1;
 358       }
 359     }
 360     if (idvol == -1) return;
 361   //
 362   // Get current particle id (ipart), track position (pos)  and momentum (mom)
 363   gMC->TrackPosition(pos);
 364   gMC->TrackMomentum(mom);
 365
 366   ipart  = gMC->TrackPid();
 367   //Int_t ipart1 = gMC->IdFromPDG(ipart);
 368   //printf("ich, ipart %d %d \n",vol[0],ipart1);
 369
 370   //
 371   // momentum loss and steplength in last step
 372   destep = gMC->Edep();
 373   step   = gMC->TrackStep();
 374
 375   //
 376   // record hits when track enters ...
 377   if( gMC->IsTrackEntering()) {
 378       gMC->SetMaxStep(fMaxStepGas);
 379       Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
 380       Double_t rt = TMath::Sqrt(tc);
 381       Double_t pmom = TMath::Sqrt(tc+mom[2]*mom[2]);
 382       Double_t tx=mom[0]/pmom;
 383       Double_t ty=mom[1]/pmom;
 384       Double_t tz=mom[2]/pmom;
 385       Double_t s=((AliMUONChamber*)(*fChambers)[idvol])
 386           ->ResponseModel()
 387           ->Pitch()/tz;
 388       theta   = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
 389       phi     = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
 390       hits[0] = Float_t(ipart);         // Geant3 particle type
 391       hits[1] = pos[0]+s*tx;                 // X-position for hit
 392       hits[2] = pos[1]+s*ty;                 // Y-position for hit
 393       hits[3] = pos[2]+s*tz;                 // Z-position for hit
 394       hits[4] = theta;                  // theta angle of incidence
 395       hits[5] = phi;                    // phi angle of incidence
 396       hits[8] = (Float_t) fNPadHits;   // first padhit
 397       hits[9] = -1;                     // last pad hit
 398
 399       // modifs perso
 400       hits[10] = mom[3]; // hit momentum P
 401       hits[11] = mom[0]; // Px/P
 402       hits[12] = mom[1]; // Py/P
 403       hits[13] = mom[2]; // Pz/P
 404       // fin modifs perso
 405
 406       // phi angle of incidence
 407       tlength = 0;
 408       eloss   = 0;
 409       eloss2  = 0;
 410       xhit    = pos[0];
 411       yhit    = pos[1];
 412       // Only if not trigger chamber
 413       if(idvol<10) {
 414           //
 415           //  Initialize hit position (cursor) in the segmentation model
 416           ((AliMUONChamber*) (*fChambers)[idvol])
 417               ->SigGenInit(pos[0], pos[1], pos[2]);
 418       } else {
 419           //geant3->Gpcxyz();
 420           //printf("In the Trigger Chamber #%d\n",idvol-9);
 421       }
 422   }
 423   eloss2+=destep;
 424
 425   //
 426   // Calculate the charge induced on a pad (disintegration) in case
 427   //
 428   // Mip left chamber ...
 429   if( gMC->IsTrackExiting() || gMC->IsTrackStop() || gMC->IsTrackDisappeared()){
 430       gMC->SetMaxStep(big);
 431       eloss   += destep;
 432       tlength += step;
 433
 434       // Only if not trigger chamber
 435       if(idvol<10) {
 436           if (eloss > 0)
 437               MakePadHits(0.5*(xhit+pos[0]),0.5*(yhit+pos[1]),eloss,0.0,idvol);
 438       }
 439
 440       hits[6]=tlength;
 441       hits[7]=eloss2;
 442       if (fNPadHits > (Int_t)hits[8]) {
 443           hits[8]= hits[8]+1;
 444           hits[9]= (Float_t) fNPadHits;
 445       }
 446
 447       new(lhits[fNhits++])
 448           AliMUONHit(fIshunt,gAlice->CurrentTrack(),vol,hits);
 449       eloss = 0;
 450       //
 451       // Check additional signal generation conditions
 452       // defined by the segmentation
 453       // model (boundary crossing conditions)
 454   } else if
 455       (((AliMUONChamber*) (*fChambers)[idvol])
 456        ->SigGenCond(pos[0], pos[1], pos[2]))
 457   {
 458       ((AliMUONChamber*) (*fChambers)[idvol])
 459           ->SigGenInit(pos[0], pos[1], pos[2]);
 460 //      printf("\n-> MakePadHits, reason special %d",ipart);
 461       if (eloss > 0)
 462           MakePadHits(0.5*(xhit+pos[0]),0.5*(yhit+pos[1]),eloss,0.0,idvol);
 463       xhit     = pos[0];
 464       yhit     = pos[1];
 465       eloss    = destep;
 466       tlength += step ;
 467       //
 468       // nothing special  happened, add up energy loss
 469   } else {
 470       eloss   += destep;
 471       tlength += step ;
 472   }
 473 }
 474
 475 //___________________________________________
 476 void AliMUONproto::BuildGeometry()
 477 {
 478     TNode *Node;
 479     TNode *Top;
 480
 481     const int kColorMUON = kBlue;
 482     //
 483     Top=gAlice->GetGeometry()->GetNode("alice");
 484     //
 485     //
 486     //
 487     float dx, dy, dz, zpos;
 488
 489     const Float_t cz[1]={975.};
 490
 491     zpos=cz[0];
 492
 493     dx=11*.975/2;
 494     dy=18*.55/2;
 495     dz=0.325;
 496
 497     new TBRIK("C_MUON101","Mother volume for Proto.","void",dx*2,dy*2,dz*2);
 498     TBRIK* PROTO = new TBRIK("PROT", "Proto. Sens. Region","void",dx,dy,dz);
 499     Top->cd();
 500     Node = new TNode("MUON101","ChamberNode","C_MUON101",0,0,zpos,"");
 501     Node->SetLineColor(kColorMUON);
 502     Node->SetVisibility(0);
 503     fNodes->Add(Node);
 504     Node->cd();
 505     Node = new TNode("MUON201", "Proto. Sens. Region Node", PROTO, -dx, dy, dz);
 506     Node->SetLineColor(kColorMUON);
 507
 508 }
 509
 510 //_____________________________________________________________________________
 511 void AliMUONproto::FindClusters(Int_t nev, Int_t last_entry)
 512 {
 513
 514 //
 515 // Loop on chambers and on cathode planes
 516 //
 517   for (Int_t icat = 0; icat < 2; icat++) {
 518             gAlice->ResetDigits();
 519             gAlice->TreeD()->GetEvent(last_entry+icat);
 520
 521       for (Int_t ich = 0; ich < kNTrackingCh; ich++) {
 522           AliMUONChamber* iChamber=(AliMUONChamber*) (*fChambers)[ich];
 523           TClonesArray *MUONdigits  = this->DigitsAddress(ich);
 524           if (MUONdigits == 0) continue;
 525           //
 526           // Get ready the current chamber stuff
 527           //
 528
 529           AliMUONResponse* response = iChamber->ResponseModel();
 530           AliMUONSegmentation* seg = iChamber->SegmentationModel(icat+1);
 531           AliMUONClusterFinder* rec = iChamber->ReconstructionModel();
 532
 533           if (seg) {
 534               rec->SetSegmentation(seg);
 535               rec->SetResponse(response);
 536               rec->SetDigits(MUONdigits);
 537               rec->SetChamber(ich);
 538               rec->FindRawClusters();
 539           }
 540
 541           TClonesArray *fRch;
 542           fRch=RawClustAddress(ich);
 543           fRch->Sort();
 544           // it seems to work
 545
 546       } // for ich
 547       // fill the tree
 548       gAlice->TreeR()->Fill();
 549
 550       ResetRawClusters();
 551   } // for icat
 552
 553   char hname[30];
 554   sprintf(hname,"TreeR%d",nev);
 555   gAlice->TreeR()->Write(hname);
 556   gAlice->TreeR()->Reset();
 557 }
 558
 559 #ifdef WE_FORGRET_THIS_SHIT
 560 //_____________________________________________________________________________
 561 void AliMUONproto::SetThreshold()
 562 {
 563
 564     ifstream inputFile("/home/alice/guernane/aliroot/pro/MUON/crped190.dat", ios::in);
 565
 566     if (inputFile.fail()) {
 567         cout << "Error opening file" << endl;
 568         exit(2);
 569     }
 570
 571     char buff[32];
 572     Int_t Serial;
 573     Int_t Ntrigger;
 574     Int_t Nchannel;
 575     Int_t i1;
 576     Int_t i2;
 577
 578     inputFile >> buff;
 579
 580     inputFile >> Serial;
 581     inputFile >> Ntrigger;
 582     inputFile >> Nchannel;
 583     inputFile >> i1;
 584     inputFile >> i2;
 585
 586     Float_t ped0[Nchannel];
 587     Float_t sig0[Nchannel];
 588     Float_t ped1[Nchannel];
 589     Float_t sig1[Nchannel];
 590     Int_t ichannel;
 591
 592     for (Int_t i = 0; i < Nchannel-1; i++) {
 593     ped0[i] = 0;
 594     sig0[i] = 0;
 595     ped1[i] = 0;
 596     sig1[i] = 0;
 597     }
 598
 599     for (Int_t i = 0; i < Nchannel-1; i++) {
 600         inputFile >> ichannel;
 601         inputFile >> ped0[i];
 602         inputFile >> sig0[i];
 603         inputFile >> ped1[i];
 604         inputFile >> sig1[i];
 605         fThreshold[i] = fNsigma*sig1[i];
 606     }
 607
 608     inputFile.close();
 609 }
 610 #endif